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

fftwlib.cpp

00001 /* Copyright (c) 2001 C. Grigorescu */
00002 
00003 #include <stdlib.h>
00004 #include <stdio.h>
00005 #include <string.h>
00006 #include <math.h>
00007 #include <config.h>
00008 
00009 #ifdef HAVE_LIBFFTW
00010 
00011 #include <fftw.h>
00012 #include <rfftw.h>
00013 #include "fftwlib.h"
00014 
00015 /*   Wrap code for performing 2D Fourier transforms and convolution
00016      with in single precision floating point by using FFTW
00017 */
00018 
00019 void _convolve(float *input_image_1, 
00020               float *input_image_2, 
00021               int size_x, int size_y, 
00022               float *convolved_image)
00023 {
00024 
00025   fftw_real a[size_x][2*(size_y/2+1)], b[size_x][2*(size_y/2+1)],c[size_x][size_y];
00026   fftw_complex *A, *B, C[size_x][size_y/2+1];
00027   fftwnd_plan p, pinv;
00028   fftw_real scale = 1.0 / (size_x * size_y);
00029   int i, j, ij;
00030 
00031   p = rfftw2d_create_plan(size_x, size_y, 
00032           FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
00033   pinv = rfftw2d_create_plan(size_x, size_y, 
00034           FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00035 
00036   half_shift(input_image_1,size_x,size_y); 
00037 
00038 
00039   A = (fftw_complex *) &a[0][0];
00040   B = (fftw_complex *) &b[0][0];
00041 
00042   for(i=0; i<size_x; ++i)
00043     for(j=0; j<size_y; ++j)
00044       {
00045       a[i][j] = input_image_1[j*size_x+i];
00046       b[i][j] = input_image_2[j*size_x+i];
00047       }
00048 
00049   rfftwnd_one_real_to_complex(p,&a[0][0],NULL);
00050   rfftwnd_one_real_to_complex(p,&b[0][0],NULL);
00051 
00052   for(i=0;i<size_x;++i)
00053     for(j=0;j<(size_y/2+1);++j)
00054       {
00055       ij = i*(size_y/2+1)+j;
00056       C[i][j].re = (A[ij].re * B[ij].re -
00057                     A[ij].im * B[ij].im)*scale;
00058 
00059       C[i][j].im = (A[ij].re * B[ij].im +
00060                     A[ij].im * B[ij].re)*scale;
00061       }
00062   
00063   rfftwnd_one_complex_to_real(pinv,&C[0][0],&c[0][0]);
00064 
00065   for(i=0; i<size_x; ++i)
00066     for(j=0; j<size_y; ++j)
00067       convolved_image[j*size_x+i] = (float)c[i][j];
00068 
00069   half_shift(input_image_1,size_x,size_y); 
00070 
00071   rfftwnd_destroy_plan(p);
00072   rfftwnd_destroy_plan(pinv);
00073 }
00074 
00075 void _fft2d(float *real_input_image, float *imag_input_image,
00076            int size_x, int size_y, 
00077            float *fourier_real_part, float *fourier_imag_part)
00078 {
00079 
00080   fftw_complex a[size_x*size_y];
00081   fftwnd_plan p;
00082   fftw_real scale;
00083   int i, j, ij; 
00084 
00085   scale = 1.0 / (size_x * size_y);
00086 
00087   p = fftw2d_create_plan(size_x, size_y, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00088 
00089   half_shift(real_input_image,size_x,size_y); 
00090   if (imag_input_image != NULL)
00091     half_shift(imag_input_image,size_x,size_y); 
00092 
00093   for(j=0; j<size_y; ++j)
00094     for(i=0; i<size_x; ++i)
00095       {
00096       ij = i*size_y+j;
00097       a[ij].re = real_input_image[j*size_x+i];
00098       if (imag_input_image != NULL)
00099         a[ij].im = imag_input_image[j*size_x+i];
00100       else
00101         a[ij].im = 0.0;
00102       }
00103 
00104   fftwnd_one(p,&a[0],NULL);
00105 
00106   for(j=0; j<size_y; ++j)
00107     for(i=0; i<size_x; ++i)
00108       {
00109       ij = i*size_y+j;
00110       fourier_real_part[j*size_x+i] = (float) a[ij].re;
00111       fourier_imag_part[j*size_x+i] = (float) a[ij].im;
00112       }
00113 
00114   half_shift(fourier_real_part,size_x,size_y);
00115   half_shift(fourier_imag_part,size_x,size_y); 
00116   fftwnd_destroy_plan(p);
00117 } 
00118 
00119 void _ifft2d(float *fourier_real_image, float *fourier_imag_image,
00120            int size_x, int size_y, float *real_image, float *imag_image)
00121 {
00122   /*  Performs a IFFT; works with FFTW package */
00123 
00124   fftw_complex a[size_x*size_y];
00125   fftwnd_plan p;
00126   fftw_real scale;
00127   int i, j, ij; 
00128 
00129   scale = 1.0 / (size_x * size_y);
00130 
00131   p = fftw2d_create_plan(size_x, size_y, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00132 
00133   half_shift(fourier_real_image,size_x,size_y);
00134   half_shift(fourier_imag_image,size_x,size_y);
00135 
00136   for(j=0; j<size_y; ++j)
00137     for(i=0; i<size_x; ++i)
00138       {
00139       ij = i*size_y+j;
00140       a[ij].re = (double) fourier_real_image[j*size_x+i];
00141       a[ij].im = (double) fourier_imag_image[j*size_x+i];
00142       }
00143 
00144   fftwnd_one(p,&a[0],NULL);
00145 
00146   for(j=0; j<size_y; ++j)
00147     for(i=0; i<size_x; ++i)
00148       {
00149       ij = i*size_y+j;
00150       real_image[j*size_x+i] = (float) a[ij].re * scale;
00151       imag_image[j*size_x+i] = (float) a[ij].im * scale;
00152       }
00153 
00154   half_shift(real_image,size_x,size_y);
00155   half_shift(imag_image,size_x,size_y);
00156   fftwnd_destroy_plan(p);
00157 } 
00158 
00159 #endif