离散傅里叶变换

目录

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

上一个教程： 改变图像的对比度和亮度！
下一个教程： 使用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] << " [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)); // 左上角 - 为每个象限创建一个感兴趣区域
    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("Input Image" , I ); // 显示结果
imshow("spectrum magnitude", magI);
cv::Mat::empty
 
    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)); // 左上角 - 为每个象限创建一个感兴趣区域
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之间）。
 
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

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] # 左上角 - 为每个象限创建一个感兴趣区域
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])

为复数值和实数值腾出空间

傅里叶变换的结果是复数。这意味着对于每个图像值，结果是两个图像值（每个分量一个）。此外，频率域的范围比其空间域对应物大得多。因此，我们通常至少以*浮点*格式存储这些值。因此，我们会将输入图像转换为此类型，并用另一个通道来存储复数值。

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)); // 左上角 - 为每个象限创建一个感兴趣区域
    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)); // 左上角 - 为每个象限创建一个感兴趣区域
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] # 左上角 - 为每个象限创建一个感兴趣区域
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之间）。

Python

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

结果

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

对于水平文本

对于旋转文本

你可以看到频率域中最具影响力的分量（幅度图像上最亮点）遵循图像上物体的几何旋转。由此我们可以计算偏移量并执行图像旋转以纠正最终的未对齐。