/*=============================================================================
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;
}