核心功能：离散傅立叶变换 OpenCV v4.8.0

上一个教程 ： 改变图像的对比度和亮度

下一个教程 ： 使用 XML 和 YAML 文件输入和输出文件

原作者	Bernát Gábor
兼容性	OpenCV >= 3.0

目标

我们将寻求以下问题的答案：

• 什么是傅立叶变换，为什么要使用它？
• 如何在 OpenCV 中进行傅立叶变换？
• 函数的使用，如：copyMakeBorder() 、merge() 、dft() 、getOptimalDFTSize() 、log() 和 normalize()。

源代码

• C++
  您可以从此处下载，或在 OpenCV 源代码库的 samples/cpp/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.cpp 中找到它。

下面是 dft() 的使用示例：

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include <iostream>
using namespace cv.NET
using namespace std；
static void help(char ** argv)
{
 cout << endl
 << "This program demonstrated the use of the discrete Fourier transform (DFT). " << endl
 << "The dft of an image is taken and it's power spectrum is displayed." << endl << endl
 << "Usage:" << endl
 << argv[0] << " [image_name -- 默认为 lena.jpg]" << endl << endl；
}
int main(int argc, char ** argv)
{
 help(argv)；
 const char* filename = argc >=2 ? argv[1] : "lena.jpg"；
 Mat I = imread( samples::findFile( filename ), IMREAD_GRAYSCALE)；
 if( I.empty()){
 cout << "Error opening image" << endl；
 return EXIT_FAILURE；
 }
 Mat padded; // 将输入图像扩展到最佳尺寸
 int m = getOptimalDFTSize( I.rows )；
 int n = getOptimalDFTSize( I.cols ); // 在边框上添加零值
 copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0))；
 Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}；
 Mat complexI；
 merge(planes, 2, complexI); // 在扩展后的平面上再添加一个零平面
 dft(complexI, complexI); // 这样得到的结果才可能适合源矩阵
 // 计算幅值并转换为对数标度
 // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
 split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
 magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
 Mat magI = planes[0]；
 magI += Scalar::all(1); // 切换到对数标度
 log(magI, magI)；
 // 如果光谱的行数或列数为奇数，则裁剪光谱
 magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2))；
 // 重新排列傅立叶图像的象限，使原点位于图像中心
 int cx = magI.cols/2；
 int cy = magI.rows/2；
 Mat q0(magI, Rect(0, 0, cx, cy)); // 左上角--每个象限创建一个 ROI
 Mat q1(magI, Rect(cx, 0, cx, cy)); // 右上角
 Mat q2(magI, Rect(0, cy, cx, cy)); // 左下角
 Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
 Mat tmp; // 交换象限（左上角和右下角）
 q0.copyTo(tmp)；
 q3.copyTo(q0)；
 tmp.copyTo(q3)；
 q1.copyTo(tmp); // 交换象限（右上角和左下角）
 q2.copyTo(q1)；
 tmp.copyTo(q2)；
 normalize(magI,magI,0,1,NORM_MINMAX);//将带有浮点数值的矩阵转换为可查看的图像形式（浮点数值介于 0 和 1 之间）。
 // 可查看的图像形式（数值在 0 和 1 之间的浮点）。
 imshow("Input Image" , I ); // 显示结果
 imshow("spectrum magnitude", magI)；
 waitKey()；
 return EXIT_SUCCESS；
}

• Java
  您可以从此处下载，或在 OpenCV 源代码库的 samples/java/tutorial_code/core/discrete_fourier_transform/DiscreteFourierTransform.java中找到它。

下面是 dft() 的使用示例：

