/*****************************************************************************/
/*									     */
/* NAME:     swarp					 		     */
/* FUNCTION: correct geometric distortion in a sequence of images	     */
/*	     according to parameters passed in command line. Params	     */
/*	     are calculated by program dots.c				     */
/* USAGE:    see below                                                       */
/* COMPILE:  cc swarp.c -lm -o swarp		       			     */
/* RUN:	     swarp -w -x 156.5 -y 113 -a .9666354 -b 9.39331e-5 -c 2.21307e-6		-o ~/lens/zzz ~/lens/spots					    */
/* AUTHOR:   velocity map & putting all together - Lee Campbell 	     */
/*	     file I/O stuff - Stephen Intille				     */
/*	     command line parsing - public domain			     */
/*	     bilinear interp - 	Sourabh Niyogi				     */
/*	     bicubic interp - John Y A Wang & Eero Simoncelli		     */
/* CHANGES:  11/29/93 polynomial changed from rnew = ar + br^2 + cr^3 	     */
/*	     to rnew = ar + br^3 + cr^5 because only odd powers are	     */
/*	     needed in radial polynomial				     */
/* steve:    fixed the part that reads multiple files (data_files)           */
/*****************************************************************************/

/* Reads a raw image file stored as a 'data' file and 'descriptor' file      */
/* and calculates center of aberration and radial polynomial for correcting  */
/* aberration.     */
/* Should work for float and char images. Works on sequences. */

/* swarp	 [-w] [-z] [-l int] [-a int] [-b int] [-c int] [-d]
                       [-e outfile] [-o outfile] [infile]
 
   -w             : Allow outfiles to overwrite existing files (if used,
                    this must be the first flag
   -l		  : use bilinear interpolation instead of default
		    bicubic interpolation
   -x float       : x component of center of aberration
   -y float       : y component of center of aberration
   -a float       : first (r^1) polynomial coeff
   -b float       : second (r^3) polynomial coeff
   -c float       : third (r^5) polynomial coeff
   -d             : assume image has been de-interlaced and is FAT

   -o outfile     : Save the file after lines have been removed in 
                    specified directory
   infile         : Name of original file
                                                                          */
#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/dir.h>
#include <sys/stat.h>

extern int optind;
extern char *optarg;
#define OPTIONS "wldx:y:a:b:c:o:"
#define NIF 1		 /* maximal number of input files. only allow one */
char *cmd;               /* used to store the name of the program         */
char *cmt;               /* used to store the entire command line         */
char *usage = "[-w] [-l] [-x float] [-y float] [-a float] [-b float] [-c float] [-o outfile] [infile]";

/*** globals ***/

/*---------------------------------------------------------------------------*/

get_input_file_name(ifile, ifile_desc, argv)
  char **ifile;
  char **ifile_desc;
  char *argv[];
{
  *ifile = (char *) malloc (sizeof(char) * (strlen(argv[optind]) + 10));
  *ifile_desc = (char *) malloc (sizeof(char) * (strlen(argv[optind]) + 12));
  
  strcpy(*ifile, argv[optind]);
  strcpy(*ifile_desc, argv[optind]);
  strcat(*ifile, "/data");
  strcat(*ifile_desc, "/descriptor");
}  

/*-----------------------------------------------------------*/

get_output_file_name(ofile, ofile_desc, overwrite_files)
  char **ofile;
  char **ofile_desc;
  int overwrite_files;
{
  FILE *fp;
  
  struct stat stbuf;         /* used to store result of stat operation */
                             /* used to check if files already exist   */
  char command[256];         /* stores unix command */
  
  *ofile = (char *) malloc (sizeof(char) * (strlen(optarg) + 11));
  *ofile_desc = (char *) malloc (sizeof(char) * (strlen(optarg) + 12));
  
  strcpy(*ofile, optarg);
  
  if (stat(*ofile, &stbuf) == -1)
	{
	      fprintf(stderr,"Couldn't access %s. Creating directory.\n",
		      *ofile);
	      sprintf(command,"mkdir %s", *ofile);
	      if ((fp = popen(command,"r")) == NULL)
		fprintf(stderr,"Couldn't run %s.", command);
	      close(fp);
	    }
	  else
	    if ((stbuf.st_mode & S_IFMT) != S_IFDIR)
	      {
		fprintf(stderr,
			"File %s exists. Aborted. No dir created.\n",
			*ofile);
		exit(1);
	      }
	/*  $$$  */
	  strcpy(*ofile_desc, *ofile);
	  strcat(*ofile_desc, "/descriptor");
	  strcat(*ofile, "/data");
	  
	  if (!overwrite_files)
	    {
	      if (stat(*ofile, &stbuf) != -1)
		{
		  fprintf(stderr,"File %s exists. Aborted.\n",*ofile);
		  exit(1);
		}
	      if (stat(*ofile_desc, &stbuf) != -1)
		{
		  fprintf(stderr,"File %s exists. Aborted.\n",*ofile_desc);
		  exit(1);
		}
	}
}	  

