积分图像的应用（二）：非局部均值去噪（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