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