+ 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);