OpenCV Tutorials 10 - 傅里叶变换

Detection
1,274 阅读10分钟

傅里叶变换

本节主要介绍：

  • OpenCV中的傅里叶变换函数

  • 使用Numpy中的FFT(快速傅里叶变换)函数

  • 傅里叶变换的应用

  • 如下函数cv.dft(), cv.idft()

  • 用法如下：

  1. cv.dft( src[, dst[, flags[, nonzeroRows]]] ) -> dst
  2. cv.idft( src[, dst[, flags[, nonzeroRows]]] ) -> dst

一、概念

傅里叶变换用于分析各种滤波器的频率特性。 对于图像，二维离散傅里叶变换 (DFT) 用于查找频域。 一种称为快速傅里叶变换 (FFT) 的快速算法用于计算 DFT。 有关这些的详细信息可以在任何图像处理或信号处理教科书中找到。 请参阅其他资源部分。

对于正弦信号，x(t)=Asin(2πft)，我们可以说 f 是信号的频率，如果取其频域，我们可以看到 f 处的尖峰。 如果对信号进行采样以形成离散信号，我们会得到相同的频域，但在 [−π,π] 或 [0,2π] 范围内是周期性的（对于 N 点 DFT，则为 [0,N]）。 您可以将图像视为在两个方向上采样的信号(二维采样)。 因此，在 X 和 Y 方向上进行傅立叶变换可以为您提供图像的频率表示。

更直观地说，对于正弦信号，如果幅度在短时间内变化很快，那么可以称其为高频信号。相反， 如果它变化缓慢，则是一个低频信号。 您可以将相同的想法扩展到图像。 图像中的幅度在哪里发生剧烈变化？在边缘点，或噪音。 所以我们可以说，边缘和噪声是图像中的高频内容。 如果幅度没有太大变化，则为低频分量。

二、Numpy中的傅里叶变换

首先，我们将了解如何使用 Numpy 找到傅立叶变换。 Numpy 有一个 FFT 包可以做到这一点。 np.fft.fft2() 为我们提供了频率变换，这将是一个复杂的数组。 它的第一个参数是输入图像，它是灰度图。 第二个参数是可选的，它决定了输出数组的大小。 如果它大于输入图像的大小，则在计算 FFT 之前用零填充输入图像。 如果小于输入图像，输入图像将被裁剪。 如果未传递任何参数，则输出数组大小将与输入相同。

现在，一旦您得到结果，零频率分量（直流分量）将位于左上角。 如果要将其置于中心，则需要将结果在两个方向上移动 N/2。 这只是通过函数 np.fft.fftshift() 完成的。 （更容易分析）。 找到频率变换后，就可以找到幅度谱。

import cv2 as cv
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

def cv_show(name, img):
    cv.imshow(name, img)
    cv.waitKey(0)
    cv.destroyAllWindows()

def compare(imgs):
  #  for i in range(len(imgs)):
 #       imgs[i][:,-3:-1,:] = [255,255,255]
    res = np.hstack(imgs)
    cv_show('Compare', res)

1. 傅里叶变换

