实验三 数字图像的空间域滤波

实验三、四数字图像的空间域滤波和频域滤波

1．实验目的

1.掌握图像滤波的基本定义及目的。

2.理解空间域滤波的基本原理及方法。

3.掌握进行图像的空域滤波的方法。

4.掌握傅立叶变换及逆变换的基本原理方法。

5.理解频域滤波的基本原理及方法。

6.掌握进行图像的频域滤波的方法。

2．实验基本原理

1.空间域增强

空间域滤波是在图像空间中借助模板对图像进行领域操作，处理图像每一个像素的取值都是根据模板对输入像素相应领域内的像素值进行计算得到的。空域滤波基本上是让图像在频域空间内某个范围的分量受到抑制，同时保证其他分量不变，从而改变输出图像的频率分布，达到增强图像的目的。

空域滤波一般分为线性滤波和非线性滤波两类。线性滤波器的设计常基于对傅立叶变换的分析，非线性空域滤波器则一般直接对领域进行操作。各种空域滤波器根据功能主要分为平滑滤波器和锐化滤波器。平滑可用低通来实现，平滑的目的可分为两类：一类是模糊，目的是在提取较大的目标前去除太小的细节或将目标内的小肩端连接起来；另一类是消除噪声。锐化可用高通滤波来实现，锐化的目的是为了增强被模糊的细节。结合这两种分类方法，可将空间滤波增强分为四类：

线性平滑滤波器（低通）

非线性平滑滤波器（低通）

线性锐化滤波器（高通）

非线性锐化滤波器（高通）

空间滤波器都是基于模板卷积，其主要工作步骤是：

1)将模板在图中移动，并将模板中心与图中某个像素位置重合；

2)将模板上的系数与模板下对应的像素相乘；

3)将所有乘积相加；

4)将和（模板的输出响应）赋给图中对应模板中心位置的像素。

2.平滑滤波器

1)线性平滑滤波器

线性低通平滑滤波器也称为均值滤波器，这种滤波器的所有系数都是正数，对3×3

的模板来说，最简单的是取所有系数为1，为了保持输出图像任然在原来图像的灰度值范围内，模板与象素邻域的乘积都要除以9。

MATLAB 提供了fspecial 函数生成滤波时所用的模板，并提供filter2 函数用指定的滤波器模板对图像进行运算。函数fspecial 的语法格式为：

h=fspecial(type);

h=fspecial(type,parameters);

其中参数type 指定滤波器的种类，parameters 是与滤波器种类有关的具体参数。

表2.1 MATLAB中预定义的滤波器种类

MATLAB 提供了一个函数imnoise 来给图像增添噪声，其语法格式为：

J=imnoise(I,type);

J=imnoise(I,type,parameters);

参数type 指定噪声的种类，parameters 是与噪声种类有关的具体参数。参数的种类见表2.2。

表2.2 噪声种类及参数说明

2)非线性平滑滤波器

中值滤波器是一种常用的非线性平滑滤波器，其滤波原理与均值滤波器方法类似，但计算的非加权求和，而是把领域中的图像的象素按灰度级进行排序，然后选择改组的中间值作为输出象素值。

MATLAB 提供了medfilt2 函数来实现中值滤波，其语法格式为：

B=medfilt2(A,[m n]);

B=medfilt2(A);

其中，A 是原图象，B 是中值滤波后输出的图像。[m n]指定滤波模板的大小，默认模板为3×3。

3.锐化滤波器

图像平滑往往使图像中的边界、轮廓变得模糊，为了减少这类不利效果的影响，需要利用图像锐化技术，使图像的边缘变得清晰。

1)线性锐化滤波器

线性高通滤波器是最常用的线性锐化滤波器。这种滤波器的中心系数都是正的，而周围的系数都是负的，所有的系数之和为0。

对3×3 的模板来说，典型的系数取值为：

[-1 -1 -1;

-1 8 -1;

-1 -1 -1]

事实上这是拉普拉斯算子。语句h=-fspecial(‘laplacian’,0.5)得到的拉普拉斯算子为：

h =-0.3333 -0.3333 -0.3333

-0.3333 2.6667 -0.3333

-0.3333 -0.3333 -0.3333

2)非线性锐化滤波

邻域平均可以模糊图像，因为平均对应积分，所以利用微分可以锐化图像。图像处理中最常用的微分方法是利用梯度。常用的空域非线性锐化滤波微分算子有sobel 算子、prewitt 算子、log 算子等。

4.频域增强

频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其他空间中，然后利用该空间的特有性质方便地进行图像处理，最后再转换回原来的图像空间中，从而得到处理后的图像。

频域增强的主要步骤是：

选择变换方法，将输入图像变换到频域空间。

在频域空间中，根据处理目的设计一个转移函数，并进行处理。

将所得结果用反变换得到增强的图像。

