数字图像处理 第四章——频率域滤波

1.1背景

1.1.1傅里叶级数和变换

• 傅里叶级数：周期函数可以用正弦和/或余弦乘以加权函数的积分来表示。
• 傅里叶变换：非周期函数（曲线下的面积有限）可以用正弦和/或余弦乘以加权函数的积分来表示，这种公式称为傅里叶变换
  用傅里叶级数或变换表示的函数特征可以通过傅里叶反变换来重建，不丢失仍和信息。

1.1.2 基本概念

• 复数：定义如下：$C=R+jI$ 其中$R$和$I$是实数，$j=\sqrt{-1}$
  复数C的共轭表示为$C^{\ast }$,定义为：$C=R-jI$
  在极坐标下复数表示为：$C=|C|(cos\theta +jsin\theta )$,其中$|C|=\sqrt{R^2+I^2}$是复平面的原点到点$(R,I)$的向量的长度。
  使用欧拉公式：$e^{j\theta }=cos\theta +jsin\theta$可以给出极坐标下复数表示：$C=|C|e^{j\theta }$

• 傅里叶级数
  具有周期T的连续变量t的周期函数$f(t)$可以被描述为乘以适当系数的正弦和余弦和，这个和是傅里叶级数，有如下形式：$f(t)={\sum }_{-\infty }^{\infty }{c}_n{e}^{j\frac{2\pi n}{T}t}$。其中$c_n=\frac{1}{T}{\int }_{-T/2}^{T/2}f(t)e^{-j\frac{2\pi n}{T}t}dt$是系数。

• 冲激及取样特性
  连续变量t在t=0出的单位冲激表示为$\delta (t)$.
  取样特性的一种更为一般的说明涉及位于任意点$t_0$的冲激，表示为：$\delta (t-t_0)$。在这种情况下，取样特性变为${\int }_{-\infty }^{\infty }f(t)\delta (t-t_0)dt=f(t_0)$

• 连续变量函数的傅里叶变换
  取连续函数$f(t)$的傅里叶变换由以下式子定义$F(\mu )={\int }_{-\infty }^{\infty }f(t)e^{-j2\pi \mu t}dt$
  $F(\mu )$可以通过傅里叶反变换获得$f(t)$,$f(t)={\int }_{-\infty }^{\infty }F(\mu )e^{j2\pi \mu t}d\mu$
  说明一个函数可以用变换来进行恢复，可以进行空间域和频率域的变换。

• 卷积的傅里叶变换
  两个连续函数的卷积公式
  $f(t)\star h(t)={\int }_{-\infty }^{\infty }f(\tau )h(t-\tau )d\tau$
  卷积的傅里叶变换公式为：$\Im [f(t)\star h(t)]={\int }_{-\infty }^{\infty }[{\int }_{-\infty }^{\infty }f(\tau )h(t-\tau )d\tau ]e^{-j2\pi \mu t}dt=H(\mu )F(\mu )$。其中$H(\mu )和F(\mu )$分别是$h(t)和f(t)$的傅里叶变换。
  得出结论：空间域中两个函数的卷积的傅里叶变换等于两个函数的傅里叶变换在频率域中的乘积。两个变换的乘积可以通过机选傅里叶反变换得到空间域的卷积。

1.1.3一维傅里叶变换

连续变量t的连续函数f(t)的傅里叶变换$F(\mu )$如式子定义：
$F(\mu )={\int }_{-\infty }^{\infty }f(t)e^{-2j\pi \mu t}dt$
对于$F(\mu )$，可以通过傅里叶反变换得到$f(t)$
$f(t)={\int }_{-\infty }^{\infty }F(\mu )e^{j2\pi \mu t}d\mu$
单变量的离散傅里叶变换(DFT)由式子给出：
$F(\mu )=\frac{1}{M}{\sum }_{x=0}^{M-1}f(x)e^{-j2\pi \mu x/M}，\mu =0,1,2,...,M-1$

$f(x)=\frac{1}{M}{\sum }_{x=0}^{M-1}F(\mu )e^{j2\pi \mu x/M}，\mu =0,1,2,...,M-1$
这是一个傅里叶变换对。

1.1.4二维傅里叶变换

两个连续变量$t和z$的来连续函数$f(t,z)$的傅里叶变换为$F(\mu ,v)$。其二维连续傅里叶变换对由以下两个表达式给出：
$F(\mu ,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)e^{-j2\pi (\mu t+vz)}dtdz$

$f(t,z)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(\mu ,v)e^{j2\pi (\mu t+vz)}d\mu dv$
其中$\mu 和v$是频率变量

