/*
to compile:  cc ppmswarp.c -lm
to run:  a.out -w -x ___ -y ___ -a ___ -b ___ -c ___ -o pork.s pork.pgm
help: a.out    (isatty is used to determine that nothing is piped in, etc.
pipe usage: cat pork.pgm | a.out > pork.s
if one unnamed argument given, that's the input file
if output is not a tty, that's pipe out
based on elwin's swarp.c, intille's i/o, sourabh's linear interp,
jywang + eero's bicubic interp 

Reads a ppm or pgm image file
and applies a coordinate transformation to it, based on a radial taylor 
series
Should later work for float, etc. images, defined by the .par format .


everything needed is in this one file, no need to link to anything but -lm
which exists on just about any platform.

*/

/* ppmswarp	 [-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
                    equivalent to usage with >! instead of just >
   -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]";

/*** 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,"ppmswarp is 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,"ppmswarp 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,"ppmswarp steve: junk=%s, frame-1 =%d\n",junk,frame-1);
          fprintf(stderr,"ppmswarp steve: read_frame_data: nframe != 1; nframe=%d\n",nframe);
          fprintf(stderr,"ppmswarp 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);
}      

/*-----------------------------------------------------------*/
/* ndf: output needs to be of same form as input.
   if input was P6
                640
                480
                255
   so should output.
   later there may be desire to have output size be just a little
   bigger or just a little smaller than input, like analogous to
   pchirp2nocrop, but don't worry that for now */


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,"ppmswarp 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 -o porks.pgm pork.pgm\n\n", cmd, usage);
     fprintf(stderr, "cat pork.pgm | %s %s > porks.pgm\n",cmd, usage);
     fprintf(stderr, "cat pork.ppm | %s %s > porks.ppm\n",cmd, usage);
     fprintf(stderr, "cat pork.ppm | %s > porks.ppm\n",cmd);
     fprintf(stderr, "cat pork.ppm | %s -o porks.ppm\n",cmd);
     exit(1);}
  
  /* get input file ...  ndf: change to read P5 or P6 instead */
  /* fscanf for "P5" or "P6";if not found, give appropriate message to stderr */
  /* then search for width and height, and 255, while ignoring lines */
  /* beginning with # which are comments.  flag error if no 255, indicating */
  /* support of other datatypes such as 65535 are forthcoming */
  
  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);
}

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

/* ndf: this part of the code loops through frames in a sequence
   and calls the subroutines that actually do the real work;
   future coding efficiency can be better by doing images ina batch
   processing, but needs to be of form pic000.ppm, pic001.ppm, pic002.ppm
   and the warp masks, etc., could be saved in memory.
   dont worry this for now, just single pictures.  later we'll
   support multiple pictgures (sequence of images) */

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. */