/* ------------------------------------------------------------------
 *
 * hfcolor.c
 *
 * copyright 1998 Mark J. Stock
 *
 * v0.0 1998-07-21	initial version
 * v0.1 1998-07-24	all in/output now same size, read river mask
 * v0.2 1998-07-26	read lake mask, compute slope, fix Linux vers.
 * v0.3 1998-07-27	speckle, include in hftools package, modified
 *                      read, write subroutines
 * v0.4 1998-07-29	add export of initial tree map
 *
 * The object of hfcolor is to be a tool for creating realistic
 * images to use as surface colormaps for heightfields for the 
 * purpose of visualization. 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 PPM file of the same size is
 * written out.
 *
 * ------------------------------------------------------------------ */


#include <stdlib.h>
#include <stdio.h>
#include <string.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 height;			/* output height */
int width;			/* output width */
int inheight;			/* returned from read_pgm_ushort */
int inwidth;


/* The heightfield array, values between -32k and 32k */
unsigned short int hf[MAX_SIZE][MAX_SIZE];

/* The river mask file (PGM) */
unsigned char rm[MAX_SIZE][MAX_SIZE];

/* The lake mask file (PGM) */
unsigned char lm[MAX_SIZE][MAX_SIZE];

/* The output imagemap array */
COLOR im[MAX_SIZE][MAX_SIZE];


/* The color table, and number of usable entries */
typedef struct color_data_at_elevation {
   int elev;     /* elevation at this color */
   COLOR c;      /* color 3:0-255 */
} COLOR_ELEV;

COLOR_ELEV ctable[30];
int cbands;

/* control variables */
int wrapEdges;		/* 1 means wrap texture around edges */

COLOR new_color;

double find_slope(int,int);

/* I need these three values global so I can convert elevations on the fly */
int min_in = 66000;		/* min-max values of input */
int min_elev = -1000;		/* model min-max elevation */
double temp_d;			/* scaling factor */

