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