非局部均值图像恢复的相关理解

        图像的恢复问题旨在从一个退化的观测Y中恢复一幅高质量的图像X。图像重建问题是一个经典的不适定问题，通常此模型可以用

Y =DHX + V来表示。这里X是用来估计的未知高质量图像。H与D都是退化操作算子，V是额外的噪声。

        当H与D都是单位矩阵时，IR问题就成了去噪问题；当D是单位矩阵而H是模糊算子时，IR问题就是去模糊问题；当D是单位矩阵而H是随意的投影矩阵，IR问题则变成压缩感知问题；当D是下采样算子而H是模糊算子，IR问题则就是单幅图像的超分辨问题。

       写这篇是想学习下非局部均值滤波的，看到了一篇很好的博文，在这里引用下(引用博文网址处):

       非局部均值是一种基于快匹配来确定滤波权值的。即先确定一个块的大小，例如7x7，然后在确定一个搜索区域，例如15x15，在15x15这个搜索区域中的每一个点，计算7x7的窗口与当前滤波点7x7窗口的绝对差值和，然后在计算一个指数函数，所有的搜索点都用指数函数计算出一个权值，当然还有权值的归一化。根据这个权值进行点的滤波操作。更具体的实现代码如下：

      

function [output]=NLmeansfilter(input,t,f,h) % 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],'symmetric');  % Used kernel kernel = make_kernel(f); kernel = kernel / sum(sum(kernel));  h=h*h;  for i=1:m for j=1:n                          i1 = i+ f;         j1 = j+ f;                         W1= input2(i1-f:i1+f , j1-f:j1+f);                  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);                  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);                                                 d = sum(sum(kernel.*(W1-W2).*(W1-W2)));                                                               w=exp(-d/h);                                                                  if w>wmax                                    wmax=w;                                   end                                sweight = sweight + w;                average = average + w*input2(r,s);                                           end          end                     average = average + wmax*input2(i1,j1);        sweight = sweight + wmax;                           if sweight > 0            output(i,j) = average / sweight;        else            output(i,j) = input(i,j);        end                 end end 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  endendkernel = kernel ./ f;

          非局部算法获得的信噪比比双边滤波略高，有时候还不如双边滤波。但是，非局部滤波是一种基于快的匹配度来计算滤波权值的，所以能获得比较好的视觉效果。然而，它的计算复杂度实在是太高了。最原始非局部均值算法是在整个图片中进行块搜索，根据块的匹配度来计算权值。实际执行过程，都会把搜索区域限定在一个局部的搜索窗口中，例如上面的代码演示的就是这个过程。

        根据快匹配的思路，人们提出了很多种有效的去噪算法，这包括BM3D去噪算法。它的去噪效果比非局部去噪效果更好。还有另外一些基于快匹配的思路，例如DCT。由于传统的非局部均值比较块的相似度时用的是时域差值比较，如果我们在DCT域去比较相似度会更合理一些。并且，网上早已经有基于DCT域非局部均值的PhotoShop插件，这个插件的速度相当的慢。对于一般720p的图像，在普通pc上都要半分钟。由此可见，这类算法，如果没有用GPU加速，要想实时，那是根本不可能的。

     于是，就有一些快速算法。主要的快速算法相关的文献如下：

     1. FAST IMAGE AND VIDEO  DENOISING VIA NON-LOCAL MEANS OF SIMILAR  NEIGHBORHOOD, Mona Mahmoudi and Guillermo Saporo

     2. FAST NON-NOCAL ALGORITHM FOR IMAGE DENOISING, Jin Wang,Yanwen Guo 等

     3. FAST NON-NOCAL ALGORITHM FOR IMAGE DENOISING, Venkateswarlu Karnati,Mithun Uliyar, Sumit Dey

     其中：

        第一篇文献的思路是用梯度和均值减少不必要的搜索点，其性能提高了不少，但还是没法处理实时视频。虽然他文献中说是拿来处理视频。

        第二篇文献的思路是用FFT来计算的，其文中声称处理30w像素只需要350ms，而相对以原始的非局部方法28.16s，加速了80倍。对这个结果，我个人持怀疑态度。首先，他文中描述的原始NL算法的时间28.16s就有问题。在我们机器上，跑上述matlab代码，需要2分钟。一种可能的解释是这篇文献的FFT代码可能是经过GPU加速过的。对于浙大的这篇文章，求实新闻网上曾津津乐道的报道过。

       第三篇文献的思路，是用积分图。单纯的积分图，明显不行。但这篇文章的作者用的是多分辨率的积分图。

     

       我觉得第三篇文献的思路是比较靠谱的，所以就实现了下。但是发现这篇文献的里面有一个公式5可能有点问题。如果按照公式5编程实现，效果相当的不好。于是，我根据文章的思路稍微改写了一下公式5，效果就可以了。今后一些列的改进和优化，我的实现最终处理结果在我的I3机器上，处理30w像素的时间是700ms。本来想移植到iphone上，但是700ms的处理时间也没法实时的处理视频，后来就放弃了。