常用的频域增强方法有低通滤波和高通滤波。

5.低通滤波

图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分。因此能降低高频成分幅度的滤波器就能减弱噪声的影响。由卷积定理，在频域实现低通滤波的数学表达式：

G(u,v) =H(u,v)F(u,v)

1)理想低通滤波器（ILPF）

2)巴特沃斯低通滤波器（BLPF）

3)指数型低通滤波器（ELPF）

6.高通滤波

由于图像中的细节部分与其高频分量相对应，所以高通滤波可以对图像进行锐化处理。高通滤波与低通滤波相反，它是高频分量顺利通过，使低频分量受到削弱。高通滤波器和低通滤波器相似，其转移函数分别为:

1)理想高通滤波器（IHPF）

2)巴特沃斯高通滤波器（BLPF）

3)指数型高通滤波器（ELPF）

图像经过高通滤波处理后，会丢失许多低频信息，所以图像的平滑区基本上会消失。所以，可以采用高频加强滤波来弥补。高频加强滤波就是在设计滤波传递函数时，加上一

个大于0小于1的常数c，即：

H′(u,v) =H(u,v)+c

3．实验内容与要求

1.平滑空间滤波：

1)读出eight.tif这幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示在同一图像窗口中。

2)对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果，要求在同一窗口中显示。

3)使用函数imfilter时，分别采用不同的填充方法（或边界选项，如零填

充、’replicate’、’symmetric’、’circular’）进行低通滤波，显示处理后的图像。

4)运用for循环，将加有椒盐噪声的图像进行10次，20次均值滤波，查看其特点,显示均值处理后的图像（提示:利用fspecial函数的’average’类型生成均值滤波器）。

5) 对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处 理，要求在同一窗口中显示结果。

6) 自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。 2. 锐化空间滤波

1) 读出blurry_moon.tif 这幅图像，采用3×3的拉普拉斯算子w = [ 1, 1, 1; 1 – 8 1; 1, 1, 1]对其进行滤波。

2) 编写函数w = genlaplacian(n)，自动产生任一奇数尺寸n 的拉普拉斯算子，如5×5 的拉普拉斯算子

w = [ 1 1 1 1 1 1 1 1 1 1 1 1 -24 1 1 1 1 1 1 1 1 1 1 1 1]

3) 分别采用5×5，9×9，15×15和25×25大小的拉普拉斯算子对blurry_moon.tif 进

行锐化滤波，并利用式2(,)(,)(,)g x y f x y f x y =-?完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。

4) 采用不同的梯度算子对blurry_moon.tif 进行锐化滤波，并比较其效果。

5) 自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像；

3. 傅立叶变换

1) 读出woman.tif 这幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位 图像。

仅对相位部分进行傅立叶反变换后查看结果图像。

2) 仅对幅度部分进行傅立叶反变换后查看结果图像。 3) 将图像的傅立叶变换F 置为其共轭后进行反变换，比较新生成图像与原始图像的差 异。

4. 平滑频域滤波

1) 设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截至频率自选，分 别给出各种滤波器的透视图。

2) 读出test_pattern.tif 这幅图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高 斯低通滤波器对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同低通滤波器得到的图像与原图像的区别，特别注意振铃效应。(提示:1)在频率域滤波同样要注意到填充问题；2）注意到(-1)x+y ;)

5. 锐化频域滤波

1) 设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截至频率自选，分 别给出各种滤波器的透视图。

2) 读出test_pattern.tif 这幅图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高 斯高通滤波器对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。

1. 实验具体实现

1.平滑空间滤波：

1)读出eight.tif这幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图

显示在同一图像窗口中。

I=imread('eight.tif');

imshow(I);

J = imnoise(I,'salt & pepper',0.05); %noise density=0.05

K= imnoise(I,'gaussian',0.01,0.01);

subplot(131),imshow(I);

subplot(132),imshow(J);

subplot(133),imshow(K);

图2.1 初始图像及椒盐噪声图像、高斯噪声污染图

2)对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果，

要求在同一窗口中显示。

I=imread('eight.tif');

H = fspecial('sobel');

Sobel = imfilter(I,H,'replicate');

H = fspecial('laplacian',0.4);

lap = imfilter(I,H,'replicate');

H = fspecial('gaussian',[3 3],0.5);

gaussian = imfilter(I,H,'replicate');

subplot(221),imshow(I);

subplot(222),imshow(Sobel);

subplot(223),imshow(lap);

subplot(224),imshow(gaussian);

图2.2 原图像及各类低通滤波处理图像

3)使用函数imfilter时，分别采用不同的填充方法（或边界选项，如零填

充、’replicate’、’symmetric’、’circular’）进行低通滤波，显示处理后的图像。originalRGB = imread('peppers.png');

h = fspecial('motion', 50, 45); %motion blurred

