hiho一下 194周（高斯消元）

题目链接

题目1 : 图像算子

时间限制: 10000ms

单点时限: 1000ms

内存限制: 256MB

描述

在图像处理的技术中，经常会用到算子与图像进行卷积运算，从而达到平滑图像或是查找边界的效果。

假设原图为H × W的矩阵A，算子矩阵为D × D的矩阵Op，则处理后的矩阵B大小为(H-D+1) × (W-D+1)。其中：

B[i][j] = ∑(A[i-1+dx][j-1+dy]*Op[dx][dy]) | (dx = 1 .. D, dy = 1 .. D), 1 ≤ i ≤ H-D+1, 1 ≤ j ≤ W-D+1

给定矩阵A和B，以及算子矩阵的边长D。你能求出算子矩阵中每个元素的值吗？

输入

第1行：3个整数，H, W, D，分别表示原图的高度和宽度，以及算子矩阵的大小。5≤H,W≤60，1≤D≤5，D一定是奇数。

第2..H+1行：每行W个整数，第i+1行第j列表示A[i][j]，0≤A[i][j]≤255

接下来H-D+1行：每行W-D+1个整数，表示B[i][j]，B[i][j]在int范围内，可能为负数。

输入保证有唯一解，并且解矩阵的每个元素都是整数。

输出

第1..D行：每行D个整数，第i行第j列表示Op[i][j]。

样例输入

5 5 3
1 6 13 10 3
13 1 5 6 15
8 2 15 0 12
19 19 17 18 18
9 18 19 5 17
22 15 6
35 -36 51
-20 3 -32

样例输出

0 1 0
1 -4 1
0 1 0

思路：开始感觉无从下手，但是你把数学方程写出来（一定要动手写！），就会发现直接用高斯消元就能解了。

AC代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

#define fbs 1e-6

using namespace std;

double a[3605][30];
double c[3605];
double ans[30];

int H,W,D;
int N,M;  //M行N列

void Swap(int i,int j)   //交换i，j行
{
    for(int k=1;k<=N;k++)
        swap(a[i][k],a[j][k]);
    swap(c[i],c[j]);
}


void Gauss()
{
    for(int i=1;i<=N;i++)  //处理出上三角矩阵
    {
        for(int j=i;j<=M;j++)
        {
            if(a[j][i]!=0)
            {
                Swap(j,i);
                break;
            }
        }

        for(int j=i+1;j<=M;j++)   //消除第i+1行到第M行的第i列
        {
            double mul=a[j][i]/a[i][i];
            for(int k=i;k<=N;k++)
                a[j][k]-=mul*a[i][k];
            c[j]-=c[i]*mul;
        }
    }

    //求出唯一解
    for(int i=N;i>0;i--)
    {
        for(int j=i+1;j<=N;j++)
        {
            c[i]-=a[i][j]*ans[j];
            a[i][j]=0;
        }
        ans[i]=c[i]/a[i][i];
    }
}

int main()
{
    scanf("%d %d %d",&H,&W,&D);
    double e[61][61];
    for(int i=1;i<=H;i++)
        for(int j=1;j<=W;j++)
            scanf("%lf",&e[i][j]);
    int x=1,y;
    for(int i=1;i<=H-D+1;i++)
        for(int j=1;j<=W-D+1;j++)
        {
            y=1;
            for(int m=0;m<D;m++)
                for(int n=0;n<D;n++)
                    a[x][y++]=e[i+m][j+n];
            x++;
        }   //寻找系数矩阵

    M=(W-D+1)*(H-D+1);
    N=D*D;

    for(int i=1;i<=D;i++)
        for(int j=1;j<=D;j++)
            scanf("%lf",&c[D*(i-1)+j]);

    Gauss();
    int n;
    for(int i=1;i<=D;i++)
    {
        n=D*(i-1)+1;
        if(ans[n]>fbs)
            printf("%d",int(ans[n]+0.5));
        else
            printf("%d",int(ans[n]-0.5));
        for(int j=2;j<=D;j++)
        {
            n=D*(i-1)+j;
            if(ans[n]>fbs)    //最后取整方式注意一下
                printf(" %d",int(ans[n]+0.5));
            else
                printf(" %d",int(ans[n]-0.5));
        }
        printf("\n");
    }

    return 0;
}