二維離散傅立葉變換圖像頻域處理

vc++實現的傅立葉變換，參考算法導論中的快速傅立葉算法，使用openCV作圖像的讀入與顯示，加入了高斯函數低通濾波；

 

  1 /************************************************************
  2 *二維離散傅立葉變換和濾波函數（高斯低通）的實現
  3 *用於計算圖像離散序列的傅立葉頻譜圖，並用於頻域圖像處理
  4 *算法：多項式快速傅立葉算法（蝶形） 理論基礎：算法導論,圖像處理
  5 *時間複雜度：一維O(NlogN)，二維O(N^2)
  6 *工具：openCV讀取圖像與展現效果
  7 *版本：測試
  8 *@auto:Bruce mu
  9 ************************************************************/
 10 #include <stdio.h>
 11 #include <iostream>
 12 #include <complex>
 13 #include <bitset>
 14 #include <fstream>
 15 #include <algorithm>
 16 #include <iterator>
 17 #include <math.h>
 18 #include "cv.h"
 19 #include "highgui.h"
 20 #include "CImg.h"
 21 #define pi 3.1415926535
 22 #define DELTA 5.0
 23 using std::iostream;
 24 using std::bitset;
 25 using std::complex;
 26 using namespace std;
 27 using namespace cimg_library;
 28 
 29 /*******爲自底向上的迭代重排序列計算位置************/
 30 int rev(int k,int n)
 31 {
 32      bitset<32> tmp(k);
 33      bitset<32> dist;
 34      for(int m = 0;m<n;m++)
 35      {
 36         if(tmp.test(m))
 37         {
 38             dist.set(n-1-m);
 39         }
 40      }
 41      int revk = (int)dist.to_ulong();
 42      return revk;
 43 }
 44 //重載形式
 45 complex<double>* FFT(const complex<double> *srcimg,int n)
 46 {
 47     double flag = log10((double)n)/log10(2.0);
 48     int N = n;
 49     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
 50         cout <<"the length of srcimg is wrong"<< endl;
 51         /*填充成2的指數項*/
 52         cout <<"need padding"<<endl;
 53         N = pow(2,(double)((int)flag+1));
 54         flag = log10((double)N)/log10(2.0);
 55     }
 56     /*改變成奇偶順序*/
 57     complex<double> *arr= new complex<double>[N];
 58     int sub;
 59     for(int k =0;k<N;k++)
 60     {
 61         sub =rev(k,(int)flag); 
 62         if(sub <= n-1){
 63             arr[k] = *(srcimg + sub);
 64         }else{
 65             complex<double> t = complex<double>(0,0);
 66             arr[k] = t;
 67         }
 68     }
 69      for(int s =1;s <= flag;s++)
 70     {
 71         int m = pow(2,(double)s);
 72         complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：從W21開始,週期變換
 73         for(int p = 0;p<N;p+=m)
 74         {
 75             complex<double> w(1,0);
 76             for(int j = 0;j<=m/2-1;j++)
 77             {
 78                 complex<double> t = w * arr[p+j+m/2];
 79                 complex<double> u = arr[p+j];
 80                 arr[p+j] = u+t;
 81                 arr[p+j+m/2] = u-t;
 82                 w = w*wm;
 83             }
 84         }
 85     }
 86     return arr;
 87 }
 88 /***********一維快速傅立葉變換********************
 89 *srcimg : 原始一維序列                          *
 90 *n     ：一維序列的長度
 91 *************************************************/
 92 complex<double>* FFT(const double *srcimg,int n)
 93 {
 94     double flag = log10((double)n)/log10(2.0);
 95     int N = n;
 96     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
 97 //        cout <<"the length of srcimg is wrong"<< endl;
 98         /*填充成2的指數項*/
 99 //        cout <<"need padding"<<endl;
100         N = pow(2,(double)((int)flag+1));
101         flag = log10((double)N)/log10(2.0);
102     }
103     /*改變成奇偶順序*/
104     complex<double> *arr= new complex<double>[N];
105     int sub;
106     for(int k =0;k<N;k++)
107     {
108         sub =rev(k,(int)flag); 
109         if(sub <= n-1){
110             arr[k] = complex<double>(*(srcimg + sub),0);
111         }else{
112             complex<double> t = complex<double>(0,0);
113             arr[k] = t;
114         }
115     }
116 //    cout<<"------------after padding and retrival-----------------"<<endl;
117 //    for(size_t y=0;y<N;y++)
118 //    {
119 //        cout << arr[y].real()<<"  "<<arr[y].imag()<<endl;
120 //    }
121     /*基於迭代的蝶形快速傅立葉變換，自底向上*/
122      for(int s =1;s <= flag;s++)
123     {
124         int m = pow(2,(double)s);
125         complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：從W21開始,週期變換
126         for(int p = 0;p<N;p+=m)
127         {
128             complex<double> w(1,0);
129             for(int j = 0;j<=m/2-1;j++)
130             {
131                 complex<double> t = w * arr[p+j+m/2];
132                 complex<double> u = arr[p+j];
133                 arr[p+j] = u+t;
134                 arr[p+j+m/2] = u-t;
135                 w = w*wm;
136             }
137         }
138     }
139     return arr;
140 }
141 
142 int countPadding(int n);
143 /*****************一維快速傅立葉逆變換********************/
144 /*fftimg:原始一維傅立葉序列
145   n     : 序列長度
146 *******************************************************/
147 complex<double>* IFFT(const complex<double> *fftimg,int n)
148 {
149     n = countPadding(n);
150     double flag = log10((double)n)/log10(2.0);
151     int N = n;
152     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
153         cout <<"the length of srcimg is wrong"<< endl;
154         /*填充成2的指數項*/
155         cout <<"need padding"<<endl;
156         N = pow(2,(double)((int)flag+1));
157         flag = log10((double)N)/log10(2.0);
158     }
159     /*改變成奇偶順序*/
160     complex<double> * spit = new complex<double>[N];
161     int sub=0;
162     for(int k =0;k<N;k++)
163     {
164         sub =rev(k,(int)flag); 
165         if(sub < n){
166             spit[k] = complex<double>(*(fftimg + sub));
167         }else{
168             spit[k] = complex<double>(0,0);
169         }
170     }
171 
172     for(int s =1;s <= flag;s++)
173     {
174         int m = pow(2,(double)s);
175         complex<double> wm = complex<double>(cos(-2*pi/m),sin(-2*pi/m));//wm1：從W2(-1)開始
176         for(int p = 0;p<N;p+=m)
177         {
178             complex<double> w(1,0);
179             for(int j = 0;j<=m/2-1;j++)
180             {
181                 complex<double> t = w * spit[p+j+m/2];
182                 complex<double> u = spit[p+j];
183                 spit[p+j] = u+t;
184                 spit[p+j+m/2] = u-t;
185                 w = w*wm;
186             }
187         }
188     }
189 
190     for(size_t p =0;p<n;p++)
191     {
192         spit[p] = spit[p]/complex<double>(N,0);
193     }
194     return spit;
195 }
196   
197 /*******使用共軛的快速傅立葉逆變換**************/
198 complex<double>* IFFT2(const complex<double> *fftimg,int n)
199 {
200     n = countPadding(n);
201     complex<double> *gfftimg = new complex<double>[n];
202     for(size_t i = 0;i<n;i++){
203         gfftimg[i] = complex<double>(fftimg[i].real(),-fftimg[i].imag());
204     }
205     complex<double> *ifft = FFT(gfftimg,n); 
206     for(size_t j = 0;j<n;j++)
207     {
208         ifft[j] = ifft[j]/complex<double>(n,0);
209     }
210     delete gfftimg;
211     return ifft;
212 }
213 
214 /*************二維快速傅立葉變換**************************
215 *srcimg: 用一維表示的二維原始序列
216 *width ：寬度
217 *height：高度
218 ********************************************************/
219 complex<double>* twoDFFT(const double *srcimg,int width,int height)
220 {
221     int w = countPadding(width);
222     int h = countPadding(height);
223     int pixes = w * h;
224     complex<double> *hdirection = new complex<double>[w];
225     complex<double> *vdirection = new complex<double>[h];
226     complex<double> *fourier = new complex<double>[pixes];
227     /*二維水平方向*/
228     for(size_t i = 0;i<h;i++){
229         for(size_t j = 0;j<w;j++){
230             if(i>=height || j >=width){
231                 hdirection[j] = complex<double>(0,0);
232             }else{
233                 hdirection[j] = complex<double>(srcimg[i*width + j],0);
234             }
235     //        cout << hdirection[j] << " ";
236         }
237     //    cout <<""<<endl;
238         complex<double> *hfourier = FFT(hdirection,w);
239         for(size_t m = 0;m<w;m++){
240             fourier[i*w+m] = hfourier[m];
241         }
242         delete hfourier;
243     }
244     /*二維垂直方向*/
245     for(size_t ii = 0;ii<w;ii++){
246         for(size_t jj = 0;jj<h;jj++){
247             vdirection[jj] = fourier[jj*w + ii];
248         }
249         complex<double> *vfourier = FFT(vdirection,h);
250         for(size_t mm = 0;mm < h;mm++){
251             fourier[mm*w +ii] = vfourier[mm];
252         }
253         delete vfourier;
254     }
255     delete hdirection;
256     delete vdirection;
257     return fourier;
258 
259 }
260 /**************二維快速傅立葉逆變換*************************
261 *fourier : 一維表示的二維傅立葉變換序列                     *
262 *width   :寬                                             *
263 *height  :高                                             *
264 ***********************************************************/
265 complex<double>* twoDIFFT(const complex<double> *fourier,int width,int height)
266 {
267     width = countPadding(width);
268     height = countPadding(height);
269     int fpoints = width * height;
270     complex<double> *hdirection = new complex<double>[width];
271     complex<double> *vdirection = new complex<double>[height];
272     complex<double> *ifourier = new complex<double>[fpoints];
273 
274     for(size_t ii = 0;ii<height;ii++)
275     {
276         for(size_t jj = 0;jj<width;jj++){
277             hdirection[jj] = fourier[ii*width+jj];
278         }
279         complex<double> *hifour = IFFT(hdirection,width);//臨時變量
280         for(size_t mm = 0;mm<width;mm++){
281             ifourier[ii*width+mm] = hifour[mm];
282         }
283         delete hifour;
284     }
285     for(size_t i = 0;i<width;i++){
286         for(size_t j = 0;j<height;j++){
287             vdirection[j] = ifourier[j*width+i];
288         }
289         complex<double> *vifour = IFFT(vdirection,height);
290         for(size_t m = 0;m<height;m++){
291             ifourier[m*width+i] = vifour[m];
292         }
293         delete vifour;
294     }    
295     delete hdirection;
296     delete vdirection;
297     return ifourier;
298 }
299 /******************計算填充數********************************************
300 *蝶形傅立葉變換算法只計算2的指數次的離散序列，對於非2的指數次的序列，使用零填充*
301 ***********************************************************************/
302 inline int countPadding(int n)
303 {
304     double lg = log10((double)n)/log10(2.0);
305     if((lg - (int)lg) == 0){
306         return n;
307     }
308     int N = pow(2.0,((int)lg+1));
309     return N;
310 }
311 
312 /*****高斯低通濾波函數*************************
313 *src：   輸入頻譜
314 *width： 輸入頻譜寬度
315 *height：輸入頻譜高度
316 *D：     高斯函數方差，即濾波閾值
317 *只對於移位居中後的傅立葉頻譜進行高斯低通濾波
318 */
319 void guass(complex<double> *src,int width,int height,double D)
320 {
321     //find centre point
322     int orgx = floor(width/2.0);
323     int orgy = floor(height/2.0);
324     double distence = 0;
325     for(size_t i = 0;i<height;i++)
326     {
327         for(size_t j = 0;j<width;j++)
328         {
329             distence = sqrt(pow(abs((int)(i-orgy)),2.0)+pow(abs((int)(j-orgx)),2.0));
330             src[i*width+j] = src[i*width+j] * exp(-distence/(2*pow(D,2.0)));
331             
332         }
333         //cout<<i<<" "<<distence<<" "<<endl;
334         
335     }
336 //    cout<<orgx<<" "<<orgy<<endl;
337 }
338 /************複數傅立葉頻譜數組轉換成256級灰度數組*****************/
339 void toGray(complex<double> *twodfourier,int pwid,int phei,char* pixval)
340 {
341         double *vals = new double[pwid*phei];
342         double max = 0;
343         double min = 255;
344         for(size_t p = 0;p<pwid*phei;p++){
345             vals[p]=log(1+sqrt(pow(twodfourier[p].real(),2.0)+pow(twodfourier[p].imag(),2.0)));//對數級的幅度鋪
346             if(vals[p] > max){
347                 max = vals[p];
348             }
349             if(vals[p] < min){
350                 min = vals[p];
351             }
352         }
353         cout<<min<< " "<<max<<endl;
354         cout<<pwid<<" "<<phei<<endl;
355         for(size_t s = 0;s<pwid*phei;s++){
356             pixval[s] = (char)((vals[s]-min)/(max-min)*255);
357         }
358         delete vals;
359 }
360 /******************opencv 讀取圖像與展現效果***********************/
361 int main(int argc,char **argv)
362 {
363     IplImage *img;
364     if((img = cvLoadImage("F:\\Fig0222(a)(face).tif",0))!=0){
365         int dim = img->imageSize;
366         //從圖像矩陣複製出單維數組，並作頻譜居中操做，對應改寫接口，返回單維數組；
367         double * src = new double[dim];
368         size_t j =0;
369         for(int y =0;y<img->height;y++)
370         {
371             uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
372             for(int x =0;x<img->width;x++){
373                 /***輸入函數乘以（-1）的（x+y）次方，頻譜將居中*/
374                 src[j] = (double)ptr[x]*pow(-1,(double)(x+y))/256;
375                 j++;
376             }
377         }
378         int w = img->width;
379         int h = img->height;
380         int pwid = countPadding(w);
381         int phei = countPadding(h);
382 
383         complex<double> *twodfourier = twoDFFT(src,w,h);
384         char * pixval = new char[pwid*phei];
385         CvMat freq;
386         toGray(twodfourier,pwid,phei,pixval);//複數頻譜轉256灰度
387         cvInitMatHeader(&freq,pwid,phei,CV_8UC1,pixval);
388         IplImage *ipl = cvCreateImage(cvGetSize(&freq),8,1);
389         cvGetImage(&freq,ipl);//獲取頻譜圖像
390 
391         guass(twodfourier,pwid,phei,DELTA);//方差DELTA的高斯平滑（高斯低通濾波）
392         CvMat gaufre;
393         char *pixvals = new char[pwid*phei];
394         toGray(twodfourier,pwid,phei,pixvals);
395         cvInitMatHeader(&gaufre,pwid,phei,CV_8UC1,pixvals);
396         IplImage *gausf = cvCreateImage(cvGetSize(&gaufre),8,1);
397         cvGetImage(&gaufre,gausf);//高斯平滑後的頻譜圖像
398 
399         complex<double> *twodifourier = twoDIFFT(twodfourier,w,h);
400         double ipix = 0;
401         size_t po = 0;
402         for(int y =0;y<img->height;y++)
403         {
404             uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
405             for(int x =0;x<img->width;x++){
406                 ipix = twodifourier[po*pwid+x].real();
407                 ptr[x] = (uchar)(ipix * 256)*pow(-1,(double)(x+y));
408             }
409             po++;
410         }
411         cvNamedWindow("frequence_domain",CV_WINDOW_AUTOSIZE);
412         cvShowImage("frequence_domain",ipl);
413         cvNamedWindow("gauss_fre",CV_WINDOW_AUTOSIZE);
414         cvShowImage("gauss_fre",gausf);
415         cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE);
416         cvShowImage("fourier_t",img);
417         cvWaitKey(0);
418         cvReleaseImage(&ipl);
419         cvReleaseImage(&gausf);
420         cvReleaseImage(&img);
421         cvDestroyWindow("fourier_t");
422         cvDestroyWindow("frequence_domain");
423         cvDestroyWindow("gauss_fre");
424         delete pixval;
425         delete pixvals;
426         delete src;
427         delete twodfourier;
428         delete twodifourier;
429         return 1;
430     }
431     return 0;
432 }

顯示輸出：

原圖

A	B
C	D

A爲原圖，B爲使用均值爲0，方差爲5.0的高斯函數作低通濾波後的效果圖。C爲原圖的傅立葉頻譜圖（已居中顯示），D爲使用方差5.0的高斯濾波後的頻譜圖；