OpenCV学习笔记（7）：图像频域（DFT）

最新推荐文章于 2024-09-02 00:14:19 发布

转载最新推荐文章于 2024-09-02 00:14:19 发布 · 1.1k 阅读

OpenCV学习笔记专栏收录该内容

28 篇文章

订阅专栏

本文详细介绍了如何使用离散傅里叶变换（DFT）将图像从空间域转换到频率域，并绘制频域图像，探讨了频域在图像处理中的作用，包括图像压缩和去噪，并解释了DFT算法原理、图像频域表示含义及DFT的性质。

#include "stdafx.h"
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

#include <iostream>

using namespace cv;
using namespace std;

static void help(char* progName)
{
	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
		<< "Usage:" << endl
		<< progName << " [image_name -- default lena.jpg] " << endl << endl;
}

int main(int argc, char ** argv)
{
	help(argv[0]);

	const char* filename = argc >= 2 ? argv[1] : "lena.png";

	Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
	if (I.empty())
		return -1;

	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)));
	//Mat _magI = magI.clone();
	//normalize(_magI, _magI, 0, 0, CV_MINMAX);
	// 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, CV_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 before shift frequency", _magI);
	imshow("spectrum magnitude", magI);
	waitKey(0);

	return 0;
}

本程序的作用是：将图像从空间域转换到频率域，并绘制频域图像。

二维图像的DFT（离散傅里叶变换），

$F(k,l)={\sum_{i=0}^{i=N-1}}{\sum_{j=0}^{j=N-1}}f(i,j)e^{-j2\pi{(\frac{ki}{N}+\frac{lj}{N}})}$

图像的频域表示的是什么含义呢？又有什么用途呢？图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量；图像的噪声大部分情况下是高频部分；图像大部分平缓的灰度变化部分则为低频分量。也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。

频域在图像处理中，就我所知的用途主要在两方面：图像压缩和图像去噪。关于这两点将在下面给出图片DFT的变换结果后说明。

有关DFT的更多性质请参考胡广书教授的《数字信号处理》教材。
请注意读图片的函数与之前有所不同：

```
Mat image = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE);
```
CV_LOAD_IMAGE_GRAYSCALE参数表示将原图像转换为灰度图后读入，这是因为后面的DFT变换都是基于二维信号的，而彩色图像是三维信号。当然，也可以对RGB每一通道都进行DFT运算。
DFT算法的原理要求输入信号的长度最好为2^n，这样可以使用快速傅里叶变换算法（FFT算法）进行加速。所以程序中使用
```
copyMakeBorder(image, padded, 0, m-image.rows, 0, n-image.cols, BORDER_CONSTANT, Scalar::all(0));
```
填充0使横纵长度都为2^n。

对于一维信号，原DFT直接运算的复杂度是O(N^2)，而快速傅里叶变换的复杂度降低到O(Nlog2(N))，假设N为512，足足提高了512/9≈57倍。
由DFT的性质知，输入为实信号（图像）的时候，频域输出为复数，因此将频域信息分为幅值和相位。频域的幅值高的代表高频分量，幅值低的地方代表低频分量，因此程序中使用
```
// => log(1+sqrt(Re(DFT(I))^2+Im(DFT(I))^2))
magI += Scalar::all(1);
log(magI, magI);

// crop the spectrum
magI = magI(Rect(0, 0, magI.cols & (-2), magI.rows & (-2)));
Mat _magI = magI.clone();
normalize(_magI, _magI, 0, 1, CV_MINMAX);
```
进行log幅值计算及归一化幅值（归一化目的主要是方便将频域通过图像的形式进行显示）。

关于频域中心平移：将图像的高频分量平移到图像的中心，便于观测。

int cx = magI.cols/2;
int cy = magI.rows/2;

Mat q0(magI, Rect(0,0,cx,cy));    // Top-Left
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

// exchange Top-Left and Bottom-Right
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);

// exchange Top-Right and Bottom-Left
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);

其原理就是将左上角的频域和右下角的互换，右上角和左下角互换。

请注意：频域点和空域点的坐标没有一一对应的关系，两者的关系只是上面的DFT公式所见到的。