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

該文章2005由Buades等人發表在CVPR上，對於single-image denoise來說，當時基本上是state-of-the-art。

去噪屬於影象復原的範疇，通常使用濾波來實現，並且往往是低通（平滑噪聲）濾波器。對於單幀影象去噪，使用空間鄰域畫素來處理，對於多幀影象去噪，則可以考慮時空域相結合的方法，即時間+空間的3DNR方法。

簡單的平滑濾波器有均值濾波器、高斯濾波器，演算法複雜度低，但會導致影象模糊，雙邊濾波器是效能較好的非線性濾波器，在去噪的同時，保留了較強的紋理細節，缺點是弱的紋理

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倍噪聲標準差

參考：