/*-----------------------------------------------------------*/

get_frame_params(ifile_desc, nframe, sets, bytes_pixel, xdim, ydim)
  /* Read input file descriptor to find size and file type */
  char *ifile_desc;
  int *nframe, *sets, *bytes_pixel;
  int *xdim, *ydim;
{
  char label[20], junk[20], junk2[20];
  char line[80];
  FILE *ifs_desc = stdin;

  if ((ifs_desc = fopen(ifile_desc, "r")) == NULL)
    {fprintf(stderr,
	     "%s: Cannot open x descriptor input file\n", cmd);
     exit(1);}

  while (fgets(line,80,ifs_desc) != NULL)
    {
      fprintf(stderr,"swarp steve: in while fgets\n");
      sscanf(line,"%s",label);
      if (!strcmp (label, "(_data_files"))
	sscanf(line,"%s%d%s", junk, nframe, junk2);
      else  if (!strcmp (label, "(_dimensions"))
	sscanf(line,"%s%d%d%s", junk, ydim, xdim, junk2);
      else  if (!strcmp (label, "(_data_sets"))
	sscanf(line,"%s%d%s", junk, sets, junk2);
      else  if (!strcmp (label, "(_data_type"))
	{
	  sscanf(line,"%s%s", junk, junk2);
	  if ((!strcmp (junk2, "\"unsigned_1\")")) ||

	      (!strcmp (junk2, "unsigned_1")))
	    *bytes_pixel = 1;
	  if ((!strcmp (junk2, "\"float\")")) ||
	      (!strcmp (junk2, "float")))
	    *bytes_pixel = 4;
	}
    }
  fclose (ifs_desc);
}

/*-----------------------------------------------------------*/

read_frame_data(ifile, frame, nframe, xdim, ydim, bytes_pixel, data)
     /* get the image */
  char *ifile;
  int frame, nframe;
  int *data;
{
  char *current_ifile;
  char junk[20];
  int bytes;
  FILE *ifs; 
  struct stat stbuf;         /* used to store result of stat operation */
                             /* used to check if files already exist   */
  
      current_ifile = (char *) malloc (sizeof(char) * (strlen(ifile) + 6));
      strcpy(current_ifile, ifile);
      fprintf(stderr,"swarp steve: current_ifile=%s\n",current_ifile);

      if (nframe != 1) {  /* if only one frame, don't tag number on end */
	  strcat(current_ifile, sprintf(junk, "%d",frame - 1));
	  /*steve*/ /*strcat(current_ifile, "steve was here");*/
	  /*steve*/ strcat(current_ifile, junk);
          fprintf(stderr,"swarp steve: junk=%s, frame-1 =%d\n",junk,frame-1);
          fprintf(stderr,"swarp steve: read_frame_data: nframe != 1; nframe=%d\n",nframe);
          fprintf(stderr,"swarp steve: current_ifile=%s\n",current_ifile);
        } 
      fprintf(stderr,"Working on %s", current_ifile);
      fprintf(stderr,".\n");

      /* open input file and read data */

      if (stat(current_ifile, &stbuf) == -1)
	{
	  fprintf(stderr,"Error: Couldn't access %s. \n",
		  current_ifile);
	  exit(1);
	}
      
      if ((ifs = fopen(current_ifile, "r")) == NULL)
	{
	  fprintf(stderr, "%s: Cannot open input file\n", cmd);
	  exit(1);
	}
      
      bytes = fread(data, bytes_pixel, xdim * ydim, ifs);
      fprintf(stderr,"Read %d.\n",bytes);
      
      fclose(ifs);
      free(current_ifile);
}      

