目录
- 二维DFT和IDFT的一些性质
- 傅里叶频谱和相角
二维DFT和IDFT的一些性质
傅里叶频谱和相角
F(u,v)=R(u,v)+jI(u,v)=∣F(u,v)∣ejϕ(u,v)(4.86)F(u, v) = R(u, v) + jI(u, v) = |F(u, v)|e^{j\phi(u,v)} \tag{4.86}F(u,v)=R(u,v)+jI(u,v)=∣F(u,v)∣ejϕ(u,v)(4.86)
- 幅度,称为傅里叶频谱(或频谱)
∣F(u,v)∣=[R2(u,v)+I2(u,v)]1/2(4.87)|F(u, v)| = \Big[ R^2(u, v) + I^2(u, v)\Big]^{1/2} \tag{4.87}∣F(u,v)∣=[R2(u,v)+I2(u,v)]1/2(4.87)
- 相角或相位谱
ϕ(u,v)=arctan[I(u,v)R(u,v)](4.88)\phi(u,v) = \text{arctan}\Bigg[ \frac{I(u, v)}{R(u, v)}\Bigg] \tag{4.88}ϕ(u,v)=arctan[R(u,v)I(u,v)](4.88)
- 功率谱
P(u,v)=∣F(u,v)∣2=R2(u,v)+I2(u,v)(4.89)P(u, v) = |F(u, v)|^2 = R^2(u, v) + I^2(u, v) \tag{4.89}P(u,v)=∣F(u,v)∣2=R2(u,v)+I2(u,v)(4.89)
实函数的傅里叶变换是共轭对称的,则有
频谱是关于原点偶对称的
∣F(u,v)∣=∣F(−u,−v)∣(4.90)|F(u, v)| = |F(-u, -v)| \tag{4.90}∣F(u,v)∣=∣F(−u,−v)∣(4.90)
相角是关于原点奇对称的:
ϕ(u,v)=−ϕ(−u,−v)(4.91)\phi(u,v) = - \phi(-u, -v) \tag{4.91}ϕ(u,v)=−ϕ(−u,−v)(4.91)
∣F(0,0)∣=MN∣fˉ∣,fˉ是f(x,y)的平均值(4.93)|F(0, 0)| = MN |\bar f|, \quad \bar f 是f(x, y)的平均值\tag{4.93}∣F(0,0)∣=MN∣fˉ∣,fˉ是f(x,y)的平均值(4.93)
def spectrum_fft(fft):"""return FFT spectrum"""return np.sqrt(np.power(fft.real, 2) + np.power(fft.imag, 2))def phase_fft(fft):"""return FFT phase angle"""return np.angle(fft)
# 矩形函数的频谱
img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0424(a)(rectangle).tif', -1)# FFT
img_fft = np.fft.fft2(img.astype(np.float32))# 非中心化的频谱
spectrum = spectrum_fft(img_fft)
spectrum = np.uint8(normalize(spectrum) * 255)# 中心化
fshift = np.fft.fftshift(img_fft) # 将变换的频率图像四角移动到中心# 中心化后的频谱
spectrum_fshift = spectrum_fft(fshift)
spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum_fshift)plt.figure(figsize=(15, 15))
plt.subplot(2, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 2), plt.imshow(spectrum, cmap='gray', vmin=0, vmax=255), plt.title('FFT Not Shift to Center'),
plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 3), plt.imshow(spectrum_fshift_n, cmap='gray', vmin=0, vmax=255),
plt.title('FFT Shift but not log'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 4), plt.imshow(spectrum_log, cmap='gray'), plt.title('FFT Shift log'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()
# 矩形函数的频谱,平移,旋转
# shift,平移后的频谱与原图像的频谱相同
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0424(a)(rectangle).tif', -1)#[[1, 0, 250], [0, 1, -250]] 表示平移变换:其中250表示水平方向上向右的平移距离,-250表示竖直方向上的向上平移距离。
mat_translation = np.float32([[1, 0, 250],[0, 1, -250]]) #变换矩阵:设置平移变换所需的计算矩阵:2行3列
img_1 = cv2.warpAffine(img_ori, mat_translation, (img_ori.shape[0], img_ori.shape[1])) #变换函数fft_shift = np.fft.fft2(img_ori.astype(np.float32))
fshift_shift = np.fft.fftshift(fft_shift) # 将变换的频率图像四角移动到中心
res_shift = spectrum_fft(fshift_shift)
res_log_shift = np.log(1 + res_shift)# rotate,旋转后的频谱也旋转了
matRotate = cv2.getRotationMatrix2D((img_ori.shape[0]*0.5, img_ori.shape[1]*0.5), -30, 1) # mat rotate: 1 center, 2 angle, 3 scale
img_rotate = cv2.warpAffine(img_ori, matRotate, (img_ori.shape[0], img_ori.shape[1]), flags=cv2.INTER_CUBIC)fft_rotate = np.fft.fft2(img_rotate.astype(np.float32))
fshift_rotate = np.fft.fftshift(fft_rotate) # 将变换的频率图像四角移动到中心
res_rotate = spectrum_fft(fshift_rotate)
res_log_rotate = np.log(1 + res_rotate)plt.figure(figsize=(15, 15))
plt.subplot(2, 2, 1), plt.imshow(img_1, cmap='gray'), plt.title('Shift Rectangle'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 2), plt.imshow(res_log_shift, cmap='gray'), plt.title('Shift rectangle FFT'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 3), plt.imshow(img_rotate, cmap='gray'), plt.title('Rotate Rectangle'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 4), plt.imshow(res_log_rotate, cmap='gray'), plt.title('Rotate Rectangle FFT'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()
中心化矩形的相角图,平移后的相角图, 旋转后的相角图
# 中心化矩形的相角图,平移后的相角图, 旋转后的相角图
phi = phase_fft(fshift)
phi_shift = phase_fft(fft_shift)
phi_rotate = phase_fft(fshift_rotate)# 下面实现是一样的功能
# phi = np.angle(fshift)
# phi_shift = np.angle(fft_shift)
# phi_rotate = np.angle(fshift_rotate)plt.figure(figsize=(15, 15))
plt.subplot(1, 3, 1), plt.imshow(phi, cmap='gray'), plt.title('Shift Rectangle'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(phi_shift, cmap='gray'), plt.title('Shift rectangle FFT'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(phi_rotate, cmap='gray'), plt.title('Rotate Rectangle'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()
频谱和相角对图像信息的贡献
# 频谱和相角对图像信息的贡献
img_woman = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0427(a)(woman).tif', -1)
img_rectangle = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0424(a)(rectangle).tif', -1)
# 原图是1024x1024,下采样为512x512
img_rectangle = img_rectangle[::2, ::2]# FFT
fft_woman = np.fft.fft2(img_woman.astype(np.float32))
fft_rectangle = np.fft.fft2(img_rectangle.astype(np.float32))# 中心化
fshift_woman = np.fft.fftshift(fft_woman)
fshift_rectangle = np.fft.fftshift(fft_rectangle)# 频谱
spectrum_woman = spectrum_fft(fshift_woman)
spectrum_rectangle = spectrum_fft(fshift_rectangle)# 相角
phi_woman = phase_fft(fshift_woman)
phi_rectangle = phase_fft(fshift_rectangle)# 只用相角重建,未做增强的情况
img_phase_recon = np.abs(np.fft.ifft2(np.fft.ifftshift(np.exp(phi_woman * 1j))))
# 下面这种与课本说明不符,但结果更接近增强的效果
# img_phase_recon = np.abs(np.fft.ifft2(np.fft.ifftshift(phi_woman)))# 只用频谱重建
img_spec_recon = np.abs(np.fft.ifft2(np.fft.ifftshift(spectrum_woman)))# 用女人的相角与矩形的频谱重建
f_woman_rectangle = spectrum_rectangle * np.exp(phi_woman * 1j)
img_woman_rectangle = np.abs(np.fft.ifft2(np.fft.ifftshift(f_woman_rectangle)))# 用女人的频谱与矩形的相角重建
f_woman_rectangle = spectrum_woman * np.exp(phi_rectangle * 1j)
img_rectangle_woman = np.abs(np.fft.ifft2(np.fft.ifftshift(f_woman_rectangle)))plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 2), plt.imshow(phi_woman, cmap='gray'), plt.title('Phase'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 3), plt.imshow(img_phase_recon, cmap='gray'), plt.title('Reconstruct with Phase only'),
plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 4), plt.imshow(img_spec_recon, cmap='gray'), plt.title('Reconstruct with Spectrume only'),
plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 5), plt.imshow(img_woman_rectangle, cmap='gray'), plt.title('Reconstruct with Woman Phase and Rectangle Spectrume'),
plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 6), plt.imshow(img_rectangle_woman, cmap='gray'), plt.title('Reconstruct with Woman Spectrum and Rectangle Phase'),
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()