第4章 Python 数字图像处理(DIP) - 频率域滤波7 - 二维DFT和IDFT的一些性质 - 傅里叶频谱和相角

二维DFT和IDFT的一些性质

傅里叶频谱和相角

$$F(u,v)=R(u,v)+jI(u,v)=|F(u,v)|e^{j\varphi (u,v)}\text{\hspace{1em}(4.86)}$$

• 幅度，称为傅里叶频谱（或频谱）

$$|F(u,v)|=[R^2(u,v)+I^2(u,v){]}^{1/2}\text{\hspace{1em}(4.87)}$$

• 相角或相位谱

$$\varphi (u,v)=\text{arctan}\left[\frac{I(u,v)}{R(u,v)}\right]\text{\hspace{1em}(4.88)}$$

• 功率谱
  $$P(u,v)=|F(u,v){|}^2=R^2(u,v)+I^2(u,v)\text{\hspace{1em}(4.89)}$$

实函数的傅里叶变换是共轭对称的，则有

频谱是关于原点偶对称的
$$|F(u,v)|=|F(-u,-v)|\text{\hspace{1em}(4.90)}$$

相角是关于原点奇对称的：
$$\varphi (u,v)=-\varphi (-u,-v)\text{\hspace{1em}(4.91)}$$

$$|F(0,0)|=MN|\bar{f}|,\bar{f}是f(x,y)的平均值\text{\hspace{1em}(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()