/* 
   ptrwarp.c: Pencigraphic Toolkit's Radial WARP about specified center
   ptwarp.c is part of pt-1.0
   see http://wearcam.org/pencigraphy
   or
   see http://www-white.media.mit.edu/~steve/pencigraphy

   ptrwarp      [-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
                    (defualt is bicubic interpolation)
   -m float       : how far down center of aberration is located
   -n float       : how far across center of aberration is located
   -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
                    (anisotropic axes)
   -o outfile     : Save the file after lines have been removed in
                    specified directory
   infile         : Name of original file

   other examples of file i/o
   -i3 pic.ppm    : specify pic000.ppm, pic001.ppm... (pic%03d.ppm)
   -i170:4:290 pic.ppm: specifies pic170.ppm, pic174.ppm,... pic290.ppm
   -i3 pic000.ppm : only numeric globbing (??? must be 3 digit NUMBER)
   -i3 pic???.ppm : 
   -i3 "170:4:290" pic.jpg ... need to determine which will work across
                   largest variety architectures
   (may simplify/robustify code and only support some of above)

   -o3 pic        : similar to above
  
   Lee Campbell,                     velocity map
   Sourabh Niyogi,                   bilinear
   John Y.A. Wang, Eero Simoncelli,  bicubic
   Stephen Intille,                  file i/o (code no longer used)
   Steve Mann,                       rewrote for HP, channels, time, opt.
   Nat Friedman                      pnm support, savewarp, loadwarp
*/

#include <math.h>
#include <stdio.h>
#include <ctype.h>
#include <fcntl.h>
#include <errno.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/dir.h>

/* Interpolation methods */
#define BILINEAR_INTERP 1
#define BICUBIC_INTERP  2

/* strdup() isn't sufficiently ubiquitous */
char *
my_strdup(char * s)
{
  char * dup_s;

  dup_s=(char *)malloc(sizeof(char)*(strlen(s)+1));
  strcpy(dup_s, s);
  return dup_s;
}

int
my_strcasecmp(char * s1, char * s2)
{
  while (1)
    {
      if (*s1!=*s2)
	{
	  if (tolower(*s1)!=tolower(*s2))
	    return *s1-*s2;
	}
      else if (*s1=='\0') break;
      s1++; s2++;
    }
  return 0;
}

void
usage(char *programname)
{
  fprintf(stderr, "use: %s [{inputfile,-i3 pic.ppm, -i pic???.ppm}] "
          "[{-o outputfile,-03 picw.ppm, -i pic???.ppm}] "
	  "-m center_m -n center_n -a float -b float -c float [-w] [-l] "
	  "[-p]\n",programname);
  exit(0); /* EXIT_SUCCESS = 0 (could also use) */
}

