实现图像的二维傅里叶变换并分析其幅度谱和相位谱的性质。

实验1 二维离散傅立叶变换及其谱分析

一、导入图片

import matplotlib.pyplot as plt
import numpy as np
import cv2

path = '.\Images\\Lab1\\'

cameraman = plt.imread(path+'cameraman.jpg')
block  = plt.imread(path+'block.bmp')
car    = plt.imread(path+'car.jpg')
desert = plt.imread(path+'desert.jpg')
face1  = plt.imread(path+'face1.jpg')
face2  = plt.imread(path+'face2.jpg')
fuse   = plt.imread(path+'fuse.jpg')

plt.rcParams['figure.figsize'] = (12.0, 7.0) 
plt.subplot(2,4,1)
plt.imshow(block)
plt.title('block')
plt.subplot(2,4,2)
plt.imshow(car)
plt.title('car')
plt.subplot(2,4,3)
plt.imshow(desert)
plt.title('desert')
plt.subplot(2,4,4)
plt.imshow(fuse, 'gray')
plt.title('fuse')
plt.subplot(2,4,5)
plt.imshow(face1)
plt.title('face1')
plt.subplot(2,4,6)
plt.imshow(face2)
plt.title('face2')
plt.subplot(2,4,7)
plt.imshow(cameraman)
plt.title('cameraman')

Text(0.5, 1.0, 'cameraman')

png

二、幅度谱与相位谱

1. 幅度谱的物理含义

block_gray = cv2.cvtColor(block, cv2.COLOR_BGR2GRAY)
block_ft = np.fft.fft2(block_gray)

plt.rcParams['figure.figsize'] = (10.0, 8.0) 
plt.subplot(1,3,1)
plt.imshow(block)
plt.title('block')
block_shift = np.fft.fftshift(block_ft)
plt.subplot(1,3,2)
plt.imshow(np.log(np.abs(block_shift)))
plt.title('amplitude spectrum')
plt.subplot(1,3,3)
plt.imshow(np.angle(block_shift)*180/np.pi)
plt.title('phase spectrum')

Text(0.5, 1.0, 'phase spectrum')

png

由上图可以看出，幅度谱的水平和竖直方向存在明显的两条谱线，这对应着源图像block中的竖直和水平边缘。

cameraman_gray = cv2.cvtColor(cameraman, cv2.COLOR_BGR2GRAY)
cameraman_ft = np.fft.fft2(cameraman_gray)
plt.rcParams['figure.figsize'] = (10.0, 8.0) 
plt.subplot(1,3,1)
plt.imshow(cameraman)
plt.title('cameraman')
cameraman_shift = np.fft.fftshift(cameraman_ft)
plt.subplot(1,3,2)
plt.imshow(np.log(np.abs(cameraman_shift)))
plt.title('amplitude spectrum')
plt.subplot(1,3,3)
plt.imshow(np.angle(cameraman_shift)*180/np.pi)
plt.title('phase spectrum')

Text(0.5, 1.0, 'phase spectrum')

png

可以看到cameraman中存在各种方向的边缘，而这也与cameraman中各个方向的谱线一一对应，例如：风衣外套是斜右上的边缘，对应着幅度谱中一条明显的斜右下的谱线。

diff = np.log(np.abs(cameraman_shift))-np.log(np.abs(block_shift))
diff_min = np.min(diff)
diff_max = np.max(diff)
plt.rcParams['figure.figsize'] = (4.0, 4.0) 
plt.imshow((diff+diff_min)/(diff_max+diff_min)*255)
plt.title('amplitude spectrum difference')

Text(0.5, 1.0, 'amplitude spectrum difference')

png

上图是cameraman的幅度谱与block的幅度谱差值图像，越亮的表明cameraman在该频点的幅值更大。由于cameraman图像的细节更多，可见其幅度谱的高频分量比block也更多更大。

以上结果均表明：幅度谱的频点相对于原点的方向指示图像边缘的法向，频点距原点的距离指示频率，细节越多的图像，高频分量越多。

2. 相位谱的重要性

face1_gray = cv2.cvtColor(face1, cv2.COLOR_BGR2GRAY)
face1_ft = np.fft.fft2(face1_gray)
face1_AM = np.abs(face1_ft)
face1_PH = np.angle(face1_ft)
face1_AM_Restructure = np.fft.ifft2(face1_AM)
face1_PH_Restructure = np.fft.ifft2(np.exp(1j*face1_PH))
face1_Restructure    = np.fft.ifft2(face1_ft)

