MATLAB中实现图像的空间域滤波和频率域滤波

一）空间域滤波与频率域滤波

 1）空间域滤波

      空间域滤波是指在图像空间中借助模板对图像领域进行操作，处理图像每一个像素值。主要分为线性滤波和非线性滤波两类，根据功能可分为平滑滤波器和锐化滤波器。平滑可通过低通来实现，平滑的目的有两类，一是模糊，目的是在提取较大的目标前去除太小的细节或将目标内的小尖端连接起来；二是去噪。锐化则可用高通滤波来实现，锐化的目的是为了增强被模糊的细节。

      在matlab中实现空间域滤波，有很多类型，如均值、中值、索贝尔、高斯、拉普拉斯、高斯-拉普拉斯等，但各有差异。下面是用matlab实现的代码：

[html]  view plain  copy

1. %空间域滤波  
2. clc;close all;  
3. I=imread('1.tif');  
4. w1=fspecial('average',[3 3]);  
5. w2=fspecial('sobel');  
6. w3=fspecial('gaussian',[3 3],0.5);  
7. w4=fspecial('laplacian',0.2);  
8. w5=fspecial('log',[5 5],0.5);  
9. g1=imfilter(I,w1,'replicate');  
10. g2=imfilter(I,w2,'replicate');  
11. g3=imfilter(I,w3,'replicate');  
12. g4=imfilter(I,w4,'replicate');  
13. g5=imfilter(I,w5,'replicate');  
14. g6=medfilt2(I);  
15. subplot(3,3,1);imshow(I);title('原图');  
16. subplot(3,3,2);imshow(g1);title('均值滤波');  
17. subplot(3,3,3);imshow(g2);title('索贝尔滤波');  
18. subplot(3,3,4);imshow(g3);title('高斯滤波');  
19. subplot(3,3,5);imshow(g4);title('拉普拉斯滤波');  
20. subplot(3,3,6);imshow(g5);title('高斯-拉普拉斯滤波');  
21. subplot(3,3,7);imshow(g6);title('中值滤波');  

其运行结果如下：

    空间域滤波均采用matlab自带的函数，参数也采用默认的值，可见中值滤波的效果最好。

    其中各滤波器的原理和优缺点如下：

   A)均值滤波：由fspecial函数生成的w1是一个大小为3*3的矩形平均滤波器，再用imfilter这个函数使这个掩模的中心逐个滑过图像的每个像素，输出为模板限定的相应领域像素与滤波器系数乘积结果的累加和。由处理结果可见均值滤波器的效果使每个点的像素都平均到它的领域去了，噪声明显减少了很多，效果较好。

   B)索贝尔滤波：w2是一个大小为3*3的sobel滤波器sv，用来近似计算垂直梯度，在图像的任何一点使用此算子，将会产生对应的梯度矢量或是其法矢量。但是Sobel导数并不是真正的导数，这是因为Sobel算子定义于一个离散空间之上，它真正表示的是多项式拟合，用较大的核的话会在更多像素上进行拟合，会更加正确。而较小的核对噪声会更加敏感，此时用sobel算子近似计算导数的缺点精度比较低，这种不精确性在试图估计图像的方向导数 (使用y/x滤波器响应的反正切得到的图像梯度的方向)。比如对于3*3的Sobel滤波器，梯度角度接近水平或者垂直方向的时候， 这样的不准确性会比较明显。由滤波效果可见到图像的边缘凸显了出来，sobel算子主要用于边缘检测。

   C）高斯滤波：高斯滤波器是平滑线性滤波器的一种，线性滤波器很适合于去除高斯噪声。而非线性滤波则很适合用于去除脉冲噪声，中值滤波就是非线性滤波的一种。高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波器是带有权重的平均值，即加权平均，中心的权重比邻近像素的权重更大，这样就可以克服边界效应。高斯滤波如果采用3×3掩模的具体公式如下:

g(x,y)={f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)+f(x+1,y+1)+[f(x-1,y)+f(x,y-1)+f(x+1,y)+f(x,y+1)]*2+f(x,y)*4}/16

其中，f(x,y)为原图像中（x,y)像素点的灰度值，g(x,y)为经过高斯滤波和的值。由处理效果可看出高斯滤波的减噪能力较好。

 

   D）拉普拉斯滤波：拉普拉斯算子是n维欧式空间的一个二阶微分算子。拉普拉斯算子会突出像素值快速变化的区域，因此常用于边缘检测。由效果可见图像的边界得到了增强。

   E）中值滤波：中值滤波法是一种非线性平滑技术，它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值。中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。方法是用某种结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为g（x,y）=med{f(x-k,y-l),(k,l∈W)} ，其中，f(x,y)，g(x,y)分别为原始图像和处理后图像。W为二维模板，通常为3*3，5*5区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。中值滤波对于斑点噪声和椒盐噪声来说尤其有用。保存边缘的特性使它在不希望出现边缘模糊的场合也很有用。由上图效果可见中值滤波的效果最好。

     由以上分析可知，各种滤波器各有优劣，适用情况也不尽相同，线性滤波器很适合于去除高斯噪声，而非线性滤波则很适合用于去除脉冲噪声，如中值滤波很适合去除椒盐噪声。使用起来要视具体实际情况而定。

