/* ------------------------------------------------------------------
 *
 * hfscale.c
 *
 * copyright 1998 Mark J. Stock
 *
 * v0.0 1998-07-24	initial version
 * v0.1 1998-07-27	working, correct scaling with interpolation
 * v0.2 1998-07-28	added rescaling elevations using bounds
 * v0.3 1998-08-20	now writes raw PGM if maxval is 255 or less
 *
 * The object of hfscale is to be a tool for rescaling heightfields
 * A heightfield, in the format of a binary of ASCII, 8 or 16-bit, 
 * PGM file, is read in, calculations are performed, and a binary 
 * PGM file of the same h/w ratio is written out.
 *
 * This is better than resizing in xv (no 16-bit gray support) and
 * better than using pnmscale (no interpolation)
 *
 * BUT, for scales less than 1.0, use pnmscale, because this prog
 * does not support averaging!
 *
 * ------------------------------------------------------------------ */


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "const.h"

/*
 * All 2D arrays are indexed from top left corner, 1st index
 * is row number, 2nd index is column
 */

/* Important sizes */
int inheight;			/* input height */
int inwidth;			/* input width */
int outheight;			/* output height */
int outwidth;			/* output width */

/* The input array, values between 0 and 65k */
unsigned short int in[MAX_SIZE][MAX_SIZE];

/* The output-size elevation array */
unsigned short int out[MAX_SIZE][MAX_SIZE];
unsigned char outc[MAX_SIZE][MAX_SIZE];

double double_interpolate(double,double,double,double,double,double);
void Usage(char [80],int);

