/* Last edited on 2007-09-07 19:55:36 by stolfi */

/* ====================================================================== */

image_vec_t *cons_image(char *fname, double ccol, double crow, image_vec_t *rest);
  /* Prepends image filename {fname} with nominal center position
    coordinates {ccol,crow} to the list {rest}.  */

int imin(int x, int y);
int imax(int x, int y);
  
void compute_optimal_shifts
  ( int ni, 
    image_vec_t *iml, 
    double dmax, 
    double radius, 
    double *dcol, 
    double *drow
  );

image_vec_t *reduce_images(image_vec_t *iml);
  /* Reduces all imagaes in {iml} by a linear factor of 2;
    returns a list of the results.  Also computes the coordinates
    of the centers of the reduced images. */

void average_images(int ni, image_vec_t *iml, int iex, float_image_t *avg, r2_t d[]);
  /* Computes average of all images, excluding image {iex}, shifting each one 
    so that its center displaced by {d[j]} matches the center of {avg}.  */
    
interval_t compute_mismatch
  ( float_image_t *tim, 
    r2_t d, 
    float_image_t *sim, 
    double scale, 
    int nwt, 
    double *gwt
  );
  /* Computes the mismatch between the image {tim}, shifted by {(dcol,drow)},
    against the square image {sim}.  The vector {wt} is a unidimensional weight
    distribution.  Assumes that the length of {wt} is the same as the number
    of rows and columns of {sim}.  

    The mismatch is the sum of the squared pixel differences, each multiplied
    by the {scale} factor and the bidimensional weight distribution {wt[col]*wt[row]}.  */

void recompute_pixel_variance(int ni, image_vec_t *iml, alignment_t *al, float_image_t *avg, int nwt, double *wt);
  /* Recomputes the mismatches {al->mism[i]} and the total pixel variance {al->tvar},
    from scratch. */
    
int max_chns(image_vec_t *iml);
  /* Returns the maximum {chns} among all images in {iml}.  */

double *gaussian_distr(int width, double sigma);
  /* Gaussian distribution {g[i] = A*exp(-((i+0.5)-m)^2/sigma^2)}
    where {m = width/2} and {A} is such that sum is 1.  */
    
void write_optimum_alignments(int ni, image_vec_t *iml, double *dcol, double *drow);
void debug_displacement(char *head, int i, double col, double row, interval_t s, char *tail);

void eval_at_point(float_image_t *im, double col, double row, interval_t *fv);
void accum_pixel(int chns, interval_t *el, double wt, interval_t *fv, double *twtP);
void make_pixel_indet(int chns, interval_t *fv);
void scale_pixel(int chns, double s, interval_t *fv);
void subtract_pixel(int chns, interval_t *av, interval_t *fv);
interval_t pixel_norm(int chns, interval_t *fv);
void convert_pixel(int old_chns, interval_t *fv, int new_chns);
interval_t interval_join(interval_t *u, interval_t *v);

void debug_floatized_pixel(char *head, int col, int row, int chns, interval_t *fv, char *tail);
void debug_quantized_pixel(char *head, int col, int row, int chns, int *fv, char *tail);
void debug_point(char *head, double w, double x, double y, char *tail);
void debug_image(char *head, char *indent, float_image_t *im);


image_vec_t *reduce_images(image_vec_t *iml)
  {
    fprintf(stderr, "r");
    if (iml == NULL)
      { return NULL; }
    else
      { float_image_t *im = iml->im;
        double rcol = iml->ccol;
        double rrow = iml->crow;
        float_image_t *rd = reduce_image(im, &rcol, &rrow);
        return cons_image(rd, rcol, rrow, reduce_images(iml->next));
      }
  }