filteredRGB = imfilter(originalRGB, h);

boundaryReplicateRGB = imfilter(originalRGB, h, 'replicate');

boundary0RGB = imfilter(originalRGB, h, 'x');

boundary0RGB = imfilter(originalRGB, h, 0);

boundarysymmetricRGB = imfilter(originalRGB, h, 'symmetric'); boundarycircularRGB = imfilter(originalRGB, h, 'circular');

图2.3 原图像及运动模糊图像

图2.4 函数imfilter各填充方式处理图像

4)运用for循环，将加有椒盐噪声的图像进行10次，20次均值滤波，查看其特点,显

示均值处理后的图像。

J = imnoise(I,'salt & pepper',0.05);

h=fspecial('average'); %Averaging Filtering

J1=imfilter(J,h);

for i=1:10

J1=imfilter(J,h);

end

for i=1:20

J2=imfilter(J,h);

图2.5 椒盐噪声污染图像经10次、20次均值滤波图像

由图2.5可得，20次滤波后的效果明显好于10次滤波，但模糊程度也更强。

5)对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处

理，要求在同一窗口中显示结果。

h=fspecial('average'); %Averaging Filtering

J1=imfilter(J,h);

J2=medfilt2(J); %Median Filtering

图2.6 椒盐噪声污染图像及均值、中值滤波图像

从图2.6中可以看出，对于椒盐噪声污染的图像处理，中值滤波效果要明显好于均值滤波。经均值滤波器处理后的图像比均值滤波器中结果图像更加模糊。

6)设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。domain=[0 0 8 0 0;

0 0 8 0 0;

8 8 8 8 8;

0 0 8 0 0;

0 0 8 0 0];

K1= ordfilt2(J,5,domain);

图2.7 椒盐噪声污染图像及5*5平滑滤波器掩模

掩模值为w=1/25*[1 1 1 1 1;1 1 1 1 1;1 1 1 1 1;1 1 1 1 1;1 1 1 1 1]

图2.8 椒盐噪声污染图像及5*5平滑滤波器掩模

掩模值为w= [0 0 8 0 0;0 0 8 0 0;8 8 8 8 8; 0 0 8 0 0;0 0 8 0 0]

2.锐化空间滤波

1)采用3×3的拉普拉斯算子w = [ 1, 1, 1; 1 – 8 1; 1, 1, 1]滤波

I=imread('blurry_moon.tif');

T=double(I);

subplot(1,2,1),imshow(T,[]);title('Original Image');

w =[1,1,1;

1,-8,1;

1,1,1];

K=conv2(T,w,'same');

图2.9 初始图像与拉普拉斯算子锐化图像

2)编写函数w = genlaplacian(n)，自动产生任一奇数尺寸n的拉普拉斯算子，如5×5的拉

普拉斯算子：

w = [ 1 1 1 1 1

1 1 1 1 1

1 1 -24 1 1

1 1 1 1 1

1 1 1 1 1]

function w = genlaplacian(n)

%Computes the Laplacian operator

w = ones(n);

x = ceil(n/2);

w(x, x) = -1 * (n * n - 1);

3)分别采用5×5，9×9，15×15和25×25大小的拉普拉斯算子对blurry_moon.tif进行锐

化滤波，并利用式

2

(,)(,)(,)

g x y f x y f x y

=-?完成图像的锐化增强，观察其有何不

同，要求在同一窗口中显示。

不同尺寸拉普拉斯算子滤波以及图像增强

w1 = genlaplacian(5);

I=imread('blurry_moon.tif');

T=double(I);

K=conv2(T,w1,'same');

J=T-K;

图2.10 初始图像与不同拉普拉斯算子锐化图像

图像锐化的实质是将原图像与梯度信息叠加，相当于对目标物的边缘进行了增强。

图2.11 拉普拉斯算子锐化与锐化增强图像

4)采用不同的梯度算子对blurry_moon.tif进行锐化滤波，并比较其效果[I,map]=imread('blurry_moon.tif');

I=double(I);

[Gx,Gy]=gradient(I); % gradient calculation

G=sqrt(Gx.*Gx+Gy.*Gy); % matrix

J1=G; % gradient1

imshow(J1,map);

J2=I; % gradient2

K=find(G>=7);

J2(K)=G(K);

imshow(J2,map);

J3=I; % gradient3

K=find(G>=7);

J3(K)=255;

imshow(J3,map);

J4=I; % gradient4

K=find(G<=7);

J4(K)=255;

imshow(J4,map);

J5=I; % gradient5

K=find(G<=7);

J5(K)=0;

Q=find(G>=7);

J5(Q)=255;

imshow(J5,map);

图2.12 原始图像与不同梯度子锐化图像

作为二阶微分算子，拉普拉斯变换在图像细节的增强处理上有明显的优点，但会产生更多的噪声。梯度变换在灰度变化区域的响应更强，但对噪声和细节的响应比拉普拉斯变换弱。

