【opencv 450 core】Discrete Fourier Transform

十年一梦实验室

已于 2022-03-22 10:24:03 修改

阅读量1.7k

点赞数

分类专栏： opencv 文章标签： opencv c++ 计算机视觉

于 2022-03-22 10:22:40 首次发布

本文链接：https://blog.csdn.net/cxyhjl/article/details/123653865

版权

opencv 专栏收录该内容

87 篇文章 12 订阅

订阅专栏

Discrete Fourier Transform

OpenCV: Discrete Fourier Transform

Goal

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

什么是傅里叶变换，为什么要使用它？

如何在 OpenCV 中做到这一点？

使用以下函数：copyMakeBorder()、merge()、dft()、getOptimalDFTSize()、log() 和 normalize()。

Source code

您可以从这里下载或在 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;

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 -- default 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; //expand input image to optimal size

int m = getOptimalDFTSize( I.rows );

int n = getOptimalDFTSize( I.cols ); // on the border add zero values

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); // Add to the expanded another plane with zeros

dft(complexI, complexI); // this way the result may fit in the source matrix

// compute the magnitude and switch to logarithmic scale

// => 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); // switch to logarithmic scale

log(magI, magI);

// crop the spectrum, if it has an odd number of rows or columns

magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

// rearrange the quadrants of Fourier image so that the origin is at the image center

int cx = magI.cols/2;

int cy = magI.rows/2;

Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant

Mat q1(magI, Rect(cx, 0, cx, cy)); // Top-Right

Mat q2(magI, Rect(0, cy, cx, cy)); // Bottom-Left

Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

Mat tmp; // swap quadrants (Top-Left with Bottom-Right)

q0.copyTo(tmp);

q3.copyTo(q0);

tmp.copyTo(q3);

q1.copyTo(tmp); // swap quadrant (Top-Right with Bottom-Left)

q2.copyTo(q1);

tmp.copyTo(q2);

normalize(magI, magI, 0, 1, NORM_MINMAX); // Transform the matrix with float values into a

// viewable image form (float between values 0 and 1).

imshow("Input Image" , I ); // Show the result

imshow("spectrum magnitude", magI);

waitKey();

return EXIT_SUCCESS;

}

Explanation

傅立叶变换将图像分解为其正弦和余弦分量。换句话说，它将图像从其空间域转换到其频域。这个想法是任何函数都可以用无限正弦和余弦函数之和来精确逼近。傅里叶变换是一种方法。数学上二维图像傅里叶变换是：

这里 f 是其空间域中的图像值， F 是其频域中的图像值。转换的结果是复数。可以通过真实图像和复数图像或通过幅度和相位图像来显示这一点。然而，在整个图像处理算法中，只有幅度图像是有趣的，因为它包含了我们需要的关于图像几何结构的所有信息。不过，如果您打算以这些形式对图像进行一些修改，然后需要重新转换它，则需要保留这两种形式。

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

Expand the image to an optimal size

将图像扩大到最佳尺寸

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

Make place for both the complex and the real values

为复杂的和真实的值分配空间

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

The result of a Fourier Transform is complex. This implies that for each image value the result is two image values (one per component). Moreover, the frequency domains range is much larger than its spatial counterpart. Therefore, we store these usually at least in a float format. Therefore we'll convert our input image to this type and expand it with another channel to hold the complex values:

Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};

Mat complexI;

merge(planes, 2, complexI); // 用零添加到扩展的另一个平面Add to the expanded another plane with zeros

merge操作：

Mat m1 = (Mat_<uchar>(2,2) << 1,4,7,10);

Mat m2 = (Mat_<uchar>(2,2) << 2,5,8,11);

Mat m3 = (Mat_<uchar>(2,2) << 3,6,9,12);

Mat channels[3] = {m1, m2, m3};

Mat m;

merge(channels, 3, m);

m =

[ 1, 2, 3, 4, 5, 6;

7, 8, 9, 10, 11, 12]

m.channels() = 3

opencv中的merge函数_alickr的博客-CSDN博客_cv::merge

void merge(const Mat* mv, size_t count, OutputArray dst);

第一个参数是图像矩阵数组，第二个参数是需要合并矩阵的个数，第三个参数是输出

void merge(const vector& mv, OutputArray dst );

第一个参数是图像矩阵向量容器，第二个参数是输出，这种方法无需说明需要合并的矩阵个数，vector对象自带说明

Make the Discrete Fourier Transform

进行离散傅里叶变换

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

dft(complexI, complexI); // this way the result may fit in the source matrix

Transform the real and complex values to magnitude

将实数和复数转换为幅度

A complex number has a real (Re) and a complex (imaginary - Im) part. The results of a DFT are complex numbers. The magnitude of a DFT is:

复数具有实数 (Re) 和复数 (虚数 - Im) 部分。 DFT 的结果是复数。 DFT 的大小为：

翻译成 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];

opencv中的split函数_alickr的博客-CSDN博客_opencv split

void split(const Mat& src,Mat *mvBegin)

void split(InputArray m, OutputArrayOfArrays mv);

用法很显然，第一个参数为要进行分离的图像矩阵，第二个参数可以是Mat数组的首地址，或者一个vector<Mat>对象

char d[] = {1,2,3,4,5,6,7,8,9,10,11,12};

Mat m(2, 2, CV_8UC3, d);

Mat channels[3];

split(m, channels);

channels[0] =

[ 1, 4;

7, 10]

channels[1] =

[ 2, 5;

8, 11]

channels[2] =

[ 3, 6;

9, 12]

#include <opencv2/core.hpp>

magnitude()
void cv::magnitude	(	InputArray 
x,
		InputArray 
y,
		OutputArray 
magnitude 
	)		
Python:
	cv.magnitude(	x, y[, magnitude]	) ->	magnitude

Calculates the magnitude of 2D vectors.

函数 cv::magnitude 计算由 x 和 y 数组的相应元素形成的二维向量的大小：

参数

x：向量的 x 坐标的浮点数组。

y：向量y坐标的浮点数组；它的大小必须与 x 相同。

幅度：与 x 大小和类型相同的输出数组。

Switch to a logarithmic scale

切换到对数刻度

事实证明，傅立叶系数的动态范围太大而无法在屏幕上显示。我们有一些小的和一些高变化的值，我们不能像这样观察到。因此，高值将全部变为白点，而小值将变为黑色。为了使用灰度值进行可视化，我们可以将线性比例转换为对数比例：

It turns out that the dynamic range of the Fourier coefficients is too large to be displayed on the screen. We have some small and some high changing values that we can't observe like this. Therefore the high values will all turn out as white points, while the small ones as black. To use the gray scale values to for visualization we can transform our linear scale to a logarithmic one:

翻译成 OpenCV 代码：

magI += Scalar::all(1); //切换到对数刻度 switch to logarithmic scale

log(magI, magI);

Crop and rearrange

裁剪和重新排列

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

// 裁剪光谱，如果它有奇数的行或列crop the spectrum, if it has an odd number of rows or columns

magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

// 重新排列傅里叶图像的象限，使原点位于图像中心rearrange the quadrants of Fourier image so that the origin is at the image center

int cx = magI.cols/2;

int cy = magI.rows/2;

Mat q0(magI, Rect(0, 0, cx, cy)); //左上角 - 每个象限创建一个 ROI Top-Left - Create a ROI per quadrant