int scale(int);
void Usage(char [80],int);

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

   /* initialize variables */
   int i,j;
   int temp_i, temp_j;
   char progname[80];			/* name of executable */

   char hfinfile[80];			/* name of the input PGM heightfield file */
   char outfileroot[80];		/* root name of output files */
   char imoutfile[80];			/* name of the output PPM imagemap file */
   int useElevation = 0;		/* ==1 means colors vary with elevation */
   int useSlope = 0;			/* ==1 means colors vary with slope */
   double slope_thresh = 400;		/* if altitudes vary by more than this, its rocky */
   int useFacing = 0;			/* ==1 means colors vary with normal direction */
   int useTemp = 0;			/* ==1 means colors vary with temperature */
   double tempF = 50;			/* average temperature, degrees Fahrenheit */
   int useRivers = 0;			/* ==1 means colors vary with rivers/flow volume */
   char riverMask[80];			/* name of the input PGM river mask file */
   int useLakes = 0;			/* ==1 means colors vary in lakes */
   char lakeMask[80];			/* name of the input PGM river mask file */
   int useSpeckle = 0;			/* ==1 means to speckle the colors */
   int write_com = 0;			/* write sample Radiance command */
   int write_tree = 0;			/* write initial tree map */
   char tmoutfile[80];			/* name of tree mask file */

   int max_in = -1;
   int use_min_max = 0;			/* ==1 means to use the min-max values to scale the input */
   int max_elev = 10000;		/* these are default values */
   double slope;
   double ocean_level = 0;

   /*
    * 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(hfinfile,argv[++i]); 
      else if (strncmp(argv[i], "-outfile", 2) == 0) {
         (void) strcpy(outfileroot,argv[++i]); 
         (void) sprintf(imoutfile,"%s.ppm",outfileroot); }
      else if (strncmp(argv[i], "-bounds",2) == 0) {
         use_min_max = 1;
         min_elev = atoi(argv[++i]);
         max_elev = atoi(argv[++i]); }
      else if (strncmp(argv[i], "-commands", 3) == 0)
         write_com = 1;			/* write sample commands */
      else if (strncmp(argv[i], "-treemap", 3) == 0) {
         write_tree = 1;		/* write initial tree map */
         (void) sprintf(tmoutfile,"%s.pgm",outfileroot); }
      else if (strncmp(argv[i], "-e", 2) == 0)
         useElevation = 1;		/* colors vary with elevation */
      else if (strncmp(argv[i], "-s", 2) == 0) {
         useSlope = 1;			/* colors vary with slope */
         if (strncmp(argv[i+1], "-", 1) != 0)
            slope_thresh = slope_thresh*atof(argv[++i]); }
      else if (strncmp(argv[i], "-f", 2) == 0)
         useFacing = 1;			/* colors vary with normal direction */
      else if (strncmp(argv[i], "-t", 2) == 0) {
         useTemp = 1;			/* colors vary with temperature */
         tempF = atof(argv[++i]); }
      else if (strncmp(argv[i], "-r", 2) == 0) {
         useRivers = 1;			/* colors vary with rivers/flow volume */
         (void) strcpy(riverMask,argv[++i]); }
      else if (strncmp(argv[i], "-l", 2) == 0) {
         useLakes = 1;			/* colors vary in lakes */
         (void) strcpy(lakeMask,argv[++i]); }
      else if (strncmp(argv[i], "-k", 2) == 0)
         useSpeckle = 1;		/* speckle the colors */
      else
         (void) Usage(progname,0);
   }

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

   /* check new elevation bounds, hf is unsigned short int */
   /* honestly, this does not make any difference */
   /* if (max_elev > 32768 || min_elev < -32768) */
      /* write_error("New elevation bounds too extreme, keep within +-32768",1); */

   /* make sure seed is not too big */
   /* if (seed > 100000) seed = 12345;
   for (i=0; i<seed; i++) (void) rand(); */

   /* create the default color table */
   if (useElevation) (void)generate_colormap(15);


   /*
    * Read in heightfield and other data
    */

   (void) read_pgm_ushort(hfinfile,hf);

   /* assign output sizes based on the size of the hf PGM */
   height = inheight;
   width = inwidth;

   /* scale the output values using min-max bounds */
   if (use_min_max) {
      /* then, find min/max of the input */
      for (i=0; i<height; i++) {
         for (j=0; j<width; j++) {
            temp_i = hf[i][j];
            if (temp_i > max_in) max_in = temp_i;
            if (temp_i < min_in) min_in = temp_i;
         }
      }
      fprintf(stdout,"scaling elevations from %d %d to %d %d for internal use\n"
                    ,min_in,max_in,min_elev,max_elev);
      /* then, compute scaling factor */
      temp_d = (max_elev - min_elev) / (double)(max_in - min_in);
      /* lastly, apply scaling to all pixels */
      /* for (i=0; i<height; i++) {
         for (j=0; j<width; j++) {
            hf[i][j] = min_elev + temp_d*(hf[i][j]-min_in);
         }
      }*/
      ocean_level = min_in + (max_in-min_in) * min_elev/(min_elev-max_elev);
      if (ocean_level > min_in) fprintf(stdout,"ocean_level is %lf units\n",ocean_level);
   } else {
      ocean_level = 0;
   }

   /* read in the river mask PGM file */
   if (useRivers) {
      (void) read_pgm_char(riverMask,rm);
      /* check the size returned */
      if (inheight != height || inwidth != width)
         write_error("River mask is not same size as heightfield.");
   }

   /* read in the lake mask PGM file */
   if (useLakes) {
      (void) read_pgm_char(lakeMask,lm);
      /* check the size returned */
      if (inheight != height || inwidth != width)
         write_error("Lake mask is not same size as heightfield.");
   }


   /*
    * Calculate colors
    */

   /* base color on elevation of output pixel */
   for (i=0; i<height; i++) {
      for (j=0; j<width; j++) {
         if (use_min_max) temp_i = scale(hf[i][j]);
         else temp_i = hf[i][j];
         /* fprintf(stderr,"orig elev= %d, new elev= %d\n",hf[i][j],temp_i); */
         if (useElevation) {
            /* find elevation of output pixel location */
            /* (void) interp_color(hf[i][j]); */
            (void) interp_color(temp_i);
            im[i][j] = new_color;
         } else {
            /* base the color on the temperature and humidity alone */
            /* (void) choose_color();
            im[i][j] = new_color; */
         }
         /* if (useSlope && hf[i][j] > -1) { */
         if (useSlope && temp_i > -1) {
            /* use 2nd derivative to find local depressions, raise
             * their effective elevations (or lower them) to determine
             * snow chance */
            if (temp_i > 5000) {
               temp_j = scale(snow_collection_elevation(i,j));
               if (temp_j == -1) {
                  /* its steep, make it rock */
                  im[i][j].r = 160;
                  im[i][j].g = 190;
                  im[i][j].b = 180;
               } else {
                  /* use color lookup, change it */
                  (void) interp_color(temp_j);
                  im[i][j] = new_color;
               }
            }

            /* calculate slope using normal */
            /* find max 1st derivative */
            slope = find_slope(i,j);
            if (slope > slope_thresh) {
               im[i][j].r = 160;
               im[i][j].g = 190;
               im[i][j].b = 180;
            }
         }
         if (useSpeckle && temp_i > 0) {
            (void) speckleColor(&im[i][j]);
         }
         if (useRivers) {
            if (rm[i][j] > 50 && temp_i > -1) {
               im[i][j].r = (im[i][j].r*(255-rm[i][j]))/256.0;
               im[i][j].g = (im[i][j].g*(255-rm[i][j]))/256.0;
               im[i][j].b = rm[i][j] + im[i][j].b*(255-rm[i][j])/256.0;
            }
         }
         if (useLakes) {
            if (lm[i][j] > 50 && temp_i > -1) {
               im[i][j].r = 4;
               im[i][j].g = 40;
               im[i][j].b = 165;
            }
         }
      }
   }


   /*
    * Write output image
    */

   (void) write_ppm(imoutfile,im,height,width);
   fprintf(stdout,"imagemap is %d wide by %d high\n",width,height);

   /* Write an initial tree mask file, use lake and river mask data */
   if (write_tree) {
      /* write over the river mask data */
      for (i=0; i<height; i++) {
         for (j=0; j<width; j++) {
            temp_i = scale(hf[i][j]);
            rm[i][j] = 128 + rm[i][j]/2;
            if (lm[i][j] > 10) rm[i][j] = 0;
            if (temp_i > 5000 && temp_i < 7000) rm[i][j] = 128 - (temp_i-5000)/16;
            else if (temp_i < 1 || temp_i >= 7000) rm[i][j] = 0;
         }
      }
      (void) write_pgm_char(tmoutfile,rm,height,width);
      fprintf(stderr,"Printed tree mask as %s\n",tmoutfile);
   }

   /*
    * Write output statistics and quit
    */
   if (write_com) write_commands(hfinfile,imoutfile,ocean_level);
   (void) write_finish();
}


