精华内容
下载资源
问答
  • 非局部均值滤波不局限于像素邻域,抗噪性能较好。本代码经过简单修改即可运行,是学习非局部均值滤波的基础,其它改进算法可以以其为基础进行修改。希望对大家有所帮助。
  • NLM非局部均值滤波

    2018-12-03 17:15:18
    NLM的源代码。使用NLM,对高斯噪声,乘性噪声,椒盐噪声进行性能测试
  • 非局部均值滤波代码,可在MATLAB上运行。已经调测。内含测试图片,方便快捷,下载即可用。非局部均值滤波代码。
  • 非局部均值滤波

    2015-03-30 13:31:16
    非局部均值滤波的效果很好, 应用也越来越广泛, 本文档时用C语言编写, 效果同MATLAB编写效果相同。
  • 基于matlab实现的非局部均值滤波器,基于matlab实现的非局部均值滤波器
  • 结合NSST和快速非局部均值滤波的图像去噪
  • 分别采用非局部均值算法、同态非局部均值算法对不同载噪比下的单目标仿真强度像进行了去噪处理,同时对比了Lee滤波和均值滤波处理效果。最后用真实的相干激光雷达强度像进行了验证。实验结果表明,采用同态两级...
  • 非局部均值滤波(NLM)的matlab程序,内含完美注释及相关文档连接。 (学习NLM自己整理所得,内含完美注释及相关文档连接)
  • 双边非局部均值滤波图像去噪算法
  • 基于开关型非局部均值滤波的指静脉图像去噪.pdf
  • 种采用3维块匹配小波变换的非局部均值滤波算法NL-3DWT (Nonlocal Filter based on 3-D Patch Matching Wavelet Transform)。该算法使用块匹配的3维抽样小波变换对极化总功率图进行预滤波,在此基础上使用边界对齐...
  • 非局部均值滤波去噪算法在目前所有滤波去噪算法中效果最为显著,但其运行速度较慢,针对该缺点提出了一种基于快速非局部均值滤波算法,运用插值算法将邻域权值计算出来,速度有了明显提升。
  • 其次对改进非局部均值滤波算法的先验信息获取方法进行了改进:对噪声图像进行提升小波变换,采用一种新型阈值函数选择低频分解系数,对高于阈值的系数进行重构得到参考图像,计算参考图像的相似度权值并将其作为改进...
  • 快速非局部均值滤波算法

    热门讨论 2012-05-24 00:31:40
    很好的非局部滤波算法,速度有较大提升 操作简单,一看就懂
  • 基于自适应非局部均值滤波的图像去块算法,王慈,刘书,JPEG和MPEG等压缩标准导致的块效应是图像降质的主要表现。以往研究提出了各种降低量化噪声的方法,但都基于噪声大小已知的假设。这�
  • 详解非局部均值滤波原理以及用MATLAB源码实现 序言 均值滤波、中值滤波、高斯滤波在滤除噪声的过程中,无可避免的使图像的边缘细节和纹理信息所被滤除。针对此问题,Buades[1]等人提出了非局部均值滤波,它使图像中...

    详解非局部均值滤波原理以及用MATLAB源码实现


    序言

    均值滤波、中值滤波、高斯滤波在滤除噪声的过程中,无可避免的使图像的边缘细节和纹理信息所被滤除。针对此问题,Buades[1]等人提出了非局部均值滤波,它使图像中的冗余信息得到了有效地利用,在滤除噪声的同时将图像的细节特征进行了最大程度的保留。

    非局部均值滤波使图像中的冗余信息得到了有效地利用,在滤除噪声的同时将图像的细节特征进行了最大程度的保留。

    非局部均值滤波利用了自然图像中的每个小块都存在关联性,与均值滤波是对邻域内的所有像素求和再平均的方法不同,它先在整幅图像中寻找相似的图像块,再根据图像块的相似度大小来赋予其不同的权值,以此来实现图像去噪。

    原理介绍

    1 借助图片

    可用N(p)来表示以待处理像素点p为中心、大小为N×N的图像块;同理,N(q1)、N(q2)、N(q3)分别代表以q1、q2、q3为中心、大小为N×N的图像块。

    由于N(q1)、N(q2)、N(q3)与N(p)的相似度均不相同,它们会通过与图像块N(p)的相似度比较后,根据其与N(p)各自的相似性被赋予不同的相似性权值,N(q1)、N(q2)、N(q3)的相似性权值分别可用w(p*,* q1)、w(p*,* q2)、w(p*,* q3)来表示。

    值得说明的是,从上图可看出,N(q1)、N(q2)与N(p)具有相同的属性,因此相似性权值w(p, q1)、w(p, q2)会大一些;而N(q3)与N(p)之间存在较大差异,相似性权值w(p, q3)会比较小。


    2 进一步用公式阐述非局部均值滤波的原理

    假定一幅离散的噪声图像为g=g{g(i)|i∈Ω},其大小为N×N,Ω代表图像邻域,i为像素索引,对噪声图像进行非局部均值滤波处理后可表达为:

    这是一个加权平均的过程。上式中分母的存在是为了使其归一化,Ωi表示中心像素为i、大小为q×q的搜索区域;w(i,j)代表赋予噪声图像g(i)的权值,数学形式为:

    h代表控制滤波平滑程度的滤波系数;N(i)和N(j)分别代表以像素中心点为i、像素中心点为j的图像块,图像块的大小都是p×p;d(i,j)表示图像块N(i)和N(j)之间的欧式加权距离,α表示高斯核的标准差,α>0。在非局部均值滤波算法中,搜索区域中将q的设置一般不会超过21,并且图像块中p的取值会比q的值小的多。

    非局部均值滤波中的一些参数选择会直接影响着滤波效果的好坏,因此是至关重要的。传统的非局部均值滤波中,用来计算两个图像块之间的欧式距离一般采用高斯加权欧式距离,上面d(i,j)可进一步表达为:

    上式中“⊗”用来表示两个矩阵之间的点乘;Gα是高斯矩阵。使用高斯加权欧式距离的一个好处是:噪声可能会影响两个图像块之间的相似性计算,采用高斯核函数卷积类似于先采用高斯滤波来滤除噪声,这样后期计算相似度会更准确。

    • 滤波参数h大小说明

    非局部均值滤波采用的是指数函数的减函数,那么滤波系数h就会对加权核函数的衰减作用有着决定性的作用。

    如果h取值较大的话,求取的权值w(i,j)的值就会趋于一致使滤波作用会等效于均值滤波,去噪的效果由于边缘细节丢失会变得比较模糊;

    若h取值较小,对噪声的滤除作用就不好从未导致图像中会有哦大量噪声存在,去噪效果不好。一般h的大小是与噪声标准差σ成正比关系的,通常将其取为噪声的标准方差,即h=σ。

    • NL-means的搜索窗口和相似窗口说明

    从理论上将,非局部均值滤波的搜索区域大小为整幅图像的大小,即以每个待处理点为中心的图像块需要和整幅图像的所有图像块都要做相似度计算来确定权值,这样计算量是巨大的。与待处理像点为中心的图像块需要进行计算相似度的图像块,我们将其叫做相似窗口。在实际应用中,为了在减少算法复杂度的同时又要尽可能找到比较多的相似度高的图像块,我们会对搜索区域的大小以及相似窗口的大小有固定的参数选择。

    在Buades等人的实验中,将搜索区域的窗口大小设为21×21,将相似窗口的大小设为7×7。由于我们的噪声图像大小为N×N,这样将参数固定后的非局部均值滤波算法的时间复杂度是49×441×N^2。事实上,尽管将搜索窗口从整幅图像的大小锁定为21×21的大小,但由于处理图像中每个像素点时,都需要计算以此像素点为中心的图像块与搜索区域中每个图像块之间的欧式距离,用来确定权值。这个过程会存在许多重复的计算,运行此滤波的效率也是相当低的。


    MATLAB源码实现

    函数

    function [output]=NLmeansfilter(input,t,f,h)
     %  input: image to be filtered
     %  t: radio of search window 搜索窗口大小,下面实际相似窗口大小为10
     %  f: radio of similarity window 相似窗口半径,下面实际相似窗口大小为3
     %  h: degree of filtering 滤波参数,控制滤波程度
    
     % 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');  %padarray扩展图像的函数,扩展m+2*f行,n+2*f列
     
     % 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)%第一次循环是r=3,s=3;满足if语句,continue;=>第二次循环r=3.s=4=>不满足,到下面W2
                        continue
                    end
                                    
                    W2= input2(r-f:r+f , s-f:s+f);                
                     
                    d = sum(sum(kernel.*(W1-W2).*(W1-W2)));%就是N(i)与N(j)的高斯加权距离
                                                   
                    w=exp(-d/h);%根据加权距离衡量相似度,确定W2窗口与W1窗口的权值w                 
                                     
                    if w>wmax                
                        wmax=w;                   
                    end
                    
                    sweight = sweight + w;
                    average = average + w*input2(r,s); % 将W2的中心点像素值乘以权值w                               
             end 
             end
                 
            average = average + wmax*input2(i1,j1); %最大权值wmax实际就是待处理像素点(i1,j1)所在的图像块,将其累计到average中
            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
      end
    end
    kernel = kernel ./ f;
    
            
    

    调用

    % demo
    clear;
    x0=double(imread('Cameraman256.png'));
    % add noise
    randn('state',0)
    sigma=25; % variance of noise 噪声方差
    x=x0+5*randn(size(x0));
    figure(1), imshow(x,[]),colorbar
    disp(['psnr of noisy image=' num2str(psnr(x,x0)) 'dB'])
    
    %--------parameters-------------
    % delta for convergence control
    f=2;   % radius of square patch in kernelized ar regresssion
    t=10;   %searching range in each direction
    
    tic  %秒表计时器
    [output]=NLmeansfilter(x,t,f,sigma);%调用函数
    toc
    figure(2), imshow(output,[])
    disp(['psnr of denoising output=' num2str(psnr(output,x0)) 'dB'])
    
    
    运行效果图
    添加了均值为0、方差为25的噪声图

    采用相似窗口t=2、搜索窗口f=10、h=25大小非局部均值滤波后效果图

    本人对代码关键点的解释(以本代码参数为例)
    • f=2,t=2,input:256×256大小; input2:260×260(将边缘对称折叠上去,扩展原input再进行相似度计算);output:256×256

    • 便于理解代码可设置断点,观察参数变化,现将第一次执行参数的值以及示意图如下展示:

    • 值得注意的一点是:当i,j确定时循环会获得wmax,也就意味着将最大权值赋予待处理像素点所在的图像块,累计到average中得到最终待处理像素点的像素值。

    参考文献

    [1] Buades A, Coll B, Morel J, et al. A non-local algorithm for image denoising[C]. computer vision and pattern recognition,2005,2(2):60-65.

    展开全文
  • 介绍了经典非局部均值滤波算法与Manjón非局部均值滤波算法,改进了非局部均值滤波方法的相似度权值,使算法在具有旋转平移不变性,保持时间复杂度的同时优化了视觉效果与信噪比。实验通过添加噪声标准差从10~100不等的...
  • 在前面的文章中,我们分别讲过均值滤波、高斯滤波、中值滤波:数字图像处理之均值滤波数字图像处理之高斯滤波数字图像处理之高斯滤波加速优化中值滤波原理及其C++实现与CUDA优化其中,均值滤波的...

    在前面的文章中,我们分别讲过均值滤波、高斯滤波、中值滤波:

    数字图像处理之均值滤波

    数字图像处理之高斯滤波

    数字图像处理之高斯滤波加速优化

    中值滤波原理及其C++实现与CUDA优化

    其中,均值滤波的核心思路是取每一个像素点邻域的矩形窗口,计算矩形窗口内所有像素点的像素平均值,作为该点滤波之后的像素值。高斯滤波与均值滤波类似,都是计算矩形窗口内所有像素点的像素值加权和,只不过其权重与均值滤波不一样,高斯滤波的权重服从二维正态分布,越靠近窗口中心点(也即当前滤波点),权重越大。

    本文我们主要讲非局部均值(NL-means)滤波算法的原理与实现。其核心思路与高斯滤波很相似:计算矩形窗口内所有像素点的像素值加权和,权重服从高斯分布。区别在于:高斯滤波使用当前滤波点与矩形窗口内其它点的空间欧式距离来计算权重,距离越近权重越大;而非局部均值滤波则使用当前滤波点的邻域块与矩形窗口内其它点的邻域块的相似度来计算权重,相似度越大则权重越大

    高斯滤波使用的点与点的空间距离d如下图所示:

    假设A的坐标为(x1,y1),B的坐标为(x2,y2),那么d可按照下式计算:

    非局部均值滤波使用的点与点的邻域块相似度如下图所示:

    要衡量两个邻域块的相似度,有多种指标,均方误差(MSE)是最常用的相似度衡量指标之一。非局部均值滤波算法就是使用MSE来计算两个邻域块的相似度。邻域块的行列相同,假设均为m行n列,那么MSE的计算如下式,其中A(i,j)为点A邻域块中的点(i,j)的像素值,B(i,j)为点B邻域块中的相同位置点的像素值:

    由以上可知,非局部均值滤波算法有两个重要的参数:ksize和ssize。假设点A为当前待滤波点,邻域块分别取以点A和点B为中心的ksize*ksize区域,矩形窗口(也称为搜索窗口)取以点A为中心ssize*ssize大小。由于ksize和ssize都是奇数,通常传入它们的半值:half_ksize和half_ssize。所以,邻域块的大小为(2*half_ksize+1)*(2*half_ksize+1),搜索窗口的大小为(2*half_ssize+1)*(2*half_ssize+1)。如下图所示:

    点A的滤波值由搜索窗口内所有点的像素值加权平均得到

    上式中,w(A,B)为点A与搜索窗口中任意一点B的高斯权重,由两点各自邻域块的MSE相似度计算而得。如下式,w(A,B)需要除以搜索窗口内所有点的高斯系数之和,以实现归一化的目的。其中h也是一个重要的参数,h越大去噪效果越好,但是图像越模糊,反之h越小去噪效果越差,但去噪之后的失真度越小

    由于原图像中每一个待滤波的点都需要取完整的一个搜索窗口和多个邻域块,因此非局部均值滤波之前,首先要对图像进行边缘扩充,上、下、左、右分别扩充half_ssize+half_ksize的行、行、列、列:

    好了,原理就讲到这里,下面上代码。

    1. 查找表的计算

    由于图像的像素值范围是确定的0~255,因此计算邻域块的MSE时可以提前计算好查找表,以节省一点时间。结合以下两个查找表可以快速计算两个数的差值平方。

    //计算0~255的平方查找表
    float table1[256];
    static void cal_lookup_table1(void)
    {
      for (int i = 0; i < 256; i++)
      {
        table1[i] = (float)(i*i);
      }
    }
    
    
    //计算两个0~255的数的绝对差值的查找表
    uchar table2[256][256];
    static void cal_lookup_table2(void)
    {
      for (int i = 0; i < 256; i++)
      {
        for (int j = i; j < 256; j++)
        {
          table2[i][j] = abs(i - j);
          table2[j][i] = table2[i][j];
        }
      }
    }
    

    2. 计算两个邻域块的MSE代码

    float MSE_block(Mat m1, Mat m2)
    {
      float sum = 0.0;
      for (int j = 0; j < m1.rows; j++)
      {
        uchar *data1 = m1.ptr<uchar>(j);
        uchar *data2 = m2.ptr<uchar>(j);
        for (int i = 0; i < m1.cols; i++)
        {
          sum += table1[table2[data1[i]][data2[i]]];
        }
      }
      sum = sum / (m1.rows*m2.cols);
      return sum;
    }
    

    3. 主体函数

    //h越大越平滑
    //halfKernelSize小框
    //halfSearchSize大框
    void NL_mean(Mat src, Mat &dst, double h, int halfKernelSize, int halfSearchSize)
    {
      Mat boardSrc;
      dst.create(src.rows, src.cols, CV_8UC1);
      int boardSize = halfKernelSize + halfSearchSize;
      copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT);   //边界扩展
      double h2 = h*h;
    
    
      int rows = src.rows;
      int cols = src.cols;
    
    
      cal_lookup_table1();
      cal_lookup_table2();
    
    
      for (int j = boardSize; j < boardSize + rows; j++)
      {
        uchar *dst_p = dst.ptr<uchar>(j - boardSize);
        for (int i = boardSize; i < boardSize + cols; i++)
        {
          Mat patchA = boardSrc(Range(j - halfKernelSize, j + halfKernelSize), Range(i - halfKernelSize, i + halfKernelSize));
          double w = 0;
          double p = 0;
          double sumw = 0;
    
    
          for (int sr = -halfSearchSize; sr <= halfSearchSize; sr++)   //在搜索框内滑动
          {
            uchar *boardSrc_p = boardSrc.ptr<uchar>(j + sr);
            for (int sc = -halfSearchSize; sc <= halfSearchSize; sc++)
            {
              Mat patchB = boardSrc(Range(j + sr - halfKernelSize, j + sr + halfKernelSize), Range(i + sc - halfKernelSize, i + sc + halfKernelSize));
              float d2 = MSE_block(patchA, patchB);
    
    
              w = exp(-d2 / h2);
              p += boardSrc_p[i + sc] * w;
              sumw += w;
            }
          }
    
    
          dst_p[i - boardSize] = saturate_cast<uchar>(p / sumw);
    
    
        }
      }
    
    
    }
    

    4. 测试代码

    分别使用上述实现的代码,以及Opencv实现的均值滤波、高斯滤波、中值滤波函数,对有很大噪声的Lene图像进行滤波。

    void nlmean_test(void)
    {
      Mat img = imread("lena_noise.jpg", CV_LOAD_IMAGE_GRAYSCALE);
      
      Mat out;
      nlmean3(img, out, 20, 3, 15);   //NL-means
      
      Mat out1;
      blur(img, out1, Size(11, 11));   //均值滤波
    
    
      Mat out2;
      GaussianBlur(img, out2, Size(11, 11), 2.5);   //高斯滤波
    
    
      Mat out3;
      medianBlur(img, out3, 9);   //中值滤波
    
    
      imshow("img", img);
      imshow("非局部", out);
      imshow("均值", out1);
      imshow("高斯", out2);
      imshow("中值", out3);
      waitKey();
    }
    

    5. 滤波结果对比

    噪声图

    均值滤波结果

    高斯滤波结果

    中值滤波结果

    非局部均值滤波结果

    由上图可知,非局部均值滤波的效果在几种滤波算法中时最好的,但是非局部均值滤波的耗时也是最大的,Lena图像的大小为496*472,这个尺寸并不算大,但是却耗时约20秒。耗时太长很影响其实际应用,因此在接下来的篇章中我们将继续讲解该算法的加速优化,敬请期待~

    本人微信公众号如下,不定期更新精彩内容,欢迎扫码关注:

    展开全文
  • 非局部均值滤波(一)

    千次阅读 2020-10-07 17:54:12
    但是相似像素并不局限于某个局部区域,自然图像含有丰富的冗余信息,所以可以采用能够描述图像结构特征的图像块在整幅图像上寻求相似块。在去噪的同时可以保留能最大程度地保持图像地细节特征。其基本思想是:当前...

    算法原理简介
    图像之中的像素点并不是孤立存在的,某一点的像素与别处的像素点存在某种关联,可以概括为灰度相关性和几何结构相似性。但是相似像素并不局限于某个局部区域,自然图像含有丰富的冗余信息,所以可以采用能够描述图像结构特征的图像块在整幅图像上寻求相似块。在去噪的同时可以保留能最大程度地保持图像地细节特征。其基本思想是:当前像素的估计值由图像中与它具有相似邻域结构的像素加权平均得到。
    以下图为例,进行详细说明:

    以上为例,假定原图像大小为7x7(图中的橙色部分),由于需要对边缘像素进行处理,所以需要将图像进行填充,这种填充采用的是边缘像素的镜像。像素x是需要处理的原图像中的第一个像素,这里把搜索区域设置为7x7,也就是说像素x的估计值由图中绿色部分的49个像素加权平均得到,而权重是通过计算邻域块之间的相似度决定(图中以x和y为中心的蓝色图像块,块大小设置为3x3),如果计算出的块与块相似度越高,那么该像素对所求像素的贡献越大。
    参考代码:

    #include <iostream>
    #include <opencv2/opencv.hpp>
    #include <cmath>
    #include <cstdlib>
    #include <ctime>
    #include <algorithm>
    #include <vector>
    using namespace cv;
    using namespace std;
    
    
    int Rows;//图像的行数
    int Cols;//图像的列数
    
    //添加高斯噪声
    void addNoise(Mat& src, Mat& dest, double sigma)
    {
    	Mat s;
    	src.convertTo(s, CV_16S);
    	Mat n(s.size(), CV_16S);
    	//randn()将n填充为高斯分布均值为0,标准差为sigma的随机数
    	randn(n, 0, sigma);
    	Mat temp = s + n;//将原图与高斯噪声叠加在一起输出
    	temp.convertTo(dest, CV_8U);
    }
    
    
    //非局部均值滤波
    void NlMeans(Mat &src, Mat &dst, double h = 10,int templateWindowSize = 7, int searchWindowSize = 21)
    {
    	//增加边界宽度
    	int r1 = searchWindowSize/ 2;
    	int r2 = templateWindowSize/ 2;
    	int r3 = -1 * r2;
    	int r = r1 + r2;
    	Mat srcCopy;
    	copyMakeBorder(src, srcCopy, r, r, r, r, BORDER_DEFAULT);
    	//使用一个矩阵用来存储图像
    	int hh = r * 2 + Rows;
    	double **matrix = new double *[hh];
    	for (int k = 0; k < hh; k++)
    	{
    		matrix[k] = new double[hh];
    	}
    	for (int x = 0; x < hh; x++)
    	{
    		for (int y = 0; y < hh; y++)
    		{
    			matrix[x][y] = srcCopy.at<uchar>(x, y);
    		}
    	}
    	//srcCopy.convertTo(srcCopy, CV_64F);
    	/*cout << srcCopy.rows << "  " << srcCopy.cols << endl;
    	imshow("image", srcCopy);
    	waitKey(0);*/
    	//非局部均值滤波
    	int i, j,m,n,p,q; //用于循环时像素的索引
    	double similarity; //图像块的相似度
    	double w;  //参考块相当于当前块的权重
    	double sum; //权重之和
    	double s;   //权重乘以像素值之和
    	double sub; //像素差值
    	//vector<double>weight(searchWindowSize*searchWindowSize);
    	//double **weight = new double *[searchWindowSize];//建立一个weight数组,用来存储权重
    	//for (int k = 0; k < searchWindowSize; k++)
    	//{
    	//	weight[k] = new double[searchWindowSize];
    	//}
    	double f = templateWindowSize * templateWindowSize;
    	double delta = 10; //高斯噪声的方差
    	double delta2 = 2 * delta*delta;
    	Mat dstImage = src.clone();
    	//遍历图像的每一个像素点
    	for (i = r; i < Rows+r; ++i)
    	{
    		for (j = r; j < Cols+r; ++j)
    		{
    			//遍历搜索域
    			sum = 0;
    			s = 0;
    			for (m = i - r1; m <= i + r1; ++m)
    			{
    				for (n = j - r1; n <= j + r1; ++n)
    				{
    					//计算每一个参考块
    					similarity = 0;
    					for (p = r3; p <= r2; ++p)
    					{
    						for (q = r3; q <= r2; ++q)
    						{
    							//计算当前块与搜索域中参考块之间的相似度
    							//similarity += pow((srcCopy.at<uchar>(m + p, n + q) - srcCopy.at<uchar>(i + p, j + q)),2)/f;
    							sub=matrix[m + p][n + q] - matrix[i + p][j + q];
    							similarity += sub*sub;
    						}
    					}
    					//计算参考块相对于当前块的权重
    					similarity=similarity/f;
    					w = exp(max((sbimilarity -delta2), 0.0) / (h*h));
    					w = 1 / w;
    					//cout << w << endl;
    					sum += w;
    					//s += w * srcCopy.at<uchar>(m, n);
    					s += w *matrix[m][n];
    				}
    			}
    			dst.at<uchar>(i - r, j - r) = s / sum;
    		}
    	}
    	//释放二维数组
    	for (i = 0; i < hh; i++)
    	{
    		delete[]matrix[i];
    	}
    	delete matrix;
    }
    
    int main()
    {
    	Mat srcImg = imread("E:/vsfiles/lena1.jpg");
    	Mat grayImg;
    	cvtColor(srcImg, grayImg, CV_RGB2GRAY);
    	resize(grayImg, grayImg, Size(256,256));
    	Rows = grayImg.rows;//行数
    	Cols = grayImg.cols;//列数
    	Mat imgWithNoise;
    	addNoise(grayImg,imgWithNoise,10);
    	Mat dstImg = imgWithNoise.clone();
    	NlMeans(imgWithNoise, dstImg);
    	//fastNlMeansDenoising(imgWithNoise, dstImg,10);
    	imshow("image1", grayImg);
    	imshow("image2", imgWithNoise);
    	imshow("image3", dstImg);
    	waitKey(0);
    	return 0;
    }
    
    展开全文
  • 局部均值滤波DEMO

    2018-06-22 14:59:55
    图像保边滤波器集锦之局部均值滤波算法博文所对应的DEMO,具体算法和代码实现请参看博文!
  • 该项目是为了研究基于深度卷积神经网络的图像去噪算法,是利用DnCNN模型,但是为了比较该算法的效果,另外实现了四种传统的图像去噪算法(均值滤波、中值滤波、非局部均值滤波NLM和三维块匹配滤波BM3D)作为对照组。...
  • 非局部均值滤波代码

    2020-02-29 06:00:28
    这是一个matlab程序:非局部均值滤波 这是一个matlab程序:非局部均值滤波 这是一个matlab程序:非局部均值滤波 这是一个matlab程序:非局部均值滤波
  • 改进的多视PolSAR非局部均值滤波算法
  • 在上一篇文章中,我们讲解了非局部均值滤波算法的原理,以及使用C++和Opencv来实现了该算法:非局部均值滤波(NL-means)算法的原理与C++实现我们知道,非局部均值滤波是非常耗时的...

    在上一篇文章中,我们讲解了非局部均值滤波算法的原理,以及使用C++和Opencv来实现了该算法:

    非局部均值滤波(NL-means)算法的原理与C++实现

    我们知道,非局部均值滤波是非常耗时的,这很影响该算法在实际场景中的应用。所以后来有研究人员提出使用积分图来加速该算法,可提升数倍的速度。本文我们将详细讲解该算法的积分图加速原理,并使用C++与Opencv来将其实现。

    积分图的原理我们之前也讲过,此处不再重复:

    数字图像处理之积分图

    数字图像处理之积分图计算优化

    积分图的一种CUDA并行运算

    首先,让我们来回顾一下上篇文章中的计算思路:

    原图像中的每一个点都是待滤波点,每一个待滤波点对应一个以它为中心的搜索窗口和一个同样以它为中心的邻域块(搜索窗口与邻域块都是矩形区域,只不过搜索窗口的尺寸要大于邻域块)。从而每一个待滤波点的滤波值由搜索窗口内所有点的像素值加权求和得到。搜索窗口中每个点的权重则由该点的邻域块与待滤波点的邻域块的MSE相似度决定。

    因此,上篇文章使用4层for循环来实现算法,伪代码如下。由此可知,对于当前的一个待滤波点,计算完其对应搜索窗口内所有点的加权和得到其滤波值之后,再开始计算下一个待滤波点的滤波值。

    for 原图像的每一行
    {
      for 原图像的每一列    //在这一层循环确定了一个待滤波点 
      {
        for 搜索窗口的每一行
        {
          for 搜索窗口的每一列   //在这一层循环确定了当前搜索窗口中的一个点
          {
              //在这里计算当前搜索窗口中确定的一个点的权重,并计算加权和
          }
        }
      }
    }
    

    为了便于使用积分图加速计算,我们把计算顺序改变一下,伪代码如下:

    for 搜索窗口的每一行
    {
      for 搜索窗口的每一列  //在这一层循环确定了所有搜索窗口中相同偏移位置的点 
      {
        for 原图像的每一行
        {
          for 原图像的每一列   //在这一层循环确定了原图像中的一个待滤波点
          {
              /*
              在这里计算当前确定的待滤波点的对应搜索窗口内的对应偏移
              位置点的权重,然后计算该点像素值与权重的乘积,并累加乘积
              */
          }
        }
      }
    }
    

    上述伪代码中讲的“所有搜索窗口中相同偏移位置的点”可能有点抽象,我们结合下图来更加详细地说明。如下图所示,点A1、点A2都是待滤波点,点B1为点A1对应搜索窗口内的一个点,点B2为点A2对应搜索窗口内的一个点,如果B1与A1的相对位置,与B2与A2的相对位置一样,那么B1与B2就是这两个搜索窗口中相同偏移位置的点。

    如果从数学角度来描述相对位置一样,那么就是下式:

    按照改变后的计算顺序,在内两层循环中计算所有待滤波点在各自的搜索窗口中相同偏移位置的点的权重,及该点像素值与权重的乘积,并累加乘积。感觉这样表达还是有点抽象,下面我们再举例子说明。

    比如,待滤波点为A1、A2、A3、A4、A5......

    A1、A2、A3、A4、A5......对应搜索窗口中所有相同偏移位置点为:

    相同偏移位置1:B1、B2、B3、B4、B5......

    相同偏移位置2:C1、C2、C3、C4、C5......

    相同偏移位置3:D1、D2、D3、D4、D5......

    .

    .

    .

    其中B1、C1、D1......组成A1对应的搜索窗口内的所有点。B2、C2、D2......组成A2对应的搜索窗口内的所有点。B3、C3、D3......组成A3对应的搜索窗口内的所有点......

    在上篇文章中我们先计算B1、C1、D1......的权重:

    WB1、WC1、WD1......

    以及像素值与权重的乘积:

    IB1*WB1、IC1*WC1、ID1*WD1......

    并将乘积与权重分别累加,得到点A1的滤波值:

    sumw1 = WB1+WB2+WB3+......

    A1的滤波值=(IB1*WB1+IC1*WC1+ID1*WD1+......)/sumw1

    改变顺序之后,我们先计算B1、B2、B3、B4、B5......的权重得到WB1、WB2、WB3、WB4、WB5......,以及像素值与权重的乘积IB1*WB1、IB2*WB2、IB3*WB3、IB4*WB4、IB5*WB5......,然后将乘积和权重分别累加到一个变量:

    sum1 = sum1+IB1*WB1

    w_sum1 = w_sum1+WB1

    sum2 = sum2+IB2*WB2

    w_sum2 = w_sum2+WB2

    sum3 = sum3+IB3*WB3

    w_sum3 = w_sum3+WB3

    sum4 = sum4+IB4*WB4

    w_sum4 = w_sum4+WB4

    sum5 = sum5+IB5*WB5

    w_sum5 = w_sum5+WB5

    .

    .
           .

    接着计算C1、C2、C3、C4、C5......的权重得到WC1、WC2、WC3、WC4、WC5......,以及像素值与权重的乘积IC1*WC1、IC2*WC2、IC3*WC3、IC4*WC4、IC5*WC5......,然后将乘积和权重分别累加到一个变量:

    sum1 = sum1+IC1*WC1

    w_sum1 = w_sum1+WC1

    sum2 = sum2+IC2*WC2

    w_sum2 = w_sum2+WC2

    sum3 = sum3+IC3*WC3

    w_sum3 = w_sum3+WC3

    sum4 = sum4+IC4*WC4

    w_sum4 = w_sum4+WC4

    sum5 = sum5+IC5*WC5

    w_sum5 = w_sum5+WC5

    .

    .
           .

    同理计算D1、D2、D3、D4、D5......,以及所有相同偏移位置的系列点。最后计算完所有偏移位置的点之后,再计算:

    A1的滤波值=sum1/w_sum1 

    A2的滤波值=sum2/w_sum2

    A3的滤波值=sum3/w_sum3

    A4的滤波值=sum4/w_sum4

    A5的滤波值=sum5/w_sum5

    .

    .
           .

    非局部均值滤波的主要耗时点在于频繁计算邻域块的MSE相似度。讲完计算顺序,接下来我们讲怎么使用积分图来加速计算邻域块的MSE相似度。由积分图的原理可知,计算矩形窗口内所有点的像素值和时,可以使用积分图来加速计算。根据两个邻域块的MSE计算公式,我们可以构造一个图像,并计算该图像的积分图,然后使用该积分图来快速计算邻域块的MSE

    首先在边缘扩充之后的图像中取图像A,相当于上、下、左、右各裁剪half_ssize行、half_ssize行、half_ssize列、half_ssize列得到(也相当于把上图中黄色虚线框内区域抠出来作为图像A)。然后根据搜索窗口内当前计算点与当前待滤波点的位置偏移(x_offset,y_offset)(-half_ssizex_offsethalf_ssize,-half_ssizey_offsethalf_ssize),获取与图像A的行、列均相同的图像B,比如偏移量为(-half_ssize,-half_ssize),那么图像B为上图中红色虚线框内区域抠出来的图像,又比如偏移量为(half_ssize,half_ssize),那么图像B为上图中黑色虚线框内区域抠出来的图像。

    得到图像A与图像B之后,计算两图像的差值图C(也即相同坐标点的像素差值):

    C = A - B

    接着,计算图C的平方图(也即对每一个坐标点的像素值计算平方):

    D = C * C

    然后,计算图D的积分图,即为位置偏移(x_offset,y_offset)的积分图。那么在不同搜索窗口内,所有相对于待滤波点的偏移位置为(x_offset,y_offset)的点,其与待滤波点的邻域块相似度均可以使用这个积分图来计算。一个搜索窗口内有(2*half_ssize+1)*(2*half_ssize+1)个点,也即有(2*half_ssize+1)*(2*half_ssize+1)个相对于其中心点的偏移位置,因此我们需要计算(2*half_ssize+1)*(2*half_ssize+1)张积分图(所有搜索窗口中相同偏移位置的点是共用一张积分图的)。

    最后,得到积分图之后,就可以根据当前待滤波点在积分图上对应的位置取得其邻域块区域,快速计算邻域块区域内所有点的像素值平均值,即为MSE相似度:

    假设待滤波点在原始图像中坐标为A(x,y),那么该点在积分图中坐标为A(x+half_ksize,y+half_ksize),从而在积分图中得到以点A为中心的邻域块,其4个顶点坐标分别为P1(x,y),P2(x+half_ksize,y-half_ksize),P3(x-half_ksize,y+half_ksize),P4(x+half_ksize,y+half_ksize)

    邻域块的尺寸是确定的:(2*half_ksize+1)*(2*half_ksize+1),最后可以根据积分图上点P1、P2、P3、P4的像素值IP1、IP2、IP3、IP4以及邻域块尺寸,计算得到点A、与点A坐标偏移为(x_offset,y_offset)的搜索窗口内点的邻域块相似度:

    讲完原理之后,下面上代码。

    1. 计算积分图代码

    void integralImgSqDiff(Mat src, Mat &dst, int Ds, int t1, int t2, int m1, int n1)
    {
      //计算图像A与图像B的差值图C
      Mat Dist2 = src(Range(Ds, src.rows - Ds), Range(Ds, src.cols - Ds)) - src(Range(Ds + t1, src.rows - Ds + t1), Range(Ds + t2, src.cols - Ds + t2));
      float *Dist2_data;
      for (int i = 0; i < m1; i++)
      {
        Dist2_data = Dist2.ptr<float>(i);
        for (int j = 0; j < n1; j++)
        {
          Dist2_data[j] *= Dist2_data[j];  //计算图像C的平方图D
        }
      }
      integral(Dist2, dst, CV_32F); //计算图像D的积分图
    }
    

    2. 快速NL-means算法主体函数代码

    void fastNLmeans(Mat src, Mat &dst, int ds, int Ds, float h)
    {
      Mat src_tmp;
      src.convertTo(src_tmp, CV_32F);
      int m = src_tmp.rows;
      int n = src_tmp.cols;
      int boardSize = Ds + ds + 1;
      Mat src_board;
      copyMakeBorder(src_tmp, src_board, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT);
      
      Mat average(m, n, CV_32FC1, 0.0);
      Mat sweight(m, n, CV_32FC1, 0.0);
    
    
      float h2 = h*h;
      int d2 = (2 * ds + 1)*(2 * ds + 1);
    
    
      int m1 = src_board.rows - 2 * Ds;   //行
      int n1 = src_board.cols - 2 * Ds;   //列
      Mat St(m1, n1, CV_32FC1, 0.0);
    
    
      for (int t1 = -Ds; t1 <= Ds; t1++)
      {
        int Dst1 = Ds + t1;
        for (int t2 = -Ds; t2 <= Ds; t2++)
        {
          int Dst2 = Ds + t2;
          integralImgSqDiff(src_board, St, Ds, t1, t2, m1, n1);
    
    
          for (int i = 0; i < m; i++)
          {
            float *sweight_p = sweight.ptr<float>(i);
            float *average_p = average.ptr<float>(i);
            float *v_p = src_board.ptr<float>(i + Ds + t1 + ds);
            int i1 = i + ds + 1;   //row
            float *St_p1 = St.ptr<float>(i1 + ds);
            float *St_p2 = St.ptr<float>(i1 - ds - 1);
    
    
            for (int j = 0; j < n; j++)
            {
    
    
              int j1 = j + ds + 1;   //col
              float Dist2 = (St_p1[j1 + ds] + St_p2[j1 - ds - 1]) - (St_p1[j1 - ds - 1] + St_p2[j1 + ds]);
    
    
              Dist2 /= (-d2*h2);
              float w = fast_exp(Dist2);
              sweight_p[j] += w;
              average_p[j] += w * v_p[j + Ds + t2 + ds];
            }
          }
    
    
        }
      }
    
    
      average = average / sweight;
      average.convertTo(dst, CV_8U);
    }
    

    运行上述代码同样对496*472的Lena图像去噪,结果如下图所示,耗时由原来的20秒左右减少为1.6秒左右,所以加速效果还是非常显著的。

    噪声图像

    快速NL-means算法去噪图像

    使用积分图加速之后,计算耗时减少了好多,不过还是秒级的。然而在实时应用场合中通常要求毫秒级的耗时,因此加速得还不够,下篇文章中我们将介绍在积分图加速得基础上,使用CUDA进一步加速优化该算法,敬请期待~

    本人微信公众号如下,不定期更新精彩内容,欢迎扫码关注:

    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 7,106
精华内容 2,842
关键字:

非局部均值滤波