dft.c

   1 /* dft.c - discrete fourier transformation program */
   2
   3 /*
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 //#define USE_FFTW3
   9
  10 #include "dft.h"
  11
  12 int main(int argc,char **argv) {
  13
  14 #ifdef USE_FFTW3
  15   fftw_plan plan;
  16   fftw_complex *in,*out;
  17 #else
  18   t_fourier fourier;
  19 #endif
  20   t_bmp src;
  21   t_bmp dst;
  22   t_bmp cut;
  23   int x,y;
  24   int offy,offt;
  25   double scale;
  26   int size;
  27   double max;
  28   double *mag;
  29
  30   bmp_init(&src,1);
  31   bmp_init(&dst,1);
  32   bmp_init(&cut,1);
  33   src.mode=READ;
  34   dst.mode=WRITE;
  35   strcpy(src.file,argv[1]);
  36   strcpy(dst.file,argv[2]);
  37
  38 #ifndef USE_FFTW3
  39   fourier_init(&fourier,1);
  40   fourier.type=DFT|FWD;
  41   fourier.dim=2;
  42 #endif
  43
  44   bmp_read_file(&src);
  45
  46   bmp_cut_grab_bottom(&cut,&src,src.info.width,GRAB);
  47
  48   dst.width=cut.width;
  49   dst.height=cut.height;
  50   bmp_alloc_map(&dst);
  51
  52   size=cut.width*cut.height;
  53
  54 #ifndef USE_FFTW3
  55   fourier.data_len[0]=cut.width;
  56   fourier.data_len[1]=cut.height;
  57   fourier_alloc_data(&fourier);
  58 #else
  59   in=fftw_malloc(sizeof(fftw_complex)*size);
  60   memset(in,0,size*sizeof(fftw_complex));
  61   out=fftw_malloc(sizeof(fftw_complex)*size);
  62   memset(out,0,size*sizeof(fftw_complex));
  63   if(in==NULL||out==NULL) {
  64     printf("fftw alloc error!\n");
  65     return -23;
  66   }
  67   plan=fftw_plan_dft_2d(cut.width,cut.height,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
  68 #endif
  69
  70   mag=(double *)malloc(size*sizeof(double));
  71   if(mag==NULL) {
  72     printf("unable to alloc mag memory\n");
  73     return -1;
  74   }
  75
  76   // copy the data
  77   for(x=0;x<size;x++) {
  78 #ifdef USE_FFTW3
  79     in[x][0]=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
  80     in[x][1]=0;
  81 #else
  82     fourier.data[x].r=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
  83 #endif
  84   }
  85
  86   // do the fourier trafo
  87 #ifdef USE_FFTW3
  88   fftw_execute(plan);
  89 #else
  90   fourier_dft_2d(&fourier);
  91 #endif
  92
  93
  94   printf("fourier done!\n");
  95
  96   // copy back the data, intensity = sqrt(r^2+i^2)
  97   max=0;
  98   offt=0;
  99   for(y=0;y<cut.height;y++) {
 100     for(x=0;x<cut.width;x++) {
 101 #ifdef USE_FFTW3
 102       mag[offt]=sqrt(out[offt][0]*out[offt][0]+out[offt][1]*out[offt][1]);
 103 #else
 104       mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
 105 #endif
 106       if(mag[offt]>max) max=mag[offt];
 107       offt+=1;
 108     }
 109   }
 110
 111   printf("found max: %f\n",max);
 112
 113   if(max!=0) scale=255.0/max;
 114   else scale=0;
 115
 116   if(scale!=0) {
 117     printf("scaling image intensity: %f\n",scale);
 118     offt=0;
 119     for(y=0;y<dst.height;y++) {
 120       for(x=0;x<dst.width;x++) {
 121         // usual image processing order
 122         offy=((y<dst.height/2)?y+dst.height/2:y-dst.height/2);
 123         offy*=dst.width;
 124         offy+=((x<dst.width/2)?x+dst.width/2:x-dst.width/2);
 125         //offy=offt;
 126         dst.map[offy].r=scale*mag[offt];
 127         dst.map[offy].g=dst.map[offy].r;
 128         dst.map[offy].b=dst.map[offy].r;
 129         offt+=1;
 130       }
 131     }
 132   }
 133
 134   free(mag);
 135
 136   bmp_write_file(&dst);
 137
 138 #ifdef USE_FFTW3
 139   fftw_destroy_plan(plan);
 140   fftw_free(in); fftw_free(out);
 141 #else
 142   fourier_shutdown(&fourier);
 143 #endif
 144
 145   bmp_shutdown(&src);
 146   bmp_shutdown(&dst);
 147
 148   printf("done ...\n");
 149
 150   return 1;
 151 }