/*
 * scale an input elevation to an internally-scaled elevation
 */
int scale(int inelev) {
   return ((int)(min_elev + temp_d*(inelev-min_in)));
}


/*
 * Generate a series of COLOR_ELEV entries, assume
 * elevations are in meters relative to sea level
 */
int generate_colormap(int avg_temp) {

   int i;

   /* do the oceanic colors first, deep to shallow */
   for (i=0; i<4; i++) {
      ctable[i].elev = -500 + 163*i;
      ctable[i].c.r = 3 + 30*i;
      ctable[i].c.g = 41 + 26*i;
      ctable[i].c.b = 163 + 20*i;
   }
   /* do sandy-like colors for near-beach */
   for (i=4; i<6; i++) {
      ctable[i].elev = -2 + 4*(i-4);
      ctable[i].c.r = 220;
      ctable[i].c.g = 210;
      ctable[i].c.b = 180;
   }
   /* do shades of green for land colors */
   for (i=6; i<25; i++) {
      ctable[i].elev = 3 + 18*(i-6)*(i-6);
      ctable[i].c.r = 55 + 4*(i-6);
      ctable[i].c.g = 88 + 6*(i-6);
      ctable[i].c.b = 50 + 4*(i-6);
   }
   /* lastly, end it with white, for mountain caps */
   i = 25;
   ctable[i].elev = 6000;
   ctable[i].c.r = 250;
   ctable[i].c.g = 250;
   ctable[i].c.b = 250;

   cbands = 26;

   /* just to be safe */
   ctable[cbands].elev = 0;

   /* for testing purposes only */
   /* for (i=-1000; i<10000; i+=100) {
      interp_color(i);
      fprintf(stderr,"elev= %d, color= %d %d %d\n",i,new_color.r,new_color.g,new_color.b);
   } */

   return(0);
}


/*
 * Interpolate a color (saved as COLOR new_color) based on
 * the elevation sent...uses data from ctable
 */
