1 arta 1.1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <time.h>
4 #include <string.h>
5 #include <math.h>
6 #include <mkl.h>
7 #include <omp.h>
|
8 arta 1.10 #include <sys/param.h>
|
9 arta 1.1 #include "cholesky_down.h"
10 #include "gapfill.h"
|
11 arta 1.10
|
12 arta 1.1
13 #define minval(x,y) (((x) < (y)) ? (x) : (y))
14 #define maxval(x,y) (((x) < (y)) ? (y) : (x))
15
|
16 arta 1.10 #define kInterpDataDir "proj/libs/interpolate/data"
17
|
18 arta 1.1 struct work_struct {
19 double *a,*rh,*a1b,*a1r,*a1x;
20 int *wgood,*wbad;
21 };
22
23 void malloc_work(
24 // Malloc work space for each thread
25 int nsample,
26 struct work_struct *work
27 )
28 {
29 int malign=32;
30
31 work->wgood=(int *)(MKL_malloc(nsample*sizeof(int),malign));
32 work->wbad=(int *)(MKL_malloc(nsample*sizeof(int),malign));
33 work->a=(double *)(MKL_malloc(nsample*nsample*sizeof(double),malign));
34 work->rh=(double *)(MKL_malloc(nsample*sizeof(double),malign));
35 work->a1b=(double *)(MKL_malloc(nsample*sizeof(double),malign));
36 work->a1r=(double *)(MKL_malloc(nsample*sizeof(double),malign));
37 work->a1x=(double *)(MKL_malloc(2*nsample*sizeof(double),malign));
38 }
39 arta 1.1
40 void free_work(
41 // Free the workspace
42 int nsample,
43 struct work_struct *work
44 )
45 {
46 MKL_free(work->wgood);
47 MKL_free(work->wbad);
48 MKL_free(work->a);
49 MKL_free(work->rh);
50 MKL_free(work->a1b);
51 MKL_free(work->a1r);
52 MKL_free(work->a1x);
53 }
54
55 int init_fill(
56 int method, // Interpolation method
57 double pnoise, // Level of photon noise for trade-off
58 int order, // Interpolation order. Generally odd.
59 int targetx, // Target point in x (normally (order-1)/2)
60 arta 1.1 int targety, // Target point in y (normally (order-1)/2)
61 struct fill_struct *pars, // Structure to save setup information etc.
|
62 arta 1.10 char **filenamep, // Pointer to name of file to read covariance from.
63 // Set to actual file used if method > 0.
64 const char *path // to data files read by this function.
|
65 arta 1.1 )
66 {
67 double *acor;
68 int nx=80;
69 int malign=32;
70 double t0,t1;
71 FILE *fileptr;
72 int i,j;
73 int *offx,*offy,ox,oy;
74 int nsample,inoise;
75 double regul,pmin;
76 double *a0,*rh0,*a0t,*rh0t;
77 int info;
78 char uplo[] = "U";
|
79 arta 1.10 char buf[PATH_MAX];
|
80 arta 1.1
81 pars->order=order;
82
83 t0=dsecnd();
84 switch (method) {
85 case 0:
86 pars->filename=strdup(*filenamep);
87 break;
88 case 8:
|
89 arta 1.10 snprintf(buf, sizeof(buf), "%s/%s/acor1_08_80.txt", path, kInterpDataDir);
90 pars->filename=strdup(buf);
|
91 arta 1.1 break;
92 case 9:
|
93 arta 1.10 snprintf(buf, sizeof(buf), "%s/%s/acor1_09_80.txt", path, kInterpDataDir);
94 pars->filename=strdup(buf);
|
95 arta 1.1 break;
96 case 10:
|
97 arta 1.10 snprintf(buf, sizeof(buf), "%s/%s/acor1_10_80.txt", path, kInterpDataDir);
98 pars->filename=strdup(buf);
|
99 arta 1.1 break;
100 case 11:
|
101 arta 1.10 snprintf(buf, sizeof(buf), "%s/%s/acor1_11_80.txt", path, kInterpDataDir);
102 pars->filename=strdup(buf);
|
103 arta 1.1 break;
104 default:
105 printf("Unknown method in init_fill.\n");
106 return 1;
107 break;
108 }
109 fileptr = fopen (pars->filename,"r");
110 if (fileptr==NULL) {
111 printf("File not found in init_fill.\n");
112 return 1;
113 }
114 if (method != 0) {
115 if (filenamep != NULL) {
116 *filenamep=strdup(pars->filename);
117 }
118 }
119 acor=(double *)(MKL_malloc(nx*nx*sizeof(double),malign));
120 //fread ((char*)acor,sizeof(double),nx*nx,fileptr);
121 for (i=0;i<nx*nx;i++) fscanf(fileptr,"%lf",acor+i);
122 fclose(fileptr);
123 t0=dsecnd()-t0;
|
124 schou 1.8 //printf("x0 %f\n",t0);
|
125 arta 1.1
126 regul=pnoise*pnoise;
127 inoise=1;
128 pmin=regul*inoise;
129
130 nsample=order*order;
131 t1=dsecnd();
132 offx=(int *)(MKL_malloc(nsample*sizeof(int),malign));
133 offy=(int *)(MKL_malloc(nsample*sizeof(int),malign));
134 for (i=0;i<nsample;i++) {
135 offx[i]=(i%order)-targetx;
136 offy[i]=(i/order)-targety;
137 }
138 // Set up complete matrix
139 pars->a0=(double *)(MKL_malloc(nsample*nsample*sizeof(double),malign));
140 pars->a00=(double *)(MKL_malloc(nsample*nsample*sizeof(double),malign));
141 pars->a0t=(double *)(MKL_malloc(nsample*nsample*sizeof(double),malign));
142 a0=pars->a0;
143 a0t=pars->a0t;
144 for (i=0;i<nsample;i++) {
145 for (j=0;j<nsample;j++) {
146 arta 1.1 ox=((offx[j]-offx[i])+nx) % nx;
147 oy=((offy[j]-offy[i])+nx) % nx;
148 a0[i+nsample*j]=acor[ox+nx*oy];
149 a0t[i+nsample*j]=acor[ox+nx*oy];
150 }
151 a0[i+nsample*i]=a0[i+nsample*i]+regul;
152 }
153 // save copy of matrix
154 memcpy(pars->a00,a0,nsample*nsample*sizeof(double));
155 pars->acort00=acor[0];
156 pars->rh0=(double *)(MKL_malloc(nsample*sizeof(double),malign));
157 pars->rh0t=(double *)(MKL_malloc(nsample*sizeof(double),malign));
158 rh0=pars->rh0;
159 rh0t=pars->rh0t;
160 for (i=0;i<nsample;i++) {
161 rh0[i]=acor[(-offx[i]+nx) % nx+nx*((-offy[i]+nx) % nx)];
162 rh0t[i]=acor[(-offx[i]+nx) % nx+nx*((-offy[i]+nx) % nx)];
163 if ((offx[i]==0) && (offy[i]==0)) rh0[i]=rh0[i]+pmin;
164 }
165 // Decompose complete matrix
166 dpotrf(uplo,&nsample,pars->a0,&nsample,&info); // Cholesky decomposition
167 arta 1.1
168 MKL_free(acor);
169 MKL_free(offx);
170 MKL_free(offy);
171
172 // Variables used in further calculations. Avoids multiple mallocs.
173 //pars->wgood=(int *)(MKL_malloc(nsample*sizeof(int),malign));
174 //pars->wbad=(int *)(MKL_malloc(nsample*sizeof(int),malign));
175
176 pars->hashmod=8191; // Hash table size
177 pars->hashcount=(int *)(MKL_malloc(pars->hashmod*sizeof(int),malign));
178 pars->hashtable=(struct fill_hash_struct **)(MKL_malloc(pars->hashmod*sizeof(struct fill_hash_struct *),malign));
179 pars->locks=(omp_lock_t *)(MKL_malloc(pars->hashmod*sizeof(omp_lock_t),malign));
180 pars->complocks=(omp_lock_t *)(MKL_malloc(pars->hashmod*sizeof(omp_lock_t),malign));
181
182 for (i=0;i<pars->hashmod;i++) {
183 pars->hashcount[i]=0;
184 // Initialize to null pointers
185 pars->hashtable[i]=(struct fill_hash_struct *)NULL;
186 omp_init_lock(&(pars->locks[i]));
187 omp_init_lock(&(pars->complocks[i]));
188 arta 1.1 }
189 pars->ndiff=0;
190
191 t1=dsecnd()-t1;
|
192 schou 1.8 //printf("x1 %f\n",t1);
|
193 arta 1.1
194 return 0;
195 }
196
197 int free_fill(
198 struct fill_struct *pars
199 )
200 {
201 int i;
202 struct fill_hash_struct *ptr,*oldptr;
203
|
204 schou 1.8 /*
|
205 arta 1.1 printf("Hashmod: %d\n",pars->hashmod);
206 printf("Ndiff: %d\n",pars->ndiff);
207 for (i=0;i<pars->hashmod;i++) {
208 printf("%d %d\n",i,pars->hashcount[i]);
209 }
|
210 schou 1.6 */
|
211 schou 1.8 /*
|
212 arta 1.1 for (i=0;i<pars->hashmod;i++) {
213 ptr=(pars->hashtable[i]);
214 while (ptr != NULL) {
|
215 schou 1.7 printf("# %d %d %d\n",i,ptr->nbad,ptr->nhit);
|
216 arta 1.1 ptr=ptr->next;
217 }
218 }
|
219 schou 1.8 */
|
220 schou 1.3
221 for (i=0;i<pars->hashmod;i++) {
222 omp_destroy_lock(&(pars->locks[i]));
223 omp_destroy_lock(&(pars->complocks[i]));
224 }
225
|
226 arta 1.1 MKL_free(pars->a0);
227 MKL_free(pars->a00);
228 MKL_free(pars->a0t);
229 MKL_free(pars->rh0);
230 MKL_free(pars->rh0t);
231
232 MKL_free(pars->locks);
233 MKL_free(pars->complocks);
234
235 // Now free hash table and contents
236
237 MKL_free(pars->hashcount);
238 for (i=0;i<pars->hashmod;i++) {
239 ptr=(pars->hashtable[i]);
240 while (ptr != NULL) {
241 oldptr=ptr;
242 MKL_free(ptr->wbad);
243 MKL_free(ptr->coeff);
244 ptr=ptr->next;
245 MKL_free(oldptr);
246 }
247 arta 1.1 }
248 MKL_free(pars->hashtable);
249 free(pars->filename);
250
251 return 0; // Currently no error conditions
252
253 }
254
255 int find_entry(
256 struct fill_struct *pars,
257 int hash0,
258 int hash,
259 int nbad,
260 int *wbad,
261 struct fill_hash_struct **ptr_out
262 )
263 {
264 int found;
265 int ident,i;
266 struct fill_hash_struct *ptr;
267
268 arta 1.1 found=0;
269 *ptr_out=NULL;
270
271 // Prevent update while starting traverse
272 omp_set_lock(&(pars->locks[hash]));
273 ptr=pars->hashtable[hash];
274 omp_unset_lock(&(pars->locks[hash]));
275
276 while (ptr!=NULL) { // Not at end of table
277 if ((hash0==ptr->hash0) && (nbad==ptr->nbad)) {
278 // Original hash and number of elements is correct, so there is hope
279 ident=1;
280 for (i=0;i<nbad;i++) {
281 if (wbad[i] != ptr->wbad[i]) {
282 ident=0;
283 break; // Out of for loop. Known not to match.
284 }
285 }
286 if (ident != 0) { // Success
|
287 schou 1.7 #pragma omp atomic
|
288 arta 1.1 ptr->nhit++;
289 found=1;
290 *ptr_out=ptr;
291 break; // Out of while loop. No need to test the rest.
292 } // if
293 } // if
294 ptr=ptr->next;
295 } // while
296
297 return found;
298 } // find_entry
299
300 double fill_point(
301 struct fill_struct *pars,
302 int *mask,
303 double *vals,
304 float *cnorm,
305 float *ierror,
306 struct work_struct *work
307 )
308 {
309 arta 1.1 double t1,t2,t3,t4,t5,t6,t7,t8/*,t9*/;
310 int nsample;
311 int ngood,nbad;
312 int *wgood,*wbad;
313 int i,info,j;
314 char uplo[] = "U";
315 int ione = 1;
316 int itwo = 2;
317 int malign=32;
318 double *a,*a0,*rh0,/* *rh,*a1b,*/bta1b,/**a1r,*/bta1r,c,*coeff,*a0t,*a00,*a1x;
319 float fillval;
320 int hash0,hash,found;
321 struct fill_hash_struct *ptr;
322 double ierror2; // interpolation error squared
323 int unlocked;
324 int nleft;
325 double tdown,tbrute;
326 double sum;
327
328 t1=0.0;t2=0.0;t3=0.0;t4=0.0;t5=0.0;t6=0.0;t7=0.0;t8=0.0;//t9=0.0;
329
330 arta 1.1 t1=dsecnd();
331
332 nsample=(pars->order)*(pars->order);
333
334 // Make work copies
335 a0=pars->a0;
336 a00=pars->a00;
337 a0t=pars->a0t;
338 rh0=pars->rh0;
339
340 wgood=work->wgood;
341 wbad=work->wbad;
342
343 ngood=0;
344 nbad=0;
345 hash0=0;
346 for (i=0;i<nsample;i++) {
347 if (mask[i]==0) {
348 wgood[ngood]=i;
349 ngood=ngood+1;
350 }
351 arta 1.1 if (mask[i]!=0) {
352 hash0=hash0+i;
353 wbad[nbad]=i;
354 nbad=nbad+1;
355 }
356 }
357 hash=hash0%pars->hashmod;
358
359 unlocked=0;
360
361 found=find_entry(pars,hash0,hash,nbad,wbad,&ptr);
362
363 if (found == 0) {
364 unlocked=omp_test_lock(&(pars->complocks[hash])); // Capture lock status
365 if (unlocked==0) omp_set_lock(&(pars->complocks[hash]));
366
367 // Rerun search. Someone may have updated while waiting for lock
368 found=find_entry(pars,hash0,hash,nbad,wbad,&ptr);
369 if (found!=0) {
370 // Clear lock
371 omp_unset_lock(&(pars->complocks[hash]));
372 arta 1.1 found=2; // Indicate collision
373 }
374 }
375
376 t1=dsecnd()-t1;
377
378 if (found == 0) { // Got to calculate new coefficients
379
380 //printf("New %d %d\n",omp_get_thread_num(),hash);
381 // malloc new entry
382 //for (i=0;i<nsample;i++) printf(" %d",mask[i]);
383 //printf("\n");
384 ptr=(struct fill_hash_struct *)(MKL_malloc(sizeof(struct fill_hash_struct),malign));
385 ptr->next=pars->hashtable[hash]; // Point to next element
386 // Set other variables
387 ptr->nbad=nbad;
388 ptr->hash0=hash0;
389 ptr->nhit=1;
390 ptr->wbad=(int *)(MKL_malloc(nbad*sizeof(int),malign));
391 for (i=0;i<nbad;i++) ptr->wbad[i]=wbad[i];
392 coeff=(double *)(MKL_malloc(ngood*sizeof(double),malign));
393 arta 1.1 pars->hashcount[hash]++; // Increment count for this hash
394 #pragma omp atomic
395 pars->ndiff++; // Increment count for total
396 ptr->coeff=coeff; // Point to new coefficients
397
398 // Point to work arrays
399 a=work->a;
400 //rh=work->rh;
401 //a1b=work->a1b;
402 //a1r=work->a1r;
403 a1x=work->a1x;
404
405 // Calculate new coefficients
406 // First estimate time for each option
407 tdown=0.0;
408 nleft=nsample;
409 for (i=nsample-1;i>=0;i--) {
410 if (mask[i]!=0) {
411 tdown=tdown+(nleft-i-1)*(nleft-i-1);
412 nleft=nleft-1;
413 }
414 arta 1.1 }
415 tdown=3.9e-9*tdown;
416 tbrute=4e-11*(ngood+0.0)*(ngood+0.0)*(ngood+0.0);
417
418 // Constants in if statement based on timings
419 // if (nbad < (1.50*sqrt(nsample)-6)) {
420 if (tdown<tbrute) {
421 t2=dsecnd();
422 // Get new decomposition by downsampling old one
423 // Copy decomposed matrix.
424 // This dominates dcholesky_down when only removing a few points.
425 // Only need to copy relevant part of matrix
426 // for (i=0;i<nsample;i++) for (j=0;j<=i;j++) a[i*nsample+j]=a0[i*nsample+j];
427 dlacpy(uplo,&nsample,&nsample,a0,&nsample,a,&nsample);
428 t2=dsecnd()-t2;
429 t3=dsecnd();
430
431 // Get decomposition of reduced matrix
432 dcholesky_down_mask(uplo, &nsample, a, &nsample, mask, &info);
433 t3=dsecnd()-t3;
434 } // if downsampling
435 arta 1.1 else {
436 t4=dsecnd();
437 // Calculate brute force from original matrix
438 // Pick good points
439 for (i=0;i<ngood;i++) for (j=0;j<=i;j++) a[i*nsample+j]=a00[wgood[i]*nsample+wgood[j]];
440 t4=dsecnd()-t4;
441 t5=dsecnd();
442
443 // Get decomposition of reduced matrix
444 dpotrf(uplo,&ngood,a,&nsample,&info); // Cholesky decomposition
445 t5=dsecnd()-t5;
446 } // Brute force
447
448 t6=dsecnd();
449
450 /*
451 for (i=0;i<ngood;i++) a1b[i]=1.;
452 for (i=0;i<ngood;i++) a1r[i]=rh0[wgood[i]];
453 dpotrs(uplo,&ngood,&ione,a,&nsample,a1b,&nsample,&info);
454 dpotrs(uplo,&ngood,&ione,a,&nsample,a1r,&nsample,&info);
455 bta1b=0.;
456 arta 1.1 for (i=0;i<ngood;i++) bta1b=bta1b+a1b[i];
457 bta1r=0.0;
458 for (i=0;i<ngood;i++) bta1r=bta1r+a1r[i];
459 c=(bta1r-1.)/bta1b;
460 for (i=0;i<ngood;i++) coeff[i]=a1r[i]-c*a1b[i];
461 ptr->cnorm=dnrm2(&ngood,coeff,&ione);
462 */
463
464 // Combine a1b and a1r and make single call to dpotrs
465 for (i=0;i<ngood;i++) a1x[i]=1.;
466 for (i=0;i<ngood;i++) a1x[i+nsample]=rh0[wgood[i]];
467 dpotrs(uplo,&ngood,&itwo,a,&nsample,a1x,&nsample,&info);
468 bta1b=0.;
469 for (i=0;i<ngood;i++) bta1b=bta1b+a1x[i];
470 bta1r=0.0;
471 for (i=0;i<ngood;i++) bta1r=bta1r+a1x[i+nsample];
472 c=(bta1r-1.)/bta1b;
473 for (i=0;i<ngood;i++) coeff[i]=a1x[i+nsample]-c*a1x[i];
474 ptr->cnorm=(float)dnrm2(&ngood,coeff,&ione);
475
476 t6=dsecnd()-t6;
477 arta 1.1 t7=dsecnd();
478
479 // Now get interpolation error
480 ierror2=0.0;
481 for (i=0;i<ngood;i++) ierror2=ierror2+coeff[i]*rh0[wgood[i]];
482 ierror2=pars->acort00-2*ierror2;
483 // Now calculate Sum_wood[i,j] coeff_i A0t_ij coeff_j
484 for (j=0;j<ngood;j++) {
485 /*
486 sum=0.0;
487 for (i=0;i<ngood;i++) {
488 sum=sum+a0t[nsample*wgood[j]+wgood[i]]*coeff[i];
489 }
490 ierror2=ierror2+sum*coeff[j];
491 */
492 // a0t is symmetric, so this is more efficient
493 sum=0.0;
494 for (i=0;i<j;i++) {
495 sum=sum+a0t[nsample*wgood[j]+wgood[i]]*coeff[i];
496 }
497 ierror2=ierror2+coeff[j]*(2*sum+coeff[j]*a0t[nsample*wgood[j]+wgood[j]]);
498 arta 1.1 }
499 ptr->ierror=(float)sqrt(ierror2/pars->acort00);
500
501 // Stick new entry at beginning of linked list
502 // Prevent others from starting to traverse
503 omp_set_lock(&(pars->locks[hash]));
504 pars->hashtable[hash]=ptr;
505 #pragma omp flush
506 omp_unset_lock(&(pars->locks[hash]));
507
508 // Done with computation
509 omp_unset_lock(&(pars->complocks[hash]));
510
511 t7=dsecnd()-t7;
512
513 // Could store coefficients in the other 7 rotations of the mask.
514 // Beware that the time for cholesky_down depends strongly on which
515 // one is calculated. See timep.pro for calculation (t2s is a good predictor).
516 } // End of calculating new coefficients
517
518 t8=dsecnd();
519 arta 1.1 fillval=(float)ddot(&ngood,ptr->coeff,&ione,vals,&ione);
520 *cnorm=ptr->cnorm;
521 *ierror=ptr->ierror;
522
523 t8=dsecnd()-t8;
524
525 //printf("%d %d %d %d %d %f %f %f %f %f %f %f %f\n",hash,unlocked,found,omp_get_thread_num(),nbad,t1,t2,t3,t4,t5,t6,t7,t8);
526 return fillval;
527 }
528
529 int fgap_fill(
530 struct fill_struct *pars,
531 float *im,
532 int nx,
533 int ny,
534 int nlead,
535 unsigned char *mask,
536 float *cnorm,
537 float *ierror
538 )
539 {
540 arta 1.1 const unsigned char maskgood=0; // Value of good entries in masks
541 const unsigned char maskfill=1; // Value of entries to be filled
542 const unsigned char masknofill=2; // Value of entries not to be filled
543 int malign=32;
544 //const int ione = 1;
545 int i,j,i1,j1,i0,j0,k;
546 int order,order2;
547 int ngood,nfill;
548 double *im1;
549 int *mask1;
550 float cnorm1,ierror1;
551 int nsample;
552 struct work_struct work;
|
553 schou 1.2 int nofill; // Number of pixels with maskfill=1 that could not be filled.
|
554 arta 1.1 //double t0;
555
556 order=pars->order;
557 order2=order/2;
558 nsample=order*order;
559
560 /*
561 int *fcount;
562 fcount=(int *)(MKL_malloc(ny*sizeof(int),malign));
563 t0=dsecnd();
564 #pragma omp parallel for private (i,j,nfill)
565 for (j=0; j<ny; j++) {
566 nfill=0;
567 for (i=0; i<nx; i++) {
568 if (mask[j*nlead+i]==maskfill) nfill++;
569 }
570 fcount[j]=nfill;
571 }
572 printf("Count: %f\n",dsecnd()-t0);
573 MKL_free(fcount);
574 */
575 arta 1.1
576 nfill=0;
|
577 schou 1.2 nofill=0;
|
578 arta 1.1 #pragma omp parallel default(none) \
579 private (im1,mask1,i,i0,i1,j,j0,j1,k,cnorm1,ierror1,ngood/*,t0*/) \
580 shared(pars,im,nx,ny,nlead,mask,cnorm,ierror) \
581 shared(malign,order,order2) \
582 private(work) shared(nsample) \
|
583 schou 1.2 reduction(+:nfill) \
584 reduction(+:nofill)
|
585 arta 1.1 { // Needed to define parallel region
586
587 im1=(double *)(MKL_malloc(order*order*sizeof(double),malign));
588 mask1=(int *)(MKL_malloc(order*order*sizeof(int),malign));
589
590 malloc_work(nsample,&work);
591 //printf("malloc %d %lx %lx\n",omp_get_thread_num(),im1,mask1);
592
593 //#pragma omp for schedule(dynamic,1)
594 //for (j=0; j<ny; j++) {
595 // for (i=0; i<nx; i++) {
596 #pragma omp for schedule(dynamic,256)
597 for (k=0;k<nx*ny;k++) {{
598 //t0=dsecnd();
599 i=k%nx;
600 j=k/nx;
601 j0=j-order2;
602 if (mask[j*nlead+i]==maskgood) {
603 cnorm1=1.0f; // Random noise is 1
604 ierror1=0.0f; // Interpolation error is 0
605 }
606 arta 1.1 if (mask[j*nlead+i]==maskfill) { // Fill if mask=1. Skip for mask>1.
607 i0=i-order2;
608 ngood=0;
609 for (j1=j0;j1<(j0+order);j1++) {
610 for (i1=i0;i1<(i0+order);i1++) {
611 if ((i1>=0) && (i1<nx) && (j1>=0) && (j1<ny)) {
612 mask1[i1-i0+(j1-j0)*order]=mask[i1+j1*nlead];
613 if (mask[j1*nlead+i1]==maskgood) {
614 im1[ngood]=im[i1+j1*nlead];
615 ngood=ngood+1;
616 }
617 } // Inside original image
618 else { // Outside original image
619 mask1[i1-i0+(j1-j0)*order]=1; // Point is missing
620 }
621 }
622 }
623 if (ngood != 0) { // Have at least one good point
624 im[j*nlead+i]=(float)fill_point(pars,mask1,im1,&cnorm1,&ierror1,&work);
625 }
626 else { // Otherwise quietly leave unchanged, except set errors
627 arta 1.1 cnorm1=0.0f;
|
628 schou 1.5 ierror1=2.0f;
|
629 schou 1.2 nofill=nofill+1;
|
630 arta 1.1 }
631 nfill=nfill+1;
632 } // end if mask==1
633 if (mask[j*nlead+i]==masknofill) { // Don't fill.
634 cnorm1=0.0f;
|
635 schou 1.5 ierror1=2.0f;
|
636 arta 1.1 }
637 cnorm[j*nlead+i]=cnorm1;
638 ierror[j*nlead+i]=ierror1;
639 // cnorm[j*nlead+i]=dsecnd()-t0;
640 // ierror[j*nlead+i]=omp_get_thread_num();
641 } // for i
642 } // for j
643 //printf("free %d %lx %lx\n",omp_get_thread_num(),im1,mask1);
644 MKL_free(im1);
645 MKL_free(mask1);
646 free_work(nsample,&work);
647 //printf("done %d %lx %lx\n",omp_get_thread_num(),im1,mask1);
648 } // End of parallel region
|
649 schou 1.8 //printf("nfill nofill %d %d\n",nfill,nofill);
|
650 arta 1.1
|
651 schou 1.2 if (nofill != 0) {
652 return 1; // Could not fill all pixels
653 }
654 else {
655 return 0; // All is good
656 }
|
657 arta 1.1 }
|
658 schou 1.9
659 char *gapfill_version() // Returns CVS version of gapfill.c
660 {
|
661 arta 1.10 return strdup("$Id: gapfill.c,v 1.9 2010/09/02 18:29:49 schou Exp $");
|
662 schou 1.9 }
663
664
|
Return to
gapfill.c
CVS log
Up to [Development] / JSOC / proj / libs / interpolate