/*-----------------------------------------------------------*/

write_frame_data(ofile, frame, nframe, xdim, ydim, bytes_pixel, vimage)
  char *ofile;
  int frame, nframe;
  int xdim, ydim;
  int bytes_pixel;
  int *vimage;
{
  char *current_ofile;
  char junk[20];
  int countout;
  FILE *ofs = stdout; 
	current_ofile = (char *) malloc (sizeof(char) * (strlen(ofile) + 6));
	strcpy(current_ofile, ofile);
	if (nframe != 1) {  /* if only one frame, don't tag number on end */
            fprintf(stderr,"swarp steve: nframe !=1; nframe = %d\n",nframe);
	    strcat(current_ofile, sprintf(junk, "%d",frame - 1));
	    /*steve*/strcat(current_ofile, junk);/*need to do it again*/
           } 
	fprintf(stderr,"assigned file to ofile %s\n", current_ofile);
	  if ((ofs = fopen(current_ofile, "w")) == NULL)
	    {
	      fprintf(stderr, "%s: Cannot open output file\n", cmd);
	      exit(1);
	    }

	countout = fwrite(vimage, bytes_pixel, xdim*ydim, ofs);
	if (countout != xdim*ydim)
	  fprintf(stderr,"output err: wrote %d pixels; should be %d \n",
		   countout, xdim*ydim);
	fclose(ofs);
	free(current_ofile);
}

/*-----------------------------------------------------------*/
  
write_output_params(ofile_desc, nframe, xdim, ydim, bytes_pixel)
  char *ofile_desc;
  int nframe;
  int xdim, ydim;
  int bytes_pixel;
{
  char *dat_type;
  FILE *ofs_desc = stdout;
      if ((ofs_desc = fopen(ofile_desc, "w")) == NULL)
	{fprintf(stderr,
		 "%s: Cannot open x descriptor output file\n", cmd);
	 exit(1);}
  
      /* Write descriptor files */

      dat_type = (char *) malloc(sizeof(char) * 30);

      if (bytes_pixel == 4)
	strcpy(dat_type, "float");
      else
	strcpy(dat_type, "unsigned_1");
  
      fprintf(ofs_desc,
       "(_data_files %d)\n(_channels 1)\n(_dimensions %d %d)\n(_data_type %s)",
	  nframe, ydim, xdim, dat_type);

      fclose(ofs_desc);
}

/*-----------------------------------------------------------*/

