作者：gnuhpc
出处：http://www.cnblogs.com/gnuhpc/
#include "cv.h"
#include "highgui.h"
#include "cxcore.h"
void cvShiftDFT(CvArr *src_arr,CvArr *dst_arr)
{
CvMat * tmp;
CvMat q1stub,q2stub;
CvMat q3stub,q4stub;
CvMat d1stub,d2stub;
CvMat d3stub,d4stub;
CvMat *q1,*q2,*q3,*q4;
CvMat *d1,*d2,*d3,*d4;
CvSize size = cvGetSize(src_arr);
CvSize dst_size = cvGetSize(dst_arr);
int cx,cy;
if ((dst_size.width!= size.width)||(dst_size.height!=size.height))
{
cvError(CV_StsUnmatchedSizes,"cvShiftDFT",
"Source and Destination arrays must have the same sizes",
__FILE__,__LINE__);
}
if (src_arr == dst_arr)
{
tmp = cvCreateMat(size.height/2,size.width/2,cvGetElemType(src_arr));
}
cx=size.width/2;//取出图像的原点
cy=size.height/2;
q1=cvGetSubRect(src_arr,&q1stub,cvRect(0,0,cx,cy));
//取出图像的第一象限，由q1指针指向它
q2=cvGetSubRect(src_arr,&q2stub,cvRect(cx,0,cx,cy));
//取出图像的第二象限，由q2指针指向它
q3=cvGetSubRect(src_arr,&q3stub,cvRect(cx,cy,cx,cy));
//取出图像的第三象限，由q3指针指向它
q4=cvGetSubRect(src_arr,&q4stub,cvRect(0,cy,cx,cy));
//取出图像的第四象限，由q4指针指向它
d1=cvGetSubRect(src_arr,&d1stub,cvRect(0,0,cx,cy));
d2=cvGetSubRect(src_arr,&d2stub,cvRect(cx,0,cx,cy));
d3=cvGetSubRect(src_arr,&d3stub,cvRect(cy,cy,cx,cy));
d4=cvGetSubRect(src_arr,&d4stub,cvRect(0,cy,cx,cy));
if (src_arr!=dst_arr)
{
if (!CV_ARE_TYPES_EQ(q1,d1))
{
cvError(CV_StsUnmatchedSizes,"cvShiftDFT",
"Source and Destination arrays must have the same sizes",
__FILE__,__LINE__);
}
//以图像中心为原点，调整傅里叶变换图像的四个象限区，
//即第一与第三象限交换，第二与第四象限交换
cvCopy(q3,d1,0);
cvCopy(q4,d2,0);
cvCopy(q1,d3,0);
cvCopy(q2,d4,0);
}
else
{//若源矩阵和目的矩阵相同则直接在源矩阵中进行操作
cvCopy(q3,tmp,0);
cvCopy(q1,q3,0);
cvCopy(tmp,q1,0);
cvCopy(q4,tmp,0);
cvCopy(q2,q4,0);
cvCopy(tmp,q2,0);
}
}
int main(int argc,char ** argv)
{
const char* filename =(argc>=2?argv[1]:"lena.jpg");
IplImage *im;
IplImage *realInput,*imaginaryInput,*complexInput;
IplImage *image_Re,*image_Im;
int dft_M,dft_N;
CvMat *dft_A;
CvMat tmp;
double m,M;
im = cvLoadImage(filename, CV_LOAD_IMAGE_GRAYSCALE );//加载图像
if (!im)
{
return -1;
}
//分配空间
realInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道
imaginaryInput =cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道
complexInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,2);//双通道
cvScale(im,realInput,1.0,0.0);
//#define cvScale cvConvertScale=>readInput=im
cvZero(imaginaryInput);
//清空这个图像的内容
cvMerge(realInput,imaginaryInput,NULL,NULL,complexInput);
//混合这两个图像作为complexInput的两个通道
/*得到最优DFT尺寸 */
dft_M = cvGetOptimalDFTSize(im->height-1);
dft_N = cvGetOptimalDFTSize(im->width-1);
dft_A = cvCreateMat(dft_M,dft_N,CV_64FC2);
image_Re = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//实部
image_Im = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//虚部
cvGetSubRect(dft_A,&tmp,cvRect(0,0,im->width,im->height));
cvCopy(complexInput,&tmp,NULL);
if( dft_A->cols > im->width )//若得到的最优DFT尺寸在宽度上大于原图，则重新取
{
cvGetSubRect(dft_A,&tmp,cvRect(im->width,0,dft_A->cols-im->width,im->height));
cvZero(&tmp);
}
cvDFT(dft_A,dft_A,CV_DXT_FORWARD,complexInput->height);
cvNamedWindow("win",0);
cvNamedWindow("magnitude",0);
cvShowImage("win",im);
//分割出实部和虚部
cvSplit(dft_A,image_Re,image_Im,0,0);
//计算功率谱 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,NULL);//image_Re<=image_Re+image_Im
cvPow(image_Re,image_Re,0.5);
//计算log(1+Mag)
cvAddS(image_Re,cvScalarAll(1.0),image_Re,NULL);
cvLog(image_Re,image_Re);
cvShiftDFT(image_Re,image_Re);
cvMinMaxLoc(image_Re,&m,&M,NULL,NULL,NULL);
cvScale(image_Re,image_Re,1.0/(M-m),1.0*(-m)/(M-m));
cvShowImage("magnitude",image_Re);
cvWaitKey(-1);
return 0;
}