第5章〓频率域图像处理 频率域图像处理是指以傅里叶变换为基础,利用傅里叶变换把空间域的数字图像变换到频率域,通过对频率域图像的变换系数进行某种修正和处理,再利用傅里叶反变换把修正或处理后的频率域图像的变换系数变换回空间域的一种图像处理方法。频率域图像是以频谱(也称为频谱图)的形式表示图像信息分布特征的一种表示方式。 本章首先介绍面向空间域图像的二维离散傅里叶变换及其频谱特性,接着介绍频率域图像处理的基本实现思路,然后介绍基于频率域的图像噪声消除(频率域低通滤波)和基于频率域的图像增强(频率域高通滤波)方法,最后介绍面向某些特殊应用的带阻滤波和带通滤波。 5.1二维离散傅里叶变换 1822年,法国工程师傅里叶(Fourier)指出: “任意”的周期函数f(x)都可以分解为无穷多个不同频率的正弦和或余弦和,即傅里叶级数。求解傅里叶级数的过程就是傅里叶变换。傅里叶级数和傅里叶变换统称为傅里叶分析或谐波分析。 离散傅里叶变换(Discrete Fourier Transform,DFT)描述的是离散信号的一维时域或二维空间域表示与频域表示之间的关系,是信号处理和图像处理中的一种最有效的数学工具之一,在频谱分析、数字滤波器设计、功率谱分析、传递函数建模、图像处理等方面具有广泛的应用。 5.1.1二维离散傅里叶变换的定义和傅里叶频谱 由于一幅数字图像可以描述成一个二维函数,所以下面仅介绍应用于图像处理的二维离散傅里叶变换。 1. 二维离散傅里叶变换的定义 设f(x,y)是在空间域上等间隔采样得到的M×N的二维离散信号,x和y是离散实变量,u和v为离散频率变量,则二维离散傅里叶变换对一般定义为 F(u,v)=1MN∑M-1x=0 ∑N-1y=0f(x,y)exp-j2πxuM+yvN u=0,1,…,M-1; v=0,1,…,N-1(5.1) f(x,y)=1MN∑M-1u=0 ∑N-1v=0F(u,v)expj2πxuM+yvN x=0,1,…,M-1; y=0,1,…,N-1(5.2) 其中,式(5.1)为傅里叶正变换,exp(-j2πxu/M-j2πyv/N)为正变换核; 式(5.2)为傅里叶反(逆)变换,exp(j2πxu/M+j2πyv/N)为反变换核。 二维离散傅里叶变换的频谱和相位角定义为 F(u,v)=R2(u,v)+I2(u,v)(5.3) (u,v)=arctan[I(u,v)/R(u,v)](5.4) 2. 图像的傅里叶频谱特性及频谱图 傅里叶变换的物理意义是将空间域图像的灰度分布函数变换为频率域图像的频率分布函数,从物理效果看,傅里叶变换是将图像从空间域转换到频率域。一幅空间域的图像f(x,y)变换到频率域F(u,v)后,其傅里叶频谱F(u,v)[即,F(u,v)的幅值|F(u,v)|称为频谱]可以通过频谱图的形式显示。例如,图5.1(a)所示的图像A的频谱图如图5.1(b)所示,图5.2(a)所示的图像B的频谱图如图5.2(b)所示。 图5.1图像A的傅里叶频谱示例 图5.2图像B的傅里叶频谱示例 在图像A、B的傅里叶频谱中,原空间域图像上的灰度突变部位、图像结构复杂的区域、图像细节及干扰噪声等信息集中在高频区,对应于图5.1(b)和图5.2(b)的傅里叶频谱的中间部位; 原空间域图像上灰度变化平缓部位(即空间域图像的平坦区域)的信息集中在低频区,对应于图5.1(b)和图5.2(b)的傅里叶频谱的4个角部分。即在傅里叶频谱图上,4个角反映的是原图像的低频特性,中间部位反映的是原图像的高频特性。 按照图像空间域和频率域的对应关系,空间域中的强相关性,即图像中一般都存在有大量的平坦区域,使得图像平坦区域中的相邻或相近像素一般趋向于取值相同或值相近的灰度值,它们在频率域中对应于低频部分,低频部分对应于傅里叶频谱的4个角部分。由于低频部分能量较集中,因而在频谱图上的视觉效果较亮。 3. 二维图像的离散傅里叶变换的实现 理论上,图像傅里叶变换的编程实现方法就是编程实现式(5.1)定义的计算功能,可用如下的MATLAB程序段描述。 clc; clear all; close all; f = imread('d:\0_matlab图像课编程\girl.jpg')% 读取灰度图像 f = im2double(f); [M,N] = sizes(f); for u=0:M-1 for v=0:N-1 sum = 0; for x=0:M-1 e1 = exp(-i*2*pi*u*x/M); for y=0:N-1 e2 = exp(-i*2*pi*v*y/N); sum = sum + f(x+1, y+1)*e1*e2 end end F(u+1,v+1) = real(sum);%取sum的实部 end end F_log=log(1+abs(F)); subplot(1,2,1); imshow(f);title('原灰度图像'); subplot(1,2,2); imshow(F_log,[]);title('原图像的傅里叶频谱图像'); 在上述程序代码中,由于exp(±iθ)=cosθ±isinθ,所以sum结果是复数,而傅里叶频谱图像F(u,v)的值应该是实数,所以需要取sum的实部,即real(sum)。 另外,由于上述程序语句中涉及指数式exp(-i*2*pi*u*x/M)和exp(-i*2*pi*v*y/N)的计算,当M×N图像阵列的M和N较大时计算量非常大,以至于受以前计算装置能力的限制而无法实现(感兴趣的读者可以进一步分析其计算量非常大的机理)。根据傅里叶变换的可分离性(见5.1.2节),一个二维的离散傅里叶变换,可以分别通过在行方向和列方向上的两次一维离散傅里叶变换来实现,所以就出现了各种各样的快速离散傅里叶变换(Fast Fourier Transformation,FFT)算法,比如按时间提取的基2一维FFT算法和按时间提取的基4一维FFT算法等。进一步介绍一维快速傅里叶变换算法已经超出了本书的内容范围,在此不再赘述。下面例子中的二维傅里叶变换用MATLAB自带的傅里叶变换函数实现。 【例5.1】图像傅里叶变换MATLAB程序,程序运行结果如图5.3所示。 图5.3例5.1傅里叶变换MATLAB程序运行结果 clc; clear all; close all; f = imread('d:\0_matlab图像课编程\lena.jpg');%读取灰度图像 F =fft2(f);%二维傅里叶正变换 F_log_img=log(1+abs(F));%对傅里叶频谱图像求对数运算 subplot(1,2,1); imshow(f);title('原灰度图像'); subplot(1,2,2); imshow(F_log_img,[]);title('傅里叶频谱图像'); 例5.1中的相关函数及功能说明如下。 (1) fft2(I): 对灰度图像I进行傅里叶正变换,变换结果的傅里叶频谱图是相对于以(0,0)为原点的显示坐标的。 (2) imshow(I,[]): 通过把图像I的数据范围线性地映射到[min(I(:)) max(I(:))],从而可以将整幅图像的信息都显示出来,一般用于显示图像频谱图。 5.1.2二维离散傅里叶变换的重要性质 二维离散傅里叶变换的性质包括线性、可分离性、平均值性质、周期性、共轭对称性、空间位置和空间频率的平移性、旋转性、尺度变换性、卷积性质等,下面介绍几种比较重要且在本书的相关内容中涉及的性质。 1. 变换系数矩阵 由二维离散傅里叶反变换式(5.2)可知,由于u有0,1,…,M-1的M个可能的取值,v有0,1,…,N-1的N个可能的取值,所以f(x,y)由M×N个频率分量组成,每个频率分量都与一个特定的(u,v)值相对应; 且对于某个特定的(u,v)值来说,当(x,y)取遍所有可能的值(x=0,1,…,M-1; y=0,1,…,N-1)时,就可得到对应于该特定的(u,v)值的一个变换系数矩阵。 exp j2π0u+0vMNexp j2π0u+1vMN…expj2π0u+(N-1)vMN expj2π1u+0vMNexpj2π1u+1vMN…expj2π1u+(N-1)vMN ……… expj2π(M-1)u+0vMNexpj2π(M-1)u+1vMN…expj2π(M-1)u+(N-1)vMN (5.5) 显然,对应于不同(u,v)值的变换系数矩阵共有M×N个,且它们与f(x,y)无关。 2. 可分离性 式(5.1)和式(5.2)的二维离散傅里叶变换对可写成如下的分离形式: F(u,v)=1MN∑M-1x=0∑N-1y=0f(x,y)exp -j2πvyNexp -j2πuxM u=0,1,…,M-1;v=0,1,…,N-1(5.6) f(x,y)=1MN∑M-1u=0∑N-1v=0F(u,v)exp-j2πyvNexp-j2πxuM x=0,1,…,M-1;y=0,1,…,N-1(5.7) 上述的可分离表示形式说明,可以连续运用两次一维DFT来实现一个二维DFT。以式(5.6)为例,可先沿y轴方向进行一维的(列)变换而求得 F(x,v)=1MN∑N-1y=0f(x,y)exp -j2πyvN v=0,1,…,N-1(5.8) 然后再对F(x,v)沿x方向进行一维的(行)变换而得到最后结果 F(u,v)=1MN∑M-1x=0F(x,v)exp -j2πuxM u=0,1,…,M-1;v=0,1,…,N-1 (5.9) 3. 平均值 一幅图像的灰度平均值可表示为 f-=1MN∑M-1x=0 ∑N-1y=0f(x,y)(5.10) 如果将u=v=0代入式(5.2),可得 F(0,0)=1MN∑M-1x=0 ∑N-1y=0f(x,y)(5.11) 所以,一幅图像的灰度平均值可由DFT在原点处的值求得,即 f-=1MNF(0,0)(5.12) 4. 周期性 对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5.1)和式(5.2),F(u,v)的周期性定义为 F(u,v)=F(u±mM,v±nN)m,n∈N(5.13) 5. 共轭对称性 设f(x,y)为实函数,则其傅里叶变换F(u,v)具有共轭对称性 F(u,v)=F(-u,-v)(5.14) |F(u,v)|=|F(-u,-v)|(5.15) 6. 平移性 对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5.1)和式(5.2),若设用符号表示函数与其傅里叶变换的对应性,则傅里叶变换的平移性可表示为 f(x,y)expj2πu0xM+v0yNF(u-u0,v-v0)(5.16) F(u,v)exp-j2πx0uM+y0vNf(x-x0,y-y0)(5.17) 式(5.16)说明,将函数f(x,y)与一个指数项相乘,就相当于把其变换后的傅里叶频谱在频率域进行平移。式(5.17)说明,将傅里叶频谱F(u,v)与一个指数项相乘,就相当于把其反变换后得到的函数在空间域进行平移。 5.1.3图像的傅里叶频谱特性分析 前面介绍的图5.1(b)和图5.2(b)是坐标原点位于(0,0)的傅里叶频谱,而在实际中进行图像的傅里叶频谱特性分析时,使用的是原点位于(M/2,N/2)的傅里叶频谱图。 1. 图像傅里叶频谱关于(M/2,N/2)的对称性 设f(x,y)是一幅大小为M×N的图像,根据离散傅里叶变换的周期性公式(5.13),有 F(u,v)=F(u-M,v-N)(5.18) 再根据式(5.18)和离散傅里叶变换的共轭对称性公式(5.15)可得 |F(u,v)|=|F(M-u,N-v)|(5.19) 下面以图5.4为参照,进一步说明图像的傅里叶变换频谱关于(M/2,N/2)的对称性。 根据式(5.19),对于u=0: 当u=0,v=0时,|F(0,0)|=|F(M,N)|; 当u=0,v=1时,|F(0,1)|=|F(M,N-1)|; 当u=0,v=2时,|F(0,2)|=|F(M,N-2)|; …… 当u=0,v=N/2时,|F(0,N/2)|=|F(M,N/2)|。 同理,对于v=0: 当u=0,v=0时,|F(0,0)|=|F(M,N)|; 当u=1,v=0时,|F(1,0)|=|F(M-1,N)|; 当u=2,v=0时,|F(2,0)|=|F(M-2,N)|; …… 当u=M/2,v=0时,|F(M/2,0)|=|F(M/2,N)|。 图5.4频谱图对称性图示 以上分析表明,如图5.4所示的傅里叶频谱的A区和D区关于点(M/2,N/2)对称。 更一般地: 当u∈[0,M/2]和v∈[0,N/2]时,有(M-u)∈[M/2,M]和(N-u)∈[N/2,N ]。 当u∈[M/2,M]和v∈[0,N/2]时,有(M-u)∈[0,M/2]和(N-u)∈[N/2,N ]。 即傅里叶频谱的A区和D区关于点(M/2,N/2)对称,傅里叶频谱的B区和C区关于点(M/2,N/2)对称。 2. 原点坐标平移到(M/2,N/2)后的傅里叶频谱 由于图像的低频分量比较集中,在频谱图的4个角所占的区域都较小,所以不利于图像的频谱特性分析。根据傅里叶频谱的周期性和平移性,如图5.5(c)所示,当把傅里叶频谱图的原点从(0,0)平移至(M/2,N/2)时,以(M/2,N/2)为原点截取大小为M×N的区间,即可得到一个低频分量位于中心的图像频谱图,如图5.5(d)所示。 图5.5将原点坐标平移到(M/2,N/2)后的傅里叶频谱特征分析示意图 对比图5.4和图5.5(d)可知,基于如图5.4所示对称性的图像的原频谱图,只要把其左上角的A区向下并向右平移到结果频谱图的D区,对应地把其右下角的D区向上并向左平移到结果频谱图的A区,同理把其右上角的B区向下并向左平移到结果频谱图的C区,对应地把其左下角的C区向上并向右平移到结果频谱图的B区,就可以得到如图5.5(d)的关于(M/2,N/2)的图像傅里叶频谱图,即中心对称的频谱图。在实际中,把通过傅里叶正变换得到的原图像的频谱图分成4块平移实现起来比较麻烦,相对比较简单的实现思路如图5.6所示,即先把原频谱图像[如图5.6(a)所示]分成上下各半,接着把由块A和块B组成的上半部分向下平移到由块C和块D组成的下半部分位置,把由块C和块D组成的下半部分向上平移到由原块A和块B组成的上半部分位置,就可得到如图5.6(b)所示的中间结果。然后把如图5.6(b)所示的中间结果分成左右各半,把由块C和块A组成的左半部分向右平移到由块D和块B组成的右半部分位置,把由块D和块B组成的右半部分向左平移到由原块C和块A组成的左半部分位置,就可得到中心对称的图像频谱图。 图5.6把原傅里叶频谱图像平移成中心对称频谱图像过程示意图 按照上述方法对图5.1(b)和图5.2(b)的原频谱图进行中心对称化平移,就可得到如图5.7(a)和图5.7(b)所示的中心对称频谱图像。 图5.7中心对称的频谱图像 显然,图5.7(a)和图5.7(b)的中心清楚地显示出了图像的低频分量的分布情况,不仅具有可视化特点,而且能简化图像的滤波过程等(详见5.3节和5.4节)。 【例5.2】对图像进行傅里叶变换,并按图5.6的原理实现频谱图中心对称的MATLAB程序,程序运行结果如图5.8所示。 clc; clear all; close all; f = imread('d:\0_matlab图像课编程\miroslava1.jpg');%读取灰度图像 [h,w] = size(f); F =fft2(f);%利用MATLAB自带函数实现二维傅里叶正变换 F1 = zeros(h,w); F2 = zeros(h,w);%暂存中间结果频谱图 %实现频谱图中心对称的平移对调 h1 = h/2;w1 = w/2; for i=1:h1 %F的上半部分平移到下半部分 for j=1:w F1(i+h1,j) = F(i,j); end end for i=h1+1:h %F的下半部分平移到上半部分 for j=1:w F1(i-h1, j) = F(i,j); end end for i=1:h %基于前面平移结果F1,再将左半部分平移到右半部分 for j=1:w1 F2(i, j+w1) = F1(i,j); end end for i=1:h %将F1的下半部分平移到上半部分 for j=1+w1:w F2(i, j-w1) = F1(i,j); end end F_log_img=log(1+abs(F)); %对原傅里叶频谱图像求对数 F2_log_img=log(1+abs(F2));%对中心对称傅里叶频谱图像求对数 subplot(1,3,1); imshow(f);title('原灰度图像'); subplot(1,3,2); imshow(F_log_img,[]);title('正变换的傅里叶频谱图像'); subplot(1,3,3); imshow(F2_log_img,[]);title('中心对称傅里叶频谱图像'); 【例5.3】利用MATLAB的自带函数,编写实现傅里叶正变换和反变换的MATLAB程序,程序运行结果如图5.9所示。 clc; clear all; close all; gray_img=imread('d:\0_matlab图像课编程\girl.jpg')% 读取灰度图像 FT_img=fft2(gray_img);% 二维傅里叶正变换 FTS_img=fftshift(FT_img);% 实现频谱图的中心对称 FTS_log_img=log(1+abs(FTS_img));% 频谱图对数运算 Inverse_I=ifftshift(FTS_img); % 反移频谱中心 Inverse_img=real(ifft2(Inverse_I)); % 傅里叶反变换,并取变换结果的实部 Inverse_FT_img=uint8(Inverse_img);% 转换为8位图像 subplot(1,3,1); imshow(gray_img);title('原灰度图像'); subplot(1,3,2); imshow(FTS_log_img,[]);title('中心对称傅里叶频谱图'); subplot(1,3,3); imshow(Inverse_FT_img); title('反变换结果图像'); 图5.8例5.2的傅里叶变换MATLAB程序运行结果 图5.9例5.3傅里叶变换与反变换MATLAB程序运行结果 例5.3中的相关函数及功能说明如下。 (1) fftshift(F): 将低频在4个角的非中心对称傅里叶频谱图,通过对其对角的平移,使低频部分移到中心,实现频谱图中心对称。 (2) ifftshift(F): 将中心对称傅里叶频谱图F分块平移成低频在4个角(非中心对称)的傅里叶频谱。 (3) ifft2(F): 对傅里叶频谱F进行傅里叶反变换。 (4) real(cI): 求傅里叶反变换结果(复数结果)cI的实部值。 3. 对图像进行傅里叶变换的意义 利用傅里叶变换把空间域的图像变换到频率域,即可得到该图像所含频率分量分布情况的傅里叶频谱,进而可以利用频率域的图像处理方法对傅里叶频谱进行处理,其意义如下。 (1) 简化计算。在空间域中处理图像时所进行的复杂的卷积运算,等同于在频率域中简单的乘积运算。 (2) 在用频谱图表示的频率域图像中,中心部位是能量集中的低频特征,反映的是图像的平滑部分。随着不断远离频谱图的中心位置,对应于空间图像中变化越来越快的细节、边缘、结构复杂区域、突变部位和噪声等高频成分逐渐加强。所以,在频率域中滤波的概念更为直观,更容易理解; 即某些在空间域中难以处理或处理起来比较复杂的问题,在频率域却比较容易处理。 (3) 将某些只能在频率域处理的特定应用需求进行转换,比如频率域图像特征提取、数据压缩、纹理分析、水印嵌入等。 5.2频率域图像处理的基本实现思路 由5.1节可知,二维离散傅里叶变换很好地描述了二维离散信号的空间域与频域之间的关系,所以对于那些在空间域中表述起来比较困难,甚至是不太可能实现的图像处理问题,可以先通过对图像进行离散傅里叶变换把图像变换到频率域,然后利用适当的频率域图像处理方式对图像进行处理,处理完后再把它转换回空间域中,从而解决那些在空间域不便于解决的图像处理问题。 5.2.1基本实现思想 由傅里叶频谱的特性可知,u和v同时为0时的频率成分对应于图像的平均灰度级。当从傅里叶频谱的原点离开时,低频对应着图像的缓慢变化分量,比如一幅图像中较平坦的区域; 当进一步离开原点时,较高的频率开始对应图像中变化越来越快的灰度级,它们反映了一幅图像中物体的边缘和灰度级突发改变(如噪声)部分的图像成分。频率域图像增强正是基于这种机理,通过对图像的傅里叶频谱进行低通滤波(使低频通过,使高频衰减)来滤除噪声; 或通过对图像的傅里叶频谱进行高通滤波(使高频通过,使低频衰减)突出图像中的边缘和轮廓来实现图像的增强。 设f(x,y)为M×N的图像,F(u,v)为原点在(0,0)的傅里叶频谱,F′(u,v)为原点平移到(M/2,N/2)的傅里叶频谱,H(u,v)为转移函数(也称为滤波器或滤波器函数),G′(u,v)为对F′(u,v)进行频率域滤波后的输出结果,G(u,v)为原点平移回(0,0)的傅里叶频谱,g(x,y)为经频率域滤波后的输出图像。则对F′(u,v)滤波有 G′(u,v)=F′(u,v)H(u,v)(5.20) 将G′(u,v)平移回原点(0,0)可得G(u,v),对G(u,v)进行傅里叶反变换有 g(x,y)=F-1[G(u,v)](5.21) 在式(5.20)中,H和F′的相乘定义为二维函数逐元素的相乘,即H的第1个元素乘以F′的第1个元素,H的第2个元素乘以F′的第2个元素,以此类推。被滤波的图像可以由G(u,v)的傅里叶反变换结果得到。通常情况下转移函数都为实函数,所以其傅里叶反变换的虚部为0。 综上所述,频率域图像处理的步骤为: (1) 对图像f(x,y)进行二维傅里叶变换,得到原点在(0,0)的傅里叶频谱F(u,v)。 (2) 把原点在(0,0)的傅里叶频谱F(u,v)按照图5.6的原理,即先把F(u,v)分成上下各半,然后把上半部分向下平移到下半部分位置,把下半部分向上平移到上半部分位置; 接着把其中间结果分成左右各半,把左半部分向右平移到右半部分位置,把右半部分向左平移到左半部分位置,就可得到原点在(M/2,N/2)的傅里叶频谱F′(u,v)。 (3) 进行频率域滤波,即用设计的滤波器函数H(u,v)乘以原点在(M/2,N/2)的傅里叶频谱F′(u,v),根据式(5.20)可得G′(u,v)。 (4) 将G′(u,v)的原点平移回(0,0)可得G(u,v)。 (5) 对G(u,v)进行傅里叶反变换,并取变换结果的实部,即计算F-1[G(u,v)]并取计算结果的实部,就可得到通过频率域滤波后的结果图像g(x,y)。 以上过程可简要地描述为图5.10。 图5.10频率域图像增强步骤 在图5.10中,前处理和后处理可能包括将输入图像向最接近的偶数维数转换(以便图像有合适的变换中心)、灰度级标定、输入向浮点数/双精度数的转换、输出向8位整数格式的转换等。 5.2.2转移函数的设计 假设原图像f(x,y)经傅里叶变换为F(u,v),频率域图像处理就是选择合适的滤波器函数H(u,v)对F(u,v)的频谱成分进行调整,然后经傅里叶反变换得到空间域的图像g(x,y)。因此,频率域图像处理的关键是设计合适的滤波转移函数H(u,v)。 关于转移函数的设计方法,一是先凭直观感觉选择一个理想的滤波器模型,然后通过反复的滤波实验和参数修正来逼近并设计出实际的滤波器; 二是利用频率成分和图像外表之间的对应关系选择频率域滤波器; 三是基于数学和统计准则设计频率域滤波器。 另外,对于大小为M×N的函数f(x,y)和h(x,y),有卷积表示形式 f(x,y)*h(x,y)=1MN∑Mm=0 ∑Nn=0f(m,n)h(x-m,y-n)(5.22) 其中,F(u,v)和H(u,v)为f(x,y)和h(x,y)的傅里叶变换。 傅里叶变换对为 f(x,y)*h(x,y)F(u,v)H(u,v)(5.23) F(u,v)*H(u,v)f(x,y)h(x,y)(5.24) 即空间域的卷积在频率域简化为相乘,频率域的卷积在空间域简化为相乘。有时也可以将频率域的滤波器函数变换到空间域,然后再在空间域对图像进行滤波运算。 一般来说,如果两个域中的滤波器尺寸(即滤波窗口大小)相同,通常在频率域中进行滤波计算更为有效,而在空间域中更适合于较小尺寸的滤波器。 5.3基于频率域的图像噪声消除——频率域低通滤波 图像中的噪声和边缘对应于傅里叶频谱的高频部分,选择能使低频通过、使高频衰减的转移函数H(u,v),就可以根据式(5.20)和式(5.21)实现低通滤波,达到滤除噪声的目的。 下面分别介绍理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器三种滤波器的有关概念和滤波原理。 5.3.1理想低通滤波器 一个理想的低通滤波器的转移函数定义为 H(u,v)=1D(u,v)≤D0 0D(u,v)>D0(5.25) 其中,D0是一个非负整数,D(u,v)为频率平面从原点到点(u,v)的距离。 设已经将傅里叶频谱的原点平移到(M/2,N/2),则点(u,v)到频率平面原点(M/2,N/2)的距离为 D(u,v)=[(u-M/2)2+(v-N/2)2]1/2(5.26) 理想低通滤波器的含义为,在半径为D0的圆内的所有频率没有衰减地通过该滤波器,而在此半径的圆外的所有频率完全被衰减掉。所以D0称为截止频率。也就是说,对于低通滤波来说,当某个频率之上的幅度统统被忽略为0时,该频率就称为截止频率。实现上,在进行低通滤波时,先要设定(或经多次实验得到)一个截止频率D0,如D0=π/3。将高于π/3 的频率乘以0,让其全部被阻止; 将低于π/3的频率乘以1,让其全部通过。此时的频率值π/3就称为截止频率。 图5.11(a)给出了理想低通滤波器的转移函数H的横截面图,其中横轴用于表示与原点的径向距离。通过将横截面绕原点旋转360°即可得到完整的理想低通滤波器转移函数,即图5.11(b)所示的转移函数H的透视图。该透视图的含义是: 只有那些位于该圆柱体内的频率范围的信号才能通过,而位于圆柱体外的频率成分都将被滤除掉。 图5.11理想低通滤波器的转移函数横截面图和透视图 在图5.11(a)及本节以下的内容中,均假设D0的频率平面原点为(M/2,N/2)。 需要说明的是,理想低通滤波器的数学意义是十分清楚的,利用计算机对其进行模拟也是可行的,但在实际中用电子元件实现直上直下的理想低通滤波器是不可能的,所以才将其称为“理想”低通滤波器。 【例5.4】频率域理想低通滤波器的滤波效果及低频特性分析。 解: 若设R为截止频率的圆周半径,EB为圆周内能量(图像功率)与原图像总能量(总功率)的百分比,根据图像信号能量在频率域上的分布有 EB=∑u∈R∑v∈RP(u,v)/∑M-1u=0 ∑N-1v=0P(u,v)×100100(5.27) 图5.12(a)是一幅包含了全部细节的原始图像; 图5.12(b)是标注半径D分别为10、20、40和80的4个截止频率圆的傅里叶频谱图; 利用其截止频率确定的理想低通滤波器对图5.12(b)的傅里叶频谱进行低通滤波,将其结果进行傅里叶反变换后得到的图像分别为图5.12(c)~图5.12(f)。根据式(5.27)计算可知,图5.12(c)~图5.12(f)分别包含了原图像中95.5%、97.9%、99.0%和99.6%的能量。 图5.12(c)~图5.12(f)的结果说明: 傅里叶频谱图的低频分量主要集中在中心。 指明了以截止频率为半径的圆内的图像功率与图像总功率的量级关系。 说明了高频成分对于表现图像的轮廓和细节是十分重要的。 图5.12(c)仅滤除掉占总能量的4.5%的高频分量,图像就变得十分模糊,并有明显的振铃效应; 图5.12(d)仅滤除掉占总能量的2.1%的高频分量,图像仍存在着一定程度的模糊和振铃效应; 当图5.12(e)仅滤除掉占总能量的1.0%的高频分量时,图像视觉效果才变得尚可接受; 而当图5.12(f)仅滤除掉占总能量的0.4%的高频分量时,图像才变得在视觉上与原图像比较接近。 图5.12频率域理想低通滤波器的滤波效果及低频特性分析 【例5.5】例5.4的频率域理想低通滤波器的滤波效果及低频特性分析MATLAB程序。 clc; clear all; close all; img=imread('d:\0_matlab图像课编程\lena.jpg'); imshow(img); title('原图像'); FT_img = fft2(double(img));%傅里叶正变换 FTS_img = fftshift(FT_img);%平移频谱图为中心对称的频谱图像数据 [h,w] = size(img); min = 255.0;%最小值赋初值 for i=1:h%求原图像傅里叶频谱的最小值,为标注截止频率做准备 for j=1:w if FTS_img < min min = FTS_img(i,j); end end end M = fix(h/2);N = fix(w/2); %将h/2和w/2分别向零方向(向下)取整 FTS_img1 = FTS_img; %标注截止频率半径的原图像傅里叶频谱图初值 for i=1:h%基于频谱图中心以10、20、40、80为半径标注截止频率 for j=1:w r=sqrt((i-M)^2+(j-N)^2); if (r>9.5&&r<10.5)||(r>19.5&&r<20.5)||(r>39.5&&r<40.5)||(r>79.5&&r<80.5) FTS_img1(i,j) = min; end end end figure; imshow(log(1+abs(FTS_img1)),[]); title('有截止频率半径的原图像傅里叶频谱图'); D0=5; for i=1:4 D0 = D0*2; %分别计算D0为10、20、40、80的理想低通滤波效果 for u=1:h forv=1:w duv=sqrt((u-M)^2+(v-N)^2);%计算频率平面到原点的距离 if duv <= D0 H = 1 %duv≤D0的理想的低通转移函数值 else H = 0 %duv>D0的理想的低通转移函数值 end huv_img(u,v)=FTS_img(u,v)*H;%对傅里叶频谱进行滤波 end end Inverse_I=ifftshift(huv_img);%反移频谱中心 Inverse_img=real(ifft2(Inverse_I));%傅里叶反变换,并取变换结果的实部 Inverse_H_img=uint8(Inverse_img); %转换为8位图像 figure; imshow(Inverse_H_img); title('理想低通滤波结果图像'); end 例5.5中的相关函数及功能说明如下。 sqrt(x): 求x的平方根的函数。 5.3.2巴特沃斯低通滤波器 一个n阶的巴特沃斯低通滤波器(Butterworth lowpass filter)的转移函数定义为 H(u,v)=11+[D(u,v)/D0]2n(5.28) 其中,D0为截止频率; D(u,v)是频率平面从原点到点(u,v)的距离,如式(5.26)所示。 图5.13(a)给出了阶数分别为1、2、3和4的巴特沃斯低通滤波器的转移函数H的横截面图。巴特沃斯低通滤波器的转移函数H的透视图如图5.13(b)所示,该透视图的含义是: 只有那些位于该“草帽体”内的频率范围的信号才能通过,而位于“草帽体”外的频率成分都将被滤除掉。由此可见,巴特沃斯低通滤波器在高低频率间的过渡比较平滑。 图5.13巴特沃斯低通滤波器的转移函数横截面图和透视图 一般情况下取当H的最大值降到某个百分比时对应的频率为截止频率。在式(5.26)中,当D(u,v)=D0时,H(u,v)=0.5(从最大值1降到50%)。一个常用的方法是,取H降到最大值的1/2时的频率为截止频率。 【例5.6】利用巴特沃斯低通滤波器进行图像去噪的MATLAB程序。 clc; clear all; close all; img0=imread('d:\0_matlab图像课编程\lena.jpg'); Noisy_img = imnoise(img0,'salt & pepper' ,0.02);%加入噪声密度为0.02的椒盐噪声 f = double(Noisy_img); FT_img = fft2(f);%傅里叶变换 FTS_img = fftshift(FT_img);% 平移频谱图为中心对称 D0=input('\n请输入非负的截止频率值(10/20/30) D0=')% 正整数 n=input('\n请输入巴特沃斯滤波器的阶数(1/2/3) n=');% 滤波器阶数1~3 [h,w] = size(img0); M=fix(h/2);N=fix(w/2);% 将h/2和w/2分别向零方向(向下)取整 for u=1:h forv=1:w duv=sqrt((u-M)^2+(v-N)^2); huv=1/(1+(duv/D0)^(2*n)); % 计算巴特沃斯低通滤波器转移函数值 huv_img(u,v)=huv*FTS_img(u,v);% 对傅里叶频谱进行滤波 end end Inverse_I=ifftshift(huv_img);% 反移频谱中心 Inverse_img=real(ifft2(Inverse_I));% 傅里叶反变换,并取变换结果的实部 Inverse_H_img=uint8(Inverse_img); % 转换为8位图像 figure; subplot(1,3,1); imshow(img0); title('原图像'); subplot(1,3,2); imshow(Noisy_img); title('加椒盐噪声图像'); subplot(1,3,3); imshow(Inverse_H_img); title('低通滤波结果图像'); 图5.14是利用巴特沃斯低通滤波器进行图像去噪的实验结果。图5.14(a)是原图像,图5.14(b)是加入了噪声密度为0.02的椒盐噪声后的图像,图5.14(c)~图5.14(h)是取不同的截止频率(10和30)和不同的滤波器阶数(1、2和3)的去噪效果图像。分析可知: 滤波器阶数n的值越大,去噪结果图像越模糊(去噪效果越差); 截止频率D0的值越大,去噪结果图像越清晰(去噪效果越好)。 图5.14例5.6利用巴特沃斯低通滤波器进行图像去噪实验结果 5.3.3高斯低通滤波器 由于高斯函数的傅里叶变换和反变换均为高斯函数,并常常用来帮助寻找空间域与频率域之间的联系,所以基于高斯函数的滤波具有特殊的重要意义。 一个二维的高斯低通滤波器(Gaussian lowpass filter)的转移函数定义为 H(u,v)=e-D2(u,v)/(2σ2)(5.29) 其中,D(u,v)是频率平面从原点到点(u,v)的距离,如式(5.26)所示。σ表示高斯曲线扩展的程度,当σ=D0时,可得到高斯低通滤波器的一种更为标准的表示形式 H(u,v)=e-D2(u,v)/(2D20)(5.30) 其中,D0是截止频率; 当D(u,v)=D0时,H下降到它的最大值的0.667处。 图5.15(a)给出了D0=10、D0=20、D0=40和D0=100的高斯低通滤波器的转移函数H的横截面图。高斯低通滤波器的转移函数H的透视图如图5.15(b)所示,该透视图的含义是: 只有那些位于该“草帽体”内的频率范围的信号才能通过,而位于“草帽体”外的频率成分都将被滤除掉。 图5.15高斯低通滤波器的转移函数横截面图和透视图 与巴特沃斯低通滤波器相比,高斯低通滤波器没有振铃现象。另外在需要严格控制低频和高频之间截止频率的过渡的情况下,选择高斯低通滤波器更合适一些。 需要说明的是,在频率域中,滤波器越窄,滤除掉的高频成分就越多,滤波后的图像就越模糊。这一特性正好对应于在空间域中,滤波器越宽(模板尺寸越大),平滑后的图像就越模糊的情况。 5.4基于频率域的图像增强——频率域高通滤波 图像中的边缘和灰度的陡峭变化对应于傅里叶频谱的高频部分,选择能使高频通过、使低频衰减的转移函数H(u,v),就可以根据式(5.20)和式(5.21)实现高通滤波,达到突出图像的高频边缘成分,实现图像增强的效果。 5.4.1理想高通滤波器 一个理想的频率域高通滤波器的转移函数定义为 H(u,v)=0D(u,v)≤D0 1D(u,v)>D0(5.31) 其中,D0为截止频率; D(u,v)为频率平面从原点到点(u,v)的距离,如式(5.26)所示。对于高通滤波来说,当某个频率之下的幅度统统被忽略为0时,该频率就称为截止频率。 理想高通滤波器的含义为,将以半径为D0的圆内的所有频率置零,而让圆外的所有频率毫不衰减地通过。 图5.16(a)给出了理想高通滤波器的转移函数H的横截面图。图5.16(b)给出了转移函数H的透视图,该透视图的含义是: 只有那些位于该圆柱体外的频率范围的信号才能通过,而位于圆柱体内的频率成分都将被滤除掉。 图5.16理想高通滤波器的转移函数横截面图和透视图 与理想低通滤波器一样,在实际中用电子元件实现直上直下的理想高通滤波器是不现实的,所以才将其称为“理想”高通滤波器。 5.4.2巴特沃斯高通滤波器 一个n阶的巴特沃斯高通滤波器(Butterworth highpass filter)的转移函数定义为 H(u,v)=11+[D0/D(u,v)]2n(5.32) 其中,D0为截止频率; D(u,v)是频率平面从原点到点(u,v)的距离,如式(5.26)所示。 图5.17(a)给出了阶数为1的巴特沃斯高通滤波器的转移函数H的横截面图。巴特沃斯高通滤波器的转移函数H的透视图如图5.17(b)所示,该透视图的含义是: 只有那些位于该“倒立形草帽体”外的频率范围的信号才能通过,而位于“倒立形草帽体”内的频率成分都将被滤除掉。与巴特沃斯低通滤波器一样,巴特沃斯高通滤波器在高低频率间的过渡比较平滑。 与巴特沃斯低通滤波器一样,一般情况下取当H的最大值降到某个百分比时对应的频率为截止频率。 图5.17巴特沃斯高通滤波器的转移函数横截面图和透视图 5.4.3高斯高通滤波器 一个截止频率距离原点为D0的高斯高通滤波器(Gaussian highpass filter)的转移函数定义为 H(u,v)=1-e-D2(u,v)/(2D20)(5.33) 其中,D(u,v)是频率平面从原点到点(u,v)的距离,如式(5.26)所示。 图5.18(a)给出了典型的高斯高通滤波器的转移函数H的横截面图。图5.18(b)给出了该高斯形高通滤波器转移函数H的透视图,该透视图的含义是: 只有那些位于该“倒立形草帽体”外的频率范围的信号才能通过,而位于“倒立形草帽体”内的频率成分都将被滤除掉。 图5.18高斯高通滤波器的转移函数横截面图和透视图 【例5.7】利用高斯高通滤波器进行图像增强的MATLAB程序。 clc; clear all; close all; img0=imread('d:\0_matlab图像课编程\girl.jpg'); f = double(img0); FT_img = fft2(f);% 傅里叶变换 FTS_img = fftshift(FT_img);% 平移频谱图为中心对称 D0=input('\n请输入非负的截止频率值(10/20/30/40) D0=')% 正整数 [h,w] = size(img0); M=fix(h/2);N=fix(w/2);% 将h/2和w/2分别向零方向(向下)取整 d=2*D0^2; for u=1:h for v=1:w duv=sqrt((u-M)^2+(v-N)^2); huv=1-exp(-duv^2/d);% 计算高斯高通滤波器转移函数值 huv_img(u,v)=huv*FTS_img(u,v);% 对傅里叶频谱进行滤波 end end Inverse_I=ifftshift(huv_img);% 反移频谱中心 Inverse_img=real(ifft2(Inverse_I));% 傅里叶反变换,并取变换结果的实部 img1=f+Inverse_img;% 形成结果图像 result_img=uint8(img1); % 转换为8位图像 figure;subplot(1,2,1);imshow(img0);title('原图像'); subplot(1,2,2);imshow(result_img);title('高斯高通滤波结果图像'); 例5.7中的相关函数及功能说明如下。 exp(x): 求以e为底的指数函数,即计算ex。 图5.19例5.7利用高斯高通滤波器进行图像增强实验结果 图5.19给出了用高斯高通滤波器进行图像增强时取不同截止频率的实验结果。分析可知: 随着截止频率D0值的增大,增强效果进一步明显,即使对于细线条的头发丝,肉眼也能看出其细微的增强效果变化。但同时,随着D0的增大,低频信息减少较多,图像中灰度剧烈变化的区域也出现了微弱的振铃现象(即图像灰度剧烈变化处产生的震荡,就好像钟被敲击后产生的空气震荡)。 5.5带阻滤波和带通滤波 带阻滤波器(bandstop filter)与带通滤波器(bandpass filter)用于某些区域的某一频率范围内的频率分量,抑制其通过或让其通过。 5.5.1带阻滤波器 在某些应用中,图像的质量可能受到带有一定规律的结构性噪声的影响。比如,图像上叠加有正弦干扰图案就是这类噪声的一个典型情况。当正弦干扰图案比较明显时,会在图像的频谱平面上出现两个比较明显的对称点(由于傅里叶变换的共轭对称性所致)。像这种用于消除以某点为对称中心的给定区域内的频率,或用于阻止以原点为对称中心的一定频率范围内信号通过的问题,就可以用带阻滤波器实现。 一个用于消除以某点为对称中心,以D为半径的圆域上的带阻滤波器,可以通过把以原点为中心的高通滤波器平移到该点得到,设该带阻滤波器的中心为点(u0,v0),半径为D0,则其传递函数定义为 H(u,v)=0D(u,v)≤D0 1D(u,v)>D0(5.34) 其中 D(u,v)=(u-u0)2+(v-v0)212(5.35) 由于傅里叶变换的共轭对称性,要求带阻滤波器必须成对出现,所以一个用于消除以(u0,v0)为中心,以D0为半径的对称区域内的所有频率的理想带阻滤波器的转移函数定义为 H(u,v)=0D1(u,v)≤D0或D2(u,v)≤D0 1其他(5.36) 其中 D1(u,v)=[(u-u0)2+(v-v0)2]1/2(5.37) D2(u,v)=[(u+u0)2+(v+v0)2]1/2(5.38) 利用上述构建带阻滤波器的思路,还可以把前面所讲的几种高通滤波器转变为带阻滤波器,如一种n阶径向对称的巴特沃斯带阻滤波器的传递函数可定义为 H(u,v)=11+D(u,v)WD2(u,v)-D202n(5.39) 其中,W为阻带带宽,D0为阻带中心半径。 图5.20给出了一个典型的带阻滤波器的转移函数H的透视图,该透视图的含义是: 只有那些位于两个圆柱体外的频率范围的信号才能通过,而位于两个圆柱体内的频率成分都将被滤除掉。 图5.20理想带阻滤波器的转移函数的透视图 5.5.2带通滤波器 带通滤波器与带阻滤波器相反,它允许以原点为对称中心的一定频率范围内的信号通过,而将其他频率衰减或抑制。理想的带通滤波器的转移函数可定义为 H(u,v)=1D1(u,v)≤D0或D2(u,v)≤D0 0其他(5.40) 带通滤波器也可以通过对相应的带阻滤波器进行“翻转”获得。若设H′(u,v)为带阻滤波器的传递函数,则对应的带通滤波器的传递函数H(u,v)可定义为 H(u,v)=1-H′(u,v)(5.41) 图5.21是一个典型的带通滤波器的转移函数H的透视图,该透视图的含义是: 只有那些位于两个圆柱体内的频率范围的信号会被通过,而位于两个圆柱体外的频率成分都将被滤除掉。 图5.21理想带通滤波器的转移函数的透视图 习题5 1. 解释下列名词概念。 (1) 低通滤波(2) 高通滤波 2. 已知有数学关系式exp(x)=ex和eix=cosx+isinx,写出下列公式的展开式。 (1) exp-j2π(xu+yv)N (2) expj2π(xu+yv)N 3. 进行图像傅里叶变换的目的是什么? 4. 设N×N的图像f(x,y)的傅里叶变换为F(u,v),请写出反映该图像的灰度平均值与傅里叶变换的直流(DC)分量的关系的公式。 5. 回答下列问题: (1) 对图像f(x,y)进行傅里叶变换后得到的频谱图,低频位于其中心处还是四角处? (2) 当把对图像f(x,y)进行傅里叶变换后,得到的频谱图的坐标原点移到中心对称时,低频位于其中心处还是四角处? (3) 频率域低通滤波实现的图像处理功能是什么? (4) 频率域高通滤波实现的图像处理功能是什么? 6. 编写利用高斯低通滤波器进行图像去噪的MATLAB程序。 7. 编写利用巴特沃斯高通滤波器进行图像增强的MATLAB程序。