理论上，该算法需要在整个图像范围内判断像素间的相似度，也就是说，每处理一个像素点时，都要计算它与图像中所有像素点间的相似度。但是考虑到效率问题，实现的时候，会设定两个固定大小的窗口：搜索窗(2*half_ssize+1)*(2*half_ssize+1)和邻域窗口(2*half_ksize+1)*(2*half_ksize+1)。邻域窗口在搜索窗口中滑动，根据邻域间的相似性确定像素的权值。
点A的滤波值由搜索窗口内所有点的像素值加权平均得到：
上式中，w(A,B)为点A与搜索窗口中任意一点B的高斯权重，由两点各自邻域块的MSE相似度计算而得。如下式，w(A,B)需要除以搜索窗口内所有点的高斯系数之和，以实现归一化的目的。其中h也是一个重要的参数，h越大去噪效果越好，但是图像越模糊，反之h越小去噪效果越差，但去噪之后的失真度越小。
代码实现类似均值+高斯滤波过程，即非局部均值滤波算法：
//sigma越大越平滑//halfKernelSize邻域窗//halfSearchSize搜索窗voidNL_mean(Matsrc, Mat&dst, doublesigma, inthalfKernelSize, inthalfSearchSize) { MatboardSrc; dst.create(src.rows, src.cols, CV_8UC1); intboardSize=halfKernelSize+halfSearchSize; copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT); //边界扩展doublesigma_square=sigma*sigma; introws=src.rows; intcols=src.cols; cal_lookup_table1(); cal_lookup_table2(); for (intj=boardSize; j<boardSize+rows; j++) { uchar*dst_p=dst.ptr<uchar>(j-boardSize); for (inti=boardSize; i<boardSize+cols; i++) { MatpatchA=boardSrc(Range(j-halfKernelSize, j+halfKernelSize), Range(i-halfKernelSize, i+halfKernelSize)); doublew=0; doublep=0; doublesumw=0; for (intsr=-halfSearchSize; sr<=halfSearchSize; sr++) //在搜索框内滑动 { uchar*boardSrc_p=boardSrc.ptr<uchar>(j+sr); for (intsc=-halfSearchSize; sc<=halfSearchSize; sc++) { MatpatchB=boardSrc(Range(j+sr-halfKernelSize, j+sr+halfKernelSize), Range(i+sc-halfKernelSize, i+sc+halfKernelSize)); floatd2=MSE_block(patchA, patchB); w=exp(-d2/sigma_square); p+=boardSrc_p[i+sc] *w; sumw+=w; } } dst_p[i-boardSize] =saturate_cast<uchar>(p/sumw); } } } ```主函数参数设置：标准差25，邻域窗2X2，搜索窗10X10```cppvoidmain() { Matsrc=imread("原图.jpg", CV_LOAD_IMAGE_GRAYSCALE); Matdst; NL_mean(src, dst, 25, 2, 10); //NL-meansimshow("原图", src); imshow("非局部均值滤波", dst); waitKey(0); }
NL_means算法效果很好，但是运行速度较慢，需要对算法进行加速，下面还看到一种积分图加速的方法：
积分图加速本质是，将图像的所有点，计算某一偏离坐标点方向的权重，一次性计算出来，而不是单次计算某一点的所有偏离点。具体，计算当前图像与一定偏移后的图像的diff并平方，继而计算与原图像同等尺寸的积分图，则每一个点与偏离它在一定偏移的坐标点的距离，通过积分图就可以计算出来了。积分图的作用体现在使用线性计算，减少重复计算，即存储换时间。
MATLAB代码实现：
clc; clearall; closeall; %---------------------------%%input%---------------------------%src=imread('lena.jpg'); src=rgb2gray(src); src=double(src); figure,imshow(src,[]),title('src image') %---------------------------%%parameter%---------------------------%[m,n]=size(src); ds=2;%blocksizeforcalculateweightDs=5;%searchblockh=10;%decayfactoroffset=ds+Ds; PaddedImg=padarray(src,[ds+Ds,ds+Ds],'symmetric','both');%扩展图像，便于处理%---------------------------%%Non-LocalMeansDenoising%---------------------------%sumimage=zeros(m,n); sumweight=zeros(m,n); maxweight=zeros(m,n); image=PaddedImg(1+Ds:Ds+m+ds,1+Ds:Ds+n+ds); [M,N]=size(image); forr=-Ds:Dsfors=-Ds:Ds%跳过当前点偏移if(r==0&&s==0) continue; end%求得差值积分图wimage=PaddedImg(1+Ds+r:Ds+m+ds+r,1+Ds+s:Ds+n+ds+s); diff=image-wimage; diff=diff.^2; J=cumsum(diff,1); J=cumsum(J,2); %计算距离distance=J(M-m+1:M,N-n+1:N)+J(1:m,1:n)-J(M-m+1:M,1:n)-J(1:m,N-n+1:N); distance=distance/((2*ds+1).^2); %计算权重并获得单个偏移下的加权图像weight=exp(-distance./(h*h)); sumimage=sumimage+weight.*wimage(ds+1:ds+m,ds+1:ds+n); sumweight=sumweight+weight; maxweight=max(maxweight,weight); endendsumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n); sumweight=sumweight+maxweight; dst=sumimage./sumweight; %---------------------------%%output%---------------------------%figure,imshow(dst,[]),title('dst');
效果图如下：
