图像算法|非局部均值滤波算法（NL-means） c++|matlab|图像处理

非局部均值滤波算法（NL-means）今天来学习一下另一类滤波算法：非局部均值滤波算法（NL-means）。非局部均值滤波算法最早于2005年由Buades等人发表在CVPR上，论文原文：A non-local algorithm for image denoising，还有一篇2011年的论文：Non-Local Means Denoising。之后还会继续介绍DCT（离散余弦变换滤波）、TV（全变分滤波）、BM3D（3维块匹配滤波）等算法。
上一篇文章介绍了均值滤波、中值滤波、高斯滤波、双边滤波和引导滤波，均为局部滤波，即对周围邻域的点加权生成当前点，加权因子反应出周围点对当前点的影响。非局部滤波就意味着它使用图像中的所有像素，这些像素根据某种相似度进行加权平均。滤波后图像清晰度高，而且不丢失细节。
与上述算法思路不同，前者根据像素点之间的距离计算权重，后者根据块之间的相似度计算权重。其核心思路与高斯滤波很相似：计算矩形窗口内所有像素点的像素值加权和，权重服从高斯分布。区别在于：高斯滤波使用当前滤波点与矩形窗口内其它点的空间欧式距离来计算权重，距离越近权重越大；而非局部均值滤波则使用当前滤波点的邻域块与矩形窗口内其它点的邻域块的相似度来计算权重，相似度越大则权重越大。
文章图片

左边是高斯滤波空间距离d，右边是非NLMeans相似度MSE。要衡量两个邻域块的相似度，有多种指标，均方误差（MSE）是最常用的相似度衡量指标之一。非局部均值滤波算法就是使用MSE来计算两个邻域块的相似度。邻域块的行列相同，假设均为m行n列，那么MSE的计算如下式，其中A(i,j)为点A邻域块中的点(i,j)的像素值，B(i,j)为点B邻域块中的相同位置点的像素值。

文章图片

代码实现计算两个区域A,B的MSE：

//计算0~255的平方查找表 float table1[256]; static void cal_lookup_table1(void) { for (int i = 0; i < 256; i++) { table1[i] = (float)(i*i); } }//计算两个0~255的数的绝对差值的查找表 uchar table2[256][256]; static void cal_lookup_table2(void) { for (int i = 0; i < 256; i++) { for (int j = i; j < 256; j++) { table2[i][j] = abs(i - j); table2[j][i] = table2[i][j]; } } }float MSE_block(Mat A, Mat B) { float sum = 0.0; for (int j = 0; j < A.rows; j++) { uchar *data1 = A.ptr(j); uchar *data2 = B.ptr(j); for (int i = 0; i < A.cols; i++) { sum += table1[table2[data1[i]][data2[i]]]; } } sum = sum / (A.rows*B.cols); return sum; }

理论上，该算法需要在整个图像范围内判断像素间的相似度，也就是说，每处理一个像素点时，都要计算它与图像中所有像素点间的相似度。但是考虑到效率问题，实现的时候，会设定两个固定大小的窗口：搜索窗(2half_ssize+1)(2half_ssize+1)和邻域窗口(2half_ksize+1)(2half_ksize+1)。邻域窗口在搜索窗口中滑动，根据邻域间的相似性确定像素的权值。

文章图片

点A的滤波值由搜索窗口内所有点的像素值加权平均得到：

文章图片

上式中，w(A,B)为点A与搜索窗口中任意一点B的高斯权重，由两点各自邻域块的MSE相似度计算而得。如下式，w(A,B)需要除以搜索窗口内所有点的高斯系数之和，以实现归一化的目的。其中h也是一个重要的参数，h越大去噪效果越好，但是图像越模糊，反之h越小去噪效果越差，但去噪之后的失真度越小。

文章图片

代码实现类似均值+高斯滤波过程，即非局部均值滤波算法：

//sigma越大越平滑 //halfKernelSize邻域窗 //halfSearchSize搜索窗 void NL_mean(Mat src, Mat &dst, double sigma, int halfKernelSize, int halfSearchSize) { Mat boardSrc; dst.create(src.rows, src.cols, CV_8UC1); int boardSize = halfKernelSize + halfSearchSize; copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT); //边界扩展 double sigma_square = sigma*sigma; int rows = src.rows; int cols = src.cols; cal_lookup_table1(); cal_lookup_table2(); for (int j = boardSize; j < boardSize + rows; j++) { uchar *dst_p = dst.ptr(j - boardSize); for (int i = boardSize; i < boardSize + cols; i++) { Mat patchA = boardSrc(Range(j - halfKernelSize, j + halfKernelSize), Range(i - halfKernelSize, i + halfKernelSize)); double w = 0; double p = 0; double sumw = 0; for (int sr = -halfSearchSize; sr <= halfSearchSize; sr++)//在搜索框内滑动 { uchar *boardSrc_p = boardSrc.ptr(j + sr); for (int sc = -halfSearchSize; sc <= halfSearchSize; sc++) { Mat patchB = boardSrc(Range(j + sr - halfKernelSize, j + sr + halfKernelSize), Range(i + sc - halfKernelSize, i + sc + halfKernelSize)); float d2 = MSE_block(patchA, patchB); w = exp(-d2 / sigma_square); p += boardSrc_p[i + sc] * w; sumw += w; } } dst_p[i - boardSize] = saturate_cast(p / sumw); } } }

主函数
参数设置：标准差25，邻域窗2X2，搜索窗10X10

void main() { Mat src = https://www.it610.com/article/imread("原图.jpg", CV_LOAD_IMAGE_GRAYSCALE); Mat dst; NL_mean(src, dst, 25, 2, 10); //NL-means imshow("原图", src); imshow("非局部均值滤波", dst); waitKey(0); }

NL_means算法效果很好，但是运行速度较慢，需要对算法进行加速，下面还看到一种积分图加速的方法：
积分图加速本质是，将图像的所有点，计算某一偏离坐标点方向的权重，一次性计算出来，而不是单次计算某一点的所有偏离点。具体，计算当前图像与一定偏移后的图像的diff并平方，继而计算与原图像同等尺寸的积分图，则每一个点与偏离它在一定偏移的坐标点的距离，通过积分图就可以计算出来了。积分图的作用体现在使用线性计算，减少重复计算，即存储换时间。
【图像算法|非局部均值滤波算法（NL-means）】MATLAB代码实现：

clc; clear all; close all; %---------------------------% %input %---------------------------% src=https://www.it610.com/article/imread('lena.jpg'); src=https://www.it610.com/article/rgb2gray(src); src=double(src); figure,imshow(src,[]),title('src image') %---------------------------% %parameter %---------------------------% [m,n]=size(src); ds=2; % block size for calculate weight Ds=5; % search block h=10; % decay factor offset=ds+Ds; PaddedImg = padarray(src,[ds+Ds,ds+Ds],'symmetric','both'); % 扩展图像，便于处理 %---------------------------% %Non-Local Means Denoising %---------------------------% 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); for r=-Ds:Ds for s=-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); end end sumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n); sumweight=sumweight+maxweight; dst=sumimage./sumweight; %---------------------------% %output %---------------------------% figure,imshow(dst,[]),title('dst');