空域高斯滤波与频域高斯滤波

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

cfqcfqcfqcfqcfq

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

阅读量3.5k

点赞数 3

分类专栏：图像处理

图像处理专栏收录该内容

9 篇文章 1 订阅

订阅专栏

卷积定理

函数空间域的卷积的傅里叶变换是函数傅里叶变换的乘积。对应地，频率域的卷积与空间域的乘积存在对应关系。

即：

$f(x,y)*g(x,y)\Leftrightarrow F(u,v)H(u,v)$

$f(x,y)g(x,y)\Leftrightarrow F(u,v)*H(u,v)$

由卷积定理可知所有频域的滤波理论上都可以转化为空域的卷积操作。

给定频率域滤波器，可对其进行傅里叶逆变换得到对应的空域滤波器；滤波在频域更为直观，但空域适合使用更

小的滤波模板以提高滤波速度。因为相同尺寸下，频域滤波器效率高于空域滤波器，故空域滤波需要一个更小尺

寸的模板近似得到需要的滤波结果。

空域卷积

将模板在图像中逐像素移动，将卷积核的每个元素分别和图像矩阵对应位置元素相乘并将结果累加，累加和作为

模板中心对应像素点的卷积结果。通俗的讲，卷积就是对整幅图像进行加权平均的过程，每一个像素点的值，都

由其本身和邻域内的其他像素值经过加权平均后得到。在像素的处理上，是先将结果暂存在于一个副本，最后统

一拷贝，故不会出现处理顺序不同而结果不同的情况。

二维连续卷积的数学定义：

$f(x,y)*g(x,y)=\int_{-\infty }^{+\infty}\int_{-\infty}^{+\infty}f(u,v)g(x-u,y-v)dudv$

离散形式：

$f(s,t)*g(s,t)=\sum_{m}\sum_{n}f(m,n)g(s-m,t-n)$

频域滤波

频率域是由傅里叶变换和频率变量 (u,v)定义的空间，频域滤波处理过程：先对图像进行傅里叶变换，转换至频率

域，在频域使用滤波函数进行滤波，最后将结果反变换至空间域。即：

$f(x,y)\overset{FFT}{\Rightarrow}F(u,v)\overset{H(u,v)}{\Rightarrow}G(u,v)\overset{IFFT}{\Rightarrow}g(x,y)$

高斯函数

公式：

$h(x,y)=\frac{1}{2\pi\sigma ^{2}}exp(-\frac{x^{2}+y^{2}}{2\sigma^{2}})$

形状：

高斯函数的特殊性：高斯函数傅里叶变换仍是高斯函数，但标准差已经变化，频域标准差越大（高斯函数越

宽），变换后空域标准差越小（高斯函数越窄）。

如频域函数 $H(u,v)=Aexp(-u^{2}/2\sigma^{2})$ 对应的空域函数为 $H(u,v)=\sqrt{2\pi} \sigma A exp(-2\pi^{2}\sigma^{2}x^{2})$

空域高斯平滑滤波

高斯模板的生成

因为图像是离散存储的，故我们需要一个高斯函数的离散近似。具体地，对高斯函数进行离散化，以离散点上的

高斯函数值作为权值，组成一定尺寸的模板，用此模板对图像进行卷积。由于高斯分布在任意点处都非零，故理
论上需要一个无穷大的模板，但根据" $3\sigma$ 准则"，即数据分布在 $(\mu-3\sigma,\mu+3\sigma]$ 的概率是0.9974，距离函数
中心超过 $3\sigma$ 数据所占权重可以忽略,因此只需要计算 $(6\sigma+1)\times(6\sigma+1)$ 的矩阵就可以保证对高斯函数的
近似了。
假设二维模板大小 $m\times n$ ,则模板上元素处的值为：
$G(x,y)=\frac{1}{2\pi\sigma^{2}}exp(-\frac{(x-m/2)^{2}+(y-n/2)^{2}}{2\sigma^{2}})$

前面的系数在实际应用中常被忽略，因为是离散取样，不能使取样和为1，最后还要做归一化操作。