void
parse_command_line(int argc, char ** argv, int * num_input_files,
		   char *** input_filename, char ** output_filename,
		   int * use_approx_warp_map, int * overwrite_files,
		   int * interp_method, int * deinterl, int * load_warp,
		   int * make_warp, char ** warp_filename, int * proc_image,
		   float * x0, float * y0, float * a1, float * a2,
		   float * a3, int * M, int * N)
{
  int oflag=0;
  int i;
 
  if (argc<2 && isatty(0)) /* STDIN_FILENO = 0 (stdin file descriptor) */
    usage(argv[0]);


  for (i=1;i<argc;i++)
    {
      if (!my_strcasecmp(argv[i], "--savewarp"))
	{
	  if (*load_warp)
	    {
	      fprintf(stderr, "--savewarp and --loadwarp cannot be "
		      "used together!\n");
	      exit(0); /* EXIT_FAILURE */
	    }
	  if (argc<=(++i)) usage(argv[0]);
	  *warp_filename=(char *)my_strdup(argv[i]);
	  *make_warp=1;
	}
      else if (!my_strcasecmp(argv[i], "--loadwarp") ||
	       !my_strcasecmp(argv[i], "--loadmap"))
	{
	  if (*make_warp)
	    {
	      fprintf(stderr, "--savewarp and --loadwarp cannot be "
		      "used together!\n");
	      exit(0);
	    }
	  if (argc<=(++i)) usage(argv[0]);
	  *warp_filename=(char *)my_strdup(argv[i]);
	  *load_warp=1;
	}
      else if (!my_strcasecmp("--noimage", argv[i]))
	{
	  if (oflag) usage(argv[0]);
	  *proc_image=0;
	}
      else if (!my_strcasecmp("--N", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *N=atoi(argv[i]);
	}
      else if (!my_strcasecmp("--M", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *M=atoi(argv[i]);
	}
      else if (!my_strcasecmp("-postfix", argv[i]))
	{
	  /* -r sets the output filename postfix for batch processing */
	  if (!*proc_image)
	    usage(argv[0]);
	  if (argc<=(++i)) usage(argv[0]);
	  if (!oflag)
	    {
	      *output_filename=(char *)my_strdup(argv[i]);
	      oflag = 1;
	    }
	  else usage(argv[0]);
	}
      else if (!my_strcasecmp("-p", argv[i]))
	*use_approx_warp_map=1;
      else if (!my_strcasecmp("-w", argv[i]))
	*overwrite_files=1;
      else if (!my_strcasecmp("-l", argv[i]))
	*interp_method=BILINEAR_INTERP;
      else if (!my_strcasecmp("-x", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *x0=atof(argv[i]);
	}
      else if (!my_strcasecmp("-y", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *y0=atof(argv[i]);
	}
      else if (!my_strcasecmp("-a", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *a1=atof(argv[i]);
	}
      else if (!my_strcasecmp("-b", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *a2=atof(argv[i]);
	}
      else if (!my_strcasecmp("-c", argv[i]))
	{
	  if (argc<=(++i)) usage(argv[0]);
	  *a3=atof(argv[i]);
	}
      else if (!my_strcasecmp("-d", argv[i]))
	*deinterl=1;
      else if (!my_strcasecmp("-o", argv[i]))
	{
	  if (!*proc_image)
	    usage(argv[0]);
	  if (argc<=(++i)) usage(argv[0]);
	  if (!oflag)
	    {
	      *output_filename=(char *)my_strdup(argv[i]);
	      oflag = 1;
	    }
	  else usage(argv[0]);
	}
      else
	{
	  if (!*proc_image) usage(argv[0]);
	  if (*argv[i]=='-') usage(argv[0]);
	  if (!*num_input_files)
	    *input_filename=(char **)malloc(sizeof(char *)*2+argc-i);
	  (*input_filename)[*num_input_files]=(char *)my_strdup(argv[i]);
	  (*num_input_files)++;
	}
    }

  if ((!*x0 || !*y0 || !*a1 || !*a2 || !*a3) && !*load_warp) usage(argv[0]);

  if (!oflag && !isatty(1) && *proc_image) /* STDOUT_FILENO = 1
                                              stdout file descriptor */
    {
      output_filename=NULL;
      oflag=1;
    }

  if (!isatty(0) && !*num_input_files)
    {
      *input_filename=NULL;
      *num_input_files=1;
    }
}

void
split_image(unsigned char * image, unsigned char ** image_r,
	    unsigned char ** image_g, unsigned char ** image_b,
	    int N, int M)
{
  int i;

  *image_r=(unsigned char *)malloc(sizeof(char)*N*M);
  *image_g=(unsigned char *)malloc(sizeof(char)*N*M);
  *image_b=(unsigned char *)malloc(sizeof(char)*N*M);
  for (i=0;i<N*M;i++)
    {
      (*image_r)[i]=image[i*3];
      (*image_g)[i]=image[(i*3)+1];
      (*image_b)[i]=image[(i*3)+2];
    }
}

void
join_image(unsigned char ** image, unsigned char * image_r,
	   unsigned char * image_g, unsigned char * image_b,
	   int N, int M)
{
  int i;

  *image=(unsigned char *)malloc(sizeof(char)*N*M*3);
  for (i=0;i<N*M;i++)
    {
	(*image)[i*3]=image_r[i];
	(*image)[(i*3)+1]=image_g[i];
	(*image)[(i*3)+2]=image_b[i];
    }
}

void
save_warp_map(char * filename, float * warp_map_x, float * warp_map_y,
	      int N, int M, float a1, float a2, float a3,
	      float x, float y)
{
  FILE *map_file;
  int count;

  if ((map_file=fopen(filename, "w"))==NULL)
    {
      fprintf(stderr, "save_warp_map: Unable to open %s\n", filename);
      exit(0);
    }

  fprintf(map_file, "# This is a pnmwarp 1.0 warp map file with the following"
	  " parameters:\n# a1: %f a2: %f a3: %f x: %f y: %f\n",
	  a1, a2, a3, x, y);
  fprintf(map_file, "%d %d\n", N, M);
  count=fwrite(warp_map_x, sizeof(float), N*M, map_file);
  if (count!=N*M)
    {
      fprintf(stderr, "save_warp_map: oh dear me, we only wrote %d "
	      "warp map points and we wanted to write %d!\n",
	      count, N*M);
      exit(0);
    }
  count=fwrite(warp_map_y, sizeof(float), N*M, map_file);
  if (count!=N*M)
    {
      fprintf(stderr, "save_warp_map: oh dear me, we only wrote %d "
	      "warp map points and we wanted to write %d!\n",
	      count, N*M);
      exit(0);
    }
  fclose(map_file);
}

void
load_warp_map(char * filename, float ** warp_map_x, float ** warp_map_y,
	      int N, int M)
{
  FILE * map_file;
  char line[80];
  int map_N, map_M;
  int count;

  if ((map_file=fopen(filename, "r"))==NULL)
    {
      fprintf(stderr, "save_warp_map: Unable to open %s\n", filename);
      exit(0);
    }

  /* Read past the initial comment lines */
  fgets(line, sizeof(line), map_file);
  while (*line=='#') fgets(line, sizeof(line), map_file);

  /* Verify that the warp map is the same size as the image that we're
     working with */
  sscanf(line, "%d %d\n", &map_N, &map_M);

  if (map_N!=N || map_M!=M)
    {
      fprintf(stderr, "load_warp_map: warp map %s is not of the same "
	      "dimensions as the image.\n", filename);
      exit(0);
    }

  *warp_map_x=(float *)malloc(sizeof(float)*N*M);
  *warp_map_y=(float *)malloc(sizeof(float)*N*M);
  count=fread(*warp_map_x, sizeof(float), N*M, map_file);
  if (count!=N*M)
    {
      fprintf(stderr, "load_warp_map: oh dear me, we only read %d "
	      "warp map points and we wanted %d!\n",
	      count, N*M);
      exit(0);
    }
  count=fread(*warp_map_y, sizeof(float), N*M, map_file);
  if (count!=N*M)
    {
      fprintf(stderr, "load_warp_map: oh dear me, we only read %d "
	      "warp map points and we wanted %d!\n",
	      count, N*M);
      exit(0);
    }
  fclose(map_file);
}

void
check_output_file_name(char * filename)
{
  
  struct stat stbuf;
  
  if (stat(filename, &stbuf) != -1)
    {
      fprintf(stderr,"File %s exists.  Aborted.\n", filename);
      exit(0);
    }
}	  

void
get_image_params(FILE * ifs, int * bytes_pixel, int * N, int * M,
		 int * max_val)
{
  char input_line[25];
  int c;
  int i;

  if (fgetc(ifs)!='P')
    {
      fprintf(stderr, "get_image_params: file is not a pnm file - format "
	      "not recognized.\n");
      exit(0);
    }
  c=fgetc(ifs);
  if (c!='5' && c!='6')
    {
      fprintf(stderr, "get_image_params: I only speak P5 and P6.\n");
      exit(0);
    }

  if (c=='5') *bytes_pixel=1;
  if (c=='6') *bytes_pixel=3;
  *N=-1; /* assume image N not -1 */
  while (*N==-1) {
    while (isspace(c=fgetc(ifs))); /* trash whitespace */
    if (c=='#') while ((c=fgetc(ifs))!='\n'); /* trash comment line */
    else if (isdigit(c)) {
      i=0;
      input_line[i]=c;
      i++;
      while (isdigit(c=fgetc(ifs))) input_line[i++]=c;
      input_line[i]=0; /* null terminate the string */
      *N=atoi(input_line);
    }
  }
  *M=-1; /* assume image N not -1 */
  while (*M==-1) {
    while (isspace(c=fgetc(ifs))); /* trash whitespace */
    if (c=='#') while ((c=fgetc(ifs))!='\n'); /* trash comment line */
    else if (isdigit(c)) {
      i=0;
      input_line[i]=c;
      i++;
      while (isdigit(c=fgetc(ifs))) input_line[i++]=c;
      input_line[i]=0; /* null terminate the string */
      *M=atoi(input_line);
    }
  }
  *max_val=-1; /* assume image N not -1 */
  while (*max_val==-1) {
    while (isspace(c=fgetc(ifs))); /* trash whitespace */
    if (c=='#') while ((c=fgetc(ifs))!='\n'); /* trash comment line */
    else if (isdigit(c)) {
      i=0;
      input_line[i]=c;
      i++;
      while (isdigit(c=fgetc(ifs))) input_line[i++]=c;
      input_line[i]=0; /* null terminate the string */
      *max_val=atoi(input_line);
    }
  }

  /* FIXME: Support max_val > 255 */
  if (*max_val>255)
    {
      fprintf(stderr, "get_image_params: Sorry, we don't support maximum "
	      "colour component values greater than 255 yet.\n");
      exit(0);
    }
}

void
read_image(FILE * ifs, int N, int M, int bytes_pixel, int max_val,
	   unsigned char ** image)
{
  int bytes;
  int i;

  *image=(unsigned char *)malloc(sizeof(char)*bytes_pixel*N*M);

  bytes=fread(*image, bytes_pixel, N*M, ifs);
  bytes*=bytes_pixel;

  if (bytes!=N*M*bytes_pixel)
    {
      if (ferror(ifs))
	{
	  fprintf(stderr, "read_image: error reading image: %s\n",
		  strerror(errno));
	  exit(0);
	}
      else if (feof(ifs))
	{
	  fprintf(stderr, "read_image: reached end of input file "
		  " prematurely!  Read %d bytes; wanted %d.\n", bytes,
		  N*M*bytes_pixel);
	  exit(0);
	}
    }

  /* Rescale the image if necessary */
  if (max_val!=255)
    for(i=0;i<N*M*bytes_pixel;i++)
      (*image)[i]=((*image)[i]*255)/max_val;
}

void
write_image(char * ofile, int N, int M, int bytes_pixel,
	    unsigned char * vimage)
{
  int countout;
  FILE *ofs = stdout;

  if (ofile!=NULL) if ((ofs = fopen(ofile, "w")) == NULL)
    {
      fprintf(stderr, "pnmwarp: Cannot open output file\n");
      exit(0);
    }

  if (bytes_pixel==1) fprintf(ofs, "P5\n");
  if (bytes_pixel==3) fprintf(ofs, "P6\n");

  fprintf(ofs, "%d %d\n255\n", N, M);
  countout = fwrite(vimage, bytes_pixel, N*M, ofs);
  if (countout != N*M)
    fprintf(stderr,"write_image: wrote %d pixels; wanted to write %d.\n",
	    countout, N*M);
  if (ofile!=NULL) fclose(ofs);
}

/* This doesn't seem to work very well  -Nat */
void
make_warp_map_approx(float x0, float y0, float a1, float a2, float a3,
		     int N, int M, float ** warp_map_x,
		     float ** warp_map_y, int 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. */
{
  int x, y, cc;
  float xx, yy, r, r2, frac;

  *warp_map_x = (float *) malloc(sizeof(float)*N*M);
  *warp_map_y = (float *) malloc(sizeof(float)*N*M);

  cc = 0;
  for (y = 0; y < M; y++)
    for (x = 0; x < N; 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 */
	(*warp_map_x)[y * N + x] = xx*frac - xx;
	(*warp_map_y)[y * N + x] = yy*frac - yy;
      }
}

void
make_warp_map(float x0, float y0, float a1, float a2, float a3, int N,
	      int M, float ** warp_map_x, float ** warp_map_y,
	      int 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	*/
{
  int x, y, i, cc;
  float xx, yy, r, r2;
  float radius, rguess, rgold;
  float f, fprime;

  *warp_map_x = (float *) malloc(sizeof(float)*N*M);
  *warp_map_y = (float *) malloc(sizeof(float)*N*M);

  cc = 0;
  for (y = 0; y < M; y++)
    {
      rguess = 0.0;	/* force first guess approx */
      for (x = 0; x < N; 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);
	    }
	  i = (y * N + x);
	  (*warp_map_x)[i] = -(xx/radius)*rguess;
	  if (deinterl != 1)
	    (*warp_map_y)[i] = -(yy/radius)*rguess;
	  else
	    (*warp_map_y)[i] = -(yy/(radius*2))*rguess;
	  cc++;
	}
    }
}

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

  *out=(unsigned char *)malloc(sizeof(char)*N*M);

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

/* John/Eero's warper, for control expt. */
void
bicube_warp(unsigned char * im, unsigned char ** res, float * wx, float *wy,
	    int N, int M)
{
  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;

  *res=(unsigned char *)malloc(sizeof(char)*N*M);
  x_off = (float)((N - 1) / 2.0);
  y_off = (float)((M - 1) / 2.0);
  
  for(j=0; j<M; j++)  
    for(i=0; i<N; i++) { 
      xfrac = i - *wx++;
      xbase = (int)(xfrac);
      yfrac = j - *wy++;
      ybase = (int)(yfrac);
      
      if((((xbase) >= 0) && ((xbase + 0) < N)) &&
	 (((ybase) >= 0) && ((ybase + 0) < M))) {
	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 < N))?1:0;
	if (x_inbounds_1) x_coeff_1 = (frac2/2.0 - frac/3.0 - frac3/6.0);
	base++;
	x_inbounds0  = ((base >= 0) && (base < N))?1:0;
	if (x_inbounds0)  x_coeff0  = (1.0 - frac/2.0 - frac2 + frac3/2.0);
	base++;
	x_inbounds1  = ((base >= 0) && (base < N))?1:0;
	if (x_inbounds1)  x_coeff1 = (frac + frac2/2.0 - frac3/2.0);
	base++;
	x_inbounds2  = ((base >= 0) && (base < N))?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*N;
	
	if ((base >= 0) && (base < M)) {
	  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 += N;
	if ((base >= 0) && (base < M)) {
	  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 += N; 
	if ((base >= 0) && (base < M)) {
	  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 += N;
	if ((base >= 0) && (base < M)) {
	  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*N + i]; 
      }
      (*res)++;
    }
}

int
main(int argc, char ** argv)
{
  /* Some options... */
  int use_approx_warp_map=0, interp_method=BILINEAR_INTERP,
    deinterl=0, overwrite_files=0, save_warp=0, load_warp=0,
    proc_image=1;

  char * output_filename = NULL;
  char * output_postfix = NULL;
  /* input_filenames is the list of input files to process */
  char ** input_filenames = NULL;
  /* input_filename is the name of the file currently being processed */
  char * input_filename=NULL;
  char * map_filename = NULL;
  int num_input_files=0;

  /* Spherical abberation correction paramaters */
  float x0=0.0, y0=0.0, a1=0.0, a2=0.0, a3=0.0;  

  /* Buffers to hold the images */
  unsigned char *input_image, *output_image;
  unsigned char *input_image_r, *input_image_g, *input_image_b;
  unsigned char *output_image_r, *output_image_b, *output_image_g;

  /* Buffers to hold the warp map */
  float * warp_map_x, * warp_map_y;

  /* Image paramaters */
  int N, M, bytes_pixel, max_val;
  int first_N, first_M, first_bytes_pixel;

  /* The input file */
  FILE * ifs=stdin;

  int input_file_num;

  parse_command_line(argc, argv, &num_input_files, &input_filenames,
		     &output_filename, &use_approx_warp_map,
		     &overwrite_files, &interp_method, &deinterl,
		     &load_warp, &save_warp, &map_filename, &proc_image,
		     &x0, &y0, &a1, &a2, &a3, &N, &M);

  if (!proc_image)
    {
      if (use_approx_warp_map)
	make_warp_map_approx(x0, y0, a1, a2, a3, N,
			     M, &warp_map_x, &warp_map_y, deinterl);
      else
	make_warp_map(x0, y0, a1, a2, a3, N, M, &warp_map_x,
		      &warp_map_y, deinterl);

      if (save_warp)
	save_warp_map(map_filename, warp_map_x, warp_map_y, N, M,
		      a1, a2, a3, x0, y0);
      exit (0);
    }

  if (num_input_files>1)
    {
      output_postfix=(char *)my_strdup(output_filename);
      free(output_filename);
    }

  for (input_file_num=0;input_file_num<num_input_files;input_file_num++)
    {

      if (input_filenames==NULL) input_filename=NULL;
      else
	input_filename=input_filenames[input_file_num];

      fprintf(stderr, "Processing %s...\n", input_filename);

      if (num_input_files>1)
	{
	  if (input_file_num) free(output_filename);
	  output_filename=(char *)malloc(sizeof(char)*
					 (strlen(output_postfix)+
					  strlen(input_filename)+1));
	  strcpy(output_filename, input_filename);
	  strcat(output_filename, output_postfix);
	}

      if (!overwrite_files) check_output_file_name(output_filename);

      if (input_filename)
	if ((ifs=fopen(input_filename, "r"))==NULL)
	  {
	    fprintf(stderr, "Cannot open %s.\n", input_filename);
	    exit(0);
	  }
	  
      fprintf(stderr, "\tGetting input paramaters for %s...\n",
	      input_filename);
      get_image_params(ifs, &bytes_pixel, &N, &M, &max_val);
      fprintf(stderr,"\t%s: successfully read %d by %d image\n",argv[0],M,N);

      if (input_file_num && ((first_N!=N) || (first_M!=M) ||
	      (first_bytes_pixel!=bytes_pixel)))
	{
	  fprintf(stderr, "%s has different dimensions/format than "
		  "the first image!\n", input_filename);
	  exit(0);
	}

      /* Only generate/load the warp map once... */
      if (input_file_num==0)
	{
	  first_N=N; first_M=M;
	  first_bytes_pixel=bytes_pixel;
	  if (!load_warp)
	    fprintf(stderr, "\tGenerating %rwarp map...\n",
		    use_approx_warp_map ? "approximate " : "");
	  else fprintf(stderr, "\tLoading warp map %s...\n", map_filename);

	  if (load_warp)
	    load_warp_map(map_filename, &warp_map_x, &warp_map_y,
			  N, M);
	  else if (use_approx_warp_map)
	    make_warp_map_approx(x0, y0, a1, a2, a3, N,
				 M, &warp_map_x, &warp_map_y, deinterl);
	  else
	    make_warp_map(x0, y0, a1, a2, a3, N, M, &warp_map_x,
			  &warp_map_y, deinterl);

	  if (save_warp)
	    save_warp_map(map_filename, warp_map_x, warp_map_y, N, M,
			  a1, a2, a3, x0, y0);
	}

      fprintf(stderr, "\tReading %s...\n", input_filename);
      read_image(ifs, N, M, bytes_pixel, max_val, &input_image);

      fprintf(stderr, "\tUsing %s interpolation to warp %s...",
	      (interp_method==BILINEAR_INTERP) ? "bilinear" : "bicubic",
	      input_filename);
      fflush(stderr);
      if (bytes_pixel==3)
	{
	  split_image(input_image, &input_image_r, &input_image_g,
		      &input_image_b, N, M);
	  free(input_image);
	  if (interp_method==BILINEAR_INTERP)
	    {
	      fprintf(stderr, "red "); fflush(stderr);
	      bilinear_warp(input_image_r, &output_image_r, warp_map_x,
			    warp_map_y, N, M);
	      fprintf(stderr, "green "); fflush(stderr);
	      bilinear_warp(input_image_g, &output_image_g, warp_map_x,
			    warp_map_y, N, M);
	      fprintf(stderr, "blue\n"); fflush(stderr);
	      bilinear_warp(input_image_b, &output_image_b, warp_map_x,
			    warp_map_y, N, M);
	    }
	  else
	    {
	      fprintf(stderr, "red "); fflush(stderr);
	      bicube_warp(input_image_r, &output_image_r, warp_map_x,
			  warp_map_y, N, M);
	      fprintf(stderr, "green "); fflush(stderr);
	      bicube_warp(input_image_g, &output_image_g, warp_map_x,
			  warp_map_y, N, M);
	      fprintf(stderr, "blue\n"); fflush(stderr);
	      bicube_warp(input_image_b, &output_image_b, warp_map_x,
			  warp_map_y, N, M);
	    }
	  join_image(&output_image, output_image_r, output_image_g,
		     output_image_b, N, M);
	}
      else if (interp_method==BILINEAR_INTERP)
	bilinear_warp(input_image, &output_image, warp_map_x, warp_map_y,
		      N, M);
      else bicube_warp(input_image, &output_image, warp_map_x, warp_map_y,
		       N, M);

      fprintf(stderr, "\tWriting output image %s...\n",
	      output_filename);
      write_image(output_filename, N, M, bytes_pixel,
		  output_image);
    }

  exit(0);
}