#ifdef ICCCOMP
#pragma warning (disable : 1572)
#endif

#pragma omp for schedule(guided)
  for (j=0; j<ny; j++) {
    xinp=xin+j*nlead;
    yinp=yin+j*nlead;

    vsFloor(nx,xinp,ixins);
    vsSub(nx,xinp,ixins,fxins);
    sscal(&nx,&pars->fover,fxins,&ione);
    vsModf(nx,fxins,ixin1s,fxin1s);

    vsFloor(nx,yinp,iyins);
    vsSub(nx,yinp,iyins,fyins);
    sscal(&nx,&pars->fover,fyins,&ione);
    vsModf(nx,fyins,iyin1s,fyin1s);

    if (edgemode==0) { // Symmetric kernels and no extrapolation
// This loop takes care of the case where the desired point
// is at the right or top edge of the image
      for (i=0; i<nx; i++) {
        if (xinp[i]==xmax) {
          ixins[i]=(float)nxin-OOO2-1;
          ixin1s[i]=pars->fover-1;
          fxin1s[i]=1.0f;
        }
        if (yinp[i]==ymax) {
          iyins[i]=(float)nyin-OOO2-1;
          iyin1s[i]=pars->fover-1;
          fyin1s[i]=1.0f;
        }
      }
    }
    else { // Allow asymmetrix kernels and extrapolation
      for (i=0; i<nx; i++) {
        x=xinp[i];
        if ((x<(OOO2-1)) && (x>=(-extrapolate))) {
          ixins[i]=(float)(OOO2-1);
          help=pars->fover*(x-ixins[i]);
          ixin1s[i]=(float)floor(help);
          fxin1s[i]=(float)(help-floor(help));
        }
        if ((x>=xmax) && (x<=(nxin-1+extrapolate))) {
          ixins[i]=(float)ixmax-1;
          help=pars->fover*(x-ixins[i]);
          ixin1s[i]=(float)floor(help);
          fxin1s[i]=(float)(help-floor(help));
        }
        y=yinp[i];
        if ((y<(OOO2-1)) && (y>=(-extrapolate))) {
          iyins[i]=(float)(OOO2-1);
          help=pars->fover*(y-iyins[i]);
          iyin1s[i]=(float)floor(help);
          fyin1s[i]=(float)(help-floor(help));
        }
        if ((y>=ymax) && (y<=(nyin-1+extrapolate))) {
          iyins[i]=(float)iymax-1;
          help=pars->fover*(y-iyins[i]);
          iyin1s[i]=(float)floor(help);
          fyin1s[i]=(float)(help-floor(help));
        }
      }
    }

    for (i=0; i<nx; i++) {
      ixin=(int)ixins[i]; // Integer pixel to interpolate to
      iyin=(int)iyins[i]; // Integer pixel to interpolate to
      if ((ixin>=(OOO2-1)) && (ixin<ixmax) && (iyin>=(OOO2-1)) && (iyin<iymax)) {
//      count=count+1;

        fxin1=fxin1s[i];
        ixin1=(int)ixin1s[i]+shift0;
        fyin1=fyin1s[i];
        iyin1=(int)iyin1s[i]+shift0;
        fxin2=1.0f-fxin1;
        fyin2=1.0f-fyin1;
  
/* Brute force addition */

        xk1=kersx+ixin1*pars->order;
        xk2=xk1+pars->order;
        yk1=kersx+iyin1*pars->order;
        yk2=yk1+pars->order;
        imp=image_in+ixin-OOO2+1+nleadin*(iyin-OOO2+1);

/*
if ((i==4095) && (j==0)) {
  printf("%d %d %d %d\n",ixin,iyin,ixin1,iyin1);
  printf("%f %f %f %f\n",xinp[i],yinp[i],ixin1s[i],iyin1s[i]);
  printf("%f %f %f %f\n",fxin1,fyin1,fxin2,fyin2);
  printf("%g %g %g %g\n",xk1[0],xk1[1],xk2[0],xk2[1]);
  printf("%g %g %g %g\n",yk1[0],yk1[1],yk2[0],yk2[1]);
}
*/

        for (i1=0; i1<OOO1; i1++) xker[i1]=fxin2*xk1[i1]+fxin1*xk2[i1];
        sum=0.0f;
        for (i1=0; i1<OOO1; i1++) {
          sum1=0.0f;
          for (j1=0; j1<OOO1; j1=j1+2) sum1=sum1+xker[j1]*imp[j1]+xker[j1+1]*imp[j1+1];
          sum=sum+sum1*(fyin2*yk1[i1]+fyin1*yk2[i1]);
          imp=imp+nleadin;
        }
        image_out[i+nlead*j]=sum;
      } // if
      else {
        image_out[i+nlead*j]=fillval;
      }
    } // i=
  } //j=


#ifdef ICCCOMP
#pragma warning (default : 1572)
#endif