二维离散傅里叶变换及反变换
二维傅里叶变换(DFT)由式子指出：
$F(\mu ,v)={\sum }_{x=0}^{M-1}{\sum }_{y=0}^{N-1}f(x,y)e^{-j2\pi (\mu x/M+vy/N)}$，其中$f(x,y)$是大小为$M\ast N$的图像。
使用傅里叶反变换可以得到$f(x,y)$，式子如下:
$\frac{1}{MN}{\sum }_{\mu =0}^{M-1}{\sum }_{v=0}^{N-1}F(\mu ,v)e^{j2\pi (\mu x/M+vy/N)}$

2.1频率域滤波

频率域滤波由修改一幅图像的傅里叶变换然后计算其反变换得到处理后的结果组成。
对于图像的处理可以在空间域和频率域进行处理。关于傅里叶变换的频率分量和衣服图像的空间特性间的关系可以做出某些一般的叙述。而高频和低频和图像有关，低频与图像中缓慢变化的灰度分量有关，高频由灰度的尖锐过度造成（边缘，噪声）。衰减高频通过低频的滤波器（低通滤波器）模糊图像。衰减低频通过高频的滤波器（高通滤波器）锐化图像。

• 频率域滤波的步骤如下：
  1.给定$M\ast N$的输入图像$f(x,y)$,使用填充参数$P=2M和Q=2N$ 。
  2.对$f(x,y)$添加必要数量的0，形成大小为$P\ast Q$的填充后图像$f_p(x,y)$。
  3.$(-1{)}^{x+y}$乘以$f_p(x,y)移到变换中心$。
  4.计算第三步图像的DFT，得到$F(\mu ,v)$。
  5.使用大小为$P\ast Q,中心在(P/2,Q/2)的$滤波函数$H(\mu ,v)$，使用阵列相乘得到乘积$G(\mu ,v)=H(\mu ,v)F(\mu ,v)$。
  6.得到处理后的图像：
  $g_p(x,y)$，将图像的左上象限提取$M\ast N$的区域。得到最终结果$g(x,y)$。

• 空间域和频率域滤波间的对应
  空间域和频率域滤波见的纽带是卷积定理，频率域中的滤波定义为滤波函数$H(\mu ,v)$与输入图像的傅里叶变换$F(\mu ,v)$的乘积。空间域滤波和频率域滤波可以通过傅里叶变换和反变换得到。
  $h(x,y)\iff H(\mu ,v)$

2.1.1频率域滤波器平滑图像

由于图像中的变换和其他尖锐的灰度转变与图像的傅里叶变换的高频内容有联系，使用在频率域平滑中可以通过衰减高频来达到，即为低通滤波。

低通滤波器：理想滤波器，布特沃斯滤波器，高斯滤波器

• 理想低通滤波器
  该滤波器以原点为圆心，以$D_0$为半径的园内，无衰减的通过所有频率，圆外所有的频率都不通过。由函数如下确定：
  $H(u,v)=\{\begin{array}{l}1,D(u,v)\le D_0\\ 0,D(u,v)>D_0\end{array}$
  其中$D(u,v)$是频率域中点(u,v)与频率矩形中间的距离。
  $D(u,v)=[(u-M/2{)}^2+(v-N/2{)}^2{]}^{1/2}$

import numpy as np
import cv2
import matplotlib.pyplot as plt
def ideal_lowpass_filter(image, cutoff_freq):
# 获取图像尺寸
height, width = image.shape[:2]
# 进行二维傅里叶变换
dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
# 构建理想低通滤波器
filter_mask = np.zeros((height, width, 2), np.float32)
cx, cy = width // 2, height // 2
filter_mask[cy - cutoff_freq:cy + cutoff_freq, cx - cutoff_freq:cx + cutoff_freq] = 1
# 将滤波器应用于频域图像
dft_shift_filtered = dft_shift * filter_mask
# 进行逆变换，将图像转换回空域
dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
img_filtered = cv2.idft(dft_filtered_shifted)
img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])
return img_filtered
# 读取图像
image = cv2.imread(‘c:\jimei.jpg’, 0) # 以灰度模式读取
#执行理想低通滤波
cutoff_frequency = 30 #截止频率
filtered_image = ideal_lowpass_filter(image, cutoff_frequency)
#显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap=‘gray’)
plt.title(‘Original Image’), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap=‘gray’)
plt.title(‘Filtered Image’), plt.xticks([]), plt.yticks([])
plt.show()

• 布特沃斯低通滤波器
  截止频率位于距原点$D_0$处n阶布特沃斯低通滤波器的传递函数定义为
  $H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}$
  该传递函数并没有在通过频率和滤除频率之间由明显截至的尖锐的不连续性。在空间域的一阶布特沃斯滤波器没有振铃现象。在二阶滤波器中，振铃现象通常很难察觉，但更高阶数的滤波器中振铃现象会很明显。

