dft.c

   1 /* dft.c - discrete fourier transformation program */
   2
   3 /*
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 #include "dft.h"
   9
  10 int main(int argc,char **argv) {
  11
  12   t_fourier fourier;
  13   t_bmp src;
  14   t_bmp dst;
  15   t_bmp cut;
  16   int x,y;
  17   int offy,offt;
  18   int scale;
  19   double max;
  20   double *mag;
  21
  22   bmp_init(&src,1);
  23   bmp_init(&dst,1);
  24   bmp_init(&cut,1);
  25   src.mode=READ;
  26   dst.mode=WRITE;
  27   strcpy(src.file,argv[1]);
  28   strcpy(dst.file,argv[2]);
  29
  30   fourier_init(&fourier,1);
  31   fourier.type=DFT|FWD;
  32   fourier.dim=2;
  33
  34   bmp_read_file(&src);
  35
  36   bmp_cut_grab_bottom(&cut,&src,src.info.width,GRAB);
  37
  38   dst.width=cut.width;
  39   dst.height=cut.height;
  40   bmp_alloc_map(&dst);
  41
  42   fourier.data_len[0]=cut.width;
  43   fourier.data_len[1]=cut.height;
  44   fourier_alloc_data(&fourier);
  45
  46   mag=(double *)malloc(fourier.data_len[0]*fourier.data_len[1]*sizeof(double));
  47   if(mag==NULL) {
  48     printf("unable to alloc mag memory\n");
  49     return -1;
  50   }
  51
  52   // copy the data
  53   offy=0;
  54   for(y=0;y<fourier.data_len[1];y++) {
  55     for(x=0;x<fourier.data_len[0];x++) {
  56       offt=offy+x;
  57       fourier.data[offt].r=(src.map[offt].r+src.map[offt].g+src.map[offt].b)/3;
  58     }
  59     offy+=fourier.data_len[0];
  60   }
  61
  62   // do the fourier trafo
  63   fourier_dft_2d(&fourier);
  64
  65   printf("fourier done!\n");
  66
  67   // copy back the data, intensity = sqrt(r^2+i^2)
  68   max=0;
  69   offt=0;
  70   for(y=0;y<fourier.data_len[1];y++) {
  71     for(x=0;x<fourier.data_len[0];x++) {
  72       mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
  73       if(mag[offt]>max) max=mag[offt];
  74       offt+=1;
  75     }
  76   }
  77
  78   printf("found max: %f\n",max);
  79
  80   if(max!=0) scale=(int)255/max;
  81   else scale=0;
  82
  83   if(scale!=0) {
  84     printf("scaling image intensity: %d\n",scale);
  85     offt=0;
  86     for(y=0;y<dst.height;y++) {
  87       for(x=0;x<dst.width;x++) {
  88         // usual image processing order
  89         offy=((y<dst.height/2)?y+dst.height/2:y-dst.height/2);
  90         offy*=dst.width;
  91         offy+=((x<dst.width/2)?x+dst.width/2:x-dst.width/2);
  92         //offy=offt;
  93         dst.map[offy].r=scale*mag[offt];
  94         dst.map[offy].g=dst.map[offy].r;
  95         dst.map[offy].b=dst.map[offy].r;
  96         offt+=1;
  97       }
  98     }
  99   }
 100
 101   free(mag);
 102
 103   bmp_write_file(&dst);
 104
 105   fourier_shutdown(&fourier);
 106
 107   bmp_shutdown(&src);
 108   bmp_shutdown(&dst);
 109
 110   printf("done ...\n");
 111
 112   return 1;
 113 }