opencv----计算互功率谱

不知道哪里有问题，出不来正确的结果

// 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 - 为每一个象限创建ROIMat q1(complexI3, Rect(cx, 0, cx, cy));  // Top-RightMat q2(complexI3, Rect(0, cy, cx, cy));  // Bottom-LeftMat q3(complexI3, Rect(cx, cy, cx, cy)); // Bottom-RightMat 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;}