import org.opencv.core.*；
import org.opencv.highgui.HighGui；
import org.opencv.imgcodecs.Imgcodecs；
import java.util.List；
import java.util.*；
 class DiscreteFourierTransformRun{
 private void help() {
 System.out.println("" +
 "This program demonstrated the use of the discrete Fourier transform (DFT). \n" +
 "The dft of an image is taken and it's power spectrum is displayed.\n" +
 "Usage:\n" +
 "./DiscreteFourierTransform [image_name -- 默认为 ../data/lena.jpg]"）；
 }
 public void run(String[] args){
 help()；
 String filename = ((args.length > 0) ? args[0] : "../data/lena.jpg")；
 Mat I = Imgcodecs.imread(filename, Imgcodecs.IMREAD_GRAYSCALE)；
 if( I.empty() ) {
 System.out.println("Error opening image")；
 System.exit(-1)；
 }
 Mat padded = new Mat(); // 将输入图像扩展到最佳尺寸
 int m = Core.getOptimalDFTSize( I.rows() )；
 int n = Core.getOptimalDFTSize( I.cols() ); // 在边框上添加零值
 Core.copyMakeBorder(I, padded, 0, m - I.rows(), 0, n - I.cols(), Core.BORDER_CONSTANT, Scalar.all(0))；
 List<Mat> planes = new ArrayList<Mat>()；
 PADEDED.CONVERT TO(PADEDED, CvType.CV_32F)；
 planes.add(padded)；
 planes.add(Mat.zeros(padded.size(), CvType.CV_32F))；
 Mat complexI = new Mat()；
 Core.merge(planes, complexI); // 在扩展后的平面中添加另一个零平面
 Core.dft(complexI, complexI); // 这样得到的结果才可能适合源矩阵
 // 计算幅值并转换为对数标度
 // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
 Core.split(complexI, planes); // planes.get(0) = Re(DFT(I)
 // planes.get(1) = Im(DFT(I))
 Core.magnitude(planes.get(0), planes.get(1), planes.get(0));// planes.get(0) = magnitude
 Mat magI = planes.get(0)；
 Mat matOfOnes = Mat.ones(magI.size(), magI.type())；
 Core.add(matOfOnes, magI, magI); // 切换到对数标度
 Core.log(magI, magI)；
 // 如果光谱的行数或列数为奇数，则裁剪光谱
 magI = magI.submat(new Rect(0, 0, magI.cols() & -2, magI.rows() & -2))；
 // 重新排列傅立叶图像的象限，使原点位于图像中心
 int cx = magI.cols()/2；
 int cy = magI.rows()/2；
 Mat q0 = new Mat(magI, new Rect(0, 0, cx, cy)); // Top-Left - 每个象限创建一个 ROI
 Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy)); // 右上角
 Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy)); // 左下角
 Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // Bottom-Right
 Mat tmp = new Mat(); // 交换象限（左上角和右下角）
 q0.copyTo(tmp)；
 q3.copyTo(q0)；
 tmp.copyTo(q3)；
 q1.copyTo(tmp); // 交换象限（右上角和左下角）
 q2.copyTo(q1)；
 tmp.copyTo(q2)；
 magI.convertTo(magI, CvType.CV_8UC1)；
 Core.normalize(magI,magI,0,255,Core.NORM_MINMAX,CvType.CV_8UC1);//将带有浮点数值的矩阵转换为可视图像形式（CV_8UC1
 // 转换为可查看的图像形式（浮点值介于
 // 值 0 和 255 之间的浮点）。
 HighGui.imshow("Input Image" , I ); // 显示结果
 HighGui.imshow("Spectrum Magnitude", magI )；
 HighGui.waitKey()；
 System.exit(0)；
 }
}
public class DiscreteFourierTransform {
 public static void main(String[] args) { 
 // 加载本地库。
 System.loadLibrary(Core.NATIVE_LIBRARY_NAME)；
 new DiscreteFourierTransformRun().run(args)；
 }
}

• Python
  您可以从此处下载，或在 OpenCV 源代码库的 samples/python/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.py 中找到它。

下面是 dft() 的用法示例：

from __future__ import print_function
import sys
import cv2 as cv
import numpy as np
def print_help()：
 print('''
 This program demonstrated the use of the discrete Fourier transform (DFT).
 The dft of an image is taken and it's power spectrum is displayed.
 Usage:
 discrete_fourier_transform.py [image_name -- default lena.jpg]''')