main(argc, argv)	/* process command line args */
     int argc;
     char *argv[];
{
  FILE *fp;
  
  int errflag = 0;     /* set to 1 on error */
  int oflag = 0;       /* infiles must be entered */
  
  int overwrite_files = 0;   /* set to 1 (allow) if user specifies -w */
  int bilinear = 0;          /* set to 1 (allow) if user specifies -l */
  int deinterl = 0;	     /* set to 1 (allow) if user specifies -d */
			     /* otherwise default is bicubic */
  struct stat stbuf;         /* used to store result of stat operation */
			     /* used to check if files already exist   */
  char command[256];         /* stores unix command */
  
  char c;              /* c is a temporary character used for breaking the  */
                       /* command line into its individual flags, names, etc*/
  char *ofile = NULL;  /* set default output file name to NULL */
  char *ofile_desc = NULL;
  char *ifile = NULL;  /* set default input  file name to NULL */
  char *ifile_desc = NULL;
  char **datafile;
  char **descfile;
  float x0, y0, a1, a2, a3;  /* aberration parameters */

  cmd = argv[0];		/* save the name of program */

  while ((c = getopt(argc, argv, OPTIONS)) != EOF)
    switch (c) {
      
    case 'w':
      overwrite_files = 1;
      break;
      
    case 'l':
      bilinear = 1;
      break;
      
    case 'x':
      x0 = atof(optarg);
      break;
      
    case 'y':
      y0 = atof(optarg);
      break;
      
    case 'a':
      a1 = atof(optarg);
      break;
      
    case 'b':
      a2 = atof(optarg);
      break;
      
    case 'c':
      a3 = atof(optarg);
      break;
      
    case 'd':
      deinterl = 1;
      break;
      
    case 'o':
      if (!oflag)                        
	{
	  datafile = &ofile;
	  descfile = &ofile_desc;
	  get_output_file_name(datafile, descfile, overwrite_files);
	  oflag = 1;
	}
      break;
      
    case '?':
      errflag = 1;
      break;
    }
  
  /* check if error was found while parsing command line */
  
  if (errflag			/* error in options */
      || !oflag
      || argc > optind + NIF)	/* we desire only NIF input files. At this */
                                /* point in the execution, getopt should   */
                                /* have parsed all the command line prompts*/
                                /* except the input file names. Therefore, */
                                /* if this relationship is not true, extra */
                                /* flags, arguments, or names were in the  */
                                /* command prompt */
    {fprintf(stderr, "\nUsage: %s %s\n\n", cmd, usage);
     fprintf(stderr, "Read ~elwin/lens/code/swarp.c for description.\n\n");
     exit(1);}
  
  /* get input file */
  
  datafile = &ifile;
  descfile = &ifile_desc;
  get_input_file_name(datafile, descfile, argv);

  /*do the processing*/
  do_it(ifile, ifile_desc, ofile, ofile_desc, x0, y0, a1, a2, a3,
	 bilinear, deinterl);

  
  free(ifile);
  free(ifile_desc);
  free(ofile);
  free(ofile_desc);
}

/*-----------------------------------------------------------*/

/* This part of the code loops through frames in a sequence
   and calls the subroutines that actually do the real work */

do_it(ifile, ifile_desc, ofile, ofile_desc, x0, y0, a1, a2, a3,
       bilinear, deinterl)
      
     char *ifile, *ifile_desc, *ofile, *ofile_desc;
     float x0, y0, a1, a2, a3;
     int bilinear, deinterl;
{
  FILE *ifs; 
  FILE *ifs_desc;
  FILE *ofs = stdout; 
  FILE *ofs_desc = stdout;
  FILE *fp;

  struct stat stbuf;         /* used to store result of stat operation */
                             /* used to check if files already exist   */
  int i, x, y;
  char label[20], junk[20], junk2[20];
  char line[80];
  char command[500];
  unsigned char *indata, *outdata;
  float *vx, *vy;
  int xdim, ydim;
  int bytes, countout;
  int type;
  int nframe;
  int frame;
  int sets;
  int bytes_pixel;
  char *current_ifile;
  char *current_ofile;
  char *dat_type;

  nframe = 1;
  bytes_pixel = 1;      /* default to picture of chars */
  sets = 0;

  get_frame_params(ifile_desc, &nframe, &sets, &bytes_pixel, &xdim, &ydim);

  fprintf(stderr,"Input image size: (%d (y), %d (x)).",ydim, xdim);
  fprintf(stderr," Frames: %d.",nframe);
  fprintf(stderr," Bytes/pixel : %d.\n",bytes_pixel);
  if (bilinear == 1)
      fprintf(stderr," bilinear interpolation method \n");
  else
      fprintf(stderr," bicubic interpolation method \n");
  
  indata = (unsigned char *) malloc(sizeof(char) * bytes_pixel * (xdim*ydim));
  outdata = (unsigned char *) malloc(sizeof(char) * bytes_pixel * (xdim*ydim));
  vx = (float *) malloc(sizeof(float) * (xdim*ydim));
  vy = (float *) malloc(sizeof(float) * (xdim*ydim));
  /* 	x0 = 156.5;	y0 = 113.0;	/* typical warp params */
  /* 	a1 = .9666354;	a2 = 9.393309e-5;	a3 = 2.213069e-6;	/* */
  
 /* make_warp_map_approx(x0, y0, a1, a2, a3, xdim, ydim, vx, vy); */
  make_warp_map(x0, y0, a1, a2, a3, xdim, ydim, vx, vy, deinterl); 

  fprintf(stderr,"Hold on, working...\n");
  for (frame = 1; frame <= nframe; frame++)
    {    
        fprintf(stderr,"swarp steve: frame = %d of %d\n",frame,nframe);
        read_frame_data(ifile, frame, nframe, xdim, ydim, bytes_pixel, indata);
	/* get the image, store it in   indata  */

/******  THIS IS WHERE THE ACTION IS  *******/

      if (bilinear == 1)
	  bilinear_warp(indata, outdata, vx, vy, xdim, ydim);
      else
	  bicube_warp(indata, outdata, vx, vy, xdim, ydim);

/*      makedisplay(xdim, ydim, outdata);   */

      /* write output file image data */
      write_frame_data(ofile, frame, nframe, xdim, ydim, bytes_pixel, outdata);

    }	/* for frame */

  /* Write descriptor files */
  write_output_params(ofile_desc, nframe, xdim, ydim, bytes_pixel);

  free(indata);	 /* "Freedom's just another word for nothin left to lose;" */
  free(outdata); /* "Nothin ain't worth nothin but its free..."		 */
}

