ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/JSOC/proj/flatfield/pzt_flat_IDL/flatfield_iter_next.c
Revision: 1.1
Committed: Fri Feb 18 00:21:48 2011 UTC (12 years, 7 months ago) by richard
Content type: text/plain
Branch: MAIN
CVS Tags: Ver_6-0, Ver_6-1, Ver_6-2, Ver_6-3, Ver_6-4, Ver_9-1, Ver_5-14, Ver_5-13, Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-2, Ver_8-8, Ver_8-2, Ver_8-3, Ver_8-0, Ver_8-1, Ver_8-6, Ver_8-7, Ver_8-4, Ver_8-5, Ver_7-1, Ver_7-0, Ver_9-5, Ver_9-4, Ver_8-10, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Log Message:
C fuction called from IDL
2011.02.17

File Contents

# Content
1
2
3 #include <stdio.h>
4 #include <math.h>
5 #include <stdlib.h>
6 #include <string.h>
7 #include <omp.h>
8
9
10
11
12
13
14
15 const int niter=10;
16 const float omega=1.2;
17 const int residuum=0;
18 const int n_foc=1;
19
20
21
22 void apod_circ(int, float, float, int, float, float, float *, int , int);
23 void printtime();
24
25 int main(int argc, char *argv[]){
26
27 int camera;
28 int nx, ny;
29 int nl;
30 float fac,facr;
31
32 if (argc != 6){fputs("Wrong number of arguments", stderr); exit(1);}
33
34 nx=atoi(argv[1]);
35 ny=atoi(argv[1]);
36 nl=atoi(argv[2]);
37 camera=atoi(argv[3]);
38 fac=(float)atoi(argv[4])/100.0;
39 facr=(float)atoi(argv[5])/100.0;
40
41 printf("nx ny nl cam %d %d %d %d\n", nx, ny, nl, camera);
42
43 float samp=4096.0/(float)nx;
44
45
46 float **imr; // image variable
47 float *imp;
48 float *imm;
49
50 float *imr_im, *imr_jm;
51
52 //float g[nx][ny]; // flatfield iteration
53 //float kap[nx][ny]; // initial value
54
55 float *vdo, *vdl; //
56 float *vdp;
57
58 vdo=(float *)(malloc(nx*ny*sizeof(float)));
59 vdl=(float *)(malloc(nx*ny*sizeof(float)));
60 vdp=(float *)(malloc(nx*ny*nl*sizeof(float)));
61
62 float xx[nl], yy[nl];
63
64 int nbad;
65 //side
66 // int nbad=32;
67 int ix_bad[32];
68
69 int ix_bad_front[32]={738657 , 1100979 , 1100980 , 1420552 , 1424648 , 1428744 , 1502539 , 1502540 , 1506635 , 1506636 , 5961670 , 6051854, 6330150 , 10171860, 10175955 , 10175956};
70
71 int ix_bad_side[32]={1800075 ,1800076 ,1804173 ,9958325 ,9958326 ,9958327 ,9962423 ,9962424 ,9962425 ,9966520 ,12176519 ,12176520 ,12176521 ,12180613 ,12180614 ,12180615 ,12180616 ,12180617 ,12184709 , 12184710 ,12184711 ,12184712 ,12184713 ,12188806 ,12188807 , 12188808 , 12188809 ,12426166 ,12426167 ,12430263 , 12430264 ,12434360};
72
73
74
75
76
77
78
79 //float n[nx][ny];
80 //float gn[nx][ny];
81 //float nn[nx][ny];
82
83
84 float *g, *kap, *n, *gn;
85 g=(float *)(malloc(nx*ny*sizeof(float)));
86 kap=(float *)(malloc(nx*ny*sizeof(float)));
87
88 n=(float *)(malloc(nx*ny*sizeof(float)));
89 gn=(float *)(malloc(nx*ny*sizeof(float)));
90
91
92 float res[niter];
93
94 // reading images
95 float datum;
96 int im_number;
97
98 float delta;
99
100 int i, j, k, im, jm, iter; // loop variables
101
102 int valfoc[7]={0,2,3,4,5,9,13}; // 0 replacing 1
103
104
105
106 //front
107 if (camera == 2){nbad=16;
108 for (k=0; k<nbad; ++k) ix_bad[k]=ix_bad_front[k];
109 }
110
111 //side camera
112 if (camera == 1){nbad=32;
113 for (k=0; k<nbad; ++k) ix_bad[k]=ix_bad_side[k];
114 }
115
116
117 printf("reading data\n");
118 int fi;
119 for (fi=0; fi < n_foc; ++fi){ //0 replacing 6
120
121 //reading data
122
123 int foc=valfoc[fi];
124 char fnumb[2]={""};
125 sprintf(fnumb, "%i", foc);
126
127 ///get input name
128 char fname[30]={""};
129
130
131 strcat(fname, "imr_hmi");
132 strcat(fname, fnumb);
133 strcat(fname, ".bin");
134
135 printf("input name %s \n", fname);
136 //////////////////////////
137 // get output flatfield name
138
139
140
141
142 char ffname[40]={""};
143
144 strcat(ffname, "flatfield_out");
145 strcat(ffname, fnumb);
146 strcat(ffname, ".bin");
147
148 printf("output name %s \n", ffname);
149 //////////////////////////////
150
151 imm=(float *)(malloc(nx*sizeof(float)));
152 imr=(float **)(malloc(nl*sizeof(float*)));
153
154 FILE *imfile;
155 size_t bytes_read;
156 imfile=fopen(fname, "rb"); // read binary data
157
158 if (imfile==NULL){fputs("File error", stderr); exit(1);};
159 {
160 for (k=0; k<nl; ++k){
161 printf("%u \n", k);
162 imp=(float *)(malloc(nx*ny*sizeof(float)));
163 *(imr+k)=imp;
164 bytes_read=fread(imp,sizeof(float), nx*ny,imfile);
165 }
166 }
167 fclose(imfile);
168
169
170 //get log of input data
171
172
173 // reading leg positions
174 FILE *legpos;
175 legpos=fopen("legpos2", "r");
176
177 if (legpos == NULL){fputs("File error", stderr); exit(1);};
178
179 {
180 for (i=0; i<nl; i++){
181 fscanf(legpos, "%f", &datum);
182 xx[i]=datum/samp;
183 printf("%f\n", xx[i]);
184 }
185 printf("\n");
186 for (i=0; i<nl; i++){
187 fscanf(legpos, "%f", &datum);
188 yy[i]=datum/samp;
189 printf("%f\n", yy[i]);
190 }
191 }
192 fclose(legpos);
193
194
195
196 // constructing apodization
197
198 int low9=0; // int(192.0/samp);
199 float nr=nx/2*fac;
200 float nb=nx/2*fac*0.0; // nb=0 !!
201
202 float nrs=nx/2*facr;
203 float nbs=nx/2*facr*0.0;
204
205 apod_circ(nx, nr, nb, 0, 0.0, 0.0, vdl,nx,ny);
206
207
208 for (k=0; k<nbad; ++k) vdl[ix_bad[k]]=0.0;
209
210 for (k=0; k<nl; ++k){
211 printf("%d\n", k);
212 apod_circ(nx, nrs, nbs, 0, xx[k], yy[k], vdo,nx,ny);
213
214 for (j=0; j<nx; ++j) for (i=0; i<ny; ++i) vdp[k*nx*ny+j*nx+i]=vdl[j*nx+i]*vdo[j*nx+i];
215 }
216
217
218
219 // zero iteration
220
221 int xf, yf, xg, yg;
222 float rx, ry;
223 float xm, ym;
224 float xn, yn;
225
226
227
228
229 int nthreads;
230 nthreads=omp_get_num_procs(); //number of threads supported by the machine where the code is running
231 omp_set_num_threads(nthreads); //set the number of threads to the maximum value
232 printf("Number of threads run in parallel by the subroutine= %d \n",nthreads);
233
234
235 printtime();
236
237 printf("zero iteration\n");
238
239 for (j=0; j<ny; ++j) for (i=0; i<nx; ++i) kap[j*nx+i]=0.0;
240
241 #pragma omp parallel for private(i,j,im,jm,xm,xf,ym,yf,yn,yg,xn,xg)
242 for (j=0; j<ny; ++j){
243 for (i=0; i<nx; ++i){
244 n[j*nx+i]=0.0;
245 for (jm=0; jm<(nl-1); ++jm){ //exchanged loops
246 im=jm+1;
247
248 imr_im=*(imr+im);
249 imr_jm=*(imr+jm);
250
251 xm=i-xx[im]+xx[jm];
252 xf=(int)(xm+0.5);
253 ym=j-yy[im]+yy[jm];
254 yf=(int)(ym+0.5);
255
256 if (xf >= 0 && xf <= nx-2 && yf >= 0 && yf <= ny-2){
257 if (vdp[im*nx*ny+j*nx+i] != 0.0 && vdp[jm*nx*ny+yf*nx+xf] != 0.0 && imr_im[j*nx+i] > 0.0 && imr_jm[yf*nx+xf] > 0.0){
258
259 kap[j*nx+i]=kap[j*nx+i]+(vdp[im*nx*ny+j*nx+i])*(vdp[jm*nx*ny+yf*nx+xf])*(log(imr_im[j*nx+i])-log(imr_jm[yf*nx+xf]));
260 n[j*nx+i]=n[j*nx+i]+(vdp[im*nx*ny+j*nx+i])*(vdp[jm*nx*ny+yf*nx+xf]);
261
262 }
263 }
264
265 yn=j+yy[im]-yy[jm];
266 yg=(int)(yn+0.5);
267 xn=i+xx[im]-xx[jm];
268 xg=(int)(xn+0.5);
269
270 if (xg >= 0 && xg <= nx-2 && yg >=0 && yg <= ny-2){
271
272 if (vdp[jm*nx*ny+j*nx+i] != 0.0 && vdp[im*nx*ny+yg*nx+xg] !=0.0 && imr_jm[j*nx+i] > 0.0 && imr_im[yg*nx+xg] > 0.0){
273
274
275 kap[j*nx+i]=kap[j*nx+i]+(vdp[jm*nx*ny+j*nx+i])*(vdp[im*nx*ny+yg*nx+xg])*(log(imr_jm[j*nx+i])-log(imr_im[yg*nx+xg]));
276
277 n[j*nx+i]=n[j*nx+i]+(vdp[jm*nx*ny+j*nx+i])*(vdp[im*nx*ny+yg*nx+xg]);
278 }
279 }
280
281 }
282
283 if (n[j*nx+i] > 0.0) kap[j*nx+i]=kap[j*nx+i]/n[j*nx+i]; else kap[j*nx+i]=0.0;
284
285 }
286 }
287
288
289
290
291
292 for (j=0; j<ny; ++j) for (i=0; i<nx; ++i) g[j*nx+i]=kap[j*nx+i];
293 for (i=0; i<niter; ++i) res[i]=0.0;
294
295
296 printtime();
297
298 // iterations
299
300 for (iter=0; iter<niter; ++iter){
301
302 printf("iteration %u \n", iter+1);
303
304 #pragma omp parallel for private(i,j,im,jm,xm,xf,ym,yf,yn,yg,xn,xg)
305 for (j=0; j<ny; ++j){
306 for (i=0; i<nx; ++i){
307 n[j*nx+i]=0.0;
308 gn[j*nx+i]=0.0;
309
310
311 for (jm=0; jm<nl; ++jm){
312 im=jm+1;
313
314 xm=i-xx[im]+xx[jm];
315 xf=(int)(xm+0.5);
316 ym=j-yy[im]+yy[jm];
317 yf=(int)(ym+0.5);
318
319 if (yf >= 0 && yf <= ny-2 && xf >= 0 && xf <= nx-2){
320 if (vdp[im*nx*ny+j*nx+i] != 0.0 && vdp[jm*nx*ny+yf*nx+xf] !=0.0){
321
322 n[j*nx+i]=n[j*nx+i]+vdp[im*nx*ny+j*nx+i]*vdp[jm*nx*ny+yf*nx+xf];
323 gn[j*nx+i]=gn[j*nx+i]+vdp[im*nx*ny+j*nx+i]*vdp[jm*nx*ny+yf*nx+xf]*g[yf*nx+xf];
324
325
326 }
327 }
328
329
330
331 xn=i+xx[im]-xx[jm];
332 xg=(int)(xn+0.5);
333 yn=j+yy[im]-yy[jm];
334 yg=(int)(yn+0.5);
335
336 if (yg >= 0 && yg <= ny-2 && xg >= 0 && xg <= nx-2){
337 if (vdp[jm*nx*ny+j*nx+i] !=0.0 && vdp[im*nx*ny+yg*nx+xg] !=0.0){
338
339 n[j*nx+i]=n[j*nx+i]+vdp[jm*nx*ny+j*nx+i]*vdp[im*nx*ny+yg*nx+xg];
340 gn[j*nx+i]=gn[j*nx+i]+vdp[jm*nx*ny+j*nx+i]*vdp[im*nx*ny+yg*nx+xg]*g[yg*nx+xg];
341
342
343 }
344 }
345
346
347 }
348
349 if (n[j*nx+i] > 0.0) gn[j*nx+i]=gn[j*nx+i]/n[j*nx+i]; else gn[j*nx+i]=0.0;
350 delta=kap[j*nx+i]-g[j*nx+i]+gn[j*nx+i];
351 g[j*nx+i]=g[j*nx+i]+omega*delta;
352
353
354 }
355 }
356
357
358
359 if (residuum){
360 for (j=0; j<ny; ++j){
361 for (i=0; i<nx; ++i){
362
363 for (jm=0; jm<(nl-1); ++jm){
364 im=jm+1;
365
366 imr_im=*(imr+im);
367 imr_jm=*(imr+jm);
368
369 xf=(int)(i+xx[im]);
370 yf=(int)(j+yy[im]);
371 xg=(int)(i+xx[jm]);
372 yg=(int)(j+yy[jm]);
373
374 if (xf >=0 && xf <nx && xg >=0 && xg <nx && yf>=0 && yf<ny && yg >=0 && yg < ny){
375 if (vdp[im*nx*ny+yf*nx+xf]*vdp[jm*nx*ny+yg*nx+xg] !=0.0 && imr_im[yf*nx+xf] > 0.0 && imr_jm[yg*nx+xg] > 0.0){
376 res[iter]=res[iter]+vdp[im*nx*ny+yf*nx+xf]*vdp[jm*nx*ny+yg*nx+xg]*pow(-g[yf*nx+xf]+g[yg*nx+xg]+log(imr_im[yf*nx+xf])-log(imr_jm[yg*nx+xg]),2);
377 }
378 }
379
380
381 }
382
383 }
384 }
385
386 printf("residuum: %f\n", res[iter]);
387 }
388 }
389
390
391 printtime();
392
393
394 FILE *outname;
395 outname = fopen (ffname, "w");
396 if (outname == NULL){fputs("File error", stderr); exit(1);};
397
398 {
399 for (j=0;j<ny;j++){
400 for (i=0; i<nx; ++i) imm[i]=exp(g[j*nx+i]);
401 fwrite ((char*)(imm),sizeof(float),nx,outname);
402 }
403 fclose(outname);
404
405
406
407 for (k=0; k<nl; ++k) free(*(imr+k));
408 free(imr);
409 }
410
411 free(imm);
412
413 return 0;
414 }
415
416
417 void apod_circ(int nn, float rad, float nb, int low9, float offx, float offy, float *vd, int nx, int ny)
418 {
419 float *rarr;
420 rarr=(float *)(malloc(nx*ny*sizeof(float)));
421 int i, j;
422
423 for (i=0; i<nn; ++i) for (j=0; j<nn; ++j) rarr[j*nx+i]=sqrt(((float)i-((float)nn/2+offx))*((float)i-((float)nn/2+offx))+((float)j-((float)nn/2+offy))*((float)j-((float)nn/2+offy)));
424
425 for (i=0; i<nn; ++i){
426 for (j=0; j<nn; ++j){
427
428 if (rarr[j*nx+i] < rad)
429 vd[j*nx+i]=1.0;
430
431 if (rarr[j*nx+i] >= rad && rarr[j*nx+i] < (rad+nb))
432 vd[j*nx+i]=0.5*cos(M_PI/nb*(rarr[j*nx+i]-rad))+0.5;
433
434 if (rarr[j*nx+i] >= (rad+nb))
435 vd[j*nx+i]=0.0;
436
437 if (low9 != 0){
438
439 if (j < low9) vd[j*nx+i]=0.0;
440 if (j >= low9 && j <= (low9+(int)nb)) vd[j*nx+i]=(-0.5*cos(M_PI/nb*(j-low9))+0.5)*vd[j*nx+i];
441 }
442
443
444
445
446 }
447 }
448
449 }
450
451
452
453
454 void printtime() // print time
455 {
456 time_t timer, timerd;
457 char *timestring;
458 int i;
459
460 timerd=time(&timer);
461 timestring=ctime(&timer);
462 for (i=0; i<24; ++i) printf("%c", *(timestring+i));
463 printf("\n");
464 }
465
466