【OpenCV学习】DFT变换

【OpenCV学习】DFT变换

作者：gnuhpc 
出处：http://www.cnblogs.com/gnuhpc/

view plain print ?

1. #include "cv.h"  
2. #include "highgui.h"  
3. #include "cxcore.h"  
4. void cvShiftDFT(CvArr *src_arr,CvArr *dst_arr)  
5. {  
6.     CvMat * tmp;  
7.     CvMat q1stub,q2stub;  
8.     CvMat q3stub,q4stub;  
9.     CvMat d1stub,d2stub;  
10.     CvMat d3stub,d4stub;  
11.     CvMat *q1,*q2,*q3,*q4;  
12.     CvMat *d1,*d2,*d3,*d4;  
13.     CvSize size = cvGetSize(src_arr);  
14.     CvSize dst_size = cvGetSize(dst_arr);  
15.     int cx,cy;  
16.     if ((dst_size.width!= size.width)||(dst_size.height!=size.height))  
17.     {  
18.         cvError(CV_StsUnmatchedSizes,"cvShiftDFT",  
19.             "Source and Destination arrays must have the same sizes",  
20.             __FILE__,__LINE__);  
21.     }  
22.     if (src_arr == dst_arr)  
23.     {  
24.         tmp = cvCreateMat(size.height/2,size.width/2,cvGetElemType(src_arr));  
25.     }  
26.     cx=size.width/2;//取出图像的原点  
27.     cy=size.height/2;  
28.     q1=cvGetSubRect(src_arr,&q1stub,cvRect(0,0,cx,cy));  
29.     //取出图像的第一象限，由q1指针指向它  
30.     q2=cvGetSubRect(src_arr,&q2stub,cvRect(cx,0,cx,cy));  
31.     //取出图像的第二象限，由q2指针指向它  
32.     q3=cvGetSubRect(src_arr,&q3stub,cvRect(cx,cy,cx,cy));  
33.     //取出图像的第三象限，由q3指针指向它  
34.     q4=cvGetSubRect(src_arr,&q4stub,cvRect(0,cy,cx,cy));  
35.     //取出图像的第四象限，由q4指针指向它  
36.      
37.     d1=cvGetSubRect(src_arr,&d1stub,cvRect(0,0,cx,cy));  
38.     d2=cvGetSubRect(src_arr,&d2stub,cvRect(cx,0,cx,cy));  
39.     d3=cvGetSubRect(src_arr,&d3stub,cvRect(cy,cy,cx,cy));  
40.     d4=cvGetSubRect(src_arr,&d4stub,cvRect(0,cy,cx,cy));  
41.     if (src_arr!=dst_arr)  
42.     {  
43.         if (!CV_ARE_TYPES_EQ(q1,d1))  
44.         {  
45.             cvError(CV_StsUnmatchedSizes,"cvShiftDFT",  
46.                 "Source and Destination arrays must have the same sizes",  
47.                 __FILE__,__LINE__);         
48.         }  
49.         //以图像中心为原点，调整傅里叶变换图像的四个象限区，  
50.         //即第一与第三象限交换，第二与第四象限交换  
51.         cvCopy(q3,d1,0);  
52.         cvCopy(q4,d2,0);  
53.         cvCopy(q1,d3,0);  
54.         cvCopy(q2,d4,0);  
55.     }  
56.     else  
57.     {//若源矩阵和目的矩阵相同则直接在源矩阵中进行操作  
58.         cvCopy(q3,tmp,0);  
59.         cvCopy(q1,q3,0);  
60.         cvCopy(tmp,q1,0);  
61.         cvCopy(q4,tmp,0);  
62.         cvCopy(q2,q4,0);  
63.         cvCopy(tmp,q2,0);  
64.     }  
65. }  
66. int main(int argc,char ** argv)  
67. {  
68.     const char* filename =(argc>=2?argv[1]:"lena.jpg");  
69.     IplImage *im;  
70.     IplImage *realInput,*imaginaryInput,*complexInput;  
71.     IplImage *image_Re,*image_Im;  
72.     int dft_M,dft_N;  
73.     CvMat *dft_A;  
74.     CvMat tmp;  
75.     double m,M;  
76.     im = cvLoadImage(filename, CV_LOAD_IMAGE_GRAYSCALE  );//加载图像  
77.     if (!im)  
78.     {  
79.         return -1;  
80.     }  
81.     //分配空间  
82.     realInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道  
83.     imaginaryInput =cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道  
84.     complexInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,2);//双通道  
85.     cvScale(im,realInput,1.0,0.0);  
86.     //#define cvScale cvConvertScale=>readInput=im  
87.     cvZero(imaginaryInput);  
88.     //清空这个图像的内容  
89.     cvMerge(realInput,imaginaryInput,NULL,NULL,complexInput);  
90.     //混合这两个图像作为complexInput的两个通道  
91.     /*得到最优DFT尺寸 */  
92.     dft_M = cvGetOptimalDFTSize(im->height-1);  
93.     dft_N = cvGetOptimalDFTSize(im->width-1);  
94.      
95.     dft_A = cvCreateMat(dft_M,dft_N,CV_64FC2);  
96.     image_Re = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//实部  
97.     image_Im = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//虚部  
98.     cvGetSubRect(dft_A,&tmp,cvRect(0,0,im->width,im->height));  
99.     cvCopy(complexInput,&tmp,NULL);  
100.     if( dft_A->cols > im->width )//若得到的最优DFT尺寸在宽度上大于原图，则重新取  
101.     {  
102.         cvGetSubRect(dft_A,&tmp,cvRect(im->width,0,dft_A->cols-im->width,im->height));  
103.         cvZero(&tmp);  
104.     }  
105.     cvDFT(dft_A,dft_A,CV_DXT_FORWARD,complexInput->height);  
106.      
107.     cvNamedWindow("win",0);  
108.     cvNamedWindow("magnitude",0);  
109.     cvShowImage("win",im);  
110.      
111.     //分割出实部和虚部  
112.     cvSplit(dft_A,image_Re,image_Im,0,0);  
113.      
114.     //计算功率谱 Mag=sqrt(Re^2+Im^2)  
115.     cvPow(image_Re,image_Re,2.0);  
116.     cvPow(image_Im,image_Im,2.0);  
117.     cvAdd(image_Re,image_Im,image_Re,NULL);//image_Re<=image_Re+image_Im  
118.     cvPow(image_Re,image_Re,0.5);  
119.     //计算log(1+Mag)  
120.     cvAddS(image_Re,cvScalarAll(1.0),image_Re,NULL);  
121.     cvLog(image_Re,image_Re);  
122.     cvShiftDFT(image_Re,image_Re);  
123.     cvMinMaxLoc(image_Re,&m,&M,NULL,NULL,NULL);  
124.     cvScale(image_Re,image_Re,1.0/(M-m),1.0*(-m)/(M-m));  
125.     cvShowImage("magnitude",image_Re);  
126.     cvWaitKey(-1);  
127.     return 0;  
128. }