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