数字图像处理(六)频域滤波器

频域滤波器简介

频域滤波器与空域滤波器相对,空域滤波器是空间卷积块与空间图像做卷积运算,而频域滤波是在将图像进行了离散傅里叶变换之后在频域上与滤波器做乘法,之后再进行傅里叶反变换得到滤波之后的图像。

使用频域滤波器对图像进行滤波的基本步骤:

  • 对图像进行DFT变换得到原图像的频域图像
  • 频域图像与频域滤波器相乘得到新的频域图像
  • 对新的频域图像进行IDFT变换到空域得到新图像。

在此之中有几个要注意的事项:

  • 在DFT运算的时候会将原图像在空间上进行上下左右周期延拓,如果直接对原图像进行DFT运算,则会以原图像四周都是原图像进行计算得到DFT的结果。因此,如果要消除图像的边缘相互影响滤波后的结果,可以将图像变换到原图的四倍大小,将原图像放在四倍图像的左上角,其余地方补零。
  • 为了方便计算,一般将频域图像的原点(即(0,0))放在图像中心,而对图像直接进行DFT得到的频域图像的(0,0)点在图像的四个角,为了移动频域图像的(0,0,)点,在DFT之前对空域图像的每个像素点乘上\((-1)^{i+j}\)即可(\(i,j\)分别为以图像左上角为原点,向右向下为正方向的横纵坐标)。当然,在变换时做了此运算,在反变换之后要记得还原。

对图像进行补零DFT并变换中心的函数如下:

1
2
3
4
5
6
7
8
def fft(img):       #对图像做二维快速傅里叶变换
shape=img.shape
new_img=np.zeros((shape[0]*2,shape[1]*2))
for i in range(shape[0]):
for j in range(shape[1]):
new_img[i][j]=img[i][j]*(-1)**(i+j)
fimg=np.fft.fft2(new_img)
return fimg

将频域图像做相应的反变换的函数如下:

