高斯滤波/高斯平滑/高斯模糊的实现及其快速算法(Gaussian Filter, Gaussian Smooth, Gaussian Blur, Fast implementation)

最新推荐文章于 2024-08-10 20:02:03 发布

无敌三角猫

最新推荐文章于 2024-08-10 20:02:03 发布

阅读量6.3k

点赞数 2

文章标签：高斯模糊

高斯滤波/高斯平滑/高斯模糊的实现及其快速算法(Gaussian Filter, Gaussian Smooth, Gaussian Blur, Fast implementation)

标签： filter 算法 fft matlab image

2012-03-07 14:11 6009人阅读评论(0) 收藏举报

 版权声明：本文为博主原创文章，未经博主允许不得转载。

网上介绍针对图像进行高斯模糊的文章不少，其原理比较简单，这里就不做过多介绍。这里简单总结一下实现高斯模糊的几种算法(假设图像大小是M*N，filter的半径是r，注意，很多文章使用的术语是filter size，指的是半径大小或者直径大小，一般情况下filter的半径r取3倍或者4倍的sigma)：

1.最原始的实现方法是，使用高斯函数生成高斯模板，然后使用该模板对图像做卷积(convolution)。该方法的算法复杂度是O(M*N*r^2)，也就是，对每个像素而言，复杂度是O(r^2)，与滤波器大小成quadratic关系。如果增大滤波器size，算法会明显减慢。

下面是Matlab实现的代码（注意，为了代码方便，以下代码中出现的filter size均指滤波器直径，而不是半径!）：

[plain]view plaincopy 
    
 function [T] = GenGaussFilterKernel(size, sigma)   
 % here the size is the diameter of the filter  
   
 if mod(size, 2) == 0  
     size = size + 1;  
 end  
   
 T = zeros(size, size);  
 distMat = zeros(size, size);  
 centerIdx = ceil(size / 2);  
   
 for i = 1:size  
     for j = 1:size  
         distMat(i,j) = sqrt(abs(i - centerIdx) ^ 2 + abs(j - centerIdx) ^ 2);  
     end  
 end  
   
 T = exp(-distMat.^2 / (2 * sigma^2));  

[plain]view plaincopy 
    
 function [img_ext] = PaddingImgByMirrorEdge(img, r)  
   
 [h,w] = size(img);  
   
 % make a larger image, mirror the edges  
 img_ext = zeros(h + 2*r, w + 2*r);  
 img_ext(r+1:h+r, r+1:w+r) = img;  
 img_ext(1+r:h+r, 1:r) = flipdim(img(1:h, 2:r+1), 2); %add left border  
 img_ext(1+r:h+r, w+r+1:w+r+r) = flipdim(img(1:h, w-r:w-1), 2); %add right border  
 img_ext(1:r, 1:w+r+r) = flipdim(img_ext(r+2:r+r+1, 1:w+r+r), 1); %add up border  
 img_ext(h+r+1:h+r+r, 1:w+r+r) = flipdim(img_ext(h+r-r:h+r-1, 1:w+r+r), 1); %add bottom border  

[plain]view plaincopy 
    
 function [I] = ImageConvolution(img, filterKernel)  
   
 [h,w] = size(img);  
 filterSize = size(filterKernel);  
 r = floor(filterSize(1) / 2);  
   
 I = zeros(h, w);  
   
 % make a larger image, mirror the border to extend it  
 imgExt = PaddingImgByMirrorEdge(img, r);  
   
 % convolution  
 sumFilter = sum(sum(filterKernel));  
 tempRegion = zeros(filterSize);  
 for i = 1 : h  
     for j = 1 : w  
         tempRegion = imgExt(i:i+2*r, j:j+2*r);  
         I(i,j) = sum(sum(double(tempRegion) .* filterKernel)) / sumFilter;  
     end  
 end  
   
 I = uint8(I);  