2）频率域滤波

  

     开始时看书看了些原理想自己用算法实现发现有困难，就去百度了一下，写了下面的一些代码：

    

[html]  view plain  copy

1. %频率域滤波  
2. clc;close all;  
3. f=imread('1.tif');  
4. f=im2double(f);  
5.   
6. F=fft2(double(f));%傅里叶变换  
7.   
8. F=fftshift(F);%将变换的原点移到频率矩形的中心  
9. [M,N]=size(f);  
10.   
11. %理想低通滤波  
12.   
13. D0=input('输入截止频率');  
14. h1=zeros(M,N);  
15. for i=1:M  
16.     for j=i:N  
17.         if(sqrt(((i-M/2)^2+(j-N/2)^2))<D0)  
18.             h1(i,j)=1;  
19.           
20.         end  
21.     end  
22. end  
23. G1=F.*h1;  
24. G1=ifftshift(G1);  
25. g1=real(ifft2(G1));  
26.   
27. %巴特沃斯低通滤波  
28. n=input('巴特沃斯滤波器的阶数 n=');    
29. n1=fix(M/2);  
30. n2=fix(N/2);  
31. h2=zeros(M,N);  
32. for i=1:M  
33.     for j=1:N  
34.         d=sqrt((i-n1)^2+(j-n2)^2);  
35.         h2=1./(1+(d./D0).^(2*n));  
36.          
37.     end  
38. end  
39. G2=F.*h2;  
40. G2=ifftshift(G2);  
41. g2=real(ifft(G2));  
42.   
43. %高斯低通滤波  
44. h3=zeros(M,N);  
45. for i=1:M  
46.     for j=1:N  
47.         h3=exp(-(d.^2)./(2*(D0^2)));  
48.     end  
49. end  
50. G3=F.*h3;  
51. G3=ifftshift(G3);  
52. g3=real(ifft(G3));  
53.   
54. subplot(2,3,1);imshow(f);title('原图');  
55. subplot(2,3,2);imshow(g1);title('理想低通滤波');  
56. subplot(2,3,3);imshow(g2);title('巴特沃斯低通滤波');  
57. subplot(2,3,4);imshow(g3);title('高斯低通滤波');  

其运行效果如下：

 

   可见出现了问题，理想低通滤波的效果并不好，巴特沃斯和高斯滤波器的输出图像成了黑色，后经反复检查，调试，并没找出问题所在，于是打算尝试另一种方法。代码和结果如下：

[html]  view plain  copy

1. %巴特沃斯低通滤波  
2. clc;close all;  
3. f=imread('1.tif');  
4. f=im2double(f);  
5. M=2*size(I,1);  
6. N=2*size(I,2);                        %滤波器的行列数    
7. u=-M/2:(M/2-1);  
8. v=-N/2:(N/2-1);  
9. [U,V]=meshgrid(u,v);  
10. D=sqrt(U.^2+V.^2);  
11. D0=50;                                %截止频率  
12. n=6;  
13. H=1./(1+(D./D0).^(2*n));              %设计巴特沃斯滤波器  
14. F=fftshift(fft2(I,size(H,1),size(H,2)));%傅里叶变换  
15. G=F.*H;  
16. L=ifft2(fftshift(G));                   %傅里叶反变换  
17. L=L(1:size(I,1),1:size(I,2));  
18. subplot(121);imshow(f);  
19. subplot(122);imshow(L);  

[html]  view plain  copy

1. %高斯低通滤波  
2. clc;close all;  
3. I=imread('1.tif');    
4. I=im2double(I);    
5. M=2*size(I,1);  
6. N=2*size(I,2);                        %滤波器的行列数    
7. u=-M/2:(M/2-1);  
8. v=-N/2:(N/2-1);  
9. [U,V]=meshgrid(u,v);  
10. D=sqrt(U.^2+V.^2);  
11. D0=20;  
12. H=exp(-(D.^2)./(2*(D0^2)));          %设计高斯滤波器  
13. J=fftshift(fft2(I,size(H,1),size(H,2)));  
14. G=J.*H;  
15. L=ifft2(fftshift(G));  
16. L=L(1:size(I,1),1:size(I,2));  
17. figure;  
18. subplot(121);imshow(I);  
19. subplot(122);imshow(L);   

    通过这种方法运行出的结果是正常的。低通滤波滤掉了图像频谱中的高频成分，仅让低频部分通过，即变化剧烈的成分减少了，结果是使图像便模糊。比较可见高斯滤波的效果最好。

   频率域高通滤波和低通滤波实现的原理差不多，只是在设计滤波器时的公式有些差异，我就没有再重复实现了。