double *gaussian_distr(int width, double sigma)
  { double m = ((double)width)/2.0;
    double *g = (double *)malloc(width*sizeof(double));
    double s = 0;
    int i;
    if (g == NULL) { demand(FALSE, "out of memory for gaussian weight"); }
    for (i=0; i < width; i++)
      { double x = (i + 0.5 - m)/sigma;
        double y = exp(-x*x); 
        g[i] = y; s += y;
      }
    for (i=0; i < width; i++)
      { g[i] /= s; }
    return g;
  }

interval_t compute_mismatch
  ( float_image_t *tim, 
    double dcol, 
    double drow, 
    float_image_t *sim, 
    double scale, 
    int nwt, 
    double *gwt
  )
  {
    int srow, scol;
    double wt2;
    interval_t d2;
    interval_t td2 = (interval_t){0.0, 0.0}; /* Total mismatch. */
    if ((sim->rows != nwt) || (sim->cols != nwt)) 
      { demand(FALSE, "ref image has wrong size"); } 
    for (srow = 0; srow < sim->rows; srow++)
      { for (scol = 0; scol < sim->cols; scol++)
          { int sindex = (srow * sim->cols + scol)*sim->chns;
            interval_t *sv = &(sim->el[sindex]);
            interval_t tv[MAX_CHANNELS];
            double tcol = scol + ((double)(tim->cols - sim->cols))/2.0 - dcol;
            double trow = srow + ((double)(tim->rows - sim->rows))/2.0 - drow;

            /* debug = ((scol % 56) == 23) & ((srow % 56) == 23); */
            eval_at_point(tim, tcol, trow, &(tv[0]));
            if (debug) 
              { debug_floatized_pixel("imi", tcol, trow, sim->chns, &(tv[0]), "\n"); } 
            if (tim->chns != sim->chns)
              { convert_pixel(tim->chns, &(tv[0]), sim->chns); }
            subtract_pixel(sim->chns, &(sv[0]), &(tv[0]));
            if (debug) 
              { debug_floatized_pixel("dif", scol, srow, sim->chns, &(tv[0]), "\n\n"); } 
            d2 = pixel_norm(sim->chns, &(tv[0]));
            wt2 = scale*wt[srow]*wt[scol];
            td2.end[0] += wt2*d2.end[0];
            td2.end[1] += wt2*d2.end[1];
          }
      }
    return td2;
  }

void average_images(int ni, image_vec_t *iml, int iex, float_image_t *avg, double *dcol, double *drow)
  {
    int row, col;
    int navg = ((iex >= 0) && (iex < ni) ? ni-1 : ni);
    double qt = 1.0/((double)navg);
    double tqt;
    for (row = 0; row < avg->rows; row++)
      { for (col = 0; col < avg->cols; col++)
          { int aindex = (row * avg->cols + col)*avg->chns;
            interval_t *fva = &(avg->el[aindex]);
            int j, k;
            image_vec_t *impj;
            for (k = 0; k < avg->chns; k++) { fva[k] = (interval_t){ 0.0, 0.0 }; }
            for (j = 0, impj = iml; j < ni; j++, impj = impj->next)
              { if (j != iex) 
                  { float_image_t *imj = impj->im;
                    interval_t fvj[MAX_CHANNELS];
                    double colj = col + ((double)(imj->cols - avg->cols))/2.0 - dcol[j];
                    double rowj = row + ((double)(imj->rows - avg->rows))/2.0 - drow[j];
                    eval_at_point(imj, colj, rowj, &(fvj[0]));
                    if (imj->chns != avg->chns)
                      { convert_pixel(imj->chns, &(fvj[0]), avg->chns); }
                    accum_pixel(avg->chns, &(fvj[0]), qt, fva, &tqt); 
                  }
              }
          }
      }
  }

int max_chns(image_vec_t *iml)
  {
    int chns = 0;
    while (iml != NULL)
      { float_image_t *im = iml->im;
        if (im->chns > chns) { chns = im->chns; }
        iml = iml->next;
      }
    return chns;
  }