/*-----------------------------------------------------------*/

make_warp_map_approx(x0, y0, a1, a2, a3, xdim, ydim,
		     xresult, yresult, deinterl)
/*	A warp map is a velocity map. It says: for the center of
	each pixel in the output picture, where did that point
	come from in the input picture? This approximate version
	does the reverse: for the center of each pixel in the
	input picture, where is that point going in the output
	picture? The approximation is to reverse that value. It
	is a pretty good approximation for small, slowly varying
	warps, but I want to be exact. This is kept around for
	historical and educational purposes only. */
     float x0, y0, a1, a2, a3;
     int   xdim, ydim;
     float *xresult, *yresult;
     int deinterl;
{
  int x, y, cc;
  float xx, yy, r, r2, rnew, frac;

  cc = 0;
  for (y = 0; y < ydim; y++)
    for (x = 0; x < xdim; x++)
      {
        xx = ((float) x) - x0;
        yy = ((float) y) - y0;
        r2 =  (xx*xx + yy*yy);
        r = ((float) sqrt((double) r2));
        /*	rnew = a1*r + a2*r2 + a3*r*r2;  */
	/*	frac = rnew/r;			*/
	/* frac = a1 + a2*r + a3*r2;	replaced by odd-powers poly */
	frac = a1 + a2*r2 + a3*r2*r2;	/* actually, never tested */
	xresult[y * xdim + x] = xx*frac - xx;
	yresult[y * xdim + x] = yy*frac - yy;
	/* if (x == y) fprintf(stderr," %4d %9.5f %9.5f  ",
			    x, xx*frac - xx, yy*frac - yy); /* */
	/* cc++; */
      }
  /* fprintf(stderr," warp map count %d \n", cc); */
}

make_warp_map(x0, y0, a1, a2, a3, xdim, ydim,
	      xresult, yresult, deinterl)
/*	For the center of each pixel in the output picture, where
	did that point come from in the input picture? rguess is
	how far along a radial line (from the center of aberr.)
	pixel p[x,y] in the output picture moved relative to the
	input picture. This code uses a neighboring pixel for a
	first approximation, and Newton's method to improve the
	approx.	*/
