/*
BABYL OPTIONS: -*- rmail -*-
Version: 5
Labels:
Note:   This is the header of an rmail file.
Note:   If you are seeing it in rmail,
Note:    it means the file has no messages in it.

1,,
From daemon Fri Jun 21 10:08:59 1996
Received: from ALEPH1.MIT.EDU by aleve.media.mit.edu; (5.65v3.2/1.1/06Jun95-8.2MPM)
	id AA14690; Fri, 21 Jun 1996 10:08:59 -0400
Received: (from ndf@localhost) by aleph1.mit.edu (8.6.12/8.6.12) id KAA17078; Fri, 21 Jun 1996 10:08:58 -0400
X-UIDL: 835366974.015
Date: Fri, 21 Jun 1996 10:08:58 -0400
Message-Id: <199606211408.KAA17078@aleph1.mit.edu>
From: Nat Friedman <ndf@aleph1.mit.edu>
To: steve@media.mit.edu
Cc: ndf@ALEPH1.MIT.EDU
Subject: pnmwarp 0.0

*** EOOH ***
From daemon Fri Jun 21 10:08:59 1996
X-UIDL: 835366974.015
Date: Fri, 21 Jun 1996 10:08:58 -0400
From: Nat Friedman <ndf@aleph1.mit.edu>
To: steve@media.mit.edu
Cc: ndf@ALEPH1.MIT.EDU
Subject: pnmwarp 0.0


Here's version 0.0 of pnmwarp.  It builds flawlessly and behaves identically
on every other machine I've tried it on.  Let me know how it works out for
you.  Here are a few sample command lines:

pnmwarp --noimage -x 100 -y 100 -a .9666 -b 9.393309e-5 -c 2.213069e-6 --width 256 --height 256 --savewarp warp.map

pnmwarp teapot.ppm -o out.ppm --loadwarp warp.map -w

pnmwarp -x 100 -y 100 -a .9666 -b 9.393309e-5 -c 2.213069e-6 < teapot.ppm > out.ppm

pnmwarp -x 100 -y 100 -a .9666 -b 9.393309e-5 -c 2.213069e-6 --savewarp warp.map < teapot.ppm > out.ppm


You get the picture....

Let me know how it does for you,
-Nat

*/

/*
 * pnmwarp.c - This program corrects spherical abberation in ppm and pgm files.
 *
 * The original authors are unknown to me.
 *
 * Converted to use pnm files, added support for multi-channel images and
 * batches of images, added support for saving and loading warp maps,
 * signifigantly rewritten by Nat Friedman, June, 1996.  <ndf@mit.edu>
 *
 */
#include <math.h>
#include <stdio.h>
#include <ctype.h>
#include <fcntl.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

/* These two externs are for getopt(3) */
extern char *optarg;
extern int optind;

#define OPTIONS "pwldx:y:a:b:c:o:"

#define BILINEAR_INTERP 1
#define BICUBIC_INTERP  2

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

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

