hdu 3359 /poj 3999 Kind of a Blur(高斯消元模板题)

本文链接：https://blog.youkuaiyun.com/acm_cxq/article/details/52207169

这是一篇关于使用高斯消元方法解决ACM竞赛题目HDU 3359 / poj 3999 Kind of a Blur的文章。题目要求逆向操作给定转换规则，通过建立方程组并进行高斯消元求解原矩阵。然而，作者在POJ上用G++编译器遇到了问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Kind of a Blur

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 2242 Accepted Submission(s): 609

Problem Description

Image blurring occurs when the object being captured is out of the camera's focus. The top two figures on the right are an example of an image and its blurred version. Restoring the original image given only the blurred version is one of the most interesting topics in image processing. This process is called deblurring, which will be your task for this problem.
In this problem, all images are in grey-scale (no colours). Images are represented as a 2 dimensional matrix of real numbers, where each cell corresponds to the brightness of the corresponding pixel. Although not mathematically accurate, one way to describe a blurred image is through averaging all the pixels that are within (less than or equal to) a certain Manhattan distance?from each pixel (including the pixel itself ). Here's an example of how to calculate the blurring of a 3x3 image with a blurring distance of 1:

Given the blurred version of an image, we are interested in reconstructing the original version assuming that the image was blurred as explained above.

Input

Input consists of several test cases. Each case is specified on H + 1 lines. The first line specifies three non negative integers specifying the width W, the height H of the blurred image and the blurring distance D respectively where (1<= W,H <= 10) and (D <= min(W/2,H/2)). The remaining H lines specify the gray-level of each pixel in the blurred image. Each line specifies W non-negative real numbers given up to the 2nd decimal place. The value of all the given real numbers will be less than 100.
Zero or more lines (made entirely of white spaces) may appear between cases. The last line of the input file consists of three zeros.

Output

For each test case, print a W * H matrix of real numbers specifying the deblurred version of the image. Each element in the matrix should be approximated to 2 decimal places and right justified in a field of width 8. Separate the output of each two consecutive test cases by an empty line. Do not print an empty line after the last test case. It is guaranteed that there is exactly one unique solution for every test case.

Sample Input

  
  
   
   2 2 1
1 1
1 1

3 3 1
19 14 20
12 15 18
13 14 16

4 4 2
14 15 14 15
14 15 14 15
14 15 14 15
14 15 14 15

0 0 0

Sample Output

  
  
   
       1.00    1.00
    1.00    1.00

    2.00   30.00   17.00
   25.00    7.00   13.00
   14.00    0.00   35.00

    1.00   27.00    2.00   28.00
   21.00   12.00   17.00    8.00
   21.00   12.00   17.00    8.00
    1.00   27.00    2.00   28.00

   
   
    
    
     
     Hint
    
    
The Manhattan Distance (sometimes called the Taxicab distance) between
two points is the sum of the (absolute) difference of their coordinates. 
The grid on the lower right illustrates the Manhattan distances from the grayed cell.

题意：跟上图的转换规则一样，求逆过程求出来的原矩阵

思路：可以写出n*m个方程，包含n*m个变元，跑一边高斯消元即可

不知为何poj的G++过不了

代码：

#include<stdio.h>
#include<algorithm>
#include<iostream>
#include<string.h>
#include<math.h>
using namespace std;
#define eps 1e-8
const int MAXN=205;
int a[MAXN][MAXN];///增广矩阵
double s[MAXN][MAXN];
double x[MAXN];///解集
bool free_x[MAXN];
inline int gcd(int a,int b)
{
    int t;
    while(b!=0)
    {
        t=b;
        b=a%b;
        a=t;
    }
    return a;
}
inline int lcm(int a,int b)
{
    return a/gcd(a,b)*b;//先除后乘防溢出
}
int Gauss(int equ,int var)
{
    int i,j,k;
    int max_r;// 当前这列绝对值最大的行.
    int col;//当前处理的列
    int ta,tb;
    int LCM;
    int temp;
    int free_x_num;
    int free_index;

    for(int i=0;i<=var;i++)
    {
        x[i]=0;
        free_x[i]=true;
    }

    //转换为阶梯阵.
    col=0; // 当前处理的列
    for(k = 0;k < equ && col < var;k++,col++)
    {// 枚举当前处理的行.
// 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
        max_r=k;
        for(i=k+1;i<equ;i++)
        {
            if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
        }
        if(max_r!=k)
        {// 与第k行交换.
            for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
        }
        if(a[k][col]==0)
        {// 说明该col列第k行以下全是0了，则处理当前行的下一列.
            k--;
            continue;
        }
        for(i=k+1;i<equ;i++)
        {// 枚举要删去的行.
            if(a[i][col]!=0)
            {
                LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                ta = LCM/abs(a[i][col]);
                tb = LCM/abs(a[k][col]);
                if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
                for(j=col;j<var+1;j++)
                {
                    a[i][j] = a[i][j]*ta-a[k][j]*tb;
                }
            }
        }
    }

  //  Debug();

    // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
    // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
    // 且出现的行数即为自由变元的个数.
    if (k < var)
    {
        // 首先，自由变元有var - k个，即不确定的变元至少有var - k个.
        for (i = k - 1; i >= 0; i--)
        {
            // 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.
            // 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.
            free_x_num = 0; // 用于判断该行中的不确定的变元的个数，如果超过1个，则无法求解，它们仍然为不确定的变元.
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
            }
            if (free_x_num > 1) continue; // 无法求解出确定的变元.
            // 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.
            temp = a[i][var];
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
            }
            x[free_index] = temp / a[i][free_index]; // 求出该变元.
            free_x[free_index] = 0; // 该变元是确定的.
        }
    }
    // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
    // 计算出Xn-1, Xn-2 ... X0.
    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (a[i][j] != 0) temp -= a[i][j] * x[j];
        }
        x[i] = temp / a[i][i];
    }
    return 0;
}
int main(void)
{
    int t=0;
    int n,m,d;
    while (~scanf("%d %d %d",&m,&n,&d)&&(n+m+d))
    {
        if(t++) printf("\n");
        memset(a,0,sizeof(a));
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
                scanf("%lf",&s[i][j]);
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
            {
                int num=0;
                for(int k=0; k<n; k++)
                    for(int l=0; l<m; l++)
                        if(abs(k-i)+abs(l-j)<=d)
                        {
                            a[i*m+j][k*m+l]=1;
                            num++;
                        }
                x[i*m+j]=num*s[i][j];
            }
        Gauss(n*m,n*m);
        for(int i=0; i<n; i++)
        {
            for(int j=0; j<m; j++)
                printf("%8.2lf",x[i*m+j]);
            printf("\n");
        }
    }
    return 0;
}