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