hihocoder图像算子(高斯消元）

描述

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

假设原图为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


高斯消元的思路如下
1、找到一个基准行，这一行左起第一个非零元素应该是 列元素中最大的
2、用该基准行去把位于该行下面的所有的列（最大列元素所在的列）化成0，最后会变成一个上三角阵
3、从最后一行开始往上求解

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

using namespace std;

const int MAXN = 1e4 + 5;
int equ, var;
double a[MAXN][100];
double x[MAXN];
double A[70][70];

void guass() {

    for (int i = 0; i < equ; i++) {
        int maxRow = i;
        double maxV = fabs(a[i][i]);
        for (int j = i + 1; j < equ; j++) {
            if (maxV < fabs(a[j][i])) {
                maxV = fabs(a[j][i]);
                maxRow = j;
            }
        }

        // 交换
        for (int j = i; j <= var; j++) {
            swap(a[i][j], a[maxRow][j]);
        }

        // 化成下三角
        for (int j = i + 1; j < equ; j++) {
            double c = a[j][i] / a[i][i];
            a[j][i] = 0;
            for (int k = i + 1; k <= var; k++) {
                a[j][k] -= c * a[i][k];
            }
        }
    }

    /*
    for (int i = equ - 1; i >= 0; i--) {
        x[i] = a[i][var] / a[i][i];
        for (int j = i - 1; j >= 0; j--) {
            a[j][var] -= x[i] * a[j][i];
        }
    }
    */
    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];
    }
}

int main() {

    int H, W, D;
    while(scanf("%d%d%d", &H, &W, &D) != EOF) {
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                scanf("%lf",  &A[i][j]);
            }
        }
        var = D * D;
        equ = (H - D + 1)  * (W - D + 1);
        memset(a, 0, sizeof(a));
        memset(x, 0, sizeof(x));
        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] = A[i / (W - D + 1) + j / D][i % (W - D + 1) + j % D];
            }
        }
        /*
        for (int i = 0; i < equ; i++) {
            for (int j = 0; j < var; j++) {
                printf("%lf ", a[i][j]);
            }
            printf("
");
        }
        */
        guass();

        for (int i = 0; i < var; i++) {
            int k = floor(x[i] + 0.5);
            if (i % D == D - 1) {
                printf("%d
", k);
            } else {
                printf("%d ", k);
            }
        }
    }
    return 0;
}