def main(argv)：
 print_help()
 filename = argv[0] if len(argv) > 0 else 'lena.jpg'.
 I = cv.imread(cv.samples.findFile(filename), cv.IMREAD_GRAYSCALE)
 if I 为空：
 print('Error opening image')
 return -1
 
 rows,cols = I.shape
 m = cv.getOptimalDFTSize( rows )
 n = cv.getOptimalDFTSize( cols )
 padded = cv.copyMakeBorder(I, 0, m - rows, 0, n - cols, cv.BORDER_CONSTANT, value=[0, 0, 0])
 
 planes = [np.float32(padded), np.zeros(padded.shape, np.float32)] 
 complexI = cv.merge(planes) # 在扩展后的平面上再添加一个零平面
 
 cv.dft(complexI, complexI) # 这样，结果就可能适合源矩阵了
 
 cv.split(complexI, planes) # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
 cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = magnitude
 magI = planes[0]
 
 matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
 cv.add(matOfOnes, magI, magI) # 切换到对数标度
 cv.log(magI, magI)
 
 magI_rows, magI_cols = magI.shape
 # 如果光谱的行数或列数为奇数，则裁剪光谱
 magI = magI[0:(magI_rows & -2)，0:(magI_cols & -2)] cx = int(magI_rows & -2)
 cx = int(magI_rows/2)
 cy = int(magI_cols/2)
 q0 = magI[0:cx, 0:cy] # 左上角 - 每个象限创建一个 ROI
 q1 = magI[cx:cx+cx, 0:cy] # 右上角
 q2 = magI[0:cx, cy:cy+cy] # 左下角
 q3 = magI[cx:cx+cx, cy:cy+cy] # 右下方
 tmp = np.copy(q0) # 交换象限（左上角和右下角）
 magI[0:cx, 0:cy] = q3
 magI[cx:cx + cx, cy:cy + cy] = tmp
 tmp = np.copy(q1) # 交换象限（从右上角到左下角）
 magI[cx:cx + cx, 0:cy] = q2
 magI[0:cx, cy:cy + cy] = tmp
 
 cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # 将包含浮点数值的矩阵转换成一个
 
 cv.imshow("Input Image" , I ) # 显示结果
 cv.imshow("spectrum magnitude", magI)
 cv.waitKey()
if __name__ == "__main__"：
 main(sys.argv[1:])

解释

傅立叶变换将图像分解为正弦和余弦成分。换句话说，它将图像从空间域转换到频率域。其原理是，任何函数都可以用无限个正弦和余弦函数之和来近似。傅立叶变换就是实现这一目的的一种方法。从数学上讲，二维图像的傅里叶变换是：

$$F\left(k,l\right)=\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}f\left(i,j\right)e^{-i2\pi \left(\frac{ki}{N}+\frac{lj}{N}\right)}$$
$$e^{ix}=\cos x+i\sin x$$
其中，f 是空间域的图像值，F 是频率域的图像值。变换的结果是复数。可以通过实数图像和复数图像，或通过幅值图像和相位图像来显示。不过，在整个图像处理算法中，只有幅值图像才是有意义的，因为它包含了我们所需的有关图像几何结构的所有信息。不过，如果您打算以这些形式对图像进行一些修改，然后需要对其进行重新变换，则需要保留这两种形式的图像。

在本示例中，我将演示如何计算和显示傅立叶变换的幅值图像。数字图像是离散的。这意味着它们可以从给定的域值中取值。例如，在基本灰度图像中，数值通常介于 0 和 255 之间。因此，傅立叶变换也需要是离散类型的，从而产生离散傅立叶变换 (DFT)。当您需要从几何角度确定图像结构时，就会用到这种方法。以下是操作步骤（在输入图像 I 为灰度图像的情况下）：

将图像扩展到最佳尺寸