image_vec_t *cons_image(float_image_t *im, double ccol, double crow, image_vec_t *rest)
  {
    image_vec_t *new = (image_t*)malloc(sizeof(image_t));

    if (new == NULL) { demand(FALSE, "out of memory for image_t"); }
    new->im = im;
    new->ccol = ccol;
    new->crow = crow;
    new->next = rest;
    return new;
  }

void debug_displacement(char *head, int i, double col, double row, interval_t s, char *tail)
  {
    fprintf(stderr, "im[%d] d=[%8.3f,%8.3f] s = [%15.6f _ %15.6f]", head, i, row, col, s.end[0], s.end[1]);
    fprintf(stderr, "%s", tail);
  }

options_t *get_options(int argc, char **argv)
  {
    argparser_t pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);

    options_t *o = (options_t *)notnull(malloc(sizeof(options_t)), "no mem");

    if (argparser_keyword_present(pp, "-radius"))
      { o->radius = argparser_get_next_double(pp, 0.0, 1.0e+6); }
    else
      { o->radius = 10.0; }
      
    if (argparser_keyword_present(pp, "-dmax"))
      { o->dmax = argparser_get_next_double(pp, 1.0e-6, 1.0e+6); }
    else
      { o->dmax = 15.0; }
      
    if (argparser_keyword_present(pp, "-epsilon"))
      { o->epsilon = argparser_get_next_double(pp, 1.0e-6, 1.0e+6);  }
    else
      { o->epsilon = 0.125; }
      
    o->hleft = FALSE;
    imgc_parse_haxis(pp, &(o->hleft));
    
    o->vdown = TRUE;
    imgc_parse_vaxis(pp, &(o->vdown));
    
    o->icenter = FALSE;
    o->iorg.c[0] = 0.0;
    o->iorg.c[1] = 0.0;
    imgc_parse_icenter_iorg(pp, &(o->icenter), &(o->iorg));
    
    argparser_get_keyword(pp, "-size");
    o->cols = argparser_get_next_int(pp, 1, MAX_SIZE);
    o->rows = argparser_get_next_int(pp, 1, MAX_SIZE);
      
    o->verbose = argparser_keyword_present(pp, "-verbose");
      
    argparser_skip_parsed(pp);

    /* Parse file arguments */
    o->ni = 0;
    o->iml = NULL;
    while (argparser_next(pp) != NULL)
      { double ccol, crow;
        char *rest; 
        char *f_name = argparser_get_next(pp);
        double ccol = argparser_get_next_double(pp, -1.0e+6, +1.0e+6);
        double crow = argparser_get_next_double(pp, -1.0e+6, +1.0e+6);
        iml = cons_image(f_name, ccol, crow, iml);
        o->ni++;
      }

    argparser_finish(pp);

    return o;
  }

/* PROTOTYPES */

int rounding_bias(int n);

#define HRADIUS 1  /* Effective radius of reconstruction kernel. */

void eval_at_point(float_image_t *im, double col, double row, Interval *fv)
/* Temporary hack: use bilinear interpolation. */
{
  if (debug) { fprintf(stderr, "evaluating im[%9.4f][%9.4f]\n", row, col); }
  if ( (col < -HRADIUS) 
    || (col > ((double)im->cols)+HRADIUS)  
    || (row < -HRADIUS) 
    || (row > ((double)im->rows)+HRADIUS) )
    { make_pixel_indet(im->chns, fv); return; }
  else
    { 
      int icol, irow, k;
      double twt = 0.0;
      for (k = 0; k < im->chns; k++) { fv[k].end[0] = 0; fv[k].end[1] = 0; }
      for (irow = floor(row); irow <= (int)ceil(row); irow++)
        { if ((irow >= 0) && (irow < im->rows))
	    { int rindex = (irow * im->cols)*im->chns;
              double sr = row - (double)irow;
              for (icol = (int)floor(col); icol <= (int)ceil(col); icol++)
                { if ((icol >= 0) && (icol < im->cols))
		    { double sc = col - (double)icol;
                      double wt = (1 - fabs(sr))*(1 - fabs(sc));
                      Interval *el = &(im->el[rindex + icol*im->chns]);
                      if (debug) 
                         { fprintf(stderr, "wt = %8.4f\n", wt);
                           debug_floatized_pixel("im", icol, irow, im->chns, el, "\n");
                         }
                      accum_pixel(im->chns, &(el[0]), wt, &(fv[0]), &twt);
                    }
                }
	    }
	}
      if (debug) 
        { fprintf(stderr, "twt = %8.4f\n", twt);
          debug_floatized_pixel("ac", col, row, im->chns, fv, "\n");
        }
      /* Normalize to unit sum of weights to account for missing pixels */
      if (twt != 1.0) 
        { if (twt == 0)
	    { make_pixel_indet(im->chns, fv); }
	  else
            { scale_pixel(im->chns, 1.0/twt, fv); }
        } 
    }
}

