HiHoCoder [Offer收割]编程练习赛6 C. 图像算子（高斯消元小数版）

传送门

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

解题思路：

其实将这个题目转化一下的话 就是一个高斯消元的问题，首先我们按照题目意思列出

$（W-D+1）*（H-D+1）$ 个方程，包含 $D*D$ 个未知数的方程，然后通过高斯消元

模板，就可以搞定了。

具体的说一下怎么构造的这些方程组，就拿样例来说吧：

$$A: \ \ \ \begin{matrix} 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 \\ \end{matrix}$$

因为 $D$ 矩阵的长度是 $3$ ，所以 $B$ 矩阵的的规模是 $(5-3+1)\ * (5-3+1)$ 的，嗯哼，那么

就是相当于 $A$ 矩阵的每个子矩阵中的数分别乘以 $D$ 矩阵的数，这个其实说不太明白 最好是在纸上

画一画基本上就出来了，分别对应相乘完之后，得到 $(H-D+1)*(W-D+1)$ 个方程和

$D*D$ 个未知数，在这里需要注意的是，构造的增广矩阵是 $double$ 的，解也是 $double$ 的。具体还是 看我的代码吧。

$My\ Code：$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
using namespace std;

typedef long long LL;
const double eps = 1e-7;
const int MAXN = 1e4+5;
int equ, var;///equ个方程 var个变量
double a[MAXN][100];///增广矩阵
double x[MAXN];///解的数目
void Gauss()
{
    int Max_r;///当前列绝对值最大的存在的行
    ///col：处理当前的列
    int row = 0;
    int free_x_num;
    int free_index;
    for(int col=0; row<equ&&col<var; row++,col++)
    {
        Max_r = row;
        for(int i=row+1; i<equ; i++)
            if(fabs(a[i][col]) > fabs(a[Max_r][col]))
                Max_r = i;

        if(Max_r != row)
            for(int i=0; i<var+1; i++)
                swap(a[row][i], a[Max_r][i]);
        for(int i=row+1; i<equ; i++)
        {
            double f = a[i][row]/a[row][col];
            for(int j=col; j<var+1; j++)
                a[i][j] -= f*a[row][j];
        }
    }
    for(int i=equ-1; i>=0; i--)
    {
        double tmp = a[i][var];
        for(int j=i+1; j<var; j++)
            tmp -= a[i][j]*x[j];
        x[i] = tmp/a[i][i];
    }
}
double aa[70][70];
int main()
{
    int H, W, D;
    while(~scanf("%d%d%d",&H,&W,&D))
    {
        for(int i=0; i<H; i++)
            for(int j=0; j<W; j++)
                scanf("%lf",&aa[i][j]);
        memset(a, 0, sizeof(a));
        memset(x, 0, sizeof(x));
        equ = (H-D+1)*(W-D+1);
        var = D*D;
        for(int i=0; i<equ; i++)
            scanf("%lf",&a[i][var]);
        for(int i=0; i<equ; i++)
            for(int j=0; j<var; j++)
                a[i][j] = aa[i/(W-D+1)+j/D][i%(W-D+1)+j%D];
        Gauss();
        for(int i=0; i<D*D; i++)
        {
            int k = floor(x[i]+0.5);
            if(i%D == D-1)
                printf("%d\n",k);
            else
                printf("%d ",k);
        }
    }
    return 0;
}