1 arta  1.1 #ifdef ICCCOMP
  2           #pragma warning (disable : 1572)
  3           #endif
  4           
  5           #pragma omp for schedule(guided)
  6             for (j=0; j<ny; j++) {
  7               xinp=xin+j*nlead;
  8               yinp=yin+j*nlead;
  9           
 10               vsFloor(nx,xinp,ixins);
 11               vsSub(nx,xinp,ixins,fxins);
 12               sscal(&nx,&pars->fover,fxins,&ione);
 13               vsModf(nx,fxins,ixin1s,fxin1s);
 14           
 15               vsFloor(nx,yinp,iyins);
 16               vsSub(nx,yinp,iyins,fyins);
 17               sscal(&nx,&pars->fover,fyins,&ione);
 18               vsModf(nx,fyins,iyin1s,fyin1s);
 19           
 20               if (edgemode==0) { // Symmetric kernels and no extrapolation
 21           // This loop takes care of the case where the desired point
 22 arta  1.1 // is at the right or top edge of the image
 23                 for (i=0; i<nx; i++) {
 24                   if (xinp[i]==xmax) {
 25                     ixins[i]=(float)nxin-OOO2-1;
 26                     ixin1s[i]=pars->fover-1;
 27                     fxin1s[i]=1.0f;
 28                   }
 29                   if (yinp[i]==ymax) {
 30                     iyins[i]=(float)nyin-OOO2-1;
 31                     iyin1s[i]=pars->fover-1;
 32                     fyin1s[i]=1.0f;
 33                   }
 34                 }
 35               }
 36               else { // Allow asymmetrix kernels and extrapolation
 37                 for (i=0; i<nx; i++) {
 38                   x=xinp[i];
 39                   if ((x<(OOO2-1)) && (x>=(-extrapolate))) {
 40                     ixins[i]=(float)(OOO2-1);
 41                     help=pars->fover*(x-ixins[i]);
 42                     ixin1s[i]=(float)floor(help);
 43 arta  1.1           fxin1s[i]=(float)(help-floor(help));
 44                   }
 45                   if ((x>=xmax) && (x<=(nxin-1+extrapolate))) {
 46                     ixins[i]=(float)ixmax-1;
 47                     help=pars->fover*(x-ixins[i]);
 48                     ixin1s[i]=(float)floor(help);
 49                     fxin1s[i]=(float)(help-floor(help));
 50                   }
 51                   y=yinp[i];
 52                   if ((y<(OOO2-1)) && (y>=(-extrapolate))) {
 53                     iyins[i]=(float)(OOO2-1);
 54                     help=pars->fover*(y-iyins[i]);
 55                     iyin1s[i]=(float)floor(help);
 56                     fyin1s[i]=(float)(help-floor(help));
 57                   }
 58                   if ((y>=ymax) && (y<=(nyin-1+extrapolate))) {
 59                     iyins[i]=(float)iymax-1;
 60                     help=pars->fover*(y-iyins[i]);
 61                     iyin1s[i]=(float)floor(help);
 62                     fyin1s[i]=(float)(help-floor(help));
 63                   }
 64 arta  1.1       }
 65               }
 66           
 67               for (i=0; i<nx; i++) {
 68                 ixin=(int)ixins[i]; // Integer pixel to interpolate to
 69                 iyin=(int)iyins[i]; // Integer pixel to interpolate to
 70                 if ((ixin>=(OOO2-1)) && (ixin<ixmax) && (iyin>=(OOO2-1)) && (iyin<iymax)) {
 71           //      count=count+1;
 72           
 73                   fxin1=fxin1s[i];
 74                   ixin1=(int)ixin1s[i]+shift0;
 75                   fyin1=fyin1s[i];
 76                   iyin1=(int)iyin1s[i]+shift0;
 77                   fxin2=1.0f-fxin1;
 78                   fyin2=1.0f-fyin1;
 79             
 80           /* Brute force addition */
 81           
 82                   xk1=kersx+ixin1*pars->order;
 83                   xk2=xk1+pars->order;
 84                   yk1=kersx+iyin1*pars->order;
 85 arta  1.1         yk2=yk1+pars->order;
 86                   imp=image_in+ixin-OOO2+1+nleadin*(iyin-OOO2+1);
 87           
 88           /*
 89           if ((i==4095) && (j==0)) {
 90             printf("%d %d %d %d\n",ixin,iyin,ixin1,iyin1);
 91             printf("%f %f %f %f\n",xinp[i],yinp[i],ixin1s[i],iyin1s[i]);
 92             printf("%f %f %f %f\n",fxin1,fyin1,fxin2,fyin2);
 93             printf("%g %g %g %g\n",xk1[0],xk1[1],xk2[0],xk2[1]);
 94             printf("%g %g %g %g\n",yk1[0],yk1[1],yk2[0],yk2[1]);
 95           }
 96           */
 97           
 98                   for (i1=0; i1<OOO1; i1++) xker[i1]=fxin2*xk1[i1]+fxin1*xk2[i1];
 99                   sum=0.0f;
100                   for (i1=0; i1<OOO1; i1++) {
101                     sum1=0.0f;
102                     for (j1=0; j1<OOO1; j1=j1+2) sum1=sum1+xker[j1]*imp[j1]+xker[j1+1]*imp[j1+1];
103                     sum=sum+sum1*(fyin2*yk1[i1]+fyin1*yk2[i1]);
104                     imp=imp+nleadin;
105                   }
106 arta  1.1         image_out[i+nlead*j]=sum;
107                 } // if

108 schou 1.2       else {
109                   image_out[i+nlead*j]=fillval;
110                 }

111 arta  1.1     } // i=
112             } //j=
113           
114           
115           #ifdef ICCCOMP
116           #pragma warning (default : 1572)
117           #endif