1
2
3
4
5
6
7
8
9
10
def ifft(f_img):    #二维傅里叶反变换
img=np.real(np.fft.ifft2(f_img))
shape=img.shape
new_img=np.zeros((shape[0]//2,shape[1]//2))
for i in range(shape[0]//2):
for j in range(shape[1]//2):
new_img[i][j]=img[i][j]*(-1)**(i+j)
new_img=new_img-np.min(new_img)
new_img=new_img/np.max(new_img)*255
return new_img.astype(np.uint8)

频域低通滤波器

低通滤波器顾名思义,即为通低频,阻高频。如果把频域原点放在图像中心,则理想的低通滤波器的图像为:

但是这种滤波器由于在边界处不连续,一是无法使用物理器件实现,而来即使在计算机上进行仿真也会对图像有一些影响。因此一般不使用这种理想的滤波器

高斯低通滤波器

高斯低通滤波器(GLPF)是最常使用的高通滤波器之一,其生成为二维高斯函数,数学表达式为: \[ H(u,v)=e^{-D^2(u,v)/2D_0 ^2} \] 其中,\(D(u,v)\)为图像中的点到图像中心的距离,\(D_0\)为截止频率,当\(D(u,v)=D_0\)时,GLPF下降到其最大值的0.607处。该滤波器使用图像表示出来为(截止频率分别为10,30,60):

生成高斯滤波器的函数为:

1
2
3
4
5
6
7
8
def gauss_filter(shape,D0):    #生成截止频率为D0的高斯低通滤波器
filt=np.zeros(shape)
P,Q=shape[0],shape[1]
for i in range(shape[0]):
for j in range(shape[1]):
D2=(i-P//2)**2+(j-Q//2)**2
filt[i][j]=math.exp(-D2/(2*D0**2))
return filt

使用高斯低通滤波器的图像滤波效果如下(分别为原图,半径为10,20的高斯滤波):

使用高斯滤波器对图像进行滤波的函数为:

1
2
3
4
5
def gaussian(img,D0):  #使用n阶,截止频率为D0的高斯低通滤波器对图像进行滤波
f_img=fft(img)
shape=f_img.shape
f_result=f_img*gauss_filter(shape,D0)
return ifft(f_result)

巴特沃斯低通滤波器

截止频率位于距原点\(D_0\)处的\(n\)阶巴特沃斯低通滤波器(BLPF)的传递函数定义为: \[ H(u,v)=\frac 1{1+[D(u,v)/D_0]^{2n}} \] 巴特沃斯滤波器用图像表示出来为(n=1,2,5):

生成巴特沃斯低通滤波器的函数为:

1
2
3
4
5
6
7
8
def butter_filter(shape,n,D0):   #生成n阶,半径为D0的butterworth低通滤波器
filt=np.zeros(shape)
P,Q=shape[0],shape[1]
for i in range(shape[0]):
for j in range(shape[1]):
D=math.sqrt((i-P//2)**2+(j-Q//2)**2)
filt[i][j]=1/(1+(D/D0)**(2*n))
return filt

其滤波效果如下(原图,使用3阶半径分别为10,20的滤波器)

使用巴特沃斯低通滤波器对图像进行滤波的函数为:

1
2
3
4
5
def butterwroth(img,n,D0):  #使用n阶,半径为D0的butterworth滤波器对图像进行滤波
f_img=fft(img)
shape=f_img.shape
f_result=f_img*butter_filter(shape,n,D0)
return ifft(f_result)

频域高通滤波器

与低通滤波器相对,其作用为阻低频,通高频,理想高通滤波器在图像上的表示为:

高斯高通滤波器

截止频率为\(D_0\)的高斯高通滤波器(GHPF)的定义为: \[ H(u,v)=1-e^{-D^2(u,v)/2D_0 ^2} \] 该滤波器在图像上表示出来为(图像大小为200*200,截止频率为5,10,30):

生成该滤波器的函数为:

1
2
3
4
5
6
7
8
def gauss_filter(shape,D0):    #生成n阶,滤波半径为D0的高斯高通滤波器
filt=np.zeros(shape)
P,Q=shape[0],shape[1]
for i in range(shape[0]):
for j in range(shape[1]):
D2=(i-P//2)**2+(j-Q//2)**2
filt[i][j]=1-math.exp(-D2/(2*D0**2))
return filt

经过该滤波器滤波之后的效果(依次为原图,半径为10,30的高斯高通滤波后增强):

巴特沃斯高通滤波器

阶数为\(n\),截止频率为\(D_0\)的巴特沃斯高通滤波器的定义为: \[ H(u,v)=\frac 1{1+[D_0/D(u,v)]^{2n}} \] 该滤波器在图像上表示出来为(图像大小为200*200,截止频率为10,阶数为1,2,4)

生成该滤波器的函数为:

1
2
3
4
5
6
7
8
def butter_filter(shape,n,D0):   #生成n阶,半径为D0的butterworth高通滤波器
filt=np.zeros(shape)
P,Q=shape[0],shape[1]
for i in range(shape[0]):
for j in range(shape[1]):
D=math.sqrt((i-P//2)**2+(j-Q//2)**2)+0.01
filt[i][j]=1/(1+(D0/D)**(2*n))
return filt

使用该滤波器对图像进行滤波的效果如下(依次为原图,使用3阶,半径分别为10,30的BHPF进行滤波):

参考文献

[1]数字图像处理[M]:第三版/(美)拉斐尔·C·冈萨雷斯(Rafael C.Gonzalez),(美)理查德·E·伍兹(Richard E. Woods)著;阮秋琦等译,—北京:电子工业出版社,2017.5

[2] lcxxcl_1234.numpy中的fft和scipy中的fft,fftshift以及fftfreq[G/OL].CSDN: 2018-07-02 [2020-04-05].https://blog.csdn.net/lcxxcl_1234/article/details/80889783

[3] 方克明.功率谱和频谱的区别[G/OL].CSDN: 2017-08-10 [2020-04-06].https://blog.csdn.net/godloveyuxu/article/details/77030793