void convert_pixel(int old_chns, Interval *fv, int new_chns)
/*
  Converts `fv' from greyscale pixel to an RGB grey color,
  or from RGB color to grayscale.
*/
{
  int k;
  if ((old_chns == 1) && (new_chns > 1))
    { for (k = 1; k < new_chns; k++) fv[k] = fv[0]; }
  else if ((old_chns == 3) && (new_chns == 1))
    { double end[0] = 0.299+fv[0].end[0] + 0.587*fv[1].end[0] + 0.114*fv[2].end[0];
      double end[1] = 0.299+fv[0].end[1] + 0.587*fv[1].end[1] + 0.114*fv[2].end[1];
      fv[0].end[0] = end[0]; fv[0].end[1] = end[1];
    }
  else if (old_chns != new_chns)
    { pm_error("cannot promote pixels to common type"); }
}

void accum_pixel(int chns, Interval *el, double wt, Interval *fv, double *twtP)
/*
  Adds `*el' scaled by `wt' onto `*fv', adds `wt' to `*twtP'. 
*/
{
  if (wt != 0.0)
    { int k;
      for (k = 0; k < chns; k++)
        { Interval *ek = &(el[k]); 
          Interval *vk = &(fv[k]);
          if (wt > 0.0)
            { vk->end[0] += wt*ek->end[0]; vk->end[1] += wt*ek->end[1]; }
          else 
            { vk->end[0] += wt*ek->end[1]; vk->end[1] += wt*ek->end[0]; }
        }
      (*twtP) += wt;
    }
}

void make_pixel_indet(int chns, Interval *fv)
{
  int k;
  for (k = 0; k < chns; k++) { fv[k].end[0] = 0.0; fv[k].end[1] = 1.0; }
}

void scale_pixel(int chns, double s, Interval *fv)
{
  int k;
  for (k = 0; k < chns; k++) 
    { Interval *vk = &(fv[k]);
      if (s >= 0) 
        { vk->end[0] *= s; vk->end[1] *= s; }
      else
        { double t = vk->end[0]; vk->end[1] = s*vk->end[0]; vk->end[0] = s*t; }
    }
}

void subtract_pixel(int chns, Interval *av, Interval *fv)
/*
  Subtracts pixel `av' from pixel `fv', leaves result in `fv'.
*/
{
  int k;
  for (k = 0; k < chns; k++)
    { fv[k].end[0] -= av[k].end[1]; fv[k].end[1] -= av[k].end[0]; }
}

Interval pixel_norm(int chns, Interval *fv)
/*
  Returns the sumof the squares of the components of `fv'.
*/
{
  int k;
  Interval nrm;
  nrm.end[0] = 0.0; nrm.end[1] = 0.0;
  for (k = 0; k < chns; k++)
    { double end[0] = fv[k].end[0];
      double end[1] = fv[k].end[1];
      if (end[0] >= 0) 
        { nrm.end[0] += end[0]*end[0]; nrm.end[1] += end[1]*end[1]; }
      else if (end[1] <= 0)
        { nrm.end[0] += end[1]*end[1]; nrm.end[1] += end[0]*end[0]; }
      else 
        { nrm.end[1] += (-end[0] > end[1] ? end[0]*end[0] : end[1]*end[1]); }
    }
  return nrm;
}

