离散傅里叶变换

目录

• 目标
• 源代码
• 解释
  • 将图像扩展到最佳大小
  • 为复数和实数值腾出空间
  • 进行离散傅里叶变换
  • 将实数值和复数值转换为幅度
  • 切换到对数刻度
  • 裁剪和重新排列
  • 标准化
• 结果

上一教程： 更改图像的对比度和亮度！
下一教程： 使用 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 中找到它。

Java

你可以 从这里下载，或者在 OpenCV 源代码库的 samples/java/tutorial_code/core/discrete_fourier_transform/DiscreteFourierTransform.java 中找到它。

Python

你可以 从这里下载，或者在 OpenCV 源代码库的 samples/python/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.py 中找到它。

以下是 dft() 的示例用法

C++

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
 
#include <iostream>
 
using namespace cv;
using namespace std;
 
static void help(char ** argv)
{
cout << endl
<< "本程序演示了离散傅里叶变换 (DFT) 的使用方法。" << endl
<< "获取图像的 DFT 并显示其功率谱。" << endl << endl
<< "用法:" << 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 << "打开图像错误" << 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] = 幅度
    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)); //右下角
 
    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 之间）。
 
imshow("输入图像" , I ); //显示结果
imshow("频谱幅度", magI);
waitKey();
 
    return EXIT_SUCCESS;
}

Java

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("" +
                "本程序演示了离散傅里叶变换 (DFT) 的使用方法。\n" +
                "获取图像的 DFT 并显示其功率谱。\n" +
                "用法:\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("打开图像错误");
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>();
padded.convertTo(padded, 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) = 幅度
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)); //左上角 - 为每个象限创建一个 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)); // 右下角
 
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); // 将具有浮点值的矩阵转换为可视图像形式（浮点数介于 0 和 255 之间）。
                                                                            // 到可视图像形式（浮点数在
                                                                            // 0 和 255 之间）。
 
HighGui.imshow("输入图像" , I ); // 显示结果
HighGui.imshow("频谱幅值", 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

from __future__ import print_function
import sys
 
import cv2 as cv
import numpy as np
 
 
def print_help()
print('''
该程序演示了离散傅里叶变换 (DFT) 的用法。
对图像进行 DFT 并显示其功率谱。
用法
discrete_fourier_transform.py [image_name -- 默认 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 is None
print('打开图像时出错')
        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] = 幅值
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)
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("输入图像" , I ) // 显示结果
    cv.imshow("频谱幅值", magI)
    cv.waitKey()
 
if __name__ == "__main__"
    main(sys.argv[1:])

说明

傅里叶变换将分解图像为正弦和余弦分量。换句话说，它将图像从空间域转换为频率域。其理念是，任何函数都可以用无限个正弦和余弦函数的总和精确逼近。傅里叶变换是一种如何做到这一点的方法。在数学上，二维图像的傅里叶变换为

\[F(k,l) = \displaystyle\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1} f(i,j)e^{-i2\pi(\frac{ki}{N}+\frac{lj}{N})}\]

\[e^{ix} = \cos{x} + i\sin {x}\]

这里，f 是图像在空间域中的值，F 是在频率域中的值。变换的结果是复数。可以通过实图像和复图像，或通过幅值图像和相位图像来显示。但是，在整个图像处理算法中，只有幅值图像是有趣的，因为它包含了关于图像几何结构的所有信息。但是，如果你想以这些形式对图像进行一些修改，然后你需要重新转换它，那么你需要保留这两者。

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

将图像扩展到最佳大小

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

C++

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

Java

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

Python

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

为复数和实数值腾出空间

傅里叶变换的结果是复数。这意味着，对于每个图像值，结果是两个图像值（每个分量一个）。此外，频率域的范围远大于其空间对应域。因此，我们通常至少将它们存储在float 格式中。因此，我们将把我们的输入图像转换为这种类型，并用另一个通道扩展它以保存复数值。

C++

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
merge(planes, 2, complexI); //向扩展图像添加一个全零平面

Java

List<Mat> planes = new ArrayList<Mat>();
padded.convertTo(padded, CvType.CV_32F);
planes.add(padded);
planes.add(Mat.zeros(padded.size(), CvType.CV_32F));
Mat complexI = new Mat();
Core.merge(planes, complexI); //向扩展图像添加一个全零平面

Python

planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
complexI = cv.merge(planes) # 在扩展的另一个平面中添加零

进行离散傅里叶变换

可以使用就地计算（输入和输出相同）。

C++

dft(complexI, complexI); //这样结果就可以保存在源矩阵中

Java

Core.dft(complexI, complexI); //这样结果就可以保存在源矩阵中

Python

    cv.dft(complexI, complexI) # 以这种方式，结果可以适合源矩阵

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

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

\[M = \sqrt[2]{ {Re(DFT(I))}^2 + {Im(DFT(I))}^2}\]

翻译成 OpenCV 代码

C++

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

Java

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) = 幅度
Mat magI = planes.get(0);

Python

    cv.split(complexI, planes) # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = 幅值
magI = planes[0]

切换到对数刻度

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

\[M_1 = \log{(1 + M)}\]

翻译成 OpenCV 代码

C++

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

Java

Mat matOfOnes = Mat.ones(magI.size(), magI.type());
Core.add(matOfOnes, magI, magI); //切换到对数刻度
Core.log(magI, magI);

Python

matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
    cv.add(matOfOnes, magI, magI) # 切换到对数刻度
    cv.log(magI, magI)

裁剪和重新排列

请记住，我们是在第一步扩展了图像？现在是时候丢弃新引入的值了。为了可视化，我们还可以重新排列结果的象限，以便原点（零，零）对应于图像中心。

C++

    //如果频谱的行或列数量是奇数，则裁剪频谱
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)); //右下角
 
    Mat tmp; //交换象限（左上角与右下角）
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
 
q1.copyTo(tmp); //交换象限（右上角与左下角）
q2.copyTo(q1);
tmp.copyTo(q2);

Java

        //如果频谱的行或列数量是奇数，则裁剪频谱
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)); //左上角 - 为每个象限创建一个 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)); // 右下角
 
Mat tmp = new Mat(); // 交换象限（左上角与右下角）
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
 
q1.copyTo(tmp); //交换象限（右上角与左下角）
q2.copyTo(q1);
tmp.copyTo(q2);

Python

magI_rows, magI_cols = magI.shape
    // 如果频谱的行数或列数为奇数，则裁剪频谱
magI = magI[0:(magI_rows & -2), 0:(magI_cols & -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

规范化

这再次是为了可视化。现在我们有了幅值，但它们仍然超出了我们的图像显示范围 0 到 1。我们使用cv::normalize() 函数将我们的值规范化到此范围。

C++

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

Java

Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // 将具有浮点值的矩阵转换为可视图像形式（浮点数介于 0 和 255 之间）。
                                                                            // 到可视图像形式（浮点数在
                                                                            // 0 和 255 之间）。

Python

    cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) // 将具有浮点值的矩阵转换为
    

结果

一个应用理念是确定图像中存在的几何方向。例如，让我们找出文本是水平的还是不是水平的？观察一些文本，你会注意到文本行在某种程度上也形成了水平线，字母在某种程度上也形成了垂直线。文本片段的这两个主要组成部分也可能在傅里叶变换中看到。让我们使用这个水平 和 这个旋转 的关于文本的图像。

在水平文本的情况下

在旋转文本的情况下

您可以看到，频域中最有影响力的成分（幅度图像中最亮的点）遵循图像中物体的几何旋转。由此，我们可以计算偏移量并执行图像旋转以校正最终的错位。