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   offy=0;
  66   for(y=0;y<fourier.data_len[1];y++) {
  67     for(x=0;x<fourier.data_len[0];x++) {
  68       offt=offy+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     }
  72     offy+=fourier.data_len[0];
  73   }
  74
  75   if(max!=0) scale=(int)255/max;
  76   else scale=0;
  77
  78   if(scale!=0) {
  79     printf("scaling image intensity: %d\n",scale);
  80     for(x=0;x<dst.width*dst.height;x++) {
  81       dst.map[x].r=scale*mag[x];
  82       dst.map[x].g=dst.map[x].r;
  83       dst.map[x].b=dst.map[x].r;
  84     }
  85   }
  86
  87   bmp_write_file(&dst);
  88
  89   fourier_shutdown(&fourier);
  90
  91   bmp_shutdown(&src);
  92   bmp_shutdown(&dst);
  93
  94   printf("done ...\n");
  95
  96   return 1;
  97 }