比上一篇blog的功能更加全面一点
上一篇只能计算单独旋转后图像的旋转角度
现在能够实现单独计算平移且旋转之后的图像的旋转角度
五一节过后第一天解决了这个让我蛋疼很久的问题
果然是磨刀不误砍柴功
刚发现这个代码的逻辑其实是有错误的
只能参考,请勿copy
不好意思,误导大家了
// test.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 <cmath>
#include <vector>
#include <cxcore.h>
using namespace cv;
using namespace std;
void fft2(IplImage *src, IplImage *dst);
void fft2shift(IplImage *src, IplImage *dst);
void pingyi(IplImage *in_1,IplImage *in_2,int &x,int &y);
int main(int argc, char ** argv)
{
IplImage *src; //源图像
IplImage *Fourier; //傅里叶系数
IplImage *dst ;
IplImage *ImageRe;
IplImage *ImageIm;
IplImage *ImageIm1;
IplImage *Image;
IplImage *ImageDst1;
IplImage *ImageDst2;
IplImage *fourier1;
IplImage *fourier2;
double m,M;
double scale;
double shift;
src = cvLoadImage("lena.bmp",0); //加载源图像,第二个参数表示将输入的图片转为单信道
Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
ImageDst1 = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
fft2(src,Fourier); //傅里叶变换
fft2shift(Fourier, Image); //中心化
cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换,并对结果进行缩放
cvSplit(dst,ImageRe,ImageIm,0,0);
cvNamedWindow("源图像",0);
cvShowImage("源图像",src);
//对数组每个元素平方并存储在第二个参数中
cvPow(ImageRe,ImageRe,2);
cvPow(ImageIm,ImageIm1,2);
cvAdd(ImageRe,ImageIm1,ImageRe,NULL);
cvPow(ImageRe,ImageRe,0.5);
cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
scale = 255/(M - m);
shift = -m * scale;
//将shift加在ImageRe各元素按比例缩放的结果上,存储为ImageDst1
cvConvertScale(ImageRe,ImageDst1,scale,shift);
cvNamedWindow("傅里叶谱",0);
cvShowImage("傅里叶谱",ImageRe);
fourier1=ImageRe;
//cvNamedWindow("傅里叶逆变换",0);
//cvShowImage("傅里叶逆变换",ImageDst1);
//
//第二幅图像进行操作/
src = cvLoadImage("75du.bmp",0); //加载源图像,第二个参数表示将输入的图片转为单信道
Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
ImageDst2 = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
fft2(src,Fourier); //傅里叶变换
fft2shift(Fourier, Image); //中心化
cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换,并对结果进行缩放
cvSplit(dst,ImageRe,ImageIm,0,0);
cvNamedWindow("90du图像",0);
cvShowImage("90du图像",src);
//对数组每个元素平方并存储在第二个参数中
cvPow(ImageRe,ImageRe,2);
cvPow(ImageIm,ImageIm1,2);
cvAdd(ImageRe,ImageIm1,ImageRe,NULL);
cvPow(ImageRe,ImageRe,0.5);
cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
scale = 255/(M - m);
shift = -m * scale;
//将shift加在ImageRe各元素按比例缩放的结果上,存储为ImageDst2
cvConvertScale(ImageRe,ImageDst2,scale,shift);
cvNormalize(ImageIm,ImageIm, 0, 1, CV_MINMAX);
cvNamedWindow("90du傅里叶谱",0);
cvShowImage("90du傅里叶谱",ImageRe);
fourier2=ImageRe;
//cvNamedWindow("2傅里叶反变换",0);
//cvShowImage("2傅里叶反变换",ImageDst2);
cvLogPolar( fourier1, fourier1, cvPoint2D32f(fourier1->width/2,fourier1->height/2), 40, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS );
cvLogPolar( fourier2, fourier2, cvPoint2D32f(fourier2->width/2,fourier2->height/2), 40, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS );
int xx=0,yy=0;
double theta=0;
pingyi(fourier1,fourier2,xx,yy);
cout<<"xx "<<xx<<" "<<"yy "<<yy<<endl;
theta=(double)(yy)/fourier1->height*360;
cout<<"theta="<<(int)theta<<endl;
cvWaitKey(0);
cvReleaseImage(&src);
cvReleaseImage(&Image);
cvReleaseImage(&ImageIm);
cvReleaseImage(&ImageRe);
cvReleaseImage(&Fourier);
cvReleaseImage(&dst);
cvReleaseImage(&ImageDst1);
cvReleaseImage(&ImageDst2);
cvDestroyAllWindows();
return 0;
return 0;
}
void fft2(IplImage *src, IplImage *dst)
{
IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0; //实部、虚部
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); //实部
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); //虚部
Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);//2 channels (image_Re, image_Im)
cvConvertScale(src, image_Re, 1, 0);// Real part conversion from u8 to 64f (double)
cvZero(image_Im);// Imaginary part (zeros)
cvMerge(image_Re, image_Im, 0, 0, Fourier);// Join real and imaginary parts and stock them in Fourier image
cvDFT(Fourier, dst, CV_DXT_FORWARD);// Application of the forward Fourier transform
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
cvReleaseImage(&Fourier);
}
void fft2shift(IplImage *src, IplImage *dst)
{
IplImage *image_Re = 0, *image_Im = 0;
int nRow, nCol, i, j, cy, cx;
double scale, shift, tmp13, tmp24;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
cvSplit( src, image_Re, image_Im, 0, 0 );
//具体原理见冈萨雷斯数字图像处理p123
// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
//计算傅里叶谱
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re);
cvPow( image_Re, image_Re, 0.5 );
//对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射一宽带输出值,具体可见冈萨雷斯数字图像处理p62)
// Compute log(1 + Mag);
cvAddS( image_Re, cvScalar(1.0), image_Re ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)
//Rearrange the quadrants of Fourier image so that the origin is at the image center
nRow = src->height; nCol = src->width;
cx = nCol/2; cy = nRow/2; // image center
//CV_IMAGE_ELEM为OpenCV定义的宏,用来读取图像的像素值,这一部分就是进行中心变换
for( j = 0; j < cy; j++ ){
for( i = 0; i < cx; i++ ){
//中心化,将整体份成四块进行对角交换
tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);
CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(image_Re, double, j+cy, i+cx);
CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;
tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);
CV_IMAGE_ELEM( image_Re, double, j, i+cx) =CV_IMAGE_ELEM( image_Re, double, j+cy, i);
CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;
}
}
//归一化处理将矩阵的元素值归一为[0,255]
//[(f(x,y)-minVal)/(maxVal-minVal)]*255
double minVal = 0, maxVal = 0;
// Localize minimum and maximum values
cvMinMaxLoc( image_Re, &minVal, &maxVal );
// Normalize image (0 - 255) to be observed as an u8 image
scale = 255/(maxVal - minVal);
shift = -minVal * scale;
cvConvertScale(image_Re, dst, scale, shift);
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
}
void pingyi(IplImage *in_1,IplImage *in_2,int &x,int &y)
{
cv::Mat padded(in_1);
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 padded2(in_2);
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(padded.size(), CV_32F), Mat::zeros(padded.size(), CV_32F)};
for(int i=0;i<padded.rows-1;i++)
{
for(int j=0;j<padded.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;
}
}
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];
complexI3 += Scalar::all(1);
log(complexI3, complexI3);
cv::normalize(complexI3,complexI3, 0, 1, CV_MINMAX);
double max=0,min=0;
cv::Point minLoc,maxLoc;
cv::minMaxLoc(complexI3,&min,&max,&minLoc,&maxLoc);
//cv::imwrite("complexI3.bmp",complexI3);
x=maxLoc.x;
y=maxLoc.y;
cout<<min<<" "<<max<<endl;
cout<<complexI3.at<float>(maxLoc)<<endl;
}