plt.rcParams['figure.figsize'] = (10.0, 8.0) 
plt.subplot(2,3,1)
plt.imshow(face1_gray, 'gray')
plt.title('face1 gray')
plt.subplot(2,3,2)
plt.imshow(np.log(np.fft.fftshift(face1_AM)))
plt.title('amplitude spectrum')
plt.subplot(2,3,3)
plt.imshow(np.fft.fftshift(face1_PH*180/np.pi))
plt.title('phase spectrum')
plt.subplot(2,3,4)
plt.imshow(np.abs(face1_AM_Restructure))
plt.title('AM Restructure')
plt.subplot(2,3,5)
plt.imshow(np.abs(face1_PH_Restructure))
plt.title('PH Restructure')
plt.subplot(2,3,6)
plt.imshow(np.abs(face1_Restructure), 'gray')
plt.title('AM+PH Restructure')

Text(0.5, 1.0, 'AM+PH Restructure')

png

由上图可以看出，幅度谱与相位谱一起的逆变换能够很好地恢复原图像的信息；而相位谱的逆变换仍有眼睛和头的大致轮廓，表明相位谱涵盖的是的位置信息，幅度谱涵盖强度信息。因此在逆变换后，幅度谱的结果基本看不出原图像的轮廓，而相位谱的结果能隐约看见。

face2_gray = cv2.cvtColor(face2, cv2.COLOR_BGR2GRAY)
face2_ft = np.fft.fft2(face2_gray)
face2_AM = np.abs(face2_ft)
face2_PH = np.angle(face2_ft)
mix_1p2a = np.fft.ifft2(face2_AM*np.exp(face1_PH*1j))
mix_2p1a = np.fft.ifft2(face1_AM*np.exp(face2_PH*1j))
face1_oDr = face1_gray-np.abs(mix_1p2a)
face2_oDr = face2_gray-np.abs(mix_2p1a)
face1_diff_min = np.min(face1_oDr); face1_diff_max = np.max(face1_oDr)
face2_diff_min = np.min(face2_oDr); face2_diff_max = np.max(face2_oDr)
face1_diff = (face1_oDr+face1_diff_min)/(face1_diff_max+face1_diff_min)*255
face2_diff = (face2_oDr+face2_diff_min)/(face2_diff_max+face2_diff_min)*255

plt.rcParams['figure.figsize'] = (14.0, 14.0) 
plt.subplot(1,4,1)
plt.imshow(np.abs(mix_1p2a), 'gray')
plt.title('girl phase + boy amplitude')
plt.subplot(1,4,2)
plt.imshow(np.abs(mix_2p1a), 'gray')
plt.title('boy phase + girl amplitude')
plt.subplot(1,4,3)
plt.imshow(face1_diff, 'gray')
plt.title('Orignal gril - (gril PH + boy AM)')
plt.subplot(1,4,4)
plt.imshow(face2_diff, 'gray')
plt.title('Orignal boy - (boy PH + gril AM)')

Text(0.5, 1.0, 'Orignal boy - (boy PH + gril AM)')

png

上图显示：女生的相位谱加男生的幅度谱恢复出的图像与女生更为近似，男生的相位谱加女生的幅度谱恢复出的图像与男生更为近似；原女生图像与女生为相位谱恢复出的图像之间只在不同位置的灰度值不同，没有轮廓上的区别，男生与男生为相位谱恢复出的图像之间的关系亦是如此。

car_gray = cv2.cvtColor(car, cv2.COLOR_BGR2GRAY)
desert_gray = cv2.cvtColor(desert, cv2.COLOR_BGR2GRAY)
fuse_gray = fuse/np.max(fuse)*255
car_ft = np.fft.fftshift(np.fft.fft2(car_gray))
desert_ft = np.fft.fftshift(np.fft.fft2(desert_gray))
fuse_ft = np.fft.fftshift(np.fft.fft2(fuse_gray))

