Non local means影象去噪演算法及其實現
論文原文: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倍噪聲標準差
參考: