离散傅里叶变换

目录

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

上一教程： 改变图像的对比度和亮度！
下一教程： 使用 XML/YAML/JSON 文件进行文件输入和输出

原作者	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] << " [图像名称 -- 默认 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] = 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)); // 右下角
 
    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("输入图像" , 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 [图像名称 -- 默认 ../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);
q1.copyTo(tmp);
 
tmp.copyTo(q2);
magI.convertTo(magI, CvType.CV_8UC1);
                                                                            Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // 将浮点型数值矩阵转换为可视图像形式
                                                                            // （浮点值范围为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 [图像名称 -- 默认 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("Input Image" , I ) # 显示结果
    cv.imshow("spectrum magnitude", 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] = magnitude
    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);
q1.copyTo(tmp);

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

归一化

这同样是为了可视化目的。我们现在有了幅度，但是这些仍然超出了我们图像显示范围（零到一）。我们使用 cv::normalize() 函数将我们的值归一化到此范围。

C++

    normalize(magI, magI, 0, 1, NORM_MINMAX); // 将具有浮点值的矩阵转换为
                                            //  可视图像形式（浮点值范围为0到1）。

Java

magI.convertTo(magI, CvType.CV_8UC1);
                                                                            Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // 将浮点型数值矩阵转换为可视图像形式
                                                                            // （浮点值范围为0到255）。

Python

    cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # 将浮点型矩阵转换为
    

结果

一个应用想法是确定图像中存在的几何方向。例如，让我们找出文本是否是水平的？观察一些文本，你会注意到文本行也形成水平线，字母也形成垂直线。文本片段的这两个主要组成部分也可以在傅里叶变换中看到。让我们使用这张水平的图像和这张旋转的图像（关于文本）。

对于水平文本的情况

对于旋转文本的情况

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