【OpenCV-Python】教程：3-11 图像变换（频域变换）

OpenCV Python 图像变换

【目标】

• 利用OpenCV 对图像进行 傅里叶变换
• 利用NumPy的FFT函数
• 傅里叶变换的应用
• cv2.dft(), cv2.idft()

【原理】

傅里叶变换常用于频域图像分析。对于图像来说，2D DFT 常用于寻找频域特征，一个快速算法 FFT（Fast Fourier Transform）用于计算DFT。更详细的资料请查找图像处理或者信号处理和 【参考】。

对于正弦信号来说，$X(t)=A\sin (2\pi ft)$, 我们称$f$ 为信号的频率，如果用了频率，那么可以在 $f$ 获得波峰。我们可以假设图像是一个二维的信号，可以在 $X$ 方向和 $Y$ 方向进行采样。

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

【代码】

• NumPy 中的傅里叶变换

我们将了解如何用 NumPy查找傅里叶变换。NumPy中有一个FFT包来实现这一点，np.fft.fft2()给我们提供了一个复数数组，第一个参数是灰度图像，第二个参数是可选的尺寸，如果大于输入，输入填充，如果小于输入，输入裁剪，默认相同。
零频率分量（DC分量）将位于左上角，如果想让它居中，那么就在两个方向上分别移动 $N/2$，可以通过函数 np.fft.fftshift()完成。
可以看到有很多的白色区域在中心，也就是说有很多低频分量。

# NumPy 的 FFT
import cv2
import numpy as np
from matplotlib import pyplot as plt

# 读入图像
img = cv2.imread('assets/messi5.jpg', 0)

# fft
f = np.fft.fft2(img)

# 移动到中心
fshift = np.fft.fftshift(f)

# 幅度
magnitude_spectrum = 20*np.log(np.abs(fshift))

# 显示
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

• NumPy 中的傅里叶变换应用（高通滤波）

# NumPy 的 FFT
import cv2
import numpy as np
from matplotlib import pyplot as plt

# 读入图像
img = cv2.imread('assets/messi5.jpg', 0)

# fft
f = np.fft.fft2(img)

# 移动到中心
fshift = np.fft.fftshift(f)

# 低频信号较多部分置0
rows, cols = img.shape
crow,ccol = rows//2 , cols//2
fshift[crow-30:crow+31, ccol-30:ccol+31] = 0

# 逆移动
f_ishift = np.fft.ifftshift(fshift)

# 逆变换
img_back = np.fft.ifft2(f_ishift)

# img_back = np.real(img_back) # 官网此处有误
img_back = np.abs(img_back)


# 显示
plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = '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()

如果仔细看结果，尤其是JET图，可以发现一些人工痕迹，这种效果称为“振铃效应”，主要是由矩形窗口的滤波导致的。所以矩形窗口不能用于滤波，最好的选项是高斯窗口。

• OpenCV 的 dft

# OpenCV dft
import numpy as np
import cv2
from matplotlib import pyplot as plt

# 读入图像
img = cv2.imread('assets/messi5.jpg', 0)

# 进行 dft 运算
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]))

# 显示
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

• OpenCV 的 dft 应用（低通滤波）

# OpenCV dft
import numpy as np
import cv2
from matplotlib import pyplot as plt

# 读入图像
img = cv2.imread('assets/messi5.jpg', 0)

# 进行 dft 运算
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)

# 移动到中心
dft_shift = np.fft.fftshift(dft)

rows, cols = img.shape
# crow,ccol = rows/2 , cols/2 # 官网此处代码也有误, 宽高除2不一定是整数
crow, ccol = rows//2, cols//2

# 建立一个mask
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow-30:crow+31, ccol-30:ccol+31] = 1

# 高频信号置零
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(131), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_back, cmap='gray')
plt.title('Image After LPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()

• DFT 的优化

DFT 计算对某些维度的数组效果更好，当尺寸为2的幂时，速度最快。当尺寸为2、3、5的乘积依然处理更有效。所以，可以在计算DFT之前获得比较好的尺寸（图像可以补零）进行DFT计算。对于OpenCV来说，需要手动补，对于NumPy来说，自动补零。

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

# 读入图像
img = cv2.imread('assets/messi5.jpg', 0)

rows, cols = img.shape
print("{} {}".format(rows, cols))
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print("{} {}".format(nrows, ncols))

342 548
360 576

【接口】

• dft

cv2.dft(	src[, dst[, flags[, nonzeroRows]]]	) ->	dst

对浮点数组执行正和逆离散傅里叶变换

• src: 输入数组，实数或复数
• dst: 输出数组，尺寸和类型取决于flag
• flags: 参数，变换的标识
• nonzeroRows: 当参数不为零时，该函数假设只有输入阵列的第一个非零行（未设置DFT_INVERSE）或只有输出阵列的第一非零行（设置DFT_REVERSE）包含非零，因此，该函数可以更有效地处理其余行并节省一些时间；该技术对于使用DFT计算阵列互相关或卷积非常有用。

• dft flags
  

• idft

cv2.dft(	src[, dst[, flags[, nonzeroRows]]]	) ->	dst

对浮点数组执行逆离散傅里叶变换

• src: 输入数组，实数或复数
• dst: 输出数组，尺寸和类型取决于flag
• flags: 参数，变换的标识
• nonzeroRows: 当参数不为零时，该函数假设只有输入阵列的第一个非零行（未设置DFT_INVERSE）或只有输出阵列的第一非零行（设置DFT_REVERSE）包含非零，因此，该函数可以更有效地处理其余行并节省一些时间；该技术对于使用DFT计算阵列互相关或卷积非常有用。

• magnitude

cv2.magnitude(	x, y[, magnitude]	) ->	magnitude

$$dst(I)=\sqrt{x(I{)}^2+y(I{)}^2}$$

计算 幅度

• x: x坐标的浮点数组
• y: y坐标的浮点数组
• magnitude: 输出的幅度

【参考】

1. OpenCV官方文档
2. An Intuitive Explanation of Fourier Theory by Steven Lehar
3. Fourier Transform at HIPR
4. What does frequency domain denote in case of images?