fixed nlsop for packaging ..

[physik/nlsop.git] / dft.c

diff --git a/dft.c b/dft.c
index 385c22d..5c1e906 100644 (file)
--- a/dft.c
+++ b/dft.c
@@ -5,37 +5,147 @@
  *
  */
 
+#define USE_FFTW3
+
 #include "dft.h"
 
 int main(int argc,char **argv) {
 
+#ifdef USE_FFTW3
+  fftw_plan plan;
+  fftw_complex *in,*out;
+#else
   t_fourier fourier;
+#endif
   t_bmp src;
   t_bmp dst;
+  t_bmp cut;
+  int x,y;
+  int offy,offt;
+  double scale;
+  int size;
+  double max;
+  double *mag;
 
   bmp_init(&src,1);
   bmp_init(&dst,1);
+  bmp_init(&cut,1);
   src.mode=READ;
   dst.mode=WRITE;
   strcpy(src.file,argv[1]);
   strcpy(dst.file,argv[2]);
 
+#ifndef USE_FFTW3
+  fourier_init(&fourier,1);
+  fourier.type=DFT|FWD;
+  fourier.dim=2;
+#endif
+
   bmp_read_file(&src);
 
-  dst.width=src.info.width;
-  dst.height=src.info.height;
+  bmp_cut_grab_bottom(&cut,&src,src.info.width,GRAB);
+
+  dst.width=cut.width;
+  dst.height=cut.height;
   bmp_alloc_map(&dst);
 
-  /* simple copy for testing by now ! */
-  memcpy(dst.map,src.map,src.info.imagesize);
+  size=cut.width*cut.height;
+
+#ifndef USE_FFTW3
+  fourier.data_len[0]=cut.width;
+  fourier.data_len[1]=cut.height;
+  fourier_alloc_data(&fourier);
+#else
+  in=fftw_malloc(sizeof(fftw_complex)*size);
+  memset(in,0,size*sizeof(fftw_complex));
+  out=fftw_malloc(sizeof(fftw_complex)*size);
+  memset(out,0,size*sizeof(fftw_complex));
+  if(in==NULL||out==NULL) {
+    printf("fftw alloc error!\n");
+    return -23;
+  }
+  plan=fftw_plan_dft_2d(cut.width,cut.height,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
+#endif
+
+  mag=(double *)malloc(size*sizeof(double));
+  if(mag==NULL) {
+    printf("unable to alloc mag memory\n");
+    return -1;
+  }
+
+  // copy the data
+  for(x=0;x<size;x++) {
+#ifdef USE_FFTW3
+    in[x][0]=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
+    in[x][1]=0;
+#else
+    fourier.data[x].r=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
+#endif
+  }
+
+  // do the fourier trafo
+#ifdef USE_FFTW3
+  fftw_execute(plan);
+#else
+  fourier_dft_2d(&fourier);
+#endif
+
+
+  printf("fourier done!\n");
+
+  // copy back the data, intensity = sqrt(r^2+i^2)
+  max=0;
+  offt=0;
+  for(y=0;y<cut.height;y++) {
+    for(x=0;x<cut.width;x++) {
+#ifdef USE_FFTW3
+      mag[offt]=sqrt(out[offt][0]*out[offt][0]+out[offt][1]*out[offt][1]);
+#else
+      mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
+#endif
+      if(mag[offt]>max) max=mag[offt];
+      offt+=1;
+    }
+  }
+
+  printf("found max: %f\n",max);
+
+  if(max!=0) scale=255.0/max;
+  else scale=0;
+
+  if(scale!=0) {
+    printf("scaling image intensity: %f\n",scale);
+    offt=0;
+    for(y=0;y<dst.height;y++) {
+      for(x=0;x<dst.width;x++) {
+       // usual image processing order
+        offy=((y<dst.height/2)?y+dst.height/2:y-dst.height/2);
+        offy*=dst.width;
+        offy+=((x<dst.width/2)?x+dst.width/2:x-dst.width/2);
+       //offy=offt;
+        dst.map[offy].r=scale*mag[offt];
+        dst.map[offy].g=dst.map[offy].r;
+        dst.map[offy].b=dst.map[offy].r;
+        offt+=1;
+      }
+    }
+  }
+
+  free(mag);
 
   bmp_write_file(&dst);
 
+#ifdef USE_FFTW3
+  fftw_destroy_plan(plan);
+  fftw_free(in); fftw_free(out);
+#else
+  fourier_shutdown(&fourier);
+#endif
+
   bmp_shutdown(&src);
   bmp_shutdown(&dst);
 
   printf("done ...\n");
 
   return 1;
-
 }