# 读入灰度图
img = cv.imread('lena.png',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
# 将结果的0频率点移动至图像中心
fshift = np.fft.fftshift(f)
# 对频谱谱进行放大处理，np.log进行对数处理
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()

output_13_0.png

从上面可以看出，中心的区域有更多更白的部分，表明低频内容更多。

得到频率变换之后，可以在频域做一些操作，比如高通滤波和重建图像(找到逆DFT)。 为此，您只需通过使用大小为 60x60 的矩形窗口进行掩蔽来去除低频。 然后使用 np.fft.ifftshift() 应用反向移位，使直流分量((0频率点)再次出现在左上角。 然后使用 np.ifft2() 函数找到逆 FFT。 结果再次将是一个复数。 你可以取它的绝对值。

  • 不能直接使用空域变换对频域结果进行处理

2. 使用高通滤波对图像进行频域处理

row,col = img.shape
fshift[row//2-50:row//2+50, col//2-50:col//2+50] = 0
magnitude_spectrum = 20*np.log(np.abs(fshift))
ffshift = np.fft.ifftshift(fshift)
f_img = np.fft.ifft2(ffshift)
res = np.real(f_img)

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

cv_show('Res', res)

F:\Temp1/ipykernel_8564/3936633447.py:3: RuntimeWarning: divide by zero encountered in log
  magnitude_spectrum = 20*np.log(np.abs(fshift))
F:\Temp1/ipykernel_8564/3936633447.py:8: RuntimeWarning: divide by zero encountered in log
  magnitude_spectrum = 20*np.log(np.abs(fshift))

output_18_1.png

傅里叶反变换结果.PNG

  • np.real作用：将傅里叶反变换的复数结果直接取实部,如下所示：
a = np.array([1+2j, 3+4j, 5+6j])
a.real
array([1.,  3.,  5.])
a.real = 9
a
array([9.+2.j,  9.+4.j,  9.+6.j])
a.real = np.array([9, 8, 7])
a
array([9.+2.j,  8.+4.j,  7.+6.j])
np.real(1 + 1j)
1.0

结果表明高通滤波是一种边缘检测操作。 这就是我们在图像梯度章节中看到的。 这也表明大部分图像数据存在于频谱的低频区域。 无论如何，我们已经看到了如何在 Numpy 中找到 DFT、IDFT 等。 现在让我们看看如何在 OpenCV 中做到这一点。

三、OpenCV中的傅里叶变换

1. 傅里叶变换

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

img = cv.imread('lena.png',0)
# 第二个参数是操作码，表明输出的是傅里叶变换的复矩阵
dft = cv.dft(np.float32(img),flags = cv.DFT_COMPLEX_OUTPUT)
# 将0频移动至图像中心
dft_shift = np.fft.fftshift(dft)
# 0维是实部， 1维是虚部，magnitude计算实部虚部的欧氏长度
magnitude_spectrum = 20*np.log(cv.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()

output_26_0.png

  • cv.dft 的 flags 是操作码：
  1. 如果设置了 DFT_ROWS 或输入数组具有单行或单列，则该函数对矩阵的每一行执行一维正向或逆变换。否则，它将执行 2D 变换。
  2. 如果输入数组为实数且未设置 DFT_INVERSE，则该函数执行一维或二维正向变换：
    1. 当设置 DFT_COMPLEX_OUTPUT 时，输出是与输入大小相同的复矩阵。
    2. 未设置 DFT_COMPLEX_OUTPUT 时，输出是与输入大小相同的实矩阵。在 2D 变换的情况下，它使用如上所示的打包格式。在单个 1D 变换的情况下，它看起来像上面矩阵的第一行。在多个 1D 变换的情况下（使用 DFT_ROWS 标志时），输出矩阵的每一行看起来像上面矩阵的第一行。
  3. 如果输入数组是复数并且未设置 DFT_INVERSE 或 DFT_REAL_OUTPUT，则输出是与输入大小相同的复数数组。该函数根据标志 DFT_INVERSE 和 DFT_ROWS 独立地对整个输入数组或输入数组的每一行执行正向或反向 1D 或 2D 变换。
  4. 当设置了 DFT_INVERSE 并且输入数组是实数，或者它是复数但设置了 DFT_REAL_OUTPUT 时，输出是与输入大小相同的实数数组。该函数根据标志 DFT_INVERSE 和 DFT_ROWS 对整个输入数组或每个单独的行执行 1D 或 2D 逆变换。如果设置了 DFT_SCALE，则在转换后进行缩放。
  • 二维矩阵的打包格式：

download.png

一般我们都会设置 DFT_COMPLEX_OUTPUT，输出复矩阵

  • 注意： 您还可以使用 cv.cartToPolar() 一次返回振幅和相位

  • 用法如下：cv.cartToPolar( x, y[, magnitude[, angle[, angleInDegrees]]] ) -> magnitude, angle

所以，现在我们必须做逆 DFT。 在上一节中，我们创建了一个 HPF，这次我们将看到如何去除图像中的高频内容，即我们将 LPF 应用于图像。 它实际上使图像模糊。 为此，我们首先在低频创建一个具有高值 (1) 的掩码，即我们通过 LF 内容，并且在 HF 区域传入 0。

2. 傅里叶逆变换

rows, cols = img.shape
# 行列必须保持为整数，使用地板除
crow,ccol = rows//2 , 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 = cv.idft(f_ishift)
# 将复矩阵转化为振幅矩阵
img_back = cv.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()

download.png

  • 注意：像往常一样，OpenCV 函数 cv.dft() 和 cv.idft() 比 Numpy 更快。 但是 Numpy 的函数更加人性化。

3. 使用OpenCV进行高通滤波操作

四、 DFT的性能优化

img = cv.imread('xy.png', 0)
row, col = img.shape
img = np.float32(img)
# 获取傅里叶变换复矩阵
f = cv.dft(img, flags = cv.DFT_COMPLEX_OUTPUT)
# 将 0 频域点转至图像中间
ffshift = np.fft.fftshift(f)

# 转化为幅度矩阵
magnitude =20 *  np.log(cv.magnitude(f[:,:,0], f[:,:,1]))


# 定义mask,mask要定义为两层，因为要对复矩阵的两层进行过滤
mask = np.ones((row, col, 2), np.uint8)
mask[row//2 - 100: row//2 + 100,col//2 - 100: col//2 + 100,:] = 0
# 掩码处理
iffshift = ffshift * mask
# log函数最好不要传入0,要加上一个较小的偏移量
magnitude_iffshift = 20*np.log(cv.magnitude(iffshift[:,:,0], iffshift[:,:,1]) + 1e-6)
# 在逆变换之前，需要将0频域点还原到左上角,否则会得到倒立图像
iffshift = np.fft.ifftshift(iffshift)

# 逆变换，得到复矩阵
iff = cv.idft(iffshift)
# 转化为振幅矩阵
iff = cv.magnitude(iff[:,:,0], iff[:,:,1])

# 注册整幅画布的大小
fig = plt.figure(figsize = (20, 8))
plt.subplot(131)
plt.imshow(magnitude, cmap = 'gray')
plt.title('Origin ')
plt.xticks([]),plt.yticks([])

plt.subplot(132)
plt.imshow(magnitude_iffshift, cmap = 'gray')
plt.title('Transforms')
plt.xticks([]),plt.yticks([])

plt.subplot(133)
plt.imshow(iff, cmap = 'gray')
plt.title('Result JET')
plt.xticks([]),plt.yticks([])
plt.show()

download.png

对于某些数组大小，DFT 计算的性能更好。 数组大小为 2 的幂(2^k)时速度最快。 大小为 2、3 和 5 的乘积的数组也得到了相当有效的处理。 因此，如果您担心代码的性能，可以在查找 DFT 之前将数组的大小修改为任何最佳大小（通过填充零）。 对于 OpenCV，您必须手动填充零。 但是对于 Numpy，您指定 FFT 计算的新大小，它会自动为您填充零。

那么我们如何找到这个最佳尺寸呢？ OpenCV 为此提供了一个函数 cv.getOptimalDFTSize()。 它适用于 cv.dft() 和 np.fft.fft2()。

需要传入tuple或者list里面是图像的尺寸

  • 用法如下：cv.getOptimalDFTSize( vecsize ) -> retval
img = cv.imread('xy.png',0)
img.shape

(673, 684)

row, col = img.shape
# 转化为整数
row = int(row)
col = int(col)
row_opt =  cv.getOptimalDFTSize(row)
col_opt =  cv.getOptimalDFTSize(col)
(row_opt, col_opt)

(675, 720)

可以看出，大小 (673, 684) 修改为 (675, 720)。 现在让我们用零填充它（对于 OpenCV）并找到它们的 DFT 计算性能。 您可以通过创建一个新的大零数组并将数据复制到其中，或使用 cv.copyMakeBorder()来实现。

  • 用法如下： cv.copyMakeBorder( src, top, bottom, left, right, borderType[, dst[, value]] ) -> dst
# 直接创建一个大0阵列，将原阵列复制到左上角
nimg = np.zeros((row_opt, col_opt))
nimg[:row,:col] = img
cv_show('test', nimg)

right = col_opt - col
bottom = row_opt - row
bordertype = cv.BORDER_CONSTANT #just to avoid line breakup in PDF file
nimg = cv.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)
cv_show('test', nimg)

下面来测试二者的Numpy傅里叶变换效率

%timeit fft1 = np.fft.fft2(img)

34.6 ms ± 1.06 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit fft2 = np.fft.fft2(img,[row_opt, col_opt])

18.3 ms ± 238 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

接下来测试OpenCV傅里叶变换效率

# 使用OpenCV中的dft函数时：第一个参数类型要为数据类型为 np.float32 的图像阵列，第二个参数要指定函数返回矩阵类型
%timeit dft1= cv.dft(np.float32(img),flags=cv.DFT_COMPLEX_OUTPUT)

6.38 ms ± 285 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit dft2= cv.dft(np.float32(nimg),flags=cv.DFT_COMPLEX_OUTPUT)

3.64 ms ± 86.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

五、为什么Laplacian是一种高通滤波？

一个论坛上也有人问过类似的问题。 问题是，为什么拉普拉斯算子是高通滤波器？ 为什么 Sobel 是 HPF？ 等等。给出的第一个答案是傅立叶变换。 只需对更大尺寸的 FFT 进行拉普拉斯算子的傅立叶变换即可。 分析一下：

# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))
# creating a gaussian filter
x = cv.getGaussianKernel(5,10)
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()

download.png

从图像中，您可以看到每个内核阻塞的频率区域，以及它通过的区域。 从这些信息中，我们可以说明为什么每个内核都是 HPF 或 LPF

avatar
文章
阅读
粉丝