數字影象處理(二)傅立葉變換
數字影象處理(二)之傅立葉變換
使用python實現數字影象處理中如下功能:
- 二維傅立葉變換
- 影象二維傅立葉逆變換
- 影象中心化傅立葉變換和譜影象
- 影象2的整數次冪填充
程式碼連結:xiaohuiduan/image_process (github.com)
二維傅立葉變換和逆變換
影象二維傅立葉變換的公式為,\(M,N\)為影象的大小,\(x,y\)為影象的座標。
\[\begin{aligned} &F(u, v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)} \\ &f(x, y)=\frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)} \end{aligned} \]其中,影象的傅立葉變換公式可以等價於:
也就是說,計算二維傅立葉變換可以先計算行的一維DFT然後再計算列的一維DFT。而一維的傅立葉變換可以使用快速傅立葉變換的方式提高運算速度。因此二維傅立葉變換這樣做可以提高計算的速度。
在逆變換中,可以進行如下的變換,將逆變換變成傅立葉正變換的形式。
\[\begin{aligned} &f(x, y)=\frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)} \\ &M N f^{*}(x, y)=\sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F^{*}(u, v) e^{-j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)} \end{aligned} \]其中\(F^*\)
也就是說,可以將傅立葉變換後的結果進行共軛,然後再次進行傅立葉變換併除以\(MN\),便可以得到IDFT後的結果。
進行一維傅立葉變換,如果使用公式\(F(u)=\sum_{x=0}^{M-1} f(x) e^{-j 2 \pi\left(\frac{u x}{M}\right)}\)進行結算,則時間複雜度為\(M^2\),而使用快速傅立葉變換,則時間複雜度為\(log_2M\)。
一維快速傅立葉變換可以藉助numpy中的函式進行實現。下面定義dft2D
函式,對影象順序進行x軸y軸進行一維傅立葉變換,便可以得到二維傅立葉變換。
在影象的傅立葉變換中,如果對傅立葉變換結果進行展示,會發現影象的亮度部分分散在邊緣。於是可以對傅立葉變換進行中心化,這樣人眼看起來就比較明顯。(影象傅立葉的中心化等價於對影象進行image_copy[i][j] = image_copy[i][j]*((-1)**(i+j))
)
def dft2D(self, image: np.array, shift: bool = False):
"""二維傅立葉變換,先對X軸進行傅立葉變換,然後再對Y軸進行傅立葉變換。
Args:
image (np.array): [影象ORkernel]
shift (bool, optional): [是否中心化]. Defaults to False.
Returns:
[type]: [返回傅立葉變換的值:a+bj]
"""
image_copy = image.copy()
if shift:
for i in range(image_copy.shape[0]):
for j in range(image_copy.shape[1]):
image_copy[i][j] = image_copy[i][j]*((-1)**(i+j))
dft_row = np.fft.fft(image_copy, axis=0)
dft_2 = np.fft.fft(dft_row, axis=1)
return dft_2
同時定義IDFT變換函式如下所示,對於IDFT變換我們只需要按照\(f(x, y)=\frac{\sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F^{*}(u, v) e^{-j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)}}{M N}\)來寫程式碼,就很簡單了。
def idft2D(self, dft_2: np.array, shift: bool = False):
"""使用dft進行idft變換,F -> F*(共軛變換)
Args:
dft_2 (np.array): [傅立葉變換]
shift (bool,optional): [是否反中心化]. Defaults to False.
Returns:
[type]: [image進行反傅立葉變換,可能會產生j虛值。返回幅值]
"""
dft_2_copy = dft_2.copy()
idft_row = np.fft.fft(dft_2_copy.conjugate(), axis=0)
image = np.fft.fft(idft_row, axis=1)
image = image/(image.shape[0]*image.shape[1])
if shift:
for i in range(image.shape[0]):
for j in range(image.shape[1]):
image[i][j] = image[i][j]*(-1)**(i+j)
return abs(image)
下面簡單的定義兩個畫圖函式,分別對傅立葉變換後的影象進行畫圖。
def show_dft2D_in_2D(self, title: str, dft2: np.array, set_log: bool = True):
"""在2維平面上展示傅立葉變換,幅值
Args:
title (str): [標題]
dtf2 (np.array): [傅立葉變換的影象]
set_log (bool): [對傅立葉變換後的結果取log]
"""
dft2_copy = dft2.copy()
dft2_copy = abs(dft2_copy)
if set_log:
dft2_copy = np.log2(dft2_copy+1)
self.show_img(title, [dft2_copy], ['gray'])
return dft2_copy
def show_dft2D_in_3D(self, title: str, image: np.array, set_log: bool = True):
"""在3維平面上展示傅立葉變換
Args:
title (str): [標題]
image (np.array): [傅立葉變換的影象]
set_log (bool): [對傅立葉變換後的結果取log]
"""
image = abs(image)
if set_log:
image = np.log10(image+1)
fig = plt.figure()
plt.title(title)
ax3 = plt.axes(projection='3d')
xx = np.arange(0, image.shape[0])
yy = np.arange(0, image.shape[1])
X, Y = np.meshgrid(xx, yy)
ax3.plot_surface(X, Y, image.astype('uint8'), cmap='rainbow')
plt.show()
比如說畫出\(\sigma=1\)的高斯傅立葉變化的圖:
畫出rose原圖,idft還原的圖以及兩者差值
影象2的整數次冪填充
在快速傅立葉變換中,假定\(N\)為2的整數次方。因此如果影象的shape不為2的整數次冪,就需要對其進行填充。
下面是numpy fft函式的一些說明
FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, by using symmetries in the calculated terms. The symmetry is highest when
n
is a power of 2, and the transform is therefore most efficient for these sizes.
對於某個數,找到最近的2的整數次冪可以用如下的演算法(演算法來源於Java HashMap):
def find2power(self,cap):
"""找到距離cap最近的2的整數次冪,activate by the hashmap in the Java
Args:
cap ([type]): [cap > 0]
Returns:
[type]: [2的整數次冪]
"""
n = cap - 1
n |= n >> 1
n |= n >> 2
n |= n >> 4
n |= n >> 8
n |= n >> 1
return n+1
因此,對某張圖片(灰度圖片)進行2的整數次冪的填充操作可以使用如下操作。
image_w = image.shape[0]
image_h = image.shape[1]
fit_w = ImageP.find2power(image_w)
fit_h = ImageP.find2power(image_h)
padding_w = int((fit_w-image_w)/2)
padding_h = int((fit_h - image_h)/2)
image_pad = np.zeros((fit_w,fit_h,1))
image_pad[padding_w:image_w+padding_w,padding_h:image_h+padding_h,:] = image
參考
- 數字影象處理(第三版)
作者: 段小輝