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   cut.mode=WRITE;
  28   strcpy(src.file,argv[1]);
  29   strcpy(dst.file,argv[2]);
  30   strcpy(cut.file,argv[3]);
  31
  32   fourier_init(&fourier,1);
  33   fourier.type=DFT|FWD;
  34   fourier.dim=2;
  35
  36   bmp_read_file(&src);
  37   bmp_cut_bottom(&cut,&src,150);
  38   bmp_write_file(&cut);
  39
  40   dst.width=src.info.width;
  41   dst.height=src.info.height;
  42   bmp_alloc_map(&dst);
  43
  44   fourier.data_len[0]=src.info.width;
  45   fourier.data_len[1]=src.info.height;
  46   fourier_alloc_data(&fourier);
  47
  48   mag=(double *)malloc(fourier.data_len[0]*fourier.data_len[1]*sizeof(double));
  49   if(mag==NULL) {
  50     printf("unable to alloc mag memory\n");
  51     return -1;
  52   }
  53
  54   // copy the data
  55   offy=0;
  56   for(y=0;y<fourier.data_len[1];y++) {
  57     for(x=0;x<fourier.data_len[0];x++) {
  58       offt=offy+x;
  59       fourier.data[offt].r=(src.map[offt].r+src.map[offt].g+src.map[offt].b)/3;
  60     }
  61     offy+=fourier.data_len[0];
  62   }
  63
  64   // do the fourier trafo
  65   fourier_dft_2d(&fourier);
  66
  67   printf("fourier done!\n");
  68
  69   // copy back the data, intensity = sqrt(r^2+i^2)
  70   max=0;
  71   offt=0;
  72   for(y=0;y<fourier.data_len[1];y++) {
  73     for(x=0;x<fourier.data_len[0];x++) {
  74       mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
  75       if(mag[offt]>max) max=mag[offt];
  76       offt+=1;
  77     }
  78   }
  79
  80   printf("found max: %f\n",max);
  81
  82   if(max!=0) scale=(int)255/max;
  83   else scale=0;
  84
  85   if(scale!=0) {
  86     printf("scaling image intensity: %d\n",scale);
  87     offt=0;
  88     for(y=0;y<dst.height;y++) {
  89       for(x=0;x<dst.width;x++) {
  90         offy=((y<dst.height/2)?y+dst.height/2:y-dst.height/2);
  91         offy*=dst.width;
  92         offy+=((x<dst.width/2)?x+dst.width/2:x-dst.width/2);
  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 }