Non local means图像去噪算法及其实现

论文原文：A non-local algorithm for image denoising

该文章2005由Buades等人发表在CVPR上，对于single-image denoise来说，当时基本上是state-of-the-art。

去噪属于图像复原的范畴，通常使用滤波来实现，并且往往是低通（平滑噪声）滤波器。对于单帧图像去噪，使用空间邻域像素来处理，对于多帧图像去噪，则可以考虑时空域相结合的方法，即时间+空间的3DNR方法。

简单的平滑滤波器有均值滤波器、高斯滤波器，算法复杂度低，但会导致图像模糊，双边滤波器是性能较好的非线性滤波器，在去噪的同时，保留了较强的纹理细节，缺点是弱的纹理被滤掉了。非局部均值（Non Local Means）方法，不仅仅考虑了像素的较小邻域，并且邻域点的权重由该点与滤波点相似度计算得到。

1. NLM方法定义

————————————————

NLM去噪后输出图像定义如下：

其中I为搜索区域，w(i,j)为权重，由匹配块的相似度决定。

块的相似度定义如下：

该值由表示点i和j邻域差值平方卷积高斯核，表征邻域相似度，Z(i)表示权重归一化系数。

2. 算法质量评估

使用该方法时，往往会选择一个搜索窗口和匹配窗口，典型值分别为21x21和7x7。

对于图像去噪、复原类问题，经常使用PSNR（峰值信噪比）来评价质量，单位为dB，峰值信号定义如下，L为图像最高灰度值，对于不同位宽的图像，峰值信号大小不一样。

MSE（mean squared error）为均方误差，定义如下：

 

3. Matlab实现

Matlab代码实现如下：

(1) 顶层代码

clc;
clear all;
close all;
tic;

imgSrc = imread('E:\00_lemonHe\01_code\matlab\04_digitalImageProcessing\learn\lena.tif');
sigma = 15;
imgRandNoise = imgSrc + uint8(sigma * randn(size(imgSrc)));     %add rand noise
% h = 10 * sqrt(sigma);
h = sigma;
imgDst = NLmeansfilter(double(imgRandNoise), 10, 3, h);
imgDst = uint8(imgDst);
figure;
subplot(2,2,1), imshow(imgSrc), title('original');
subplot(2,2,2), imshow(imgRandNoise), title('noisyImage');
subplot(2,2,3), imshow(imgDst), title('denoising');
subplot(2,2,4), imshow(imgRandNoise - imgDst), title('noisy');

filterGaussian = fspecial('gaussian');      %3x3 gaussian filter
filterAverage = fspecial('average');        %3x3 average filter
imgGaussian = imfilter(imgRandNoise, filterGaussian, 'replicate');
imgAverage = imfilter(imgRandNoise, filterAverage, 'replicate');

figure;
subplot(2,2,1),imshow(imgSrc);
subplot(2,2,2),imshow(imgDst),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgDst, 8))]);
subplot(2,2,3),imshow(imgGaussian),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgGaussian, 8))]);
subplot(2,2,4),imshow(imgAverage),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgAverage, 8))]);

toc;

(2). NLM处理函数

function [output] = NLmeansfilter(input,t,f,h)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %  https://cn.mathworks.com/matlabcentral/fileexchange/13176-non-local-means-filter
 %  input: image to be filtered
 %  t: radius of search window
 %  f: radius of similarity window
 %  h: degree of filtering
 %
 %  Author: Jose Vicente Manjon Herrera & Antoni Buades
 %  Date: 09-03-2006
 %  Modify: lemonHe
 %  Modify Date: 09-15-2017
 %
 %  Implementation of the Non local filter proposed for A. Buades, B. Coll and J.M. Morel in
 %  "A non-local algorithm for image denoising"
 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Size of the image
[m n]=size(input);
% Memory for the output
output=zeros(m,n);

% Replicate the boundaries of the input image
input2 = padarray(input,[f f],'replicate');

