OpenCV-Python学习笔记（十四）：图像变换：傅里叶变换与反变换

最新推荐文章于 2024-09-18 21:42:55 发布

GiffordY

最新推荐文章于 2024-09-18 21:42:55 发布

阅读量1w

点赞数 2

分类专栏： OpenCV 机器视觉文章标签： OpenCV 机器视觉 Python

OpenCV 同时被 2 个专栏收录

21 篇文章 10 订阅

订阅专栏

机器视觉

20 篇文章 1 订阅

订阅专栏

OpenCV-Python官方文档关于图像傅里叶变换和反变换的教程网址：

https://docs.opencv.org/4.1.0/de/dbc/tutorial_py_fourier_transform.html

目标
我们将要学习：
• 使用 OpenCV 对图像进行傅里叶变换（DFT）：cv2.dft()，cv2.idft()
• 使用 Numpy 中 FFT（快速傅里叶变换）函数：
• 傅里叶变换的一些用处
• 我们将要学习的函数有： cv2.dft()，cv2.idft() 等

原理

把图像想象成沿着两个方向采集的信号。所以对图像同时进行 X 方向和 Y 方向的傅里叶变换，我们就会得到这幅图像的频域表示（频谱图）。更直观一点，对于一个正弦信号，如果它的幅值变化非常快，我们可以说它是高频信号，如果变化非常慢，我们称之为低频信号。把这种想法应用到图像中，图像哪里的幅度变化非常大呢？边界点或者噪声。所以我们说边界和噪声是图像中的高频分量。

1.1 Numpy 中的傅里叶变换

Numpy 中的 FFT 包可以帮助我们实现快速傅里叶变换。函数 np.fft.fft2() 可以对信号进行频率转换，输出结果是一个复杂的组。本函数的第一个参数是输入图像，要求是灰度格式。第二个参数是可选的, 决定输出数组的大小。输出数组的大小和输入图像大小一样。如果输出结果比输入图像大，输入图像就需要在进行 FFT 前补0。如果输出结果比输入图像小的话，输入图像就会被切割。

现在我们得到了结果，频率为 0 的部分（直流分量）在输出图像的左上角。如果想让它（直流分量）在输出图像的中心，我们还需要将结果沿两个方向平移 N/2 。函数 np.fft.fftshift() 可以帮助我们实现这一步。（这样更容易分析）。进行完频率变换之后，我们就可以构建振幅谱了。

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

