引言
“任何”周期信號都可以用一系列成諧波關系的正弦曲線來表示。
—— 讓·巴普蒂斯·約瑟夫·傅里葉
由於頻率和與時間變化率直接相關，因此容易想到 傅里葉變換中的頻率與圖像中的灰度值變化模式 也可以關聯起來。
空間域 英文： spatial domain。 釋義： 又稱圖像空間(image space)。由圖像像元組成的空間。在圖像空間中以長度(距離)為自變量直接對像元值進行處理稱為空間域處理。
頻率域 英文： spatial frequency domain。 釋義： 以頻率(即波數)為自變量描述圖像的特征,可以將一幅圖像像元值在空間上的變化分解為具有不同振幅、空間頻率和相位的簡振函數的線性疊加,圖像中各種頻率成分的組成和分布稱為空間頻譜。這種對圖像的頻率特征進行分解、處理和分析稱為頻率域處理或波數域處理。
二者關系：
空間域與頻率域可互相轉換。在頻率域中可以引用已經很成熟的頻率域技術,處理的一般步驟為：
①對圖像施行二維離散傅立葉變換或小波變換,將圖像由圖像空間轉換到頻域空間。
②在頻率域中對圖像的頻譜作分析處理,以改變圖像的頻率特征。即設計不同的數字濾波器,對圖像的頻譜進行濾波。
③通過反變換將圖像從頻率域變換到空間域。
頻率域處理主要用於與圖像空間頻率有關的處理中。如圖像恢復、圖像重建、輻射變換、邊緣增強、圖像銳化、圖像平滑、噪聲壓制、頻譜分析、紋理分析等處理和分析中。
1.圖像傅里葉變換的物理意義
頻譜圖(點擊查看代碼)
import cv2
import numpy as np
import matplotlib.pyplot as plt
# plt 顯示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False
img = cv2.imread("640_400.jpg", 0)
# 傅里葉變換
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
f_magnitude_spectrum = np.log(np.abs(f))
fshift_magnitude_spectrum = np.log(np.abs(fshift))
# opencv 顯示
# magnitude_spectrum = magnitude_spectrum.astype("int8")
# cv2.imshow("test", magnitude_spectrum)
# cv2.waitKey(0)
# plt 顯示
plt.subplot(131)
plt.imshow(img, cmap='gray')
plt.title("原圖"), plt.xticks([]), plt.yticks([])
plt.subplot(132)
plt.imshow(f_magnitude_spectrum, cmap='gray')
plt.title("二維FFT"), plt.xticks([]), plt.yticks([])
plt.subplot(133)
plt.imshow(fshift_magnitude_spectrum, cmap='gray')
plt.title("零頻率分量移至頻譜中心"), plt.xticks([]), plt.yticks([])
plt.show()
# 1. np.fft.fft2() # Numpy 提供的實現傅里葉變換的函數
# - 這里需要注意的是，參數“原始圖像”的類型是灰度圖像，函數的返回值是一個復數數組（complex ndarray）。
# - 經過該函數的處理，就能得到圖像的頻譜信息，此時，圖像頻譜中的零頻率分量位於頻譜圖像（頻域圖像）的左上角
# 2. np.fft.fftshift() # 將零頻率成分移動到頻域圖像的中心位置
# 3. np.log(np.abs(頻譜值)) # 對圖像進行傅里葉變換后，得到的是一個復數數組。為了顯示為圖像，需要將它們的值調整到 [0,255] 的灰度空間內
說明：
- 圖像的頻率是表征圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。
- 如：大面積的沙漠在圖像中是一片灰度變化緩慢的區域，對應的頻率值很低；而對於地表屬性變換劇烈的邊緣區域在圖像中是一片灰度變化劇烈的區域，對應的頻率值較高。
- 傅里葉變換在實際中有非常明顯的物理意義，從物理效果看，傅里葉變換是將圖像從空間域轉換到頻率域，其逆變換是將圖像從頻率域轉換到空間域。換句話說，傅里葉變換的物理意義是將圖像的灰度分布函數變換為圖像的頻率分布函數。
- 為了便於頻域和頻譜分析，常在變換后進行頻譜中心化，即對調頻譜的四個象限。
- 頻譜中心化后，中間最亮的點是低頻率，屬於直流分量，越往外頻率越高。
- 低頻 對應 圖像中變化緩慢的灰度分量
- 高頻 對應 圖像中灰度急劇變化的物體邊緣或其他分量
低通濾波器(點擊查看代碼)
import cv2
import numpy as np
import matplotlib.pyplot as plt
# plt 顯示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False
img = cv2.imread("640_400.jpg", 0)
# 傅里葉變換
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 獲取原圖長寬和中心點
(h, w) = img.shape
ch, cw = h//2, w//2
# 低通濾波提取高頻分量
mask = np.zeros((h, w), np.uint8)
mask[ch-30:ch+30, cw-30:cw+30] = 1 # 低通濾波器
fshift = fshift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.subplot(131)
plt.imshow(img, cmap= "gray")
plt.title("原圖"), plt.xticks([]), plt.yticks([])
plt.subplot(132)
plt.imshow(mask * 255, cmap= "gray")
plt.title("低通濾波器"), plt.xticks([]), plt.yticks([])
plt.subplot(133)
plt.imshow(img_back, cmap="gray")
plt.title("低通濾波后（低頻分量）"), plt.xticks([]), plt.yticks([])
plt.show()
高通濾波器（點擊查看代碼）
import cv2
import numpy as np
import matplotlib.pyplot as plt
# plt 顯示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False
img = cv2.imread("640_400.jpg", 0)
# 傅里葉變換
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 獲取原圖長寬和中心點
(h, w) = img.shape
ch, cw = h//2, w//2
# 低通濾波提取高頻分量
mask = np.zeros((h, w), np.uint8)
mask[ch-30:ch+30, cw-30:cw+30] = 1 # 低通濾波器
fshift = fshift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.subplot(131)
plt.imshow(img, cmap= "gray")
plt.title("原圖"), plt.xticks([]), plt.yticks([])
plt.subplot(132)
plt.imshow(mask * 255, cmap= "gray")
plt.title("低通濾波器"), plt.xticks([]), plt.yticks([])
plt.subplot(133)
plt.imshow(img_back, cmap="gray")
plt.title("低通濾波后（低頻分量）"), plt.xticks([]), plt.yticks([])
plt.show()
2. 頻率濾波步驟
頻率濾波步驟
- 輸入圖像f(x,y)
- 做二維FFT得到 F(u,v)
- 頻譜中心化
- 構建一個實對稱濾波器傳遞函數 H(u,v)
- 采用對應像素相乘得到G(u,v)=H(u,v)F(u,v)
- 對 G(u,v) 進行 IDFT 得到濾波后的圖像
3. 常見濾波器
3.1 低通濾波器
常見低通濾波器（點擊查看代碼）
import matplotlib.pyplot as plt
import numpy as np
def cal_distance(pa, pb):
from math import sqrt
dis = sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
return dis
# 濾波器模板
class Filter:
@classmethod
def generate_filter(cls, d, shape, *args, **kwargs):
# 這里要反過來 shape 的兩個維度
transfer_matrix = np.zeros((shape[0], shape[1]))
center_point = tuple(map(lambda x: (x - 1) // 2, shape))
for i in range(transfer_matrix.shape[0]):
for j in range(transfer_matrix.shape[1]):
dist = cal_distance(center_point, (i, j))
transfer_matrix[i, j] = cls.get_one(d, dist, *args, **kwargs)
return transfer_matrix
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return 1
# 理想低通濾波器
class ILPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
if dist <= d:
return 1
else:
return 0
# 高斯低通濾波器
class GLPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return np.exp(-(dist ** 2) / (2 * (d ** 2)))
# 巴特沃斯低通濾波器
class BLPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
n = kwargs["n"]
return 1 / ((1 + dist / d) ** (2 * n))
if __name__ == "__main__":
# plt 顯示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False
ILPFFilter = ILPFFilter.generate_filter(10, (100, 100), w=5)
GLPFFilter = GLPFFilter.generate_filter(10, (100, 100), w=5)
BLPFFilter = BLPFFilter.generate_filter(10, (100, 100), w=5, n=1)
plt.subplot(131)
plt.imshow(ILPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
plt.title("理想低通濾波器")
plt.subplot(132)
plt.imshow(GLPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
plt.title("高斯低通濾波器")
plt.subplot(133)
plt.imshow(BLPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
plt.title("巴特沃斯低通濾波器 n=1")
plt.show()
理想低通濾波器
傳遞函數為:
\[H(u,v)= \left\{\begin{matrix} 1, D(u,v)\leqslant D_0 \\ 0, D(u,v)> D_0 \end{matrix}\right. \]
式中，\(D_0\)是一個正常數，\(D(u,v)\)是頻率域中點\((u,v)\)到\(P×Q\)頻率矩形中心的距離，即
\[D(u,v) = [(u-\frac{P}{2})^{2}+(v-\frac{Q}{2})^{2}]^{\frac{1}{2}} \]
高斯低通濾波器
傳遞函數為:
\[H(u,v)=e^{-D^{2}(u,v)/2D_{0}^{2}} \]
巴特沃斯低通濾波器
傳遞函數為:
\[H(u,v)=\frac{1}{1+[D(u,v)/D_0)]^{2n}]} \]
巴特沃斯濾波器的形狀由濾波器階數 \(n\)控制，當階數取值較大時，巴特沃斯濾波器接近於理想濾波器，取值較小時，接近於高斯濾波器。
3.2 高通濾波器
常見高通濾波器（點擊查看代碼）
# 理想高通濾波器
class IHPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
if dist <= d:
return 0
else:
return 1
# 高斯高通濾波器
class GHPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return 1 - np.exp(-(dist ** 2) / (2 * (d ** 2)))
# 巴特沃斯高通濾波器
class BHPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
n = kwargs["n"]
return 1 - 1 / ((1 + dist / d) ** (2 * n))
# 拉普拉斯
class LPFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return 4 * math.pi**2 * dist**2
在頻率域中用1減去低通濾波器傳遞函數，會得到對應的高通濾波器傳遞函數。
3.3 選擇性濾波
處理特定頻帶的濾波器稱為頻帶濾波器。
- 若某個頻帶中的頻率被濾除，則稱該頻帶濾波器為帶阻濾波器；
- 若某個頻帶中的頻率被通過，則稱該頻帶濾波器為帶通濾波器；
- 可處理小頻率區域的濾波器稱為陷波濾波器，根據陷波區域中的頻率被濾除還是被通過，可將陷波濾波器進一步分為帶阻陷波濾波器和帶通陷波濾波器。
常見選擇性濾波（點擊查看代碼）
# 巴特沃斯帶通濾波
class BBPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
n = kwargs["n"]
if dist == d:
return 1
else:
return 1- 1 / (1 + ((dist * w) / (dist ** 2 - d ** 2)) ** (2 * n))
# 巴特沃斯帶阻濾波
class BBRFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
n = kwargs["n"]
if dist == d:
return 0
else:
return 1 / (1 + ((dist * w) / (dist ** 2 - d ** 2)) ** (2 * n))
# 高斯帶通濾波
class GBPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
return np.exp(-((dist ** 2 - d ** 2) / (d * w)) ** 2)
# 高斯帶阻濾波
class GBRFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
return 1-np.exp(-((dist ** 2 - d ** 2) / (d * w)) ** 2)
# 理想帶通濾波器
class IBPFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
if d - w / 2 <= dist <= d + w / 2:
return 1
else:
return 0
# 理想帶阻濾波器
class IBRFFilter(Filter):
@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
w = kwargs["w"]
if d - w / 2 <= dist <= d + w / 2:
return 0
else:
return 1
空間域和頻率域濾波之間的對應關系
空間卷積運算非常適合采用硬件和/或固件的小型核。但在使用普通計算機時，使用快速傅里葉變換（FFT）算法計算DFT的頻率域方法，要比空間卷積運算的速度快幾百倍，具體取決於所用核的大小。
頻率域中的濾波概念更加直觀，且頻率域中的濾波器設計也更容易。充分利用頻率域和空間域性質的一種方法是，在頻率域中規定一個濾波器，計算其IDFT，然后利用生成的全尺寸空間核的性質，指導構建較小的核。
\[h(x,y) = F ^{-1}[H(u,v)] \]
驗證1：使用低通頻率域濾波器平滑圖像
- 理想低通濾波器
- 高斯低通濾波器
- 巴特沃斯低通濾波器
驗證2：使用高通濾波器銳化圖像
- 理想高通濾波器
- 高斯高通濾波器
- 巴特沃斯高通濾波器
- 拉普拉斯
- 同態濾波
總結
- 使用低通頻率域濾波器平滑圖像，常見的低通濾波器有理想低通濾波器、高斯低通濾波器、巴特沃斯低通濾波器。
- 使用高通濾波器銳化圖像，常見的高通濾波器有理想高通濾波器、高斯高通濾波器、巴特沃斯高通濾波。
- 理想濾波器非常尖銳；高斯濾波器非常平滑；巴特沃斯濾波器的形狀由濾波器階數控制，當階數取值較大時，巴特沃斯濾波器接近於理想濾波器，取值較小時，接近於高斯濾波器。
- 理想濾波器因為具有振鈴效應，並不實用；頻率域高斯函數的傅里葉反變化也是高斯函數，說明通過反傅里葉變換后得到的空間高斯濾波器核沒有振鈴效應；針對巴特沃斯濾波器，濾波器階數為1時，沒有振鈴效應，階數越大，振鈴效應越明顯。
- 在頻率域中用1減去低通濾波器傳遞函數，會得到對應的高通濾波器傳遞函數。
- 陷波濾波器是最有用的選擇性濾波器，可去除周期干擾、印刷物圖像中的莫爾模式。
- 空間域濾波核頻率域濾波之間的紐帶是卷積定理。
- 頻率域卷積比空間域卷積快幾百倍（具體取決於所用核的大小）。