2.高斯滤波器的kernel是可分离的(separable)，也就是说，可以将2D的高斯kernel分解为两个1D的kernel，先沿x方向对图像进行1D高斯kernel的卷积，然后沿y方向对图像进行1D的高斯kernel卷积，最后的结果和使用一个2D高斯kernel对图像卷积效果是一样的。这样一来，针对每个像素，滤波器的算法复杂度降为O(r)。

3.使用FFT-based convolution。因为图像的convolution等价于FFT转换后频域的乘积，所以可以将高斯kernel扩展至与图像同等大小，然后分别对图像和高斯kernel进行FFT变换，然后直接将两者的FFT结果图像按照per pixel的方式相乘，最终结果再逆FFT变换回去。该算法总体复杂度为O(M*N*log(M*N))，针对每个像素，滤波器的复杂度为O(log(M*N))。可以看到，当滤波器size比较小时，该算法不如使用可分离的1D高斯卷积，而当滤波器size比较大(r >log(M*N))时，FFT-based卷积速度会比较快。使用Matlab实现基于FFT的高斯模糊，这里有一篇讲解非常详细的文章。

Matlab代码如下（注意，为了代码方便，以下代码中出现的filter size均指滤波器直径，而不是半径!）：

[plain]view plaincopy 
    
 function [I] = GaussianSmooth(img, filterSize, sigma, is_use_fft)  
 % here the size is the diameter of the filter  
   
 if nargin < 4  
     is_use_fft = 1;  
 end  
   
 T = GenGaussFilterKernel(filterSize, sigma); %generate gaussian kernel  
 [h, w] = size(img);  
 [d1, d2] = size(T);  
 r = floor(d1 / 2);  
 r1 = floor(h / 2);  
 r2 = floor(w / 2);  
   
 if is_use_fft == 0  
     % use convolution to implement  
     I = ImageConvolution(img, T);  
 else  
     % use fft to implement Gaussian smooth  
     extImg = PaddingImgByMirrorEdge(img, r);  
     extT = zeros(size(extImg)); %extend filter kernel to the same size  
     extT(1:end-h+1, 1:end-w+1) = T; %just place the original kernel on the topleft  
       
     fftImg = fft2(extImg);  
     fftT = fft2(extT);  
     resFFT = fftImg .* fftT; %multiply the FFT result  
     resI = ifft2(resFFT); %inverse FFT  
           
     % calculate coefficients of Gaussian smooth  
     extCoef = ones(size(extImg));  
     fftCoef = fft2(extCoef);  
     resCoef = fftCoef .* fftT;  
     resC = ifft2(resCoef);  
       
     % output image  
     I = uint8(resI ./ resC);  
     I = I(2*r+1:2*r+h, 2*r+1:2*r+w);   
 end  

4.递归高斯滤波(Recursive Gaussian Filter)算法。比上面的算法更进一步，有人提出一些approximate高斯滤波的算法，可以使算法速度不依赖于滤波器size，针对每个像素而言，滤波器的算法复杂度是O(1)级别的。这方面的算法原理介绍比较复杂，要理解信号处理相关的知识才能理解(如z变换等，参见如下参考文献)。所幸的是，算法本身其实相当简单，大概思路是正向扫描一遍并累积像素值（累积时使用特定的系数可以approximate高斯kernel），然后再反向扫描一遍并累积像素值，done！而且算法也是separable的，可以x和y方向分别做，所以速度可以达到非常之快！有相关的代码在此，NVIDIA也有CUDA代码在此，最新的GPU快速算法见这里。原理可以参考这里和这里。以下是几篇这方面的参考文献：

[1] Rachid Deriche - "Recursively implementing the Gaussian and its derivatives", 1993.
[2] Lucas J. van Vliet, Ian T. Young and Piet W. Verbeek - "Recursive Gaussian derivative Filters", 1998
[3] Dave Hale, "Recursive Gaussian Filters", CWP-546
[4] Luis Alvarez, Rachid Deriche, Francisco Santana - "Recursivity and PDE's in Image Processing", 2000