Interval interval_join(Interval a, Interval b)
{
  Interval res;
  res.end[0] = (a.end[0] < b.end[0] ? a.end[0] : b.end[0]);
  res.end[1] = (a.end[1] > b.end[1] ? a.end[1] : b.end[1]);
  return res;
}

void debug_floatized_pixel(char *head, int col, int row, int chns, Interval *fv, char *tail)
{ 
  int k;
  fprintf(stderr, "%s[%3d][%3d] = (", head, row, col);
  for (k = 0; k < chns; k++) 
    { fprintf(stderr, " [%7.4f _ %7.4f]", fv[k].end[0], fv[k].end[1]); }
  fprintf(stderr, " )%s", tail);
}

void debug_quantized_pixel(char *head, int col, int row, int chns, int *fv, char *tail)
{ 
  int k;
  fprintf(stderr, "%s[%3d][%3d] = (", head, row, col);
  for (k = 0; k < chns; k++) 
    { fprintf(stderr, " %3d", fv[k]); }
  fprintf(stderr, " )%s", tail);
}

void debug_point(char *head, double w, double x, double y, char *tail)
{
  fprintf(stderr, "%s = [%9.5f %9.5f %9.5f]", head, w, x, y);
  if (w != 0.0)
    {  fprintf(stderr, " = (%12.5f %12.5f)", x/w, y/w); }
  fprintf(stderr, "%s , tail);
}

void debug_image(char *head, char *indent, float_image_t *im)
{
  int row, col;
  fprintf(stderr, "%sbegin image %s (%d × %d × %d)\n",
    indent, head, im->rows, im->cols, im->chns
  );
  for (row = 0; row < im->rows; ++row)
    { for (col = 0;  col < im->cols; ++col)
        { Interval *el = &(im->el[(row*im->cols + col)*im->chns]);
          debug_floatized_pixel(indent, col, row, im->chns, el, "\n");
        }
    }
  fprintf(stderr, "%send image\n", indent);
}


/*
** AUTHORSHIP, RIGHTS AND CONDITIONS NOTICE.  This program was written
** in aug/2002 by Jorge Stolfi at the University of Campinas, Brazil.
** 
** Permission to use, copy, modify, and redistribute this software and
** its documentation for any purpose and without fee is hereby
** granted, provided that: (1) the copyright notice at the top of this
** file and this authorship, rights and conditions notice is retained
** in all derived source files and documentation; (2) no executable
** code derived from this file is published or distributed without the
** corresponding source code; and (3) these same rights are granted to
** any recipient of such code, under the same conditions.  This
** software is provided "as is" without express or implied warranty.
** Neither the author nor his employer are liable for any damages or
** losses that may result from its use. END OF NOTICE.
**
*/

/* COPYRIGHT, AUTHORSHIP, AND WARRANTY NOTICE:
**
** Permission to use, copy, modify, and redistribute this software and
** its documentation for any purpose and without fee is hereby
** granted, provided that: (1) the copyright notice at the top of this
** file and this copyright, authorship, and warranty notice is retained
** in all derived source files and documentation; (2) no executable
** code derived from this file is published or distributed without the
** corresponding source code; and (3) these same rights are granted to
** any recipient of such code, under the same conditions.
** This software is provided "as is", WITHOUT ANY EXPLICIT OR IMPLICIT
** WARRANTIES, not even the implied warranties of merchantibility and
** fitness for a particular purpose. END OF NOTICE.
*/

/* ====================================================================== */

#define HELP \
  "  -precision {XPREC} {YPREC} \\\n"
  
#define HELP_INFO \
  "  -precision {XPREC} {YPREC}\n" \
  "    This required argument specifies the desired accuracy" \
  " of the adjustments.\n" \
  "\n" \

typedef struct foo_t
  {
    r2_t precision;     /* Desired precision of alignments. */
  } foo_t;

void read_pnm_file_header(FILE *f, PNMHeader *hd);
float_image_t *allocate_input_image(PNMHeader *hd);
void read_pnm_file_pixels(FILE *f, PNMHeader *hd, double zero, double unit, float_image_t *im);
Interval floatize(xelval ival, xelval maxval, double zero, double scale);
float_image_t *read_image(char *name);
float_image_t* allocate_image(int cols, int rows, int chns);

float_image_t *read_image(char *name)
/* 
  Reads an input image from the named pgm/ppm/pnm file.
*/
{
  double zero, unit;
  float_image_t *im;
  FILE *f;
  PNMHeader hd;  /* Specs from input file header. */
  
  pm_message("reading %s...", name);
  f = pm_openr(name); 
  read_pnm_file_header(f, &hd);
  zero = 0;
  unit = (double)hd.maxval;
    
  /* Allocate input image and read input pixels: */
  im = allocate_input_image(&hd);
  read_pnm_file_pixels(f, &hd, zero, unit, im);
  return im;
}

void read_pnm_file_header(FILE *f, PNMHeader *hd)
{
  /* Read header: */
  pnm_readpnminit(f, &(hd->cols), &(hd->rows), &(hd->maxval), &(hd->format));
  pm_message("input format = %c%c maxval = %d cols = %d rows = %d", 
    (hd->format >> 8)&255, hd->format&255, hd->maxval, hd->rows, hd->cols);
}

float_image_t *allocate_input_image(PNMHeader *hd)
{
  int chns = 0;
  /* Decide number of channels: */
  switch (hd->format)
    { case PPM_FORMAT:
      case RPPM_FORMAT:
        chns = 3;
      break;
      case PGM_FORMAT:
      case RPGM_FORMAT:
      case PBM_FORMAT:
      case RPBM_FORMAT:
        chns = 1;
      break;
    default:
      pm_error("bad input file format");
    }
  return allocate_image(hd->cols, hd->rows, chns);
} 

float_image_t* allocate_image(int cols, int rows, int chns)
{
  float_image_t *im = (float_image_t*)malloc(sizeof(float_image_t));
  if (im == NULL) { pm_error("out of memory for image descriptor"); }
  im->cols = cols; im->rows = rows; im->chns = chns;
  if (MAX_PIXELS / cols < rows) { pm_error("image too big"); }
  im->el = (Interval *)malloc(rows*cols*chns*sizeof(Interval));
  if (im->el == NULL) { pm_error("out of memory for image pixels"); }
  return im;
}

void read_pnm_file_pixels(FILE *f, PNMHeader *hd, double zero, double unit, float_image_t *im)
{
  register xel* xelrow;
  register xel* xP;
  int row, col;
  double scale = unit - zero;

  pm_message("reading pixels: zero = %.4f unit = %.4f scale = %4f", zero, unit, scale);

  /* Read pixels: */
  xelrow = pnm_allocrow(hd->cols);
  for (row = 0; row < hd->rows; ++row)
    { pnm_readpnmrow(f, xelrow, hd->cols, hd->maxval, hd->format);
      for (col = 0, xP = xelrow;  col < hd->cols; ++col, ++xP)
        { Interval *el = &(im->el[(row*hd->cols + col)*im->chns]);
          debug = ((col % 56) == 23) & ((row % 56) == 23);
          switch (hd->format)
            { case PPM_FORMAT:
              case RPPM_FORMAT:
                el[0] = floatize(PPM_GETR(*xP), hd->maxval, zero, scale);
                el[1] = floatize(PPM_GETG(*xP), hd->maxval, zero, scale);
                el[2] = floatize(PPM_GETB(*xP), hd->maxval, zero, scale);
              break;
              case PGM_FORMAT:
              case RPGM_FORMAT:
              case PBM_FORMAT:
              case RPBM_FORMAT:
                el[0] = floatize(PNM_GET1(*xP), hd->maxval, zero, scale);
              break;
              default:
                pm_error("bad input file format");
            }
          if (debug) { debug_floatized_pixel("im", col, row, im->chns, el, "\n"); }
        }
    }
  pm_close(f);
}

Interval floatize(xelval ival, xelval maxval, double zero, double scale)
/* Maps `zero' to 0, `zero+scale' to 1, with proper uncertainty. */
{ 
  double end[0] = ((double)ival - zero - 0.5)/scale;
  double end[1] = ((double)ival - zero + 0.5)/scale;
  Interval res;
  if (scale > 0.0) 
    { res.end[0] = end[0]; res.end[1] = end[1]; }
  else
    { res.end[0] = end[1]; res.end[1] = end[0]; }
  return res;
}

void foo
  (
    r2_t precision 
  )
  {
    o->precision;
    /* Compute the maximum number {maxsols} of solutions: */
    int nsx = (int)(floor(dmax->c[0]/step));  /* Max number of steps in {x}. */
    int nsy = (int)(floor(dmax->c[1]/step));  /* Max number of steps in {y}. */
    double maxdx = num_rel_shifts(ni, nsx, 0); /* Max number of alignments in {x} */
    double maxdy = num_rel_shifts(ni, nsy, 0); /* Max number of alignments in {y} */
    double maxsols = maxdx*maxdy;
    fprintf(stderr, "  maxsols = %g\n", maxsols);
    
    for (i = 0; i < n; i++)
      { 

        /* Check channel count of {img[i]}: */
        int NCi = img[i]->sz[0];
        if (i == 0) 
          { nc = NCi; }
        else 
          { demand(nc == NCi, "images must have the same channel counts"); }
        /* Check whether search domain lies inside the image's domain: */
        for (k = 0; k < 2; k++)
          { int NPi = img[i]->sz[k+1]; /* Number of pixels along axis {k}. */
            demand(p[i].c[k] - dmax.c[k] >= 0 - fudge, "point too close to low edge"); 
            demand(p[i].c[k] + dmax.c[k] <= NPi + fudge, "point too close to high edge"); 
          }
      }

  }

typedef struct pnm_header_t { /* Data for/from a PBM/PPM/PGM file header */ 
  int format, forceplain;
  int rows, cols;
  pnm_sample_t maxval;
} pnm_header_t;


#define HELP1 \
  "\\\n" \
  "  -outprefix {NAME}; namely," \
  " a projective transformation {T[i]} for each image {IM[i]} so that" \
  " the transformed images become as similar as possible.  It then " \
  " writes the transformation matrices" \
  " {T[i]} to standard output, and the transformed" \
  " images to files named \"{PREFIX}-{NN}.{EXT}\", where {PREFIX} is" \
  " a given string, {NN} is the image number (\"01\", \"02\", etc.), and" \
  " {EXT} is the appropriate extension (\"pgm\" or \"ppm\").\n" \
  "\n" \
  "  Each transformation {T[i]} is computed by finding a set" \
  " of four points {P[i,0..3]} in each image {IM[i]} such that, for" \
  " each {j}, the neighborhoods of points {P[0..N-1,j]} in images" \
  " {IM[0..N-1]} are maximally similar.\n" \
  "\n" \
  "  The domain of image {IM[i]} is assumed to be a rectangle of" \
  " the plane, such that the origin {(0,0)} is located {XC[i]} pixels" \
  " from the image's left edge, and {YC[i]} pixels from its top" \
  " edge.  The same convention is assumed for the average image, which" \
  " has {NH} by {NV} pixels. "