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.6 #include <sys/param.h>
|
9 arta 1.1 #include "finterpolate.h"
|
10 arta 1.6
|
11 arta 1.1
12 #define minval(x,y) (((x) < (y)) ? (x) : (y))
13 #define maxval(x,y) (((x) < (y)) ? (y) : (x))
|
14 arta 1.6
15 #define kInterpDataDir "proj/libs/interpolate/data"
16
|
17 arta 1.1 #define fint_wiener 1
18 #define fint_linear 2
19 #define fint_cubic_conv 3
20
|
21 schou 1.3 static double sinc(double x)
|
22 arta 1.1 {
23 if (fabs(x) < (1.e-10))
24 return 1.;
25 else
26 return sin(M_PI*x)/(M_PI*x);
27 }
28
29 int init_finterpolate_wiener_old(
30 struct fint_struct *pars,
31 int order,
32 int edgemode, // 0 to go as far as you can with symmetric kernel
33 // Otherwise go further (as set by extrapolate)
34 float extrapolate, // How far to extrapolate
35 int minorder, // Minimum order to use when approaching edge or beyond
|
36 arta 1.6 int nconst, // Number of polynomial constraints
37 // 0 None, 1 for norm, 2 for linear preserved, etc.
38 const char *path // to data files read by this function.
|
39 arta 1.1 )
40 {
41 int status;
42 int table=0;
|
43 arta 1.6 char filename[PATH_MAX];
|
44 arta 1.1
|
45 arta 1.6 snprintf(filename, sizeof(filename), "%s/%s/acor1d_80x100_double.txt", path, kInterpDataDir);
46 status=init_finterpolate_wiener(pars,order,edgemode,extrapolate,minorder,nconst,table,&filename, path);
|
47 arta 1.1
48 return status;
49 }
50
51 int init_finterpolate_wiener(
52 struct fint_struct *pars,
53 int order,
54 int edgemode, // 0 to go as far as you can with symmetric kernel
55 // Otherwise go further (as set by extrapolate)
56 float extrapolate, // How far to extrapolate
57 int minorder, // Minimum order to use when approaching edge or beyond
58 int nconst, // Number of polynomial constraints
59 // 0 None, 1 for norm, 2 for linear preserved, etc.
60 int cortable, // Which of the hardcoded tables to use.
61 // 0 To use table pointed to by filenamep
|
62 arta 1.6 char **filenamep, // Pointer to name of file to read covariance from.
|
63 arta 1.1 // Set to actual file used if cortable>0
|
64 arta 1.6 const char *path // to data files read by this function.
|
65 arta 1.1 )
66 {
67 const int malign=32;
68 const int nsubin0=100;
69 const int maxoff=40;
70 int i,j,k;
71 int ishift,offset,order2,nextra,nshifts,shift0;
72 FILE *fileptr;
73 double *acor,*a,*rh,*a1b,*a1r,*coeffd;
74 double bta1b,bta1r,c;
75 double *b,*rhc,*bta1b1,*help;
76 double xc,xp,sum;
77 int info;
78 char uplo[] = "U";
79 int ione = 1;
80 double pmin,regul;
81 int minorder1,curorder,imin,imax,icent,is0,i0,nrh;
82
83 pars->method=fint_wiener;
84 pars->kersx=NULL;
85
86 arta 1.1 if ((order%2) != 0) {
87 printf("Order must be even!\n");
88 return 2;
89 }
90 if ((minorder%2) != 0) {
91 printf("Minorder must be even!\n");
92 return 2;
93 }
94 if (nconst < 0) {
95 printf("Number of constraints must be non-negative!\n");
96 return 2;
97 }
98 order2=order/2;
99 if (edgemode==0) { // Only do area with symmetric set of points available
100 nshifts=nsubin0+1;
101 shift0=0;
102 if (extrapolate!=0) {
103 printf("Warning: Can't extrapolate for edgemode==0.\n");
104 extrapolate=0.0f;
105 }
106 }
107 arta 1.1 nextra=maxval(0,extrapolate)*nsubin0+2;
108 shift0=(order2-1)*nsubin0+nextra;
109 nshifts=(order-1)*nsubin0+2*nextra+1;
110 icent=(nshifts-1)/2; // Center of range. icent=shift0+nsubin0/2
111 pars->edgemode=edgemode;
112 pars->extrapolate=extrapolate;
113 pars->order=order;
114 pars->nshifts=nshifts;
115 pars->shift0=shift0;
116 pars->fover=(float)nsubin0;
117 pars->nsub=nsubin0+1;
118 pars->kersx=(float *)(MKL_malloc(nshifts*order*sizeof(float),malign));
119
120 regul=1/400./400.;
121 pmin=regul;
122
123 if ((cortable < 0) || (cortable > 1)) {
124 printf("Illegal cortable passed to finterpolate!\n");
125 return 1;
126 }
127 if (cortable == 0) {
128 arta 1.1 pars->filename=strdup(*filenamep);
129 }
130 if (cortable == 1) {
|
131 arta 1.6 char dfile[PATH_MAX];
132 snprintf(dfile, sizeof(dfile), "%s/%s/acor1d_80x100_double.txt", path, kInterpDataDir);
133 pars->filename=strdup(dfile);
|
134 arta 1.1 }
135 fileptr = fopen (pars->filename, "r");
136 if (fileptr == NULL) {
137 printf("File not found in finterpolate!\n");
138 return 2;
139 }
140 if (cortable != 0) {
141 if (filenamep != NULL) {
142 *filenamep=strdup(pars->filename);
143 }
144 }
145 acor=(double *)(MKL_malloc((nsubin0*maxoff+1)*sizeof(double),malign));
146 for (i=0;i<nsubin0*maxoff+1;i++) fscanf(fileptr,"%lf",acor+i);
147 fclose(fileptr);
148
149 nrh=maxval(1,nconst);
150 a=(double *)(MKL_malloc(order*order*sizeof(double),malign));
151 rh=(double *)(MKL_malloc(nrh*order*sizeof(double),malign));
152 a1r=(double *)(MKL_malloc(nrh*order*sizeof(double),malign));
153 a1b=(double *)(MKL_malloc(nrh*order*sizeof(double),malign));
154 coeffd=(double *)(MKL_malloc(order*sizeof(double),malign));
155 arta 1.1 b=(double *)(MKL_malloc(nrh*order*sizeof(double),malign));
156 rhc=(double *)(MKL_malloc(nrh*sizeof(double),malign));
157 help=(double *)(MKL_malloc(nrh*sizeof(double),malign));
158 bta1b1=(double *)(MKL_malloc(nrh*nrh*sizeof(double),malign));
159 minorder1=minorder;
160 if (edgemode == 0) minorder1=order; // Effectively ignore minorder
161 for (curorder=minorder1;curorder<=order;curorder=curorder+2) {
162 for (i=0;i<curorder;i++) {
163 for (j=i;j<curorder;j++) {
164 offset=nsubin0*(j-i);
165 a[i*curorder+j]=acor[offset];
166 a[j*curorder+i]=acor[offset];
167 }
168 a[i*curorder+i]=a[i*curorder+i]+regul;
169 }
170 dpotrf(uplo,&curorder,a,&curorder,&info); // Cholesky decomposition
171
172 imin=icent-(order-curorder+1)*nsubin0/2;
173 imax=icent+(order-curorder+1)*nsubin0/2;
174 if (curorder == minorder) { // Extrapolate using minorder
175 imin=0;
176 arta 1.1 imax=nshifts-1;
177 }
178 for (ishift=imin;ishift<=imax;ishift++) {
179 // Find shift equivalent of first pixel to be used for interpolation
180 // First find remainder
181 is0=(ishift-shift0+nsubin0*2*order) % nsubin0;
182 // Then find pixel
183 is0=ishift-is0-nsubin0*(curorder/2-1);
184 // Then limit to range of valid pixels
185 is0=maxval(is0,icent-(order-1)*nsubin0/2);
186 is0=minval(is0,icent+(order-1)*nsubin0/2-(curorder-1)*nsubin0);
187 i0=(is0-icent+(order-1)*nsubin0/2)/nsubin0;
188 //printf("%d %d %d %d %d %d\n",order,imin,imax,ishift,is0,i0);
189 for (i=0;i<curorder;i++) {
190 offset=ishift-(shift0+nsubin0*(i+i0-order/2+1));
191 if (offset<0) offset=-offset;
192 rh[i]=acor[offset]+pmin*sinc((offset+0.0)/nsubin0);
193 } // i
194
195 switch (nconst) {
196 case 0: // No constraints
197 arta 1.1 for (i=0;i<curorder;i++) coeffd[i]=rh[i];
198 dpotrs(uplo,&curorder,&ione,a,&curorder,coeffd,&curorder,&info);
199 break;
200 case 1: // Unit norm. Some tricks applied.
201 for (i=0;i<curorder;i++) a1b[i]=1.;
202 dpotrs(uplo,&curorder,&ione,a,&curorder,a1b,&curorder,&info);
203 bta1b=0.;
204 for (i=0;i<curorder;i++) bta1b=bta1b+a1b[i];
205 for (i=0;i<curorder;i++) a1r[i]=rh[i];
206 dpotrs(uplo,&curorder,&ione,a,&curorder,a1r,&curorder,&info);
207 bta1r=0.0;
208 for (i=0;i<curorder;i++) bta1r=bta1r+a1r[i];
209 c=(bta1r-1.)/bta1b;
210 for (i=0;i<curorder;i++) coeffd[i]=a1r[i]-c*a1b[i];
211 break;
212 default: // General case
213 for (i=0;i<curorder;i++) {
214 offset=ishift-(shift0+nsubin0*(i+i0-order/2+1));
215 xc=(offset+0.0)/nsubin0;
216 xp=1.0;
217 for (j=0;j<nconst;j++) {
218 arta 1.1 b[j*curorder+i]=xp;
219 xp=xp*xc;
220 }
221 }
222 rhc[0]=1;
223 for (j=1;j<nconst;j++) rhc[j]=0;
224 for (i=0;i<curorder*nconst;i++) a1b[i]=b[i];
225 dpotrs(uplo,&curorder,&nconst,a,&curorder,a1b,&curorder,&info);
226 for (i=0;i<nconst;i++) {
227 for (j=0;j<nconst;j++) {
228 sum=0.0;
229 for (k=0;k<curorder;k++)
230 sum=sum+a1b[i*curorder+k]*b[j*curorder+k];
231 bta1b1[i*nconst+j]=sum;
232 }
233 }
234 dpotrf(uplo,&nconst,bta1b1,&nconst,&info); // Cholesky decomposition
235 for (i=0;i<nconst;i++) {
236 sum=0.0;
237 for (k=0;k<curorder;k++) sum=sum+a1b[i*curorder+k]*rh[k];
238 help[i]=sum-rhc[i];
239 arta 1.1 }
240 dpotrs(uplo,&nconst,&ione,bta1b1,&nconst,help,&nconst,&info);
241 for (i=0;i<curorder;i++) a1r[i]=rh[i];
242 dpotrs(uplo,&curorder,&ione,a,&curorder,a1r,&curorder,&info);
243 for (i=0;i<curorder;i++) {
244 sum=0.0;
245 for (k=0;k<nconst;k++) sum=sum+a1b[k*curorder+i]*help[k];
246 coeffd[i]=a1r[i]-sum;
247 }
248 break;
249 }
250
251 // Copy to output array
252 // First zero unused entries
253 for (i=0;i<order;i++) pars->kersx[ishift*order+i]=0.0f;
254 for (i=0;i<curorder;i++) pars->kersx[ishift*order+i+i0]=(float)coeffd[i];
255 // for (i=0;i<curorder;i++) printf(" %f",coeffd[i]);printf("\n");
256 } // ishift
257 } // order
258 //for (ishift=0;ishift<nshifts;ishift++) {
259 // for (i=0;i<order;i++) printf(" %f",pars->kersx[ishift*order+i]);
260 arta 1.1 // printf("\n");
261 //}
262
263 MKL_free(acor);
264 MKL_free(a);
265 MKL_free(a1r);
266 MKL_free(a1b);
267 MKL_free(rh);
268 MKL_free(coeffd);
269 MKL_free(b);
270 MKL_free(rhc);
271 MKL_free(help);
272 MKL_free(bta1b1);
273
274 return 0;
275 }
276
277 int init_finterpolate_linear(
278 struct fint_struct *pars,
279 float extrapolate // How far to extrapolate
280 )
281 arta 1.1 {
282
283 pars->method=fint_linear;
284 pars->extrapolate=extrapolate;
285
286 return 0;
287 }
288
289 int init_finterpolate_cubic_conv(
290 struct fint_struct *pars,
291 int edgemode, // 0 to go as far as you can with symmetric kernel
292 // Otherwise go further (as set by extrapolate)
293 float extrapolate // How far to extrapolate
294 )
295 {
296
297 pars->method=fint_cubic_conv;
298 pars->edgemode=edgemode;
299 pars->extrapolate=extrapolate;
300
301 return 0;
302 arta 1.1 }
303
304 int free_finterpolate(
305 struct fint_struct *pars
306 )
307 {
308 if (pars->method==fint_wiener) {
309 // Check for null if called after error in init.
310 if (pars->kersx!=NULL) {
311 MKL_free (pars->kersx);
312 free (pars->filename);
313 }
314 }
315
316 return 0;
317 }
318
319 int winterpolate(
|
320 schou 1.2 struct fint_struct *pars,
321 float *image_in,
322 float *xin,
323 float *yin,
324 float *image_out,
325 int nxin,
326 int nyin,
327 int nleadin,
328 int nx,
329 int ny,
330 int nlead,
331 float fillval
|
332 arta 1.1 )
333
334 {
335 int i,j,i1,j1;
336 int order2;
337 float *xker;
338 const int ione = 1;
339 int malign=32;
340 int ixin,iyin,ixin1,iyin1;
341 float *ixins,*iyins,*ixin1s,*iyin1s;
342 float fxin1,fyin1,fxin2,fyin2;
343 float *fxins,*fyins,*xinp,*yinp,*fxin1s,*fyin1s;
344 float *xk1,*xk2,*yk1,*yk2;
345 float *imp;
346 float sum,sum1;
347 float *kersx;
348 int ixmax,iymax;
349 float xmax,ymax;
350 int order,shift0,edgemode;
351 float extrapolate;
352 float x,y,help;
353 arta 1.1
354 order=pars->order;
355 kersx=pars->kersx;
356 shift0=pars->shift0;
357 edgemode=pars->edgemode;
358 extrapolate=pars->extrapolate;
359
360 order2=order/2;
361 ixmax=nxin-order2; // Max index for first element of kernel
362 iymax=nyin-order2; // Max index for first element of kernel
363 xmax=(float)ixmax;
364 ymax=(float)iymax;
365
366 #pragma omp parallel default(none)\
367 private(ixins,iyins,fxins,fyins,ixin1s,iyin1s,fxin1s,fyin1s,xker) \
368 private (i,j,xinp,yinp,ixin,iyin,ixin1,iyin1,fxin1,fyin1,fxin2,fyin2,imp) \
369 private (xk1,xk2,yk1,yk2,sum,sum1,i1,j1,x,y,help) \
370 shared (pars,nlead,nx,ny,xin,yin,nleadin,nxin,nyin,image_in,image_out,order) \
|
371 schou 1.2 shared (malign,order2,kersx,ixmax,iymax,xmax,ymax,shift0,edgemode,extrapolate,fillval)
|
372 arta 1.1 { // Needed to define parallel region
|
373 schou 1.4 ixins=(float *)(MKL_malloc(nx*sizeof(float),malign));
374 iyins=(float *)(MKL_malloc(nx*sizeof(float),malign));
|
375 arta 1.1 fxins=(float *)(MKL_malloc(nx*sizeof(float),malign));
376 fyins=(float *)(MKL_malloc(nx*sizeof(float),malign));
377 ixin1s=(float *)(MKL_malloc(nx*sizeof(float),malign));
378 iyin1s=(float *)(MKL_malloc(nx*sizeof(float),malign));
379 fxin1s=(float *)(MKL_malloc(nx*sizeof(float),malign));
380 fyin1s=(float *)(MKL_malloc(nx*sizeof(float),malign));
381
382 xker=(float *)(MKL_malloc(order*sizeof(float),malign));
383
384 // This awful code basically hardcodes the common cases and allows the
385 // compiler to optimize the code for a significant speed improvement
386 switch(order) {
387 case 2:
388 #define OOO1 2
389 #define OOO2 1
390 #include "finterpolate.include"
391 #undef OOO1
392 #undef OOO2
393 break;
394 case 4:
395 #define OOO1 4
396 arta 1.1 #define OOO2 2
397 #include "finterpolate.include"
398 #undef OOO1
399 #undef OOO2
400 break;
401 case 6:
402 #define OOO1 6
403 #define OOO2 3
404 #include "finterpolate.include"
405 #undef OOO1
406 #undef OOO2
407 break;
408 case 8:
409 #define OOO1 8
410 #define OOO2 4
411 #include "finterpolate.include"
412 #undef OOO1
413 #undef OOO2
414 break;
415 case 10:
416 #define OOO1 10
417 arta 1.1 #define OOO2 5
418 #include "finterpolate.include"
419 #undef OOO1
420 #undef OOO2
421 break;
422 case 12:
423 #define OOO1 12
424 #define OOO2 6
425 #include "finterpolate.include"
426 #undef OOO1
427 #undef OOO2
428 break;
429 case 14:
430 #define OOO1 14
431 #define OOO2 7
432 #include "finterpolate.include"
433 #undef OOO1
434 #undef OOO2
435 break;
436 case 16:
437 #define OOO1 16
438 arta 1.1 #define OOO2 8
439 #include "finterpolate.include"
440 #undef OOO1
441 #undef OOO2
442 break;
443 case 18:
444 #define OOO1 18
445 #define OOO2 9
446 #include "finterpolate.include"
447 #undef OOO1
448 #undef OOO2
449 break;
450 case 20:
451 #define OOO1 20
452 #define OOO2 10
453 #include "finterpolate.include"
454 #undef OOO1
455 #undef OOO2
456 break;
457 case 22:
458 #define OOO1 22
459 arta 1.1 #define OOO2 11
460 #include "finterpolate.include"
461 #undef OOO1
462 #undef OOO2
463 break;
464 case 24:
465 #define OOO1 24
466 #define OOO2 12
467 #include "finterpolate.include"
468 #undef OOO1
469 #undef OOO2
470 break;
471 case 26:
472 #define OOO1 26
473 #define OOO2 13
474 #include "finterpolate.include"
475 #undef OOO1
476 #undef OOO2
477 break;
478 case 28:
479 #define OOO1 28
480 arta 1.1 #define OOO2 14
481 #include "finterpolate.include"
482 #undef OOO1
483 #undef OOO2
484 break;
485 case 30:
486 #define OOO1 30
487 #define OOO2 15
488 #include "finterpolate.include"
489 #undef OOO1
490 #undef OOO2
491 break;
492 case 32:
493 #define OOO1 32
494 #define OOO2 16
495 #include "finterpolate.include"
496 #undef OOO1
497 #undef OOO2
498 break;
499 // If not hardcoded do general calculation
500 default:
501 arta 1.1 #define OOO1 order
502 #define OOO2 order2
503 #include "finterpolate.include"
504 #undef OOO1
505 #undef OOO2
506 break;
507 }
508
509 MKL_free(xker);
510 MKL_free(ixins);
511 MKL_free(iyins);
512 MKL_free(fxins);
513 MKL_free(fyins);
514 MKL_free(ixin1s);
515 MKL_free(iyin1s);
516 MKL_free(fxin1s);
517 MKL_free(fyin1s);
518 } // End of parallel region
519
520 return 0;
521 }
522 arta 1.1
523 int linterpolate(
524 float *image_in,
525 float *xin,
526 float *yin,
527 float *image_out,
528 int nxin,
529 int nyin,
530 int nleadin,
531 int nx,
532 int ny,
533 int nlead,
534 float extrapolate,
535 float fillval
536 )
537
538 {
539 int i,j;
540 int malign=32;
541 int ixin,iyin;
542 float *ixins,*iyins;
543 arta 1.1 float fxin1,fyin1,fxin2,fyin2;
544 float *xinp,*yinp;
545 float *imp;
546 float xmin,xmax,ymin,ymax;
547
548 xmin=-extrapolate;
549 xmax=nxin-1.0f+extrapolate;
550 ymin=-extrapolate;
551 ymax=nyin-1.0f+extrapolate;
552
553 #pragma omp parallel default(none) \
554 private(ixins,iyins,i,j,xinp,yinp,ixin,iyin,fxin1,fyin1,fxin2,fyin2,imp) \
555 shared(image_in,xin,yin,image_out,extrapolate,fillval) \
556 shared(nx,ny,nlead,nxin,nyin,nleadin) \
557 shared(malign,xmin,xmax,ymin,ymax)
558 { // Needed to define parallel region
|
559 schou 1.4 ixins=(float *)(MKL_malloc(nx*sizeof(float),malign));
560 iyins=(float *)(MKL_malloc(nx*sizeof(float),malign));
|
561 arta 1.1
562 #pragma omp for
563 for (j=0; j<ny; j++) {
564 xinp=xin+j*nlead;
565 vsFloor(nx,xinp,ixins);
566 yinp=yin+j*nlead;
567 vsFloor(nx,yinp,iyins);
568 for (i=0; i<nx; i++) {
569 ixin=(int)ixins[i]; // Integer pixel to interpolate to
570 iyin=(int)iyins[i]; // Integer pixel to interpolate to
571 if ((ixin>=0) && (ixin <(nxin-1)) && (iyin>=0) && (iyin <(nyin-1))) {
572 // Normal case
573
574 fxin1=xinp[i]-ixin;
575 fyin1=yinp[i]-iyin;
576 fxin2=1.0f-fxin1;
577 fyin2=1.0f-fyin1;
578
579 imp=image_in+ixin+nleadin*iyin;
580
581 /* Brute force addition */
582 arta 1.1 image_out[i+nlead*j]=fyin2*(fxin2*imp[0]+fxin1*imp[1])+
583 fyin1*(fxin2*imp[nleadin]+fxin1*imp[nleadin+1]);
584 } // if
585 else {
586 if ((xinp[i]>=xmin) && (xinp[i]<=xmax) && (yinp[i]>=ymin) && (yinp[i]<=ymax)) {
587 // Extrapolating or point exactly at edge
588 // Move ixin and iyin inside of array
589 ixin=maxval(0,minval(ixin,nxin-2));
590 iyin=maxval(0,minval(iyin,nyin-2));
591
592 fxin1=xinp[i]-ixin;
593 fyin1=yinp[i]-iyin;
594 fxin2=1.0f-fxin1;
595 fyin2=1.0f-fyin1;
596
597 imp=image_in+ixin+nleadin*iyin;
598
599 /* Brute force addition */
600 image_out[i+nlead*j]=fyin2*(fxin2*imp[0]+fxin1*imp[1])+
601 fyin1*(fxin2*imp[nleadin]+fxin1*imp[nleadin+1]);
602 }
603 arta 1.1 else {
604 // Outside of array
605 image_out[i+nlead*j]=fillval;
606 }
607 }
608 } // i=
609 } //j=
610
611 MKL_free(ixins);
612 MKL_free(iyins);
613 } // End of parallel region
614
615 return 0;
616 }
617
618 int ccinterpolate( // Cubic convolution interpolation
619 float *image_in,
620 float *xin,
621 float *yin,
622 float *image_out,
623 int nxin,
624 arta 1.1 int nyin,
625 int nleadin,
626 int nx,
627 int ny,
628 int nlead,
629 int edgemode,
630 float extrapolate,
631 float fillval
632 )
633
634 {
635 int i,j,i1;
636 float xker[4],yker[4];
637 int malign=32;
638 int ixin,iyin;
639 float *ixins,*iyins;
640 float *xinp,*yinp;
641 float *imp;
642 float sum;
643 float s,s2;
644 float xmin,xmax,ymin,ymax;
645 arta 1.1
646 if (edgemode == 0) {
647 xmin=1.0f;
648 xmax=nxin-2.0f;
649 ymin=1.0f;
650 ymax=nyin-2.0f;
651 }
652 else {
653 xmin=-extrapolate;
654 xmax=nxin-1.0f+extrapolate;
655 ymin=-extrapolate;
656 ymax=nyin-1.0f+extrapolate;
657 }
658
659 #pragma omp parallel default(none) \
660 private(ixins,iyins,i,j,xinp,yinp,ixin,iyin,imp,s,s2,xker,yker,sum,i1) \
661 shared(image_in,xin,yin,image_out) \
662 shared(nxin,nyin,nleadin,nx,ny,nlead,extrapolate,fillval) \
663 shared(malign,xmin,xmax,ymin,ymax)
664 { // Needed to define parallel region
|
665 schou 1.4 ixins=(float *)(MKL_malloc(nx*sizeof(float),malign));
666 iyins=(float *)(MKL_malloc(nx*sizeof(float),malign));
|
667 arta 1.1
668 #pragma omp for
669 for (j=0; j<ny; j++) {
670 xinp=xin+j*nlead;
671 vsFloor(nx,xinp,ixins);
672 yinp=yin+j*nlead;
673 vsFloor(nx,yinp,iyins);
674 for (i=0; i<nx; i++) {
675 ixin=(int)ixins[i]; // Integer pixel to interpolate to
676 iyin=(int)iyins[i]; // Integer pixel to interpolate to
677 if ((ixin>=1) && (ixin <(nxin-2)) && (iyin>=1) && (iyin <(nyin-2))) {
678
679 // RML kernel calculation based on fc3kernel in winterpolate_old.c
680 // Kernels are unnormalized and need to be divided by 2
681
682 s=xinp[i]-ixin;
683 s2 = s*s;
684 xker[0] = -s*(1.0f - s*(2.0f - s));
685 xker[1] = 2.0f - s2*(5.0f - 3.0f*s);
686 xker[2] = s*(1.0f + s*(4.0f - 3.0f*s));
687 xker[3] = -s2*(1.0f-s);
688 arta 1.1
689 s=yinp[i]-iyin;
690 s2 = s*s;
691 yker[0] = -s*(1.0f - s*(2.0f - s));
692 yker[1] = 2.0f - s2*(5.0f - 3.0f*s);
693 yker[2] = s*(1.0f + s*(4.0f - 3.0f*s));
694 yker[3] = -s2*(1.0f-s);
695
696 imp=image_in+ixin-1+nleadin*(iyin-1);
697
698 /* Brute force addition */
699 sum=0.0f;
700 for (i1=0; i1<4; i1++) {
701 sum=sum+yker[i1]*(xker[0]*imp[0]+xker[1]*imp[1]+xker[2]*imp[2]+xker[3]*imp[3]);
702 imp=imp+nleadin;
703 }
704 image_out[i+nlead*j]=sum/4;
705 } // if
706 else {
707 if ((xinp[i]>=xmin) && (xinp[i]<=xmax) && (yinp[i]>=ymin) && (yinp[i]<=ymax)) {
708 // Extrapolating or point exactly at edge
709 arta 1.1 // Move ixin and iyin inside of array
710 ixin=maxval(1,minval(ixin,nxin-3));
711 iyin=maxval(1,minval(iyin,nyin-3));
712 // RML kernel calculation based on fc3kernel in winterpolate_old.c
713 // Kernels are unnormalized and need to be divided by 2
714
715 s=xinp[i]-ixin;
716 s2 = s*s;
717 xker[0] = -s*(1.0f - s*(2.0f - s));
718 xker[1] = 2.0f - s2*(5.0f - 3.0f*s);
719 xker[2] = s*(1.0f + s*(4.0f - 3.0f*s));
720 xker[3] = -s2*(1.0f-s);
721
722 s=yinp[i]-iyin;
723 s2 = s*s;
724 yker[0] = -s*(1.0f - s*(2.0f - s));
725 yker[1] = 2.0f - s2*(5.0f - 3.0f*s);
726 yker[2] = s*(1.0f + s*(4.0f - 3.0f*s));
727 yker[3] = -s2*(1.0f-s);
728
729 imp=image_in+ixin-1+nleadin*(iyin-1);
730 arta 1.1
731 /* Brute force addition */
732 sum=0.0f;
733 for (i1=0; i1<4; i1++) {
734 sum=sum+yker[i1]*(xker[0]*imp[0]+xker[1]*imp[1]+xker[2]*imp[2]+xker[3]*imp[3]);
735 imp=imp+nleadin;
736 }
737 image_out[i+nlead*j]=sum/4;
738 }
739 else {
740 // Outside of array
741 image_out[i+nlead*j]=fillval;
742 }
743 }
744 } // i=
745 } //j=
746
747 MKL_free(ixins);
748 MKL_free(iyins);
749 } // End of parallel region
750
751 arta 1.1 return 0;
752 }
753
754 int finterpolate(
755 struct fint_struct *pars,
756 float *image_in,
757 float *xin,
758 float *yin,
759 float *image_out,
760 int nxin,
761 int nyin,
762 int nleadin,
763 int nx,
764 int ny,
765 int nlead,
766 float fillval
767 )
768 {
769 int status;
770
771 switch (pars->method) {
772 arta 1.1 case fint_wiener:
|
773 schou 1.2 status=winterpolate(pars,image_in,xin,yin,image_out,nxin,nyin,nleadin,nx,ny,nlead,fillval);
|
774 arta 1.1 break;
775
776 case fint_linear:
777 status=linterpolate(image_in,xin,yin,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->extrapolate,fillval);
778 break;
779
780 case fint_cubic_conv:
781 status=ccinterpolate(image_in,xin,yin,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->edgemode,pars->extrapolate,fillval);
782 break;
783
784 default:
785 status=1;
786 break;
787 }
788
789 return status;
790
791 }
792
|
793 schou 1.5 char *finterpolate_version() // Returns CVS version of finterpolate.c
794 {
|
795 arta 1.6 return strdup("$Id: finterpolate.c,v 1.5 2010/09/02 18:29:49 schou Exp $");
|
796 schou 1.5 }
797
798
|
Return to
finterpolate.c
CVS log
Up to [Development] / JSOC / proj / libs / interpolate