import numpy as np
import cv2
import matplotlib.pyplot as plt
def butterworth_lowpass_filter(image, cutoff_freq, order):
# 获取图像尺寸
height, width = image.shape[:2]
# 进行二维傅里叶变换
dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
# 构建布特沃斯低通滤波器
filter_mask = np.zeros((height, width, 2), np.float32)
cx, cy = width // 2, height // 2
for i in range(height):
for j in range(width):
dist = np.sqrt((i - cy) ** 2 + (j - cx) ** 2)
filter_mask[i, j] = 1 / (1 + (dist / cutoff_freq) ** (2 * order))
# 将滤波器应用于频域图像
dft_shift_filtered = dft_shift * filter_mask
# 进行逆变换，将图像转换回空域
dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
img_filtered = cv2.idft(dft_filtered_shifted)
img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])
return img_filtered
#读取图像
image = cv2.imread(‘c:\jimei.jpg’, 0) # 以灰度模式读取
#执行布特沃斯低通滤波
cutoff_frequency = 30 # 截止频率
order = 2 # 阶数
filtered_image = butterworth_lowpass_filter(image, cutoff_frequency, order)
#显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap=‘gray’)
plt.title(‘Original Image’), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap=‘gray’)
plt.title(‘Filtered Image’), plt.xticks([]), plt.yticks([])
plt.show()

• 高斯低通滤波器
  该滤波器的二维形式公式如下：
  $H(u,v)=e^{-D^2(u,v)/2D_0^2}$
  其中，$D_0$是截至频率

import numpy as np
import cv2
import matplotlib.pyplot as plt
def gaussian_lowpass_filter(image, cutoff_freq):
# 获取图像尺寸
height, width = image.shape[:2]
# 进行二维傅里叶变换
dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
# 构建高斯低通滤波器
filter_mask = np.zeros((height, width, 2), np.float32)
cx, cy = width // 2, height // 2
for i in range(height):
for j in range(width):
dist = np.sqrt((i - cy) ** 2 + (j - cx) ** 2)
filter_mask[i, j] = np.exp(-0.5 * (dist / cutoff_freq) ** 2)
# 将滤波器应用于频域图像
dft_shift_filtered = dft_shift * filter_mask
# 进行逆变换，将图像转换回空域
dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
img_filtered = cv2.idft(dft_filtered_shifted)
img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])
return img_filtered
#读取图像
image = cv2.imread(‘c:\jimei.jpg’, 0) # 以灰度模式读取
#执行高斯低通滤波
cutoff_frequency = 30 # 截止频率
filtered_image = gaussian_lowpass_filter(image, cutoff_frequency)
#显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap=‘gray’)
plt.title(‘Original Image’), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap=‘gray’)
plt.title(‘Filtered Image’), plt.xticks([]), plt.yticks([])
plt.show()

2.1.2频率域滤波器锐化图像

锐化图像，因为边缘和其他灰度的急剧变化域高频分量有关，图像的锐化可在频率域通过高通滤波来实现。
高通滤波器是从给定的低通滤波器用下式得到的：
$H_{HP}(u,v)=1-H_{LP}(u,v)$。
其中，$H_{LP}(u,v)$是低通滤波器的传递函数。被低通滤波器衰减的频率能够通过高通滤波器。

• 理想高通滤波器
  二维理想高通滤波器定义为：
  $H(u,v)=\{\begin{array}{l}0,D(u,v)\le D_0\\ 1,D(u,v)>D_0\end{array}$

• 布特沃斯高通滤波器
  截至频率为$D_0$的n阶布特沃斯高通滤波器定义为：
  $H(u,v)=\frac{1}{1+[D_0/D(u,v){]}^{2n}}$
  截至频率越大，使用该滤波器得到的结果越平滑。

• 高斯高通滤波器
  截至频率处在距频率矩形中心距离为$D_0$的高斯高通滤波器的传递函数由式子给出:
  $H(u,v)=1-e^{-D^2(u,v)/2D_0^2}$

2.1.3选择性滤波

当处理指定的频段或频率矩形的小区域。用到两类滤波器，第一类叫做带阻滤波器或带通滤波器。第二类叫做陷波滤波器。

• 带阻滤波器和带通滤波器
  
  一个带通滤波器可以从带阻滤波器得到：
  $H_{BP}(u,v)=1-H_{BR}(u,v)$

• 陷波滤波器
  陷波滤波器是更有用的选择性滤波器。陷波滤波器拒绝(或通过)事先定义的关于频率矩形中心的个邻域的频率。
  $H_{NR}(u,v)={\prod }_{k=1}^Q{H}_k(u,v)H_{-k}(u,v)$