5)自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像；

图2.13 原始图像与不同边缘锐化图像

3.傅立叶变换

1)读出woman.tif这幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。domain=[8 8 0 8 8;

8 8 0 8 8;

0 0 0 0 0;

8 8 0 8 8;

8 8 0 8 8];

K1= ordfilt2(J,5,domain);

F=imread('woman.tif');

F1=fft2(F);

F2=log(1+abs(F1)); %amplitude spectrum

F3=fftshift(F1);

imshow(log(1+abs(F3)),[]);

F4=angle(F1); %phase spectrum

图2.14 原始图像与快速傅立叶变换图像

2)仅对相位部分进行傅立叶反变换后查看结果图像。

F1=fft2(F);

i=sqrt(-1);

f2=ifft2(exp(i*angle(F1)));

imshow(real(f2),[]);

图2.15 原始图像与对全部信息进行傅立叶逆变换结果

图2.16 仅对相位信息进行傅立叶逆变换结果

相位谱决定了图像信号中各频率分量的位置。

3)仅对幅度部分进行傅立叶反变换后查看结果图像。

f1=ifft2(abs(F1));

imshow(log(1+abs(f1)),[]);

图2.17 仅对幅度信息进行傅立叶逆变换结果

4)将图像的傅立叶变换F置为其共轭后进行反变换，比较新生成图像与原始图像的差异。F1=fft2(F);

F2=log(1+abs(F1)); %amplitude spectrum

F3=fftshift(F1);

F4=angle(F1); %phase spectrum

F5=-F4

F6= double(F3*exp(F4)); %the complex conjugate of the fourier transform F7=ifft2(F6); %inverse fourier transform

imshow(real(F7),[]);

图2.18 共轭傅立叶逆变换结果

傅立叶变换的相位谱为对称的，原变换与其共轭变换间仅频率谱互为相反。

4.平滑频域滤波

1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器

理想低通滤波器

I=imread('test_pattern.tif');

f=double(I); % chage into double as MATLAB doesn’t suppor calculation % of image in unsigned int type

g=fft2(f); % fourier transform

g=fftshift(g); % zero-frequency area centralized

[M,N]=size(g);

d0=100; %cutoff frequency

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

if(d<=d0)

h=1;

else h=0;

end

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

imshow(J2)

巴特沃斯低通滤波器（二阶）

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

nn=2; % 2-grade Butterworth lowpss filter

d0=100;

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

h=1/(1+0.414*(d/d0)^(2*nn)); % filter transform function

%h=1./(1+(d./d0).^(2*n))

%h=exp(-(d.^2)./(2*(d0^2)));

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

imshow(J2);

高斯低通滤波

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

d0=100;

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

h=exp(-(d.^2)./(2*(d0^2))); % gaussian filter transform

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

图2.19 理想低通滤波器透视图

图2.20 巴特沃斯低通滤波器透视图

图2.21高斯低通滤波器透视图

由各类低通滤波器透视图可见，高斯滤波器剖面线最平滑，二阶巴特沃斯低通滤波器函数剖面线较为紧凑，而理想滤波器完全为圆筒状结构，未考虑选择范围内不同信息频率的差

别化处理。

图2.22 理想低通滤波器滤波效果(d0=15,30,100)

当截止频率do = 15时，滤波后的图像模糊，难以分辨，振铃现象明显。

当do = 30时，滤波后的图像模糊减弱，能分辨出字母与图形轮廓，但由于理想低通滤波器在频率域的锐截止特性，滤波后的图像仍有较明显的振铃现象。

当do = 100时，滤波后的图像比较清晰，但高频分量损失后，图像边沿与文字变的有些模糊，在图像的边框（如条带和矩形轮廓）附近仍有轻微振铃现象。

图2.23 巴特沃斯低通滤波器滤波效果(d0=15,30,100)

图2.23中显示了3种二阶巴特沃斯低通滤波器的滤波效果，各截止频率同图2.22。二

阶的巴特沃斯低通滤波器显示了轻微的振铃和较小的负值，但远不如理想滤波器明显。

一阶巴特沃斯滤波器无振铃现象，在二阶中振铃通常很微小。阶数越高振铃现象越明显，一个20阶的巴特沃斯低通滤波器已经呈现出理想低通滤波器的特性。

图2.24 高斯低通滤波器滤波效果(d0=15,30,100)

图2.24中显示了3种高斯低通滤波器的滤波效果，各截止频率同图2.22。高斯低通滤波器无法达到有相同截止频率的二阶巴特沃斯低通滤波器的平滑效果，但此时结果图像中无振铃现象产生。

5.锐化频域滤波

设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器

理想高通滤波器

图2.25 理想高通滤波器透视图

巴特沃斯高通滤波器