不知道哪里有问题,出不来正确的结果
// fourier.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/legacy/legacy.hpp>
#include <iostream>
#include <vector>
using namespace cv;
using namespace std;
int _tmain(int argc, _TCHAR* argv[])
{
string name1="0-0.bmp";
string name2="15-10.bmp";
/
cv::Mat src1=cv::imread(name1,0);
cv::Mat padded=src1;
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);
cv::dft(complexI, complexI);
cv::split(complexI, planes);
//
cv::Mat src2=cv::imread(name2,0);//读取输入图像
cv::Mat padded2=src2;
//cv::Mat dst2=cv::Mat::zeros(src2.size(),CV_32FC2);
Mat planes2[] = {Mat_<float>(padded2), Mat::zeros(padded2.size(), CV_32F)};
Mat complexI2;
merge(planes2, 2, complexI2);
cv::dft(complexI2, complexI2);
cv::split(complexI2, planes2);
//cv::normalize(planes[1],planes[1], 0, 1, CV_MINMAX);
//src1变换后实部 planes[0] 虚部planes[1]
//src2变换后实部 planes2[0] 虚部planes2[1]
//互功率谱存放矩阵 实部 planes3[0] 虚部planes3[1]
//
Mat planes3[] = {Mat::zeros(planes[0].size(), CV_32F), Mat::zeros(planes[0].size(), CV_32F)};
for(int i=0;i<src1.rows-1;i++)
{
for(int j=0;j<src1.cols-1;j++)
{
double r1=planes[0].at<float>(j,i);//获取像素值
double i1=planes[1].at<float>(j,i);
double r2=planes2[0].at<float>(j,i);
double i2=planes2[1].at<float>(j,i);
double r3 =r1*r2+i1*i2;//计算互功率谱的值
double i3 =r1*i2-r2*i1;
double abs =sqrt((r3*r3)+(i3*i3));
double r_exp =r3/abs;
double i_exp =i3/abs;
planes3[0].at<float>(j,i)=r_exp;
planes3[1].at<float>(j,i)=i_exp;
}
}
//cv::Mat dst=cv::Mat::zeros(src1.size(),CV_32FC1);
magnitude(planes3[0], planes3[1], planes3[0]);
Mat complexI3=planes3[0];
//merge(planes3, 2, complexI3);
cv::dft(complexI3,complexI3,CV_DXT_INV_SCALE);
double max=0,min=0;
cv::Point minLoc,maxLoc;
cv::minMaxLoc(complexI3,&min,&max,&minLoc,&maxLoc);
cout<<"max "<<max<<" min "<<min<<endl;
cout<<maxLoc.x<<" "<<maxLoc.y<<endl;
cv::namedWindow("re");
cv::imshow("re",complexI3);
cv::waitKey();
return 0;
}
修改好的版本
// fourier.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/legacy/legacy.hpp>
#include <iostream>
#include <vector>
using namespace cv;
using namespace std;
int _tmain(int argc, _TCHAR* argv[])
{
string name1="0-0.bmp";
string name2="100-100.bmp";
/
cv::Mat src1=cv::imread(name1,0);
cv::Mat padded=src1;
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);
cv::dft(complexI, complexI);
cv::split(complexI, planes);
//
cv::Mat src2=cv::imread(name2,0);//读取输入图像
cv::Mat padded2=src2;
//cv::Mat dst2=cv::Mat::zeros(src2.size(),CV_32FC2);
Mat planes2[] = {Mat_<float>(padded2), Mat::zeros(padded2.size(), CV_32F)};
Mat complexI2;
merge(planes2, 2, complexI2);
cv::dft(complexI2, complexI2);
cv::split(complexI2, planes2);
//cv::normalize(planes[1],planes[1], 0, 1, CV_MINMAX);
//src1变换后实部 planes[0] 虚部planes[1]
//src2变换后实部 planes2[0] 虚部planes2[1]
//互功率谱存放矩阵 实部 planes3[0] 虚部planes3[1]
//
Mat planes3[] = {Mat::zeros(src1.size(), CV_32F), Mat::zeros(src1.size(), CV_32F)};
for(int i=0;i<src1.rows-1;i++)
{
for(int j=0;j<src1.cols-1;j++)
{
double r1=planes[0].at<float>(j,i);//获取像素值
double i1=planes[1].at<float>(j,i);
double r2=planes2[0].at<float>(j,i);
double i2=planes2[1].at<float>(j,i);
double r3 =r1*r2+i1*i2;//计算互功率谱的值
double i3 =r1*i2-r2*i1;
double abs =sqrt((r3*r3)+(i3*i3));
double r_exp =r3/abs;
double i_exp =i3/abs;
planes3[0].at<float>(j,i)=r_exp;
planes3[1].at<float>(j,i)=i_exp;
}
}
//magnitude(planes3[0], planes3[1], planes3[0]);
//Mat complexI3=planes3[0];
//cv::Mat =cv::Mat::zeros(dst.size(),dst.type());
//cv::flip(complexI3,complexI3,-1);
Mat complexI3;
merge(planes3, 2, complexI3);
cv::dft(complexI3,complexI3,CV_DXT_INV_SCALE);
cv::split(complexI3,planes3);
magnitude(planes3[0], planes3[1], planes3[0]);
complexI3=planes3[0];
//cv::add(complexI3,cv::Scalar(1.0),complexI3,NULL);
//cv::log(complexI3,complexI3);
complexI3 += Scalar::all(1);
log(complexI3, complexI3);
cv::normalize(complexI3,complexI3, 0, 1, CV_MINMAX);
/*
complexI3 = complexI3(Rect(0, 0, complexI3.cols & -2, complexI3.rows & -2));
int cx = complexI3.cols/2;
int cy = complexI3.rows/2;
Mat q0(complexI3, Rect(0, 0, cx, cy)); // Top-Left - 为每一个象限创建ROI
Mat q1(complexI3, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(complexI3, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(complexI3, Rect(cx, cy, cx, cy)); // Bottom-Right
Mat tmp; // 交换象限 (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 交换象限 (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);
*/
double max=0,min=0;
cv::Point minLoc,maxLoc;
cv::minMaxLoc(complexI3,&min,&max,&minLoc,&maxLoc);
//cv::imwrite("complexI3.bmp",complexI3);
cout<<"max "<<max<<" min "<<min<<endl;
cout<<maxLoc.x<<" "<<maxLoc.y<<endl;
cv::namedWindow("re");
cv::imshow("re",planes3[1]);
cv::waitKey();
return 0;
}