/*	11/29/93 changed rguess, f, fprime for odd-powers poly	*/
     float x0, y0, a1, a2, a3;
     int   xdim, ydim;
     float *xresult, *yresult;
     int deinterl;
{
  int x, y, i, cc, pf;
  float xx, yy, r, r2;
  float radius, rguess, rgold;
  float f, fprime;

  cc = 0;
  for (y = 0; y < ydim; y++)
    {
      rguess = 0.0;	/* force first guess approx */
      for (x = 0; x < xdim; x++)
	{
	  xx = ((float) x) - x0;
	  yy = ((float) y) - y0;
	  if (deinterl == 1) yy = ((float) y+y) - y0;
	  radius = ((float) sqrt((double) (xx*xx + yy*yy)));
	  if (fabs((double) rguess) < 1e-4)
	    {
	      /* Use the warp_map_approx() method for first guess if rguess
		 is small. For rest of row, prev sol'n is first guess.  */
	      xx = ((float) x) - x0;
	      yy = ((float) y) - y0;
	      if (deinterl == 1) yy = ((float) y+y) - y0;
	      r2 =  (xx*xx + yy*yy);
	      r = ((float) sqrt((double) r2));
	     /* rguess = r - (a1*r + a2*r2 + a3*r*r2); */
	      rguess = r - (a1 + (a2 + a3*r2)*r2)*r;
	    }
	  rgold = 1.25*rguess;	/* force newton loop to try */
	 /* if (x == y) pf = 1; else pf = 0; */
	 /* if (pf == 1) fprintf(stderr," guess  %10.6f",rguess); */
	  for (i=1; (i<30)&&(fabs(((double)(rgold-rguess)))>1e-5); i++)
	    {	/* newton loop until error is small */
	      r = radius + rguess;
	      r2 = r*r;
	    /*  f = radius - (a1*r + a2*r2 + a3*r*r2); */
	    /*  fprime = -(a1 + 2*a2*r + 3*a3*r2);     */
	      f = radius - (a1 + (a2 + a3*r2)*r2)*r;
	      fprime = -(a1 + (3*a2 + 5*a3*r2)*r2);
	      rgold = rguess;
	      rguess = .9*(rguess - f/fprime);
	    }
	 /* if (pf == 1) fprintf(stderr," %10.6f %10.6f %2d %4d %4d \n",
			       rgold, rguess, i, x, y); /* */
	  i = (y * xdim + x);
	  xresult[i] = -(xx/radius)*rguess;
	  if (deinterl != 1)
	    yresult[i] = -(yy/radius)*rguess;
	  else
	    yresult[i] = -(yy/(radius*2))*rguess;
	  cc++;
	}
    }
  /* fprintf(stderr," warp map count %d \n", cc); */
}

/* bilinear interpolation code from sourabh */
bilinear_warp(im, out, vx, vy, xdim, ydim)
unsigned char  *im, *out;
float *vx;
float *vy;
int  xdim;
int  ydim;
{
  int x, y;
  float xpf, ypf;
  int xpi, ypi;
  int xbound, ybound;
  float dx, dy, ans;

  xbound = xdim-1;
  ybound = ydim-1;  
  for (y=0; y<ydim; y++) 
    for (x=0; x<xdim; x++) {
      xpf = (float)x-*(vx+y*xdim+x);
      ypf = (float)y-*(vy+y*xdim+x);
      xpi = (int)floor(xpf);
      ypi = (int)floor(ypf);
      if ( ( xpi < 0 ) || ( ypi < 0 ) || 
	   ( (xpi+1) > xbound) ||  ( (ypi+1) > ybound) ) 
	*(out+y*xdim+x) = *(im+y*xdim+x);
      else {
	dx  = xpf-(float)xpi;
	dy  = ypf-(float)ypi;
	ans =	(1-dx) * (1-dy) * (*(im + ypi   * xdim + xpi)) + 
		  dx   * (1-dy) * (*(im + ypi   * xdim + xpi + 1)) +
		(1-dx) *   dy   * (*(im+(ypi+1) * xdim + xpi)) + 
		  dx   *   dy   * (*(im+(ypi+1) * xdim + xpi + 1));
	if (ans > 255) ans = 255;
	if (ans < 0) ans = 0;
	*(out+y*xdim+x) = ((unsigned char) (ans + 0.5));
      }	/* */
   }
}