% Used kernel
% kernel = make_kernel(f);
% kernel = kernel / sum(sum(kernel));
kernel = fspecial('gaussian', [2*f+1 2*f+1], 1);    %gaussian filter kernel
hwait=waitbar(0,'计算中...');                       %由于计算太慢，此处加入进度bar
for i=1:m
    for j=1:n
        value = 100 * i / m;
        waitbar(i/m, hwait, sprintf('计算中:%3.2f%%',value));
        i1 = i + f;
        j1 = j + f;
        W1= input2(i1-f:i1+f , j1-f:j1+f);      %输入图像的2f+1领域
        wmax = 0;
        average = 0;
        sweight = 0;
        rmin = max(i1-t,f+1);
        rmax = min(i1+t,m+f);
        smin = max(j1-t,f+1);
        smax = min(j1+t,n+f);
        %遍历21x21的search领域
        for r = rmin:1:rmax
            for s = smin:1:smax
                if(r==i1 && s==j1)
                    continue;
                end
                W2= input2(r-f:r+f , s-f:s+f);          %search window中的2f+1领域
                d = sum(sum(kernel.*((W1-W2).^2)));     %两个2f+1邻域的高斯加权欧式距离
                w = exp(-d/(h^2));                      %像素点(r-f,s-f)的权重系数Z(i)
                if w>wmax
                    wmax = w;
                end
                sweight = sweight + w;                  %除了像素点(i,j)外的所有点的权值和
                average = average + w * input2(r,s);    %除了像素点(i,j)外的所有点的加权值
            end
        end
        average = average + wmax*input2(i1,j1);         %search window的加权和
        sweight = sweight + wmax;                       %search window的权值和   
        if sweight > 0
            output(i,j) = uint8(average / sweight);
        else
            output(i,j) = input(i,j);
        end
    end
end
close(hwait);

(3). 高斯滤波窗口函数

function [kernel] = make_kernel(f)

kernel=zeros(2*f+1,2*f+1);
for d=1:f
    value= 1 / (2*d+1)^2;
    for i=-d:d
        for j=-d:d
            kernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value;
        end
    end
end
kernel = kernel ./ f;

(4) PSNR计算函数

function [imgPSNR] = calcPSNR(imgSrc, imgDst, bitWidth)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% imgSrc: original image
% imgDst: dst image
% 
% Author: lemonHe
% Date: 20170920
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(size(imgSrc,3) == 1)
    imgMSE = sum(sum((double(imgDst) - double(imgSrc)).^2)) / (size(imgSrc,1) * size(imgSrc,2));
    imgPSNR = 10 * log10((2^bitWidth-1)^2 / imgMSE);
else
    imgSrcGray = rgb2ycbcr(imgSrc);
    imgDstGray = rgb2ycbcr(imgDst);
    imgMSE = sum(sum((double(imgDstGray(:,:,1)) - double(imgSrcGray(:,:,1))).^2)) / (size(imgSrc,1) * size(imgSrc,2));
    imgPSNR = 10 * log10((2^bitWidth-1)^2 / imgMSE);
end

4. 效果

使用21x21搜索窗口，7x7匹配窗口进行测试，效果如下，从左到右分别为原图、带噪声图像、NLM去噪图像和噪声图像

添加PSNR打印，并对比高斯和均值去噪，结果如下，从左到右分别为原图、NLM去噪图像、gaussian去噪图像和均值去噪图像，NLM方法的PSNR（31.9512 dB）要比gaussian和average去噪的高。

使用5x5匹配窗口进行匹配，效果如下，PSNR仍然有30.9469 dB。

 

作者论文中有张图很清晰地显示出了算法原理，以图e为例，左图中白色亮点标记的点在该图中进行搜索，找到的权重较大的点，见右图中白色亮点，正好说明了滤波点的权重与邻域相似度相关。

该算法对纹理区域及周期性结构区域的去噪效果比较好，原因在于对这类图像，在search window中能找到多个比较相似的邻域，进而较好的抑制噪声。

5. 复杂度及总结

图像大小为MxN，搜索范围为S，匹配范围为F，那么算法复杂度为

NLM方法总结：

(1). 确定搜索半径，作者论文中建议使用21x21的搜索窗口

(2). 确定匹配块半径，论文中建议使用7x7窗口，我尝试使用5x5，效果不会差太多，速度拉升近一倍

(3). 设定好滤波参数h，论文建议使用10倍噪声标准差

 

参考：

[1] https://web.stanford.edu/class/ee367/reading/A%20non-local%20algorithm%20for%20image%20denoising.pdf

[2] https://en.wikipedia.org/wiki/Non-local_means