plt.rcParams['figure.figsize'] = (10.0, 8.0) 
plt.subplot(3,3,1)
plt.imshow(car_gray, 'gray')
plt.title('car')
plt.subplot(3,3,2)
plt.imshow(desert_gray, 'gray')
plt.title('desert')
plt.subplot(3,3,3)
plt.imshow(fuse_gray, 'gray')
plt.title('fuse')
plt.subplot(3,3,4)
plt.imshow(np.log(np.abs(car_ft)))
plt.title('car AM')
plt.subplot(3,3,5)
plt.imshow(np.log(np.abs(desert_ft)))
plt.title('desert AM')
plt.subplot(3,3,6)
plt.imshow(np.log(np.abs(fuse_ft)))
plt.title('fuse AM')
plt.subplot(3,3,7)
plt.imshow(np.angle(car_ft))
plt.title('car PH')
plt.subplot(3,3,8)
plt.imshow(np.angle(desert_ft))
plt.title('desert PH')
plt.subplot(3,3,9)
plt.imshow(np.angle(fuse_ft))
plt.title('fuse PH')

Text(0.5, 1.0, 'fuse PH')

png

有上图可以看出，fuse是desert和car两张图片的叠加，因此可见fuse的幅度谱在desert的幅度谱的基础上多出了一些低频分量，这些分量与car幅度谱中的一致；同时由于fuse的大面积轮廓与desert一致，因此其幅度谱总体上与desert更为近似，但由于多了car的局部轮廓，因此相位谱与desert的不完全相同。

以上现象均表明：幅度谱含有的是强度信息，相位谱含有的是位置信息。失去了相位谱的信息，灰度值只是平面波的直接叠加，丢失了原图像的轮廓信息；失去了幅度谱的信息，虽然能恢复出原图像的大致轮廓，但各处的灰度值与原图有较大差异；而只有在幅度谱与相位谱的叠加下才能完整地恢复出原图像。

3. 频谱滤波

3.1 高通滤波

(rows,cols,_) = cameraman.shape
r_in = 50
r_out = max(rows, cols)
mask = np.ones((rows, cols),np.uint8)
center = [rows//2, cols//2]
x, y = np.ogrid[:rows, :cols]
mask_area = np.logical_and(((x - center[0]) ** 2 + (y - center[1]) ** 2 >= r_in ** 2),
                            ((x - center[0]) ** 2 + (y - center[1]) ** 2 <= r_out ** 2))

plt.rcParams['figure.figsize'] = (14.0, 14.0)
plt.subplot(1,4,1)
plt.imshow(cameraman)

plt.subplot(1,4,2)
plt.imshow(np.log(np.abs(cameraman_shift)))

plt.subplot(1,4,3)
amp_af_f = mask_area*np.abs(cameraman_shift)
plt.imshow(np.log(amp_af_f))

plt.subplot(1,4,4)
cameraman_af_f = np.fft.ifft2(np.fft.fftshift(amp_af_f)*np.exp(np.angle(cameraman_ft)*1j))
plt.imshow(np.abs(cameraman_af_f), 'gray')

<ipython-input-51-21b805fa3f94>:19: RuntimeWarning: divide by zero encountered in log
  plt.imshow(np.log(amp_af_f))





<matplotlib.image.AxesImage at 0x2a1cfd79d90>

png

3.2 低通滤波

(rows,cols,_) = cameraman.shape
r_in = 0
r_out = 50
mask = np.ones((rows, cols),np.uint8)
center = [rows//2, cols//2]
x, y = np.ogrid[:rows, :cols]
mask_area = np.logical_and(((x - center[0]) ** 2 + (y - center[1]) ** 2 >= r_in ** 2),
                            ((x - center[0]) ** 2 + (y - center[1]) ** 2 <= r_out ** 2))

plt.rcParams['figure.figsize'] = (14.0, 14.0)
plt.subplot(1,4,1)
plt.imshow(cameraman)

plt.subplot(1,4,2)
plt.imshow(np.log(np.abs(cameraman_shift)))

plt.subplot(1,4,3)
amp_af_f = mask_area*np.abs(cameraman_shift)
plt.imshow(np.log(amp_af_f))

plt.subplot(1,4,4)
cameraman_af_f = np.fft.ifft2(np.fft.fftshift(amp_af_f)*np.exp(np.angle(cameraman_ft)*1j))
plt.imshow(np.abs(cameraman_af_f), 'gray')

<ipython-input-48-0bce9c5c0026>:19: RuntimeWarning: divide by zero encountered in log
  plt.imshow(np.log(amp_af_f))





<matplotlib.image.AxesImage at 0x2a1cc2dc100>

png

通过上面两个实验可以看出：低频分量主要对整副图像的强度的综合度量，高频分量主要是对图像边缘和轮廓的度量，图像的能量主要集中在低频分量。

实验代码见： https://github.com/Yangliangzhe/Image-Processing

二维傅里叶变换