/* John/Eero's warper, for control expt. */
bicube_warp(im, res, wx, wy, xdim, ydim)
  unsigned char *im, *res;     
  float *wx, *wy;
  int xdim, ydim;
{
  register unsigned char *im_ptr; /* accessed 38x in loop */
  register int base;		/* accessed 30x in loop */
  register float frac, temp;	/* accessed 20x in loop */
  register float frac2, frac3;  /* accessed 10x in loop */
  int x_inbounds_1, x_inbounds0, x_inbounds1, x_inbounds2; /* 6x in loop */
  register float x_coeff_1, x_coeff0, x_coeff1, x_coeff2; /* 5x in loop */
  register float res_val;
  register int i,j,trunc;
  register float xfrac,xbase,yfrac,ybase;
  float x_off, y_off;
    
  x_off = (float)((xdim - 1) / 2.0);
  y_off = (float)((ydim - 1) / 2.0);
  
  for(j=0; j<ydim; j++)  /* j is the y var */
    for(i=0; i<xdim; i++) { 
      xfrac = i - *wx++;
      xbase = (int)(xfrac);
      yfrac = j - *wy++;
      ybase = (int)(yfrac);
      
      if((((xbase) >= 0) && ((xbase + 0) < xdim)) &&
	 (((ybase) >= 0) && ((ybase + 0) < ydim))) {
	res_val = 0.0;
	
	frac = xfrac - xbase;
	frac2 = frac*frac;
	frac3 = frac*frac2;
	base  = xbase - 1;
	im_ptr = im + base;
	
	x_inbounds_1 = ((base >= 0) && (base < xdim))?1:0;
	if (x_inbounds_1) x_coeff_1 = (frac2/2.0 - frac/3.0 - frac3/6.0);
	base++;
	x_inbounds0  = ((base >= 0) && (base < xdim))?1:0;
	if (x_inbounds0)  x_coeff0  = (1.0 - frac/2.0 - frac2 + frac3/2.0);
	base++;
	x_inbounds1  = ((base >= 0) && (base < xdim))?1:0;
	if (x_inbounds1)  x_coeff1 = (frac + frac2/2.0 - frac3/2.0);
	base++;
	x_inbounds2  = ((base >= 0) && (base < xdim))?1:0;
	if (x_inbounds2)  x_coeff2 = (frac3/6.0 - frac/6.0);
	
	frac = yfrac - ybase;
	frac2 = frac*frac;
	frac3 = frac*frac2;
	base = ybase - 1 ;
	im_ptr += base*xdim;
	
	if ((base >= 0) && (base < ydim)) {
	  temp = 0.0;
	  if (x_inbounds_1) temp += x_coeff_1 * *im_ptr;
	  im_ptr++;
	  if (x_inbounds0)  temp += x_coeff0  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds1)  temp += x_coeff1  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds2)  temp += x_coeff2  * *im_ptr;
	  im_ptr -= 3;
	  res_val += temp * (frac2/2.0 - frac/3.0 - frac3/6.0);
	}
	base++;  im_ptr += xdim;
	if ((base >= 0) && (base < ydim)) {
	  temp = 0.0;
	  if (x_inbounds_1) temp += x_coeff_1 * *im_ptr;
	  im_ptr++;
	  if (x_inbounds0)  temp += x_coeff0  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds1)  temp += x_coeff1  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds2)  temp += x_coeff2  * *im_ptr;
	  im_ptr -= 3;	    
	  res_val += temp * (1 - frac/2.0 - frac2 + frac3/2.0);
	}
	base++;  im_ptr += xdim; 
	if ((base >= 0) && (base < ydim)) {
	  temp = 0.0;
	  if (x_inbounds_1) temp += x_coeff_1 * *im_ptr;
	  im_ptr++;
	  if (x_inbounds0)  temp += x_coeff0  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds1)  temp += x_coeff1  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds2)  temp += x_coeff2  * *im_ptr;
	  im_ptr -= 3;	    
	  res_val += temp * (frac + frac2/2.0 - frac3/2.0);
	}
	base++;  im_ptr += xdim;
	if ((base >= 0) && (base < ydim)) {
	  temp = 0.0;
	  if (x_inbounds_1) temp += x_coeff_1 * *im_ptr;
	  im_ptr++;
	  if (x_inbounds0)  temp += x_coeff0  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds1)  temp += x_coeff1  * *im_ptr;
	  im_ptr++;
	  if (x_inbounds2)  temp += x_coeff2  * *im_ptr;
	  im_ptr -= 3;	    
	  res_val += temp * (frac3/6.0 - frac/6.0);
	}
	trunc = ((int) (res_val + 0.5));
	if (trunc > 255) trunc = 255;
	if (trunc < 0) trunc = 0;
	*res = trunc;
      } 
      else {
	*res = im[j*xdim + i]; 
      }
      res++;
    }
}

/* */	/* This is the end. */