程序：

[objc]view plaincopy 
    
 function filt=mygaussian(varargin)  
 %参数初始化，使用varargin处理可变参数情况  
 siz=varargin{1};%模板尺寸  
 if(numel(siz)==1)  
     siz=[siz,siz];  
 end  
 std=varargin{2};%方差  
 centa = (siz(1)+1)/2;%此处不要取整  
 centb = (siz(1)+1)/2;  
   
 filt = zeros(siz(1),siz(2));  
 summ=0;  
 for i=1:siz(1)  
     for j=1:siz(2)  
         radius = ((i-centa)^2+(j-centb)^2);  
         filt(i,j) = exp(-(radius/(2*std^2)));  
         summ=summ+filt(i,j);  
     end  
 end  
 filt=filt/summ;%归一化  

测试：

执行mygaussian(4,1)得：

0.0181 0.0492 0.0492 0.0181

0.0492 0.1336 0.1336 0.0492

0.0181 0.0492 0.0492 0.0181

执行fspecial('gaussian',4,1)得：

0.0181 0.0492 0.0492 0.0181

0.0492 0.1336 0.1336 0.0492

0.0181 0.0492 0.0492 0.0181

可以看出与Matlab结果相同。

查看fspecial的Matlab源码，将gaussian部分提取出来：

[objc]view plaincopy 
    
 function h = hisfspecial(varargin)  
 siz=varargin{1};  
 if(numel(siz)==1)  
     siz=[siz,siz];  
 end  
 std   = varargin{2};  
 siz   = (siz-1)/2;  
 [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));  
 arg   = -(x.*x + y.*y)/(2*std*std);  
 h     = exp(arg);  
 h(h<eps*max(h(:))) = 0;  
 sumh = sum(h(:));  
 if sumh ~= 0,  
     h  = h/sumh;  
 end;  

与mygaussian类似，只是更具Matlab特色。

模板与图像卷积

直接利用Matlab提供的fspecial、imfilter。

[objc]view plaincopy 
     
 I=imread('lena.bmp');  
 h=fspecial('gaussian',7,8);  
 O=imfilter(I,h,'replicate','same','conv');  
 subplot(1,2,1);imshow(I,[]);title('input');  
 subplot(1,2,2);imshow(O,[]);title('output');   

具体的实现细节还有很多值得研究的，如边界处理、将二维高斯模板分解为两个一维模板等。可参考文后链接。

结果：

频域高斯低通滤波

在频率域，高斯函数表示为:

$H(u,v)=exp(-\frac{u^{2}+v^{2}}{2\sigma^{2}})$ 或写成 $H(u,v)=exp(-\frac{D(u,v)^{2}}{2D_{0}^{2}})$

其中 $D_{0}$ 为截止频率， $D(u,v)$ 为图像频率空间点(u,v)处的值， $D(u,v)=u^{2}+v^{2}$ ；可以看出距离频率空间的

(0,0)处越远， $D(u,v)$ 越大，也就是频率越高。而对高斯函数 $H(u,v)$ 距离中心点越远值越小，故利用

$G(u,v)=F(u,v)H(u,v)$ 可以滤除高频，保留低频。

程序：

[objc]view plaincopy 
    
 I=imread('lena.bmp');  
 I=double(I);  
 [M,N]=size(I);  
 m=(M+1)/2;  
 n=(N+1)/2;   
 d0=30;                       %截止频率  
 F=fftshift(fft2(I));               %空域转频域，平移中心  
 H=zeros(M,N);  
 G=zeros(M,N);  
 for i=1:M  
    for j=1:N  
        d=(i-m)^2+(j-n)^2;  
        H(i,j)=exp(-d/(2*d0^2));  
        G(i,j)=H(i,j)*F(i,j);  
     end  
 end  
 O=ifft2(ifftshift(G));                %平移中心，频域转空域  
 O=real(O);                        %取实数部分  
   
 subplot(1,2,1);imshow(I,[]);title('inpit');  
 subplot(1,2,2);imshow(O,[]);title('output');  