高斯消元_HihoCoderOffer6_03

题目3 : 图像算子

时间限制: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

高斯消元解齐次线性方程组。A * OP = B  ，这里OP矩阵为未知矩阵。作为一个卷积算子矩阵，它完成的功能是以自身矩阵（d * d）在A上扫描并运算。比如在A上的第一个位置(1,1)上，

1 6 13　　　　　　　　　　　　　　　x1 x2 x3

13 1 5   mul直接对应位置元素相乘      x4 x5 x6      =    B[1]][1]= 22。  移动d * d次刚好 d * d个方程组，求解d*d个未知数。

8 2 15                                          x7 x8 x9

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <stack>
#define PII pair<int, int>
#define INF 0x3f3f3f3f
#define eps 1e-9
using namespace std ;

const int Maxn = 5000;

int h, w, d;
int a[Maxn][Maxn], b[Maxn][Maxn];
//Mat:系数矩阵 V:解向量
double Mat[Maxn][Maxn], V[Maxn];
//Gauss消元：n是行列式个数，m是未知变量个数
void Gauss(int n, int m)
{
    int k = 0, i , j;
    for(j = 0; j < m; j ++){
        for(i = k; i < n; i ++){
            if(fabs(Mat[i][j]) > eps)
                break;
        }
        if(i == n) continue;
        for(int p = 0; p < m; p ++){
            swap(Mat[i][p], Mat[k][p]);
        }swap(V[i], V[k]);
        double tmp = Mat[k][j];
        for(int p = j; p < m; p ++)
            Mat[k][p]/=tmp;
        V[k]/=tmp;
        for(int p = 0; p < n; p ++){
            if(p != k && (fabs(Mat[p][j]) > eps)){
                tmp = Mat[p][j];
                for(int q = 0; q < m; q ++)
                    Mat[p][q] -= tmp * Mat[k][q];
                V[p] -= tmp*V[k];
            }
        }
        k ++;
    }
}

int main()
{
    scanf("%d%d%d",&h,&w,&d);
    for(int i = 0; i < h; i ++){
        for(int j = 0; j < w; j ++){
            scanf("%d",&a[i][j]);
        }
    }for(int i = 0; i < h - d + 1; i ++ ){
        for(int j = 0; j < w - d + 1; j ++){
            scanf("%d",&b[i][j]);
        }
    }
    int r = 0;
    for(int i = 0; i < h - d + 1; i ++){
        for(int j = 0; j < w - d + 1; j ++){
            for(int p = 0; p < d; p ++){
                for(int q = 0; q < d; q ++){
                    Mat[r][p*d + q] = a[i + p][j + q];
                }
            }
            V[r] = b[i][j];
            r ++;
        }
    }
    Gauss(r, d*d);
    for(int i = 0; i < d*d; i ++){
        if(V[i] > -1e-6) printf("%.0f",V[i] + 1e-6);
        else printf("%0.f",V[i] - 1e-6);
        if(i % d == d - 1) printf("
");
        else printf(" ");
    }
    return 0;
}