# 非局部均值滤波算法（NL-means）

## 非局部均值滤波算法（NL-means）

``````//计算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<uchar>(j);
uchar *data2 = B.ptr<uchar>(j);
for (int i = 0; i < A.cols; i++)
{
sum += table1[table2[data1[i]][data2[i]]];
}
}
sum = sum / (A.rows*B.cols);
return sum;
}
``````

``````//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<uchar>(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<uchar>(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<uchar>(p / sumw);
}
}
}
``````

``````void main()
{
Mat dst;
NL_mean(src, dst, 25, 2, 10);   //NL-means
imshow("原图", src);
imshow("非局部均值滤波", dst);
waitKey(0);
}
``````

NL_means算法效果很好，但是运行速度较慢，需要对算法进行加速，下面还看到一种积分图加速的方法：

MATLAB代码实现：

``````clc;
clear all;
close all;

%---------------------------%
%       input
%---------------------------%
src=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;

%---------------------------%
%  Non-Local Means Denoising
%---------------------------%
sumimage=zeros(m,n);
sumweight=zeros(m,n);
maxweight=zeros(m,n);
[M,N]=size(image);
for r=-Ds:Ds
for s=-Ds:Ds

%跳过当前点偏移
if(r==0&&s==0)
continue;
end

%求得差值积分图
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');
``````

THE END