opencv 中 傅里叶变换 FFT

[cpp]  view plain copy

1. void fft2(IplImage *src, IplImage *dst)  
2. {   //实部、虚部  
3.     IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;  
4.     //   int i, j;  
5.     image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //实部  
6.     //Imaginary part  
7.     image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虚部  
8.     //2 channels (image_Re, image_Im)  
9.     Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);  
10.     // Real part conversion from u8 to 64f (double)  
11.     cvConvertScale(src, image_Re);  
12.     // Imaginary part (zeros)  
13.     cvZero(image_Im);  
14.     // Join real and imaginary parts and stock them in Fourier image  
15.     cvMerge(image_Re, image_Im, 0, 0, Fourier);  
16.   
17.     // Application of the forward Fourier transform  
18.     cvDFT(Fourier, dst, CV_DXT_FORWARD);  
19.     cvReleaseImage(&image_Re);  
20.     cvReleaseImage(&image_Im);  
21.     cvReleaseImage(&Fourier);  
22. }  
23.   
24. void fft2shift(IplImage *src, IplImage *dst)  
25. {  
26.     IplImage *image_Re = 0, *image_Im = 0;  
27.     int nRow, nCol, i, j, cy, cx;  
28.     double scale, shift, tmp13, tmp24;  
29.     image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  
30.     //Imaginary part  
31.     image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  
32.     cvSplit( src, image_Re, image_Im, 0, 0 );  
33.     //具体原理见冈萨雷斯数字图像处理p123  
34.     // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)  
35.     //计算傅里叶谱  
36.     cvPow( image_Re, image_Re, 2.0);  
37.     cvPow( image_Im, image_Im, 2.0);  
38.     cvAdd( image_Re, image_Im, image_Re);  
39.     cvPow( image_Re, image_Re, 0.5 );  
40.     //对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射  
41.     //一宽带输出值，具体可见冈萨雷斯数字图像处理p62)  
42.     // Compute log(1 + Mag);  
43.     cvAddS( image_Re, cvScalar(1.0), image_Re ); // 1 + Mag  
44.     cvLog( image_Re, image_Re ); // log(1 + Mag)  
45.   
46.     //Rearrange the quadrants of Fourier image so that the origin is at the image center  
47.     nRow = src->height;  
48.     nCol = src->width;  
49.     cy = nRow/2; // image center  
50.     cx = nCol/2;  
51.     //CV_IMAGE_ELEM为OpenCV定义的宏，用来读取图像的像素值，这一部分就是进行中心变换  
52.     for( j = 0; j < cy; j++ ){  
53.         for( i = 0; i < cx; i++ ){  
54.             //中心化，将整体份成四块进行对角交换  
55.             tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);  
56.             CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(  
57.                 image_Re, double, j+cy, i+cx);  
58.             CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;  
59.   
60.             tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);  
61.             CV_IMAGE_ELEM( image_Re, double, j, i+cx) =  
62.                 CV_IMAGE_ELEM( image_Re, double, j+cy, i);  
63.             CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;  
64.         }  
65.     }  
66.     //归一化处理将矩阵的元素值归一为[0,255]  
67.     //[(f(x,y)-minVal)/(maxVal-minVal)]*255  
68.     double minVal = 0, maxVal = 0;  
69.     // Localize minimum and maximum values  
70.     cvMinMaxLoc( image_Re, &minVal, &maxVal );  
71.     // Normalize image (0 - 255) to be observed as an u8 image  
72.     scale = 255/(maxVal - minVal);  
73.     shift = -minVal * scale;  
74.     cvConvertScale(image_Re, dst, scale, shift);  
75.     cvReleaseImage(&image_Re);  
76.     cvReleaseImage(&image_Im);  
77. }  
78.   
79. void CCVMFCView::OnFuliyeTransform()  
80. {  
81.     IplImage *src;  
82.     IplImage *Fourier;   //傅里叶系数  
83.     IplImage *dst ;  
84.     IplImage *ImageRe;  
85.     IplImage *ImageIm;  
86.     IplImage *Image;  
87.     IplImage *ImageDst;  
88.     double m,M;  
89.     double scale;  
90.     double shift;  
91.     //src = workImg;  
92.     if(workImg->nChannels==3)  
93.         OnColorToGray();  
94.     src=cvCreateImage(cvGetSize(workImg),IPL_DEPTH_64F,workImg->nChannels);  //源图像  
95.     imageClone(workImg,&src);  
96.     cvFlip(src);  
97.   
98.     Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);  
99.     dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);  
100.     ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
101.     ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
102.     Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);  
103.     ImageDst = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);  
104.     fft2(src,Fourier);                  //傅里叶变换  
105.     fft2shift(Fourier, Image);          //中心化  
106.     cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换，并对结果进行缩放  
107.     cvSplit(dst,ImageRe,ImageIm,0,0);  
108.   
109.     cvNamedWindow("源图像",0);  
110.     cvShowImage("源图像",src);               
111.     //对数组每个元素平方并存储在第二个参数中  
112.     cvPow(ImageRe,ImageRe,2);                 
113.     cvPow(ImageIm,ImageIm,2);  
114.     cvAdd(ImageRe,ImageIm,ImageRe,NULL);  
115.     cvPow(ImageRe,ImageRe,0.5);  
116.     cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);  
117.     scale = 255/(M - m);  
118.     shift = -m * scale;  
119.     //将shift加在ImageRe各元素按比例缩放的结果上，存储为ImageDst  
120.     cvConvertScale(ImageRe,ImageDst,scale,shift);  
121.   
122.     cvNamedWindow("傅里叶谱",0);  
123.     cvShowImage("傅里叶谱",Image);  
124.     cvNamedWindow("傅里叶逆变换",0);  
125.     cvShowImage("傅里叶逆变换",ImageDst);  
126.     //释放图像  
127.     cvWaitKey(10000);  
128.     cvReleaseImage(&src);  
129.     cvReleaseImage(&Image);  
130.     cvReleaseImage(&ImageIm);  
131.     cvReleaseImage(&ImageRe);  
132.     cvReleaseImage(&Fourier);  
133.     cvReleaseImage(&dst);  
134.     cvReleaseImage(&ImageDst);  
135.     Invalidate();  
136. }