最小二乘法 图像去噪（空间域）

简单理解：从一个3*3的邻域里，使用x,y坐标，作为输入（f(x,y) = ax2 + by2 + cxt + dx + ey + f ），x,y对应像素作为输出，构建一个矩阵（9*6）乘以一个未知数矩阵（6*1）得到输出矩阵（9*1），利用高斯消元法求出未知数f（因为要求的像素为x,y坐标为0）,即可。f为消去噪声的像素值。

 

我们注意到，对于原点（0,0）, 所以很自然的用高斯消去方法得到f就是我们所求的该点应该的像素值，不需要再去求其他五个参数了(这也是为什么我们要选用参照坐标系的原因，可以反复套用)，对于每一个像素值依次处理，即可得到结果.

注:对于边界值，可以直接取0填充.

以下是处理结果对比:

此为CPU版本，但是可以用GPU cuda加速处理，先放个坑等不忙的时候把GPU版本补上放Git链接.

核心部分代码

//选取的图片大小是512X512size的，构建mask
for(j=0;j<ydim;j++){
    for(i=0;i<xdim;i++){
      for(h=0;h<3;h++){
        for(g=0;g<3;g++){
          y=j+h-1;
          x=i+g-1;
          if(x==-1 || y==-1 || x==512 || y==512){mask[h*3+g]=0;}
          else{mask[h*3+g]=image[y*xdim+x];}
        }
      }
      
//ATAx=ATb,计算矩阵A的转置矩阵和b的乘积
      for(int k=0;k<6;k++){
        B[k]=0;
        for(int p=0;p<9;p++){
          B[k]+=b[k][p]*mask[p];
        }
      }
//高斯消去
      for(int p=0;p<5;p++){
        if(!c[p][p]){
          return -1;
        }
        for(int q=p+1;q<6;q++){
            forward=c[q][p]/c[p][p];
            for(int r=p;r<6;r++){
              c[q][p]=c[q][p]-forward*c[p][r];
            }
            B[q]=B[q]-forward*B[p];
        }
      }
      result=B[5]/c[5][5];
      result=(int)result;
      if(result<0){
        image[j*xdim+i]=0;
      }
      else if(result>255){
        image[j*xdim+i]=255;
      }
      else{
        image[j*xdim+i]=result;
      }