(十九) 傅里叶变换(Fourier Transform)

1、原理

对于图像，使用2D离散傅里叶变换(DFT)来查找频域。对于一个正弦信号 ,x(t)=Asin(2pift)，我们可以说f是信号的频率，如果取其频率域，我们可以在f处看到一个尖峰。如果信号被采样以形成离散信号，我们得到相同的频域。但是在范围[-pi,pi]或[0,2*pi]或[0,N]，你可以将图像视为在两个方向上采样的信号。因此，在X方向和Y方向进行傅里叶变换将为你提供图像的频率表示。

更直观地说，对于正弦信号，如果幅度在短时间内变化得如此之快，则可以说它是高频信号。如果变化缓慢，则是低频信号。你可以将相同的想法扩展到图像。图像中振幅的变化幅度如何？在边缘点，或噪音？所以我们可以说，边缘和噪音是图像中的高频内容，如果振幅没有太大的变化，则它是低频分量。

2、Numpy中的傅里叶变换

首先，我们将看到如何使用Numpy查找傅里叶变换。Numpy有一个FFT软件包来做到这一点。np.fft.fft2()为我们提供了频率变换，这将是一个复杂的数组。它的第一个参数是灰度的输入图像，第二个参数是可选的，它决定输出数组的大小。如果它大于输入图像的大小，则在计算FFT之前，输入图像用0填充。如果它小于输入图像，输入图像将被裁剪。如果没有参数传递，则输出数组大小将与输入相同。现在一旦得到结果，零频率分量(直流分量)将位于左上角。如果要将它置于中心，则需要在两个方向上将结果移动N/2。这只是通过np.fft.fft2（）完成的(分析起来更容易）。一旦你找到频率变换，就可以找到幅度谱。

import cv2 
import numpy as np 
from matplotlib import pyplot as plt 
img = cv2.imread('5.jpg', 0) 
f = np.fft.fft2(img) # 快速傅里叶变换算法得到频率分布 
fshift = np.fft.fftshift(f) # 默认结果中心点位置是在左上角，转移到中间位置 
magnitude_spectrum = 20 * np.log(np.abs(fshift)) # 结果是复数，求绝对值才是振幅 
# 结果展示 
plt.subplot(121),plt.imshow(img, 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(122), plt.imshow(magnitude_spectrum, 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) 
plt.show()

可以在中心看到更多更白的区域，显示低频内容更多.你找到了频率变换，现在你可以在频域中做一些操作，如高通滤波和重构图像，即找到逆DFT。为此，你只需通过使用尺寸为60*60的矩形窗口进行遮盖来移除低频。然后使用np.fft.ifftshit()应用反向偏移，以便DC分量再次出现在左上角。然后使用np.ifft2()函数来找到反FFT。结果也是一个复杂的数字。你可以拿它的绝对值。

import cv2 
import numpy as np 
from matplotlib import pyplot as plt 
img = cv2.imread('5.jpg', 0) 
f = np.fft.fft2(img) # 快速傅里叶变换算法得到频率分布 
fshift = np.fft.fftshift(f) # 默认结果中心点位置是在左上角，转移到中间位置 
magnitude_spectrum = 20 * np.log(np.abs(fshift)) # 结果是复数，求绝对值才是振幅 # 找到频域变换后，可以在频域中做一些操作，如高通滤波和重构图像，即找到逆DFT 
rows, cols = img.shape 
crow, ccol = rows / 2, cols / 2 
fshift[crow - 30:crow + 30, ccol - 30:ccol + 30] = 0 
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, 'gray') plt.title("Input Image"), plt.xticks([]), plt.yticks([]) 
plt.subplot(132),plt.imshow(img_back, 'gray') plt.title('Image After HPF'), plt.xticks([]), plt.yticks([]) 
plt.subplot(133),plt.imshow(img_back) plt.title('Result in JET'), plt.xticks([]), plt.yticks([]) plt.show()

3、OpenCV中的傅里叶变换

OpenCV为此提供了cv2.dft()和cv2.idft().它返回与以前相同的结果，但具有两个通道。第一个通道将具有结果的实际部分，第二个通道将具有结果的虚部。输入图像应首先转换为np.float32.我们将看到如何去做。

现在我们必须做逆DFT。在之前的会话中，我们创建了一个HPF，这次我们将看到如何去除图像中的高频内容，即我们将LPF应用于图像。它实际上模糊了图像，为此，我们首先在低频处创建一个高值(1)的掩膜，即我们通过LF内容，在HF区域传递0.

import cv2 
import numpy as np 
from matplotlib import pyplot as plt 
img = cv2.imread('5.jpg', 0) 
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT) 
dft_shift = np.fft.fftshift(dft) 
magnitude_spectrum = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
rows, cols = img.shape 
crow, ccol = int(rows / 2), int(cols / 2) # 这里还是必须转换为整数哦 
# 首先创建一个掩膜，中心点为1，其余值为0 
mask = np.zeros((rows, cols, 2), np.uint8) 
mask[crow - 30:crow + 30, ccol - 30: ccol + 30] = 1 
# 应用掩膜和逆DFT 
fshift = dft_shift * mask 
f_ishift = np.fft.ifftshift(fshift) 
img_back = cv2.idft(f_ishift) 
img_back = cv2.magnitude(img_back[:, :, 0], img_back[:, :, 1]) 
plt.subplot(121), plt.imshow(img, 'gray') 
plt.title('Input Image'), plt.xticks([]), plt.yticks([]) 
plt.subplot(122), plt.imshow(img_back, 'gray') 
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) 
plt.show()

4、DFT的性能优化

对于某些数组大小，DFT计算的性能更好。当阵列大小是2的幂时，它是最快的。大小为2，3和5的乘积的数组也被相当有效地处理。因此，如果你担心代码的性能，可以在查找DFT之前将数组的大小修改为任何最佳大小（通过填充零）。对于OpenCV，你必须手动填充零。但对于Numpy，你指定了FFT计算的新大小，并且它会自动为你填充零。那么我们如何寻找这个最佳尺寸？OpenCV为此专门提供了一个函数cv2.getOptimalDFTSize().它适用于cv2.dft()和np.fft.fft2()。