DFT 的性能取决于图像的大小。当图像大小是 2、3 和 5 的倍数时，DFT 的速度最快。因此，要获得最佳性能，通常最好在图像中填充边框值，以获得具有这些特征的图像大小。getOptimalDFTSize() 会返回这种最佳尺寸，我们可以使用 copyMakeBorder() 函数来扩展图像的边框（添加的像素初始化为零）：

 Mat padded; // 将输入图像扩展到最佳尺寸
 int m = getOptimalDFTSize( I.rows )；
 int n = getOptimalDFTSize( I.cols ); // 在边框上添加零值
 copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0))；

同时放置复数值和实数值

傅立叶变换的结果是复值。这意味着每个图像值的结果都是两个图像值（每个分量一个）。此外，频域范围远大于空间范围。因此，我们通常至少以浮点格式存储这些值。因此，我们要将输入图像转换为这种类型，并用另一个通道来扩展它，以保存复数值：

 Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}；
 Mat complexI；
 merge(planes, 2, complexI); // 在扩展后的另一个平面上添加零点

进行离散傅立叶变换

可以进行就地计算（输入与输出相同）：

 dft(complexI, complexI); // 这样计算结果就可以放入源矩阵中了

将实数和复数值转换为幅度

复数有实部（Re）和虚部（Im）。DFT 的结果就是复数。DFT 的幅值为

$$M=\sqrt{\text{Re}{\left(DFT\left(I\right)\right)}^2+\text{Im}{\left(DFT\left(I\right)\right)}^2}$$
翻译为 OpenCV 代码：

 split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
 magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
 Mat magI = planes[0]；

切换到对数标度

事实证明，傅立叶系数的动态范围太大，无法在屏幕上显示。我们有一些变化较小的值和一些变化较大的值，无法像这样观察。因此，数值大的会显示为白点，数值小的会显示为黑点。为了将灰度值用于可视化，我们可以将线性刻度转换为对数刻度：

$$M_1=log\left(1+M\right)$$

转换为 OpenCV 代码：

 magI += Scalar::all(1); // 切换到对数标度
 log(magI, magI)；

裁剪并重新排列

还记得在第一步，我们展开了图像吗？现在是时候丢掉新引入的值了。出于可视化目的，我们还可以重新排列结果的象限，使原点（零，零）与图像中心相对应。

 // 如果光谱的行数或列数为奇数，则裁剪光谱
 magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2))；
 // 重新排列傅立叶图像的象限，使原点位于图像中心
 int cx = magI.cols/2；
 int cy = magI.rows/2；
 Mat q0(magI, Rect(0, 0, cx, cy)); // 左上角--每个象限创建一个 ROI
 Mat q1(magI, Rect(cx, 0, cx, cy)); // 右上角
 Mat q2(magI, Rect(0, cy, cx, cy)); // 左下角
 Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
 Mat tmp; // 交换象限（左上角和右下角）
 q0.copyTo(tmp)；
 q3.copyTo(q0)；
 tmp.copyTo(q3)；
 q1.copyTo(tmp); // 交换象限（右上角和左下角）
 q2.copyTo(q1)；
 tmp.copyTo(q2)；

规范化

为了实现可视化，我们再次这样做。现在我们有了幅度，但仍超出了图像显示范围（0 到 1）。我们使用 cv::normalize() 函数将数值归一化到这个范围。

 normalize(magI, magI, 0, 1, NORM_MINMAX); // 将带有浮点数值的矩阵转换为可视图像形式（浮点数值介于 0 和 1 之间）。
 // 可查看的图像形式（数值在 0 和 1 之间的浮点）。

结果

一种应用思路是确定图像中存在的几何方向。例如，让我们找出文字是否水平？通过观察一些文字，你会发现文字的线条有点像水平线，而字母则有点像垂直线。在傅立叶变换中也可以看到文本片段的这两个主要组成部分。让我们用水平图像和旋转图像来表示文字。

如果是水平文本，则

如果是旋转文字

您可以看到，频域中影响最大的分量（幅度图像上最亮的点）随着图像上物体的几何旋转而变化。由此，我们可以计算出偏移量，并进行图像旋转，以纠正最终的错位。