傅里叶变换
二维DFT的极坐标表示
幅度或频率谱为
R(u,v)和I(u,v)分别是F(u,v)的实部和虚部
相角或相位谱为
使用书中图4-24和图4-25测试,显示傅里叶谱和相位谱
代码如下:
#include "opencv2/opencv.hpp"
using namespace cv;
#define PI2 2*3.141592654
int main( int argc, char *argv[] )
{
const char* filename = argc >=2 ? argv[1] : "lena.jpg";
Mat image = imread(filename, IMREAD_GRAYSCALE);
if( image.empty())
return -1;
imshow("src",image);
image.convertTo(image,CV_32FC1);
////////////////////////////////////////二维基本傅里叶变换//////////////////////////////////////////////////
// for(int i=0; i<image.rows; i++) { //中心化
// float *p = image.ptr<float>(i);
// for(int j=0; j<image.cols; j++){
// p[j] = p[j] * pow(-1, i+j);
// }
// }
// Mat dftRe = Mat::zeros(image.size(), CV_32FC1);
// Mat dftIm = Mat::zeros(image.size(), CV_32FC1);
// for(int u=0; u<image.rows; u++)
// {
// float *pRe = dftRe.ptr<float>(u);
// float *pIm = dftIm.ptr<float>(u);
// for(int v=0; v<image.cols; v++)
// {
// float sinDft=0, cosDft=0;
// for(int i=0; i<image.rows; i++)
// {
// float *q = image.ptr<float>(i);
// for(int j=0; j<image.cols; j++)
// {
// float temp = PI2 *((float)u*i/image.rows + (float)v*j/image.cols);
// sinDft -= q[j] * sin(temp);
// cosDft += q[j] * cos(temp);
// }
// }
// pRe[v] = sinDft;
// pIm[v] = cosDft;
// }
// }
// divide(dftRe, image.rows*image.rows, dftRe);
// divide(dftIm, image.rows*image.rows, dftIm);
// multiply(dftIm, dftIm, dftIm);
// multiply(dftRe, dftRe, dftRe);
// add(dftRe, dftIm, dftRe);
// pow(dftRe, 0.5, dftRe);
// imshow("mydft", dftRe);
///////////////////////////////////////快速傅里叶变换/////////////////////////////////////////////////////
int oph = getOptimalDFTSize(image.rows);
int opw = getOptimalDFTSize(image.cols);
Mat padded;
copyMakeBorder(image, padded, 0, oph-image.rows, 0, opw-image.cols, BORDER_CONSTANT, Scalar::all(0));
Mat temp[] = {padded, Mat::zeros(image.size(),CV_32FC1)};
Mat complexI;
merge(temp, 2, complexI);
dft(complexI, complexI); //傅里叶变换
//显示频谱图
split(complexI, temp);
Mat amplitude, angle;
magnitude(temp[0], temp[1], amplitude);
phase(temp[0], temp[1], angle);
// cartToPolar(temp[0], temp[1],amplitude, angle);
int cx = amplitude.cols/2;
int cy = amplitude.rows/2;
Mat q0(amplitude, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant
Mat q1(amplitude, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(amplitude, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(amplitude, 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);
Mat amplitude_src;
divide(amplitude, oph*opw, amplitude_src );
imshow("amplitude_src",amplitude_src);
amplitude += Scalar::all(1); // switch to logarithmic scale
log(amplitude, amplitude);
normalize(amplitude, amplitude, 0, 255, NORM_MINMAX); //归一化 方便显示,和实际数据没有关系
amplitude.convertTo(amplitude, CV_8U);
imshow("amplitude",amplitude);
normalize(angle, angle, 0, 255, NORM_MINMAX); //归一化 方便显示,和实际数据没有关系
angle.convertTo(angle, CV_8U);
imshow("angle",angle);
imwrite("phase.jpg", angle)
waitKey(0);
return 1;
}
结果:
图4-24 原图及原始傅里叶图
log变换后的傅里叶图及相位图
图4-25的图像