main(int argc, char *argv[]) {

   /* initialize variables */
   int i,j;
   double temp_d;
   double x_offset, y_offset, x_exact, y_exact;
   char progname[80];			/* name of executable */

   char infile[80];			/* name of the input PGM heightfield file */
   char outfile[80];			/* name of the output PPM imagemap file */
   double scale = 1.0;			/* output size, as multiple of input size */
   int output_width = 0;		/* forced output width, pixels */

   int min_in = 66000;                  /* min-max values of input */
   int max_in = -1;
   int use_min_max = 0;                 /* ==1 means to use the min-max values to scale the input */
   int min_out = 0;                    /* model min-max elevation */
   int max_out = 255;
   double range;

   /*
    * Parse the command-line options
    */
   (void) strcpy(progname,argv[0]);
   if (argc < 4) (void) Usage(progname,0);
   for (i=1; i<argc; i++) {
      if (strncmp(argv[i], "-help", 2) == 0)
         (void) Usage(progname,0);
      else if (strncmp(argv[i], "-infile", 2) == 0)
         (void) strcpy(infile,argv[++i]); 
      else if (strncmp(argv[i], "-outfile", 2) == 0)
         (void) strcpy(outfile,argv[++i]); 
      else if (strncmp(argv[i], "-scale", 2) == 0)
         scale = atof(argv[++i]);	/* scale the output imagemap up or down */
      else if (strncmp(argv[i], "-width", 2) == 0)
         output_width = atoi(argv[++i]); /* force output width, overrides scale */
      else if (strncmp(argv[i], "-bounds",2) == 0) {
         use_min_max = 1;
         min_out = atoi(argv[++i]);
         max_out = atoi(argv[++i]); }
      else
         (void) Usage(progname,0);
   }


   /*
    * Act on user inputs, set up variables
    */

   /* warn user that this prog is for upsizing, really */
   if (scale < 0.8) fprintf(stderr,"You should really use pnmscale %lf %s > %s for this operation.\n\n",scale,infile,outfile);

   /* check output size */
   if (output_width > MAX_SIZE)
      write_error("Output size too large, reduce and try again.");


   /*
    * Read in heightfield
    */

   (void) read_pgm_ushort(infile,in);

   /* determine output size, based on input size */ 
   if (output_width != 0) {
      /* then user defined output_width, not scale */
      scale = output_width/(double)inwidth;
      outheight = scale*inheight;
      outwidth = output_width;
   } else {
      /* use scale, which defaults to 1.0 */
      outheight = scale*inheight;
      outwidth = scale*inwidth;
   }

   /* make sure output size is within bounds */
   if (outwidth > MAX_SIZE || outheight > MAX_SIZE)
      write_error("Output size too large, try a different scale.",1);
   if (outwidth < 1 || outheight < 1)
      write_error("Output size is zero, try a different scale.",1);

   /* scale heightfield (in) up to output size, save as double */
   if (outwidth != inwidth) {
      fprintf(stderr,"Scaling image from %d %d to %d %d\n",inwidth,inheight,outwidth,outheight);
      x_offset = 0.5/outwidth - 0.5/inwidth;
      y_offset = 0.5/outheight - 0.5/inheight;
      for (i=0; i<outheight; i++) {
         y_exact = i/scale + y_offset;
         if (y_exact < 0.0) y_exact = 0.0;
         if (y_exact > inheight - 1.0) y_exact = inheight - 1.0;
         for (j=0; j<outwidth; j++) {
            /* find elevation of output pixel location */
            x_exact = j/scale + x_offset;
            if (x_exact < 0.0) x_exact = 0.0;
            if (x_exact > inwidth - 1.0) x_exact = inwidth - 1.0;
            temp_d = double_interpolate(in[(int)y_exact+1][(int)x_exact],
                                   in[(int)y_exact+1][(int)x_exact+1],
                                   in[(int)y_exact][(int)x_exact+1],
                                   in[(int)y_exact][(int)x_exact],
                                   x_exact-floor(x_exact),
                                   1.0-y_exact+floor(y_exact));
            out[i][j] = (unsigned short int)temp_d;
         }
      }
   } else {
      /* just copy the values from in to out */
      fprintf(stderr,"Not scaling image, will stay %d %d\n",outwidth,outheight);
      for (i=0; i<outheight; i++) {
         for (j=0; j<outwidth; j++) {
            out[i][j] = in[i][j];
         }
      }
   }

   /* scale the elevations */
   if (use_min_max) {
      /* first, find min/max */
      for (i=0; i<outheight; i++) {
         for (j=0; j<outwidth; j++) {
            if (out[i][j] < min_in) min_in = out[i][j];
            if (out[i][j] > max_in) max_in = out[i][j];
         }
      }
      /* then, scale all of the values */
      fprintf(stderr,"Rescaling elevations from %d %d to %d %d\n",min_in,max_in,min_out,max_out);
      range = (double)(max_out-min_out)/(double)(max_in-min_in);
      for (i=0; i<outheight; i++) {
         for (j=0; j<outwidth; j++) {
            out[i][j] = min_out + range*(out[i][j]-min_in);
         }
      }
   }


   /*
    * Write output image
    */

   if (outwidth != inwidth || use_min_max) {
      if (max_out < 256) {
         for (i=0; i<outheight; i++) {
            for (j=0; j<outwidth; j++) {
               outc[i][j] = (unsigned char)out[i][j];
            }
         }
         (void) write_pgm_char(outfile,outc,outheight,outwidth);
      } else
         (void) write_pgm_ushort(outfile,out,outheight,outwidth);
   } else 
      fprintf(stderr,"No changes made, no file written\n\n");


   /*
    * Write output statistics and quit
    */
   (void) write_finish();
}


/*
 * Do a double interpolation using 4 values as the corners,
 * and a dx and dy fractions
 */
double double_interpolate(double bl, double br, double tr,
                          double tl, double dx, double dy) {

   double topx, bottomx, retval;

   topx = tl + dx*(tr-tl);
   bottomx = bl + dx*(br-bl);
   retval = bottomx + dy*(topx-bottomx);

   return(retval);
}


void Usage(char progname[80],int status) {

   /* Usage for hfscale */
   static char **cpp, *help_message[] =
   {
       "where options include:",
       "  -help          writes this help information",
       " ",
       "  -infile name   name of input PGM file, binary or raw",
       "  -out name      name of output PGM file, raw (ascii)",
       "  -scale val     rescale the heightfield this multiple",
       "  -width val     rescale the heightfield to this width",
       "  -bounds n n    minimum and maximum elevation in the input, def= 0,255",
       " ",
       "Options may be abbreviated to an unambiguous length.",
       NULL
   };

   fprintf(stderr, "usage:\n  %s [-options ...]\n\n", progname);
   for (cpp = help_message; *cpp; cpp++)
      fprintf(stderr, "%s\n", *cpp);
   exit(status);
}