int interp_color(int elev) {

   int i = 0;
   double elev_frac = 0.0;

   while (elev > ctable[i].elev && i < cbands) i++;

   if (i==0) new_color = ctable[0].c;
   else if (i==cbands) new_color = ctable[cbands-1].c;
   else {
      /* do some interpolating */
      /* fprintf(stderr,"ctable (i-1) is %d %d %d\n",ctable[i-1].c.r,ctable[i-1].c.g,ctable[i-1].c.b); */
      elev_frac = (elev - ctable[i-1].elev)/(double)(ctable[i].elev - ctable[i-1].elev);
      new_color.r = ctable[i-1].c.r + elev_frac * (ctable[i].c.r - ctable[i-1].c.r);
      new_color.g = ctable[i-1].c.g + elev_frac * (ctable[i].c.g - ctable[i-1].c.g);
      new_color.b = ctable[i-1].c.b + elev_frac * (ctable[i].c.b - ctable[i-1].c.b);
   }

   /* fprintf(stderr,"elev of %d produced color %d %d %d, frac= %lf\n",elev,new_color.r,new_color.g,new_color.b,elev_frac); */
}

double find_slope(int i, int j) {

   double top, bottom, left, right;

   if (i==0) top = hf[0][j];
   else top =  hf[i-1][j];
   if (i==height-1) bottom = hf[height-1][j];
   else bottom = hf[i+1][j];
   if (j==0) left = hf[i][0];
   else left = hf[i][j-1];
   if (j==width-1) right = hf[i][0];
   else right = hf[i][j+1];

   right = abs(left-right);
   left = abs(top-bottom);
   if (right > left) return(right);
   else return(left);

}


/*
 * If the location is in a local depression, there is a
 * greater liklihood of snow collecting there. To simulate
 * this, that location's effective elevation is changed
 * before the color table lookup. This assumes -e is 
 * used, too. Bummer.
 */
int snow_collection_elevation(int i,int j) {

   int top, bottom, left, right;

   if (i==0) top = hf[0][j];
   else top =  hf[i-1][j];
   if (i==height-1) bottom = hf[height-1][j];
   else bottom = hf[i+1][j];
   if (j==0) left = hf[i][0];
   else left = hf[i][j-1];
   if (j==width-1) right = hf[i][0];
   else right = hf[i][j+1];

   top = top+bottom+left+right - 4*hf[i][j];

   if (top > -50) {
      return(hf[i][j] + 15*top);
   } else {
      return(hf[i][j]);
   }
}


/*
 * change the color's intensity an equal %age
 */
int speckleColor(COLOR* color) {

   double mult = 0.2*rand()/(RAND_MAX+1.0) + 0.9;
   double tmp_d;

   tmp_d = color->r*mult;
   if ((int)tmp_d > 255) color->r = 255;
   else color->r = (unsigned char)tmp_d;
   tmp_d = color->g*mult;
   if ((int)tmp_d > 255) color->g = 255;
   else color->g = (unsigned char)tmp_d;
   tmp_d = color->b*mult;
   if ((int)tmp_d > 255) color->b = 255;
   else color->b = (unsigned char)tmp_d;

   return(0);
}


/*
 * Help the user out by writing sample Unix and Radiance
 * statements and commands to incorporate the imagemap into
 * a rendering
 */
int write_commands(char hffile[80],char ppmfile[80],double oceanlevel) {

   fprintf(stderr,"\nRun the following command:\n");
   fprintf(stderr,"ra_ppm -r %s > groundimage.pic\n",ppmfile);

   fprintf(stderr,"\nThen, try scaling the heightfield with one of these:\n");
   fprintf(stderr,"hfscale -in %s -out larger.pgm -scale x.0\n",hffile);
   fprintf(stderr,"pnmscale 0.x %s > smaller.pgm\n",hffile);
   fprintf(stderr,"\nLastly, make the .rad files with:\n");
   if (oceanlevel > 0.001)
      fprintf(stderr,"hf2rad -in %s -out hf -ocean %lf\n",hffile,oceanlevel);
   else
      fprintf(stderr,"hf2rad -in %s -out hf\n",hffile);

   return(0);
}


/*
 * This function writes basic usage information to stderr,
 * and then quits. Too bad.
 */
void Usage(char progname[80],int status) {

   /* Usage for hfcolor */
   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      root name of output files",
       "  -commands      write series of commands to incorporate imagemap",
       "  -tree          write an initial tree mask as a PGM image",
       " ",
       "  -e             use elevation data to determine image colors",
       "  -s cutoff      use slope to determine image colors, slope>cutoff means rocks",
       "  -r maskfile    use river mask PGM file to determine image colors",
       "  -l lakefile    use lake mask PGM file to determine image colors",
       "  -k             speckle the land colors",
       " ",
       "  -bounds n n    minimum and maximum elevation in the input, def= -1000, 1000",
       /* "  -seed val      set the seed for the psuedo-random number generator", */
       " ",
       "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);
}
