Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   Related Pages   Examples  

FourierTransform.cpp

00001 /* Copyright (c) 2001 C. Grigorescu */
00002 
00003 #include <math.h>
00004 #include <sys/wait.h>  
00005 #include <signal.h> //only for kill
00006 #include "misc.h"
00007 #include "config.h"
00008 #include "FourierTransform.h"
00009 
00010 #ifdef HAVE_LIBFFTW
00011 #include "fftwlib.h"
00012 #endif
00013 
00014 void FourierTransform::_fft(float *data, int nn, int isign)
00015 {
00016   int n, mmax, m, j, istep, i;
00017   float wtemp, wr, wpr, wpi, wi, theta;
00018   float tempr, tempi;
00019 
00020   n = nn<<1;
00021   j = 1;
00022   for (i=1;i<n;i+=2) {
00023     if (j>i) {
00024       SWAP(data[j],data[i]);
00025       SWAP(data[j+1],data[i+1]);
00026     }
00027     m = n>>1;
00028     while((m>=2) && (j>m)) {
00029       j -= m;
00030       m >>= 1;
00031     }
00032     j+=m;
00033   }
00034 
00035   mmax = 2; 
00036   while( n > mmax) {
00037     istep = mmax<<1;
00038     theta = isign*(6.28318530/mmax);
00039     wtemp = sin(0.5*theta);
00040     wpr = -2.0*wtemp*wtemp;
00041     wpi = sin(theta);
00042     wr = 1.0;
00043     wi = 0.0;
00044     for (m=1;m<mmax;m+=2) {
00045       for (i=m;i<=n;i+=istep) {
00046         j=i+mmax;
00047         tempr = wr*data[j]-wi*data[j+1];
00048         tempi = wr*data[j+1]+wi*data[j];
00049         data[j] = data[i]-tempr;
00050         data[j+1] = data[i+1]-tempi;
00051         data[i] += tempr;
00052         data[i+1] += tempi;
00053       }
00054       wr = (wtemp=wr)*wpr-wi*wpi+wr;
00055       wi = wi*wpr+wtemp*wpi+wi;
00056     }
00057     mmax=istep;
00058   }
00059 }
00060 
00061 FourierTransform::FourierTransform(FloatImage &input)
00062 {
00063 #ifdef HAVE_LIBFFTW
00064   input.makeEvenSized();
00065   real_part = input;
00066   imag_part.makeImage(input.getWidth(),input.getHeight());
00067   imag_part.clearImage(0);
00068 #else
00069   input.makeEvenSized();
00070   real_part = input;
00071   int w = real_part.getWidth();
00072   int h = real_part.getHeight();
00073   int p2_w = 1; int p2_h = 1;
00074 
00075   for(int i=0;i<16;i++) {
00076     if (w <= p2_w) break;
00077     p2_w = p2_w<<1;
00078   }
00079   for(int i=0;i<16;i++) {
00080     if (h <= p2_h) break;
00081     p2_h = p2_h<<1;
00082   }
00083   real_part.pad((p2_w - w)/2, (p2_h - h)/2, ZEROS);
00084   imag_part.makeImage(real_part.getWidth(),real_part.getHeight());
00085   imag_part.clearImage(0);
00086 #endif
00087 }
00088 
00089 FourierTransform::FourierTransform(FloatImage &input_real, FloatImage &input_imag)
00090 {
00091   if ((input_real.getWidth() != input_imag.getWidth()) ||
00092       (input_real.getHeight() != input_imag.getHeight())) {
00093     ERROR("Real and imaginary input images have different sizes");
00094   }
00095   else {
00096 #ifdef HAVE_LIBFFTW
00097     input_real.makeEvenSized();
00098     input_imag.makeEvenSized();
00099     real_part = input_real;
00100     imag_part = input_imag;
00101 #else
00102     int w = input_real.getWidth();
00103     int h = input_imag.getHeight();
00104     int p2_w = 1; int p2_h = 1;
00105     
00106     int flag_w = 0, flag_h =0;
00107     for(int i=0;i<16;i++) {
00108       if (w == p2_w) {
00109         flag_w = 1; 
00110         break;
00111       }
00112       p2_w = p2_w<<1;
00113     }
00114     if (!flag_w) 
00115       ERROR("Image width not a power of 2; results may be incorrect");
00116     for(int i=0;i<16;i++) {
00117       if (h == p2_h) {
00118         flag_h = 1;
00119         break;
00120       }
00121       p2_h = p2_h<<1;
00122     }
00123     if (!flag_h) 
00124       ERROR("Image height not a power of 2; results may be incorrect");
00125     real_part = input_real;
00126     imag_part = input_imag;
00127 #endif
00128   }
00129 }
00130 
00131 void FourierTransform::fft(void)
00132 {
00133 
00134 #ifdef HAVE_LIBFFTW
00135    FloatImage temp_real = real_part;
00136    FloatImage temp_imag = imag_part;
00137    _fft2d(temp_real.getPixels(), temp_imag.getPixels(), 
00138           real_part.getWidth() , real_part.getHeight(),
00139           real_part.getPixels(), imag_part.getPixels());
00140 #else
00141   int w =  real_part.getWidth();
00142   int h =  real_part.getHeight();
00143 
00144   float * temp = (float *) malloc(2*w* sizeof(float));
00145   for (int i=0;i<h;i++) {
00146     for(int j=0;j<w;j++) {
00147       temp[2*j] = real_part[i][j];
00148       temp[2*j+1] = imag_part[i][j];
00149     }
00150     _fft(temp-1, w, 1);
00151     for(int j=0;j<w;j++) {
00152       real_part[i][j] = temp[2*j];
00153       imag_part[i][j] = temp[2*j+1];
00154     }
00155   }
00156   free(temp);
00157   temp = (float *) malloc(2*h * sizeof(float));
00158   for(int j=0;j<w;j++) {
00159     for (int i=0;i<h;i++) {
00160       temp[2*i] = real_part[i][j];
00161       temp[2*i+1] = imag_part[i][j];
00162     }
00163     _fft(temp-1, h, 1);
00164     for(int i=0;i<h;i++) {
00165       real_part[i][j] = temp[2*i];
00166       imag_part[i][j] = temp[2*i+1];
00167     }
00168   }
00169   free(temp);
00170   for(int j=0;j<w;j++) 
00171     for(int i=0;i<h;i++) {
00172       real_part[i][j] /= w*h;
00173       imag_part[i][j] /= w*h;
00174     }
00175   HALF_SHIFT(real_part.getPixels(),w,h);
00176   HALF_SHIFT(imag_part.getPixels(),w,h);
00177 #endif
00178 }
00179 
00180 void FourierTransform::ifft(void)
00181 {
00182 #ifdef HAVE_LIBFFTW
00183   FloatImage temp_real = real_part;
00184   FloatImage temp_imag = imag_part;
00185   _ifft2d(temp_real.getPixels(), temp_imag.getPixels(), 
00186           real_part.getWidth() , real_part.getHeight(),
00187           real_part.getPixels(), imag_part.getPixels());
00188 #else
00189   int w = real_part.getWidth();
00190   int h = real_part.getHeight();
00191 
00192   HALF_SHIFT(real_part.getPixels(),w,h);
00193   HALF_SHIFT(imag_part.getPixels(),w,h);
00194   float *temp = (float *) malloc(2*w * sizeof(float));
00195   for (int i=0;i<h;i++) {
00196     for(int j=0;j<w;j++) {
00197       temp[2*j] = real_part[i][j];
00198       temp[2*j+1] = imag_part[i][j];
00199     }
00200     _fft(temp-1, w, -1);
00201     for(int j=0;j<w;j++) {
00202       real_part[i][j] = temp[2*j];
00203       imag_part[i][j] = temp[2*j+1];
00204     }
00205   }
00206   free(temp);
00207   temp = (float *) malloc(2*h * sizeof(float));
00208   for(int j=0;j<w;j++) {
00209     for (int i=0;i<h;i++) {
00210       temp[2*i] = real_part[i][j];
00211       temp[2*i+1] = imag_part[i][j];
00212     }
00213     _fft(temp-1, h, -1);
00214     for(int i=0;i<h;i++) {
00215       real_part[i][j] = temp[2*i];
00216       imag_part[i][j] = temp[2*i+1];
00217     }
00218   }
00219   free(temp);
00220 #endif
00221 }
00222 
00223 void FourierTransform::showMagnitude(int log_values)
00224 {
00225   FloatImage magnitude(real_part.getWidth(), real_part.getHeight());
00226 
00227   for (int j=0;j<real_part.getWidth();j++)
00228     for (int i=0;i<real_part.getHeight();i++) 
00229       magnitude[i][j] = sqrt(real_part[i][j]*real_part[i][j] +
00230                              imag_part[i][j]*imag_part[i][j]);
00231   if (log_values) {
00232     for (int j=0;j<real_part.getWidth();j++)
00233       for (int i=0;i<real_part.getHeight();i++) 
00234         magnitude[i][j] = log(magnitude[i][j]);
00235   } 
00236   if(pid_magnitude) {
00237     kill(pid_magnitude, SIGUSR1); 
00238     waitpid(pid_magnitude, NULL, 0);
00239     pid_magnitude = magnitude.showImage();
00240   }
00241   else
00242     pid_magnitude = magnitude.showImage();
00243 }
00244 
00245 void FourierTransform::showPhase(int log_values)
00246 {
00247   FloatImage phase(imag_part.getWidth(), imag_part.getHeight());
00248   for (int j=0;j<imag_part.getWidth();j++)
00249     for (int i=0;i<imag_part.getHeight();i++)
00250       phase[i][j] = atan(imag_part[i][j]/real_part[i][j]);
00251   if (log_values) {
00252     for (int j=0;j<imag_part.getWidth();j++)
00253       for (int i=0;i<imag_part.getHeight();i++) 
00254         phase[i][j] = log(phase[i][j]);
00255   } 
00256   if(pid_phase) {
00257     kill(pid_phase, SIGUSR1); 
00258     waitpid(pid_phase, NULL, 0);
00259     pid_phase = phase.showImage();
00260   }
00261   else
00262     pid_phase = phase.showImage();
00263 }
00264 
00265 void FourierTransform::closeWindows()
00266 {
00267   if(pid_magnitude) {
00268     kill(pid_magnitude, SIGUSR1); 
00269     waitpid(pid_magnitude, NULL, 0);
00270     pid_magnitude = 0;
00271   }
00272   if(pid_phase) {
00273     kill(pid_phase, SIGUSR1); 
00274     waitpid(pid_phase, NULL, 0);
00275     pid_phase = 0;
00276   }
00277 }
00278 
00279 void  fft(FloatImage &input, FloatImage &output_real, FloatImage &output_imag)
00280 {
00281   FourierTransform a(input);
00282   a.fft();
00283   output_real = a.getReal();
00284   output_imag = a.getImag();
00285 }
00286 
00287 void ifft(FloatImage &input_real, FloatImage &input_imag, FloatImage &output)
00288 {
00289   FourierTransform a(input_real, input_imag);
00290   a.ifft();
00291   output = a.getReal();
00292 }
00293