积分图像的应用(二):非局部均值去噪(NL-means)

非局部均值去噪(NL-means)一文介绍了NL-means基本算法,同时指出了该算法效率低的问题,本文将使用积分图像技术对该算法进行加速。

假设图像共像个素点,搜索窗口大小,领域窗口大小, 计算两个矩形邻域间相似度的时间为,对于每个像素点需要计算它与搜索窗口内个像素间的相似度,故NL-means复杂度为 。

经过分析可以发现,该算法可以提高之处只有邻域间相似度的计算,即耗时的操作。基本算法中,每次计算邻域间距离时都需要遍历两个邻域,逐对像素点求差值。

如果我们先构造一个关于像素差值的积分图像:

其中

这样在计算两个邻域和  间的距离时,就可以在常量时间内完成:

这样,整个算法复杂度将降为 。

具体的算法描述可以参考[1]中:

 

为了降低空间复杂度,上述算法将偏移量作为最外层循环,即每次只需要在一个偏移方向上求取积分图像,并对该积分图像进行处理。而不需要一次性求取出所有积分图像。

程序:

close all;
clear all;
clc
I=double(imread('lena.tif'));
I=I+10*randn(size(I));
tic
O1=NLmeans(I,2,5,10);
toc
tic
O2=fastNLmeans(I,2,5,10);
toc
figure;
imshow([I,O1,O2],[]);
function DenoisedImg=fastNLmeans(I,ds,Ds,h)
%I:含噪声图像
%ds:邻域窗口半径
%Ds:搜索窗口半径
%h:高斯函数平滑参数
%DenoisedImg:去噪图像
I=double(I);
[m,n]=size(I);
PaddedImg = padarray(I,[Ds+ds+1,Ds+ds+1],'symmetric','both');
PaddedV = padarray(I,[Ds,Ds],'symmetric','both');
average=zeros(m,n);
sweight=average;
wmax=average;
h2=h*h;
d2=(2*ds+1)^2;
for t1=-Ds:Ds
    for t2=-Ds:Ds
        if(t1==0&&t2==0)
            continue;
        end
        St=integralImgSqDiff(PaddedImg,Ds,t1,t2);
        v = PaddedV(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2);
        w=zeros(m,n);
        for i=1:m
            for j=1:n
                i1=i+ds+1;
                j1=j+ds+1;
                Dist2=St(i1+ds,j1+ds)+St(i1-ds-1,j1-ds-1)-St(i1+ds,j1-ds-1)-St(i1-ds-1,j1+ds);
                Dist2=Dist2/d2;
                w(i,j)=exp(-Dist2/h2);
                sweight(i,j)=sweight(i,j)+w(i,j);
                average(i,j)=average(i,j)+w(i,j)*v(i,j);
            end
        end
        wmax=max(wmax,w);
    end
end
average=average+wmax.*I;
sweight=sweight+wmax;
DenoisedImg=average./sweight;

function Sd = integralImgSqDiff(PaddedImg,Ds,t1,t2)
%PaddedImg:边缘填充后的图像
%Ds:搜索窗口半径
%(t1,t2):偏移量
%Sd:积分图像
[m,n]=size(PaddedImg);
m1=m-2*Ds;
n1=n-2*Ds;
Sd=zeros(m1,n1);
Dist2=(PaddedImg(1+Ds:end-Ds,1+Ds:end-Ds)-PaddedImg(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2)).^2;
for i=1:m1
    for j=1:n1
         if i==1 && j==1
             Sd(i,j)=Dist2(i,j);
         elseif i==1 && j~=1
             Sd(i,j)=Sd(i,j-1)+Dist2(i,j); 
         elseif i~=1 && j==1
             Sd(i,j)=Sd(i-1,j)+Dist2(i,j);
         else
             Sd(i,j)=Dist2(i,j)+Sd(i-1,j)+Sd(i,j-1)-Sd(i-1,j-1);
         end
     end
end

结果:

三幅图像依次是含噪声原图,原始NL-means算法去噪结果、使用积分图像加速的NL-means算法去噪结果。对于256*256的lena图,原始算法耗时 36.251389s,使用积分图像加速的算法耗时 4.647372s。

当然,对于Matlab而言,若充分利用它的函数和矩阵操作,可进一步在编程上加速:

function DenoisedImg=fastNLmeans2(I,ds,Ds,h)
I=double(I);
[m,n]=size(I);
PaddedImg = padarray(I,[Ds+ds+1,Ds+ds+1],'symmetric','both');
PaddedV = padarray(I,[Ds,Ds],'symmetric','both');
average=zeros(m,n);
wmax=average;
sweight=average;
h2=h*h;
d=(2*ds+1)^2;
for t1=-Ds:Ds
    for t2=-Ds:Ds
        if(t1==0&&t2==0)
            continue;
        end
        Sd=integralImgSqDiff(PaddedImg,Ds,t1,t2);
        SqDist2=Sd(2*ds+2:end-1,2*ds+2:end-1)+Sd(1:end-2*ds-2,1:end-2*ds-2)...
               -Sd(2*ds+2:end-1,1:end-2*ds-2)-Sd(1:end-2*ds-2,2*ds+2:end-1);
        SqDist2=SqDist2/d;
        w=exp(-SqDist2/h2);
        v = PaddedV(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2);
        average=average+w.*v;
        wmax=max(wmax,w);
        sweight=sweight+w;
    end
end
average=average+wmax.*I;
average=average./(wmax+sweight);
DenoisedImg = average;

function Sd = integralImgSqDiff(PaddedImg,Ds,t1,t2)
Dist2=(PaddedImg(1+Ds:end-Ds,1+Ds:end-Ds)-PaddedImg(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2)).^2;
Sd = cumsum(Dist2,1);
Sd = cumsum(Sd,2);

使用上述fastNLmeans2函数对该lena图处理仅耗时0.416442s。

参考:

[1]FromentJ. Parameter-Free Fast Pixelwise Non-Local Means Denoising[J]. Image ProcessingOn Line, 2014, 4: 300-326

posted @ 2015-09-06 15:36  xiaoluo91  阅读(11135)  评论(4编辑  收藏  举报