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

                               pamshaderelief

===============================================================================

  Generate a shaded relief image of terrain, given a terrain map - a two

  dimensional map of elevations.  A shaded relief image is an image of

  what terrain with the given elevations would look like illuminated by

  oblique light.

  The input array is a one-channel PAM image.  The sample values are

  elevations of terrain.

  This is derived from John Walker's 'pgmcrater' which not only does this

  shading, but first generates a terrain map of fractal craters on which to

  run it.

  The original program carried this attribution and license:

       Designed and implemented in November of 1989 by:

        John Walker

        Autodesk SA

        Avenue des Champs-Montants 14b

        CH-2074 MARIN

        Switzerland

        Usenet: kelvin@Autodesk.com

        Fax:    038/33 88 15

        Voice:  038/33 76 33

  Permission  to  use, copy, modify, and distribute this software and

  its documentation  for  any  purpose  and  without  fee  is  hereby

  granted,  without any conditions or restrictions.  This software is

  provided "as is" without express or implied warranty.

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

/* Modifications by Arjen Bax, 2001-06-21: Remove black vertical line at right

   edge.

*/

#define _XOPEN_SOURCE 500  /* get M_PI in math.h */

#include <assert.h>

#include <math.h>

#include "pm_c_util.h"

#include "mallocvar.h"

#include "nstring.h"

#include "shhopt.h"

#include "pam.h"

struct CmdlineInfo {

    /* All the information the user supplied in the command line,

       in a form easy for the program to use.

    */

    const char * inputFileName;

    float        gamma;

};

static void

parseCommandLine(int argc, const char ** const argv,

                 struct CmdlineInfo * const cmdlineP) {

/*----------------------------------------------------------------------------

   Note that the file spec array we return is stored in the storage that

   was passed to us as the argv array.

-----------------------------------------------------------------------------*/

    optEntry * option_def;

        /* Instructions to OptParseOptions3 on how to parse our options.

         */

    optStruct3 opt;

    unsigned int option_def_index;

    unsigned int gammaSpec;

    MALLOCARRAY_NOFAIL(option_def, 100);

    option_def_index = 0;   /* incremented by OPTENT3 */

    OPTENT3(0,   "gamma",    OPT_FLOAT,   &cmdlineP->gamma,

            &gammaSpec,       0);

    opt.opt_table = option_def;

    opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */

    opt.allowNegNum = FALSE;  /* We may have parms that are negative numbers */

    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);

        /* Uses and sets argc, argv, and some of *cmdlineP and others. */

    if (!gammaSpec)

        cmdlineP->gamma = 1.0;

    if (cmdlineP->gamma <= 0.0)

        pm_error("gamma correction must be greater than 0");

    if (argc-1 == 0) 

        cmdlineP->inputFileName = "-";

    else if (argc-1 != 1)

        pm_error("Program takes zero or one argument (filename).  You "

                 "specified %u", argc-1);

    else

        cmdlineP->inputFileName = argv[1];

    free(option_def);

}

/* Definitions for obtaining random numbers. */

/*  Display parameters  */

static double const ImageGamma = 0.5;     /* Inherent gamma of mapped image */

static int    const slopemin   = -52;

static int    const slopemax   = 52;

static void

generateSlopeGrayMap(sample * const slopeGrayMap,

                     double   const dgamma) {

/*----------------------------------------------------------------------------

   Map each possible slope to the brightness that terrain with that

   left-to-right slope should have in the shaded relief.

   The brightness is what would result from light incident from the left

   falling on the terrain.

-----------------------------------------------------------------------------*/

    double const gamma = dgamma * ImageGamma;

    int i;

    for (i = slopemin; i <= 0; ++i) {   /* Negative, downhill, dark */

        slopeGrayMap[i - slopemin] =

            128 - 127.0 * pow(sin((M_PI / 2) * i / slopemin), gamma);

    }

    for (i = 0; i <= slopemax; ++i) {   /* Positive, uphill, bright */

        slopeGrayMap[i - slopemin] =

            128 + 127.0 * pow(sin((M_PI / 2) * i / slopemax), gamma);

    }

    /* Confused?   OK,  we're using the  left-to-right slope to

       calculate a shade based on the sine of  the  angle  with

       respect  to the vertical (light incident from the left).

       Then, with one exponentiation, we account for  both  the

       inherent   gamma   of   the   image  (ad-hoc),  and  the

       user-specified display gamma, using the identity:

       (x^y)^z = (x^(y*z))

    */

}

static gray

brightnessOfSlope(int      const slope,

                  sample * const slopeGrayMap) {

    return slopeGrayMap[MIN(MAX(slopemin, slope), slopemax) - slopemin];

}

static void

writeShadedRelief(struct pam * const terrainPamP,

                  tuple **     const terrain,

                  double       const dgamma,

                  FILE *       const ofP) {

    unsigned int row;

    tuple * outrow;

    sample * slopeGrayMap; /* Slope to gray value map */

    struct pam outpam;

    outpam.size   = sizeof(outpam);

    outpam.len    = PAM_STRUCT_SIZE(tuple_type);

    outpam.file   = ofP;

    outpam.format = PAM_FORMAT;

    outpam.height = terrainPamP->height;

    outpam.width  = terrainPamP->width;

    outpam.depth  = 1;

    outpam.maxval = 255;

    outpam.bytes_per_sample = 1;

    STRSCPY(outpam.tuple_type, "GRAYSCALE");

    outrow = pnm_allocpamrow(&outpam);

    pnm_writepaminit(&outpam);

    MALLOCARRAY(slopeGrayMap, slopemax - slopemin + 1);

    generateSlopeGrayMap(slopeGrayMap, dgamma);

    for (row = 0; row < terrainPamP->height; ++row) {

        unsigned int col;

        for (col = 0; col < terrainPamP->width - 1; ++col) {

            int const slope = terrain[row][col+1][0] - terrain[row][col][0];

            outrow[col][0] = brightnessOfSlope(slope, slopeGrayMap);

        }

        {

            /* Wrap around to determine shade of pixel on right edge */

            int const slope = 

                terrain[row][0][0] - terrain[row][outpam.width-1][0];

            outrow[outpam.width - 1][0] =

                brightnessOfSlope(slope, slopeGrayMap);

        }

        pnm_writepamrow(&outpam, outrow);

    }

    free(slopeGrayMap);

    pnm_freepamrow(outrow);

}

static void

readTerrain(FILE *       const ifP,

            struct pam * const pamP,

            tuple ***    const tuplesP) {

    *tuplesP = pnm_readpam(ifP, pamP, PAM_STRUCT_SIZE(tuple_type));

}

int

main(int argc, const char ** argv) {

    struct CmdlineInfo cmdline;

    FILE * ifP;

    struct pam terrainPam;

    tuple ** terrain;

        /* Array of elevations */

    pm_proginit(&argc, argv);

    parseCommandLine(argc, argv, &cmdline);

    ifP = pm_openr(cmdline.inputFileName);

    readTerrain(ifP, &terrainPam, &terrain);

    writeShadedRelief(&terrainPam, terrain, cmdline.gamma, stdout);

    return 0;

}