int
my_strcasecmp(const char * s1, const 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(void)
{
  fprintf(stderr, "Usage: pnmwarp [inputfile] [-o outputfile] "
	  "-x center_x -y center_y -a float -b float -c float [-w] [-l] "
	  "[-p]\n");
  exit(EXIT_SUCCESS);
}

void
parse_command_line(int argc, char ** argv, 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 * width, int * height)
{
  int oflag=0, iflag=0;
  int i;


  if (argc<2 && isatty(STDIN_FILENO))
    usage();

  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(EXIT_FAILURE);
	    }
	  if (argc<=(++i)) usage();
	  *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(EXIT_FAILURE);
	    }
	  if (argc<=(++i)) usage();
	  *warp_filename=(char *)my_strdup(argv[i]);
	  *load_warp=1;
	}
      else if (!my_strcasecmp("--noimage", argv[i]))
	{
	  if (oflag) usage();
	  *proc_image=0;
	}
      else if (!my_strcasecmp("--width", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *width=atoi(argv[i]);
	}
      else if (!my_strcasecmp("--height", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *height=atoi(argv[i]);
	}
      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();
	  *x0=atof(argv[i]);
	}
      else if (!my_strcasecmp("-y", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *y0=atof(argv[i]);
	}
      else if (!my_strcasecmp("-a", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *a1=atof(argv[i]);
	}
      else if (!my_strcasecmp("-b", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *a2=atof(argv[i]);
	}
      else if (!my_strcasecmp("-c", argv[i]))
	{
	  if (argc<=(++i)) usage();
	  *a3=atof(argv[i]);
	}
      else if (!my_strcasecmp("-d", argv[i]))
	*deinterl=1;
      else if (!my_strcasecmp("-o", argv[i]))
	{
	  if (!*proc_image)
	    usage();
	  if (argc<=(++i)) usage();
	  if (!oflag)
	    {
	      *output_filename=(char *)my_strdup(argv[i]);
	      oflag = 1;
	    }
	  else usage();
	}
      else
	{
	  if (*argv[i]=='-') usage();
	  if (!iflag) *input_filename=(char *)my_strdup(argv[i]);
	  else usage();
	}
    }

  if ((!x0 || !y0 || !a1 || !a2 || !a3) && !*load_warp) usage();

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

  if (!isatty(STDIN_FILENO) && !iflag) *input_filename=NULL;
}

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

  *image_r=(unsigned char *)malloc(sizeof(char)*width*height);
  *image_g=(unsigned char *)malloc(sizeof(char)*width*height);
  *image_b=(unsigned char *)malloc(sizeof(char)*width*height);
  for (i=0;i<width*height;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 width, int height)
{
  int i;

  *image=(unsigned char *)malloc(sizeof(char)*width*height*3);
  for (i=0;i<width*height;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 width, int height, 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(EXIT_FAILURE);
    }

  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", width, height);
  count=fwrite(warp_map_x, sizeof(float), width*height, map_file);
  if (count!=width*height)
    {
      fprintf(stderr, "save_warp_map: oh dear me, we only wrote %d "
	      "warp map points and we wanted to write %d!\n",
	      count, width*height);
      exit(EXIT_FAILURE);
    }
  count=fwrite(warp_map_y, sizeof(float), width*height, map_file);
  if (count!=width*height)
    {
      fprintf(stderr, "save_warp_map: oh dear me, we only wrote %d "
	      "warp map points and we wanted to write %d!\n",
	      count, width*height);
      exit(EXIT_FAILURE);
    }
  fclose(map_file);
}

void
load_warp_map(char * filename, float ** warp_map_x, float ** warp_map_y,
	      int width, int height)
{
  FILE * map_file;
  char line[80];
  int map_width, map_height;
  int count;

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

  /* 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_width, &map_height);

  if (map_width!=width || map_height!=height)
    {
      fprintf(stderr, "load_warp_map: warp map %s is not of the same "
	      "dimensions as the image.\n", filename);
      exit(EXIT_FAILURE);
    }

  *warp_map_x=(float *)malloc(sizeof(float)*width*height);
  *warp_map_y=(float *)malloc(sizeof(float)*width*height);
  count=fread(*warp_map_x, sizeof(float), width*height, map_file);
  if (count!=width*height)
    {
      fprintf(stderr, "load_warp_map: oh dear me, we only read %d "
	      "warp map points and we wanted %d!\n",
	      count, width*height);
      exit(EXIT_FAILURE);
    }
  count=fread(*warp_map_y, sizeof(float), width*height, map_file);
  if (count!=width*height)
    {
      fprintf(stderr, "load_warp_map: oh dear me, we only read %d "
	      "warp map points and we wanted %d!\n",
	      count, width*height);
      exit(EXIT_FAILURE);
    }
  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(EXIT_FAILURE);
    }
}	  

void
get_image_params(FILE * ifs, int * bytes_pixel, int * width, int * height,
		 int * max_val)
{
  int c;

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

  /* FIXME: Support comments.  */
  fscanf(ifs, "%d", width);
  fscanf(ifs, "%d", height);
  fscanf(ifs, "%d\n", max_val);

  /* 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(EXIT_FAILURE);
    }

  if (c=='5') *bytes_pixel=1;
  if (c=='6') *bytes_pixel=3;
}

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

  *image=(unsigned char *)malloc(sizeof(char)*bytes_pixel*width*height);

  bytes=fread(*image, bytes_pixel, width*height, ifs);
  bytes*=bytes_pixel;

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

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

void
write_image(char * ofile, int width, int height, 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(EXIT_FAILURE);
    }

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

  fprintf(ofs, "%d %d\n255\n", width, height);
  countout = fwrite(vimage, bytes_pixel, width*height, ofs);
  if (countout != width*height)
    fprintf(stderr,"write_image: wrote %d pixels; wanted to write %d.\n",
	    countout, width*height);
  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 width, int height, 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)*width*height);
  *warp_map_y = (float *) malloc(sizeof(float)*width*height);

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

void
make_warp_map(float x0, float y0, float a1, float a2, float a3, int width,
	      int height, 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)*width*height);
  *warp_map_y = (float *) malloc(sizeof(float)*width*height);

  cc = 0;
  for (y = 0; y < height; y++)
    {
      rguess = 0.0;	/* force first guess approx */
      for (x = 0; x < width; 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 * width + 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 width, int height)
{
  int x, y;
  float xpf, ypf;
  int xpi, ypi;
  int xbound, ybound;
  float dx, dy, ans;

  *out=(unsigned char *)malloc(sizeof(char)*width*height);

  xbound = width-1;
  ybound = height-1;  
  for (y=0; y<height; y++) 
    for (x=0; x<width; x++) {
      xpf = (float)x-*(vx+y*width+x);
      ypf = (float)y-*(vy+y*width+x);
      xpi = (int)floor(xpf);
      ypi = (int)floor(ypf);
      if ( ( xpi < 0 ) || ( ypi < 0 ) || 
	   ( (xpi+1) > xbound) ||  ( (ypi+1) > ybound) ) 
	*((*out)+y*width+x) = *(im+y*width+x);
      else {
	dx  = xpf-(float)xpi;
	dy  = ypf-(float)ypi;
	ans =	(1-dx) * (1-dy) * (*(im + ypi   * width + xpi)) + 
	  dx   * (1-dy) * (*(im + ypi   * width + xpi + 1)) +
	  (1-dx) *   dy   * (*(im+(ypi+1) * width + xpi)) + 
	  dx   *   dy   * (*(im+(ypi+1) * width + xpi + 1));
	if (ans > 255) ans = 255;
	if (ans < 0) ans = 0;
	*((*out)+y*width+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 width, int height)
{
  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)*width*height);
  x_off = (float)((width - 1) / 2.0);
  y_off = (float)((height - 1) / 2.0);
  
  for(j=0; j<height; j++)  
    for(i=0; i<width; i++) { 
      xfrac = i - *wx++;
      xbase = (int)(xfrac);
      yfrac = j - *wy++;
      ybase = (int)(yfrac);
      
      if((((xbase) >= 0) && ((xbase + 0) < width)) &&
	 (((ybase) >= 0) && ((ybase + 0) < height))) {
	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 < width))?1:0;
	if (x_inbounds_1) x_coeff_1 = (frac2/2.0 - frac/3.0 - frac3/6.0);
	base++;
	x_inbounds0  = ((base >= 0) && (base < width))?1:0;
	if (x_inbounds0)  x_coeff0  = (1.0 - frac/2.0 - frac2 + frac3/2.0);
	base++;
	x_inbounds1  = ((base >= 0) && (base < width))?1:0;
	if (x_inbounds1)  x_coeff1 = (frac + frac2/2.0 - frac3/2.0);
	base++;
	x_inbounds2  = ((base >= 0) && (base < width))?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*width;
	
	if ((base >= 0) && (base < height)) {
	  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 += width;
	if ((base >= 0) && (base < height)) {
	  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 += width; 
	if ((base >= 0) && (base < height)) {
	  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 += width;
	if ((base >= 0) && (base < height)) {
	  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*width + 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 *input_filename = NULL;
  char *map_filename = NULL;

  /* 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 width, height, bytes_pixel, max_val;

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

  parse_command_line(argc, argv, &input_filename, &output_filename,
		     &use_approx_warp_map, &overwrite_files,
		     &interp_method, &deinterl, &load_warp, &save_warp,
		     &map_filename, &proc_image, &x0, &y0, &a1, &a2, &a3,
		     &width, &height);

  if (!proc_image) goto make_warp;

  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(EXIT_FAILURE);
      }

  fprintf(stderr, "Getting input paramaters...\n");
  get_image_params(ifs, &bytes_pixel, &width, &height, &max_val);

  fprintf(stderr, "Reading image...\n");
  read_image(ifs, width, height, bytes_pixel, max_val, &input_image);

  if (!load_warp)
    fprintf(stderr, "Generating %swarp map...\n",
	    use_approx_warp_map ? "approximate " : "");
  else fprintf(stderr, "Loading warp map %s...\n", map_filename);

make_warp:
  if (load_warp)
    load_warp_map(map_filename, &warp_map_x, &warp_map_y, width, height);
  else if (use_approx_warp_map)
    make_warp_map_approx(x0, y0, a1, a2, a3, width, height, &warp_map_x,
			 &warp_map_y, deinterl);
  else
    make_warp_map(x0, y0, a1, a2, a3, width, height, &warp_map_x,
		  &warp_map_y, deinterl);

  if (save_warp)
      save_warp_map(map_filename, warp_map_x, warp_map_y, width, height,
		    a1, a2, a3, x0, y0);
  if (!proc_image) exit(EXIT_SUCCESS);

  fprintf(stderr, "Using %s interpolation to warp image...",
	  (interp_method==BILINEAR_INTERP) ? "bilinear" : "bicubic");
  fflush(stderr);
  if (bytes_pixel==3)
    {
      split_image(input_image, &input_image_r, &input_image_g,
		  &input_image_b, width, height);
      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, width, height);
	  fprintf(stderr, "green "); fflush(stderr);
	  bilinear_warp(input_image_g, &output_image_g, warp_map_x,
			warp_map_y, width, height);
	  fprintf(stderr, "blue\n"); fflush(stderr);
	  bilinear_warp(input_image_b, &output_image_b, warp_map_x,
			warp_map_y, width, height);
	}
      else
	{
	  fprintf(stderr, "red "); fflush(stderr);
	  bicube_warp(input_image_r, &output_image_r, warp_map_x,
			warp_map_y, width, height);
	  fprintf(stderr, "green "); fflush(stderr);
	  bicube_warp(input_image_g, &output_image_g, warp_map_x,
			warp_map_y, width, height);
	  fprintf(stderr, "blue\n"); fflush(stderr);
	  bicube_warp(input_image_b, &output_image_b, warp_map_x,
			warp_map_y, width, height);
	}
      join_image(&output_image, output_image_r, output_image_g,
		 output_image_b, width, height);
    }
  else if (interp_method==BILINEAR_INTERP)
    bilinear_warp(input_image, &output_image, warp_map_x, warp_map_y,
		  width, height);
  else bicube_warp(input_image, &output_image, warp_map_x, warp_map_y,
		   width, height);

  fprintf(stderr, "Writing output image...\n");
  write_image(output_filename, width, height, bytes_pixel, output_image);

  exit(EXIT_SUCCESS);
}