img = cv2.imread('messi5.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, 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()

结果：

然后就可以在频域对图像进行操作了，比如进行高通滤波，再傅里叶反变换回空间域。比如我们可以使用一个60x60 的矩形窗口对图像进行掩模操作从而去除低频分量。然后再使用函数np.fft.ifftshift() 进行逆平移操作，所以现在直流分量又回到左上角了，然后使用函数 np.ifft2() 进行 FFT 逆变换。同样又得到一堆复杂的数字，我们可以对他们取绝对值：

rows, cols = img.shape
crow,ccol = int(rows/2) , int(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, 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 颜色的图像，你会看到一些不自然的东西（如我用红色箭头标出的区域）。看上图那里有些条带装的结构，这被成为振铃效应。这是由于我们使用矩形窗口做掩模造成的。这个掩模被转换成正弦形状时就会出现这个问题。所以一般我们不使用矩形窗口滤波。最好的选择是高斯窗口。

1.2 OpenCV 中的傅里叶变换

OpenCV 中相应的函数是 cv2.dft() 和 cv2.idft()。和前面输出的结果一样，但是是双通道的。第一个通道是结果的实数部分，第二个通道是结果的虚数部分。输入图像要首先转换成 np.float32 格式。

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

img = cv2.imread('messi5.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]))

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()

注意：你可以使用函数 cv2.cartToPolar()，它会同时返回幅度和相位。

现在我们来做逆 DFT。在前面的部分我们实现了一个 HPF（高通滤波），现在我们来做 LPF（低通滤波）将高频部分去除。其实就是对图像进行模糊操作。首先我们需要构建一个掩模，与低频区域对应的地方设置为 1, 与高频区域对应的地方设置为 0。

rows, cols = img.shape
crow,ccol = int(rows/2) , int(cols/2)

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse 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, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

结果：

注意： OpenCV 中的函数 cv2.dft() 和 cv2.idft() 要比 Numpy 快。但是Numpy 函数更加用户友好。关于性能的描述，请看下面的小节。

1.3 DFT 的性能优化

当数组的大小为某些值时 DFT 的性能会更好。当数组的大小是 2 的指数时 DFT 效率最高。当数组的大小是 2， 3， 5 的倍数时效率也会很高。所以如果你想提高代码的运行效率时，你可以修改输入图像的大小（补 0）。对于OpenCV 你必须自己手动补0。但是 Numpy，你只需要指定 FFT 运算的大小，它会自动补 0。

那我们怎样确定最佳大小呢？OpenCV 提供了一个函数 cv2.getOptimalDFTSize()。它可以同时被 cv2.dft() 和 np.fft.fft2() 使用。让我们一起使用 IPython的魔法命令%timeit 来测试一下吧。

In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548

In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576

看到了吧，数组的大小从（342， 548）变成了（360， 576）。现在我们为它补 0，然后看看性能有没有提升。你可以创建一个大的 0 数组，然后把我们的数据拷贝过去，或者使用函数 cv2.copyMakeBoder()。

nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

或者：

right = ncols - cols
bottom = nrows - rows
#just to avoid line breakup in PDF file
bordertype = cv2.BORDER_CONSTANT
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

现在我们看看 Numpy 的表现：

In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop

速度提高了 4 倍。我们再看看 OpenCV 的表现：

In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

也提高了 4 倍，同时我们也会发现 OpenCV 的速度是 Numpy 的 3 倍。你也可以测试一下逆 FFT 的表现。

1.4 为什么拉普拉斯算子是高通滤波器？

在论坛中遇到了一个类似的问题。为什么拉普拉斯算子是高通滤波器？为什么 Sobel 是 HPF？等等。对于第一个问题的答案我们以傅里叶变换的形式给出。我们一起来对不同的算子进行傅里叶变换并分析它们：

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

# simple averaging filter without scaling parameter
mean_filter = np.ones((3, 3))
# creating a guassian filter
x = cv2.getGaussianKernel(5, 10)
# x.T 为矩阵转置
gaussian = x * x.T

# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
                   [-10, 0, 10],
                   [-3, 0, 3]])
# sobel in x direction
sobel_x = np.array([[-1, 0, 1],
                    [-2, 0, 2],
                    [-1, 0, 1]])
# sobel in y direction
sobel_y = np.array([[-1, -2, -1],
                    [0, 0, 0],
                    [1, 2, 1]])
# laplacian
laplacian = np.array([[0, 1, 0],
                      [1, -4, 1],
                      [0, 1, 0]])

filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian', 'laplacian', 'sobel_x', \
               'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z) + 1) for z in fft_shift]

for i in range(6):
    plt.subplot(2, 3, i + 1), plt.imshow(mag_spectrum[i], cmap='gray')
    plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])
plt.show()

结果：

从图像中我们就可以看出每一个算子允许通过那些信号。从这些信息中我们就可以知道那些是 HPF 那是 LPF。

更多资源：

An Intuitive Explanation of Fourier Theory by Steven Lehar
Fourier Transform at HIPR
What does frequency domain denote in case of images?

另：为了更清楚的观察振铃效应，换用OpenCV的logo图像做测试，如下：

"""
使用 numpy 进行 FFT 和 IFFT
"""

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

img = cv2.imread('opencv_logo.jpg', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20 * np.log(np.abs(fshift))

rows, cols = img.shape
crow, ccol = int(rows / 2), int(cols / 2)

fshift[crow - 20: crow + 20, ccol - 20: ccol + 20] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

"""
使用 OpenCV 进行 FFT 和 IFFT
"""
# cv2.dft( )函数返回值是双通道的，第一个通道是实数部分，第二个是虚数部分
# cv2.dft( )函数输入图像要先转为 np.float32 格式
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
magnitude_spectrum2 = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))

# 创建低通滤波器掩模、使用掩模滤波、IDFT
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow - 30:crow + 30, ccol - 30:ccol + 30] = 1

f_dft_shift = dft_shift * mask
f_idft_shift = np.fft.ifftshift(f_dft_shift)

img_back2 = cv2.idft(f_idft_shift)
img_back2 = cv2.magnitude(img_back2[:, :, 0], img_back2[:, :, 1])

# 显示结果图像
imgList = [img, magnitude_spectrum, img_back, img, magnitude_spectrum2, img_back2]
imgName = ['img', 'magnitude_spectrum', 'img_back', 'img', 'magnitude_spectrum2', 'img_back2']

for i in range(6):
    plt.subplot(2, 3, i + 1), plt.imshow(imgList[i], cmap='gray'), plt.title(imgName[i])
    plt.xticks([]), plt.yticks([])
plt.show()

结果如下：