每日一题之 hiho194 图像算子

题目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

思路：

高斯消元求解。顺便复习复习一下列主元消元的算法。

#include <iostream>
#include <cmath>
#include <cstdio>
#include <vector>

using namespace std;
const int maxn=70;
const int maxm=90;

int H,W,D;

double A[maxn][maxn];
double B[maxn][maxn];
double X[maxm][maxm];
double b[maxm];


const double eps=1e-9;


// 列主元消元法
int gauss_elimination(int n,double a[][maxm],double b[])
{   
    double mx;
    int i,j,k,row;
    for ( k = 1; k <= n; ++k) { //先列后行

        for ( i = k,mx = 0; i <= n; ++i) { //找当前列的主元
            if (fabs(a[i][k]) > fabs(mx)) {
                mx = a[i][k];
                row = i;
            }
        }
        if (fabs(mx) < eps) return 0; // mx = 0 则无解

        if (row != k) { // 进行行变换将列主元所在的行换到当前的行

            for (j = k; j <= n; ++j) {
                swap(a[row][j],a[k][j]);
            }
            swap(b[row],b[k]);
        }

        //消元过程
        for (j = k+1; j <= n; ++j) { // col 
            for (a[k][j] /= mx, i = k+1; i <= n; ++i) { //row
                a[i][j] -= a[i][k] * a[k][j];
            }
        }

        for (b[k] /= mx, i = k+1; i <= n; ++i)
            b[i] -= b[k]*a[i][k];
    }
    // 回代求解 因为都归一了 所以 a[i][i] = 1;

    for ( i = n; i >= 1; --i) 
        for ( j = i+1; j <= n; ++j)
            b[i] -= a[i][j] * b[j];

    return 1;
}

int main(){

    cin>>H>>W>>D;

    for(int i = 1; i <= H; ++i) {
        for(int j = 1; j <= W; ++j){
            cin>>A[i][j];
        }
    }
    for(int i = 1; i <= H-D+1; ++i){
        for(int j = 1;j <= W-D+1; ++j){
            cin>>B[i][j];
        }
    }

    for (int i = 1 ; i <= D*D; ++i) {
        int xoffset = (i-1) / (W-D + 1); //x 的偏移量
        int yoffset = (i-1) % (W-D + 1); //y 的偏移量

        for (int j = 1; j <= D; ++j) {
            for (int k = 1; k <= D; ++k) {
                X[i][(j-1)*D+k] = A[xoffset + j][yoffset + k];
            }
        }
        b[i] = B[xoffset+1][yoffset+1];
    }

    gauss_elimination(D*D,X,b);

    for (int i = 1; i <= D; ++i) {
        for (int j = 1; j <= D; ++j)
            cout << int(round(b[(i-1)*D + j])) <<" ";
        cout << endl;
    }

    return 0;
}