Split Pgmcrater into Pamcrater and Pamshadedrelief

git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@2231 9d0c8265-081b-0410-96cb-a4ca84ce46f8

6 files changed, 792 insertions, 496 deletions

diff --git a/converter/other/pamtopnm.c b/converter/other/pamtopnm.c
index 9bb662b7..f043d721 100644
--- a/converter/other/pamtopnm.c
+++ b/converter/other/pamtopnm.c
@@ -17,23 +17,23 @@
 #include "shhopt.h"
 #include "mallocvar.h"
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
-    const char *inputFilespec;  /* Filespecs of input files */
-    unsigned int assume;    /* -assume option */
+    const char * inputFileName;  /* Name of input file */
+    unsigned int assume;
 };
 
 
 static void
-parseCommandLine(int argc, char ** argv,
-                 struct cmdlineInfo *cmdlineP) {
+parseCommandLine(int argc, const char ** argv,
+                 struct CmdlineInfo *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;
+    optEntry * option_def;
         /* Instructions to OptParseOptions3 on how to parse our options.
          */
     optStruct3 opt;
@@ -49,16 +49,18 @@ parseCommandLine(int argc, char ** argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    pm_optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (argc-1 == 0) 
-        cmdlineP->inputFilespec = "-";
+        cmdlineP->inputFileName = "-";
     else if (argc-1 != 1)
         pm_error("Program takes zero or one argument (filename).  You "
                  "specified %d", argc-1);
     else
-        cmdlineP->inputFilespec = argv[1];
+        cmdlineP->inputFileName = argv[1];
+
+    free(option_def);
 }
 
 
@@ -105,19 +107,19 @@ validateTupleType(struct pam const inpam,
 
 
 int
-main(int argc, char *argv[]) {
+main(int argc, const char **argv) {
 
-    struct cmdlineInfo cmdline;
-    FILE* ifP;
+    struct CmdlineInfo cmdline;
+    FILE * ifP;
     int eof;   /* no more images in input stream */
     struct pam inpam;   /* Input PAM image */
     struct pam outpam;  /* Output PNM image */
 
-    pnm_init(&argc, argv);
+    pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    ifP = pm_openr(cmdline.inputFilespec);
+    ifP = pm_openr(cmdline.inputFileName);
 
     eof = FALSE;
     while (!eof) {
diff --git a/generator/Makefile b/generator/Makefile
index 3c30cdd0..d0ea6b60 100644
--- a/generator/Makefile
+++ b/generator/Makefile
@@ -14,9 +14,10 @@ include $(BUILDDIR)/config.mk
 # This package is so big, it's useful even when some parts won't 
 # build.
 
-PORTBINARIES = pamgauss pamgradient pamseq pamstereogram \
+PORTBINARIES = pamcrater pamgauss pamgradient \
+	       pamseq pamshadedrelief pamstereogram \
 	       pbmpage pbmmake pbmtext pbmtextps pbmupc \
-	       pgmcrater pgmkernel pgmmake pgmnoise pgmramp \
+	       pgmkernel pgmmake pgmnoise pgmramp \
 	       ppmcie ppmcolors ppmforge ppmmake ppmpat ppmrough ppmwheel \
 
 # We don't include programs that have special library dependencies in the
@@ -28,7 +29,7 @@ MERGEBINARIES = $(PORTBINARIES)
 
 
 BINARIES = $(MERGEBINARIES) $(NOMERGEBINARIES)
-SCRIPTS = ppmrainbow
+SCRIPTS = pgmcrater ppmrainbow
 
 OBJECTS = $(BINARIES:%=%.o)
 
diff --git a/generator/pamcrater.c b/generator/pamcrater.c
new file mode 100644
index 00000000..d61ce548
--- /dev/null
+++ b/generator/pamcrater.c
@@ -0,0 +1,428 @@
+/*=============================================================================
+                               pamcrater
+===============================================================================
+  Fractal cratering
+
+  This is derived from John Walker's 'pgmcrater' which not only creates
+  the terrain map as this program does, but then does a relief filter to
+  convert it to a shaded visual image.
+
+  The  algorithm  used  to  determine crater size is as described on
+  pages 31 and 32 of:
+
+  Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
+      Images, New York: Springer Verlag, 1988.
+
+  The  mathematical  technique  used  to calculate crater radii that
+  obey the proper area law distribution from a uniformly distributed
+  pseudorandom sequence was developed by Rudy Rucker.
+
+  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. Make craters wrap around the image (enables tiling of image).
+ */
+
+#define _XOPEN_SOURCE   /* get M_PI in math.h */
+
+#include <assert.h>
+#include <math.h>
+
+#include "pm_c_util.h"
+#include "mallocvar.h"
+#include "shhopt.h"
+#include "nstring.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.
+    */
+    unsigned int number;
+    unsigned int height;
+    unsigned int width;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
+    unsigned int test;
+    unsigned int radius;
+};
+
+
+
+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 numberSpec, heightSpec, widthSpec, radiusSpec;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+
+    option_def_index = 0;   /* incremented by OPTENT3 */
+    OPTENT3(0,   "number",     OPT_UINT,    &cmdlineP->number,
+            &numberSpec,                 0);
+    OPTENT3(0,   "height",     OPT_UINT,    &cmdlineP->height,
+            &heightSpec,                 0);
+    OPTENT3(0,   "width",      OPT_UINT,    &cmdlineP->width,
+            &widthSpec,                  0);
+    OPTENT3(0,   "randomseed", OPT_UINT,    &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec,   0);
+    OPTENT3(0,   "test",       OPT_FLAG,    NULL,
+            &cmdlineP->test,       0);
+    OPTENT3(0,   "radius",     OPT_UINT,    &cmdlineP->radius,
+            &radiusSpec,           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 (argc-1 > 0)
+        pm_error("There are no non-option arguments.  You specified %u",
+                 argc-1);
+
+    if (!heightSpec)
+        cmdlineP->height = 256;
+
+    if (cmdlineP->height == 0)
+        pm_error("-height must be positive");
+
+    if (!widthSpec)
+        cmdlineP->width = 256;
+
+    if (cmdlineP->width == 0)
+        pm_error("-width must be positive");
+
+    if (cmdlineP->test) {
+        if (!radiusSpec)
+            pm_error("With -test, you must specify -radius");
+        else {
+            if(MAX(cmdlineP->height, cmdlineP->width) * 2 < cmdlineP->radius)
+                pm_error("Radius (%u) too large", cmdlineP->radius);
+
+            if (numberSpec)
+                pm_error("-number is meaningless with -test");
+
+            if (cmdlineP->randomseedSpec)
+                pm_error("-randomseed is meaningless with -test");
+        }
+    } else {
+        if (radiusSpec)
+            pm_error("-radius is meaningful only with -test");
+
+        if (!numberSpec)
+            cmdlineP->number = 50000;
+
+        if (cmdlineP->number == 0)
+            pm_error("-number must be positive");
+    }
+    free(option_def);
+}
+
+
+/* Definitions for obtaining random numbers. */
+
+/*  Display parameters  */
+
+static double const arand       = 32767.0;  /* Random number parameters */
+static double const CdepthPower = 1.5;      /* Crater depth power factor */
+static double const DepthBias2  = 0.5;      /* Square of depth bias */
+
+
+
+static double const
+cast(double const high) {
+
+    return high * ((rand() & 0x7FFF) / arand);
+}
+
+
+
+static unsigned int
+mod(int          const t,
+    unsigned int const n) {
+
+    /* This is used to transform coordinates beyond bounds into ones
+       within: craters "wrap around" the edges.  This enables tiling
+       of the image.
+
+       Produces strange effects when crater radius is very large compared
+       to image size.
+    */
+
+    int m;
+
+    m = t % (int)n;
+
+    if (m < 0)
+        m += n;
+
+    return m;
+}
+
+
+
+static sample *
+terrainModP(struct pam * const pamP,
+            tuple **     const terrain,
+            int          const x,
+            int          const y) {
+
+    return &terrain[mod(y, pamP->height)][mod(x, pamP->width)][0];
+}
+
+
+
+
+static sample
+terrainMod(struct pam * const pamP,
+           tuple **     const terrain,
+           int          const x,
+           int          const y) {
+
+    return *terrainModP(pamP, terrain, x, y);
+}
+
+
+
+static void
+smallCrater(struct pam * const pamP,
+            tuple **     const terrain,
+            int          const cx,
+            int          const cy,
+            double       const g) {
+/*----------------------------------------------------------------------------
+   Generate a crater with a special method for tiny craters.
+-----------------------------------------------------------------------------*/
+    int y;
+    unsigned int amptot;
+    unsigned int npatch;
+
+    /* Set pixel to the average of its Moore neighborhood. */
+
+    for (y = cy - 1, amptot = 0, npatch = 0; y <= cy + 1; ++y) {
+        int x;
+        for (x = cx - 1; x <= cx + 1; ++x) {
+            amptot += terrainMod(pamP, terrain, x, y);
+            ++npatch;
+        }
+    }
+    {
+        unsigned int const axelev = amptot / npatch;
+        
+        /* Perturb the mean elevation by a small random factor. */
+
+        int const x = g >= 1 ? ((rand() >> 8) & 3) - 1 : 0;
+        *terrainModP(pamP, terrain, cx, cy) = axelev + x;
+    }
+}
+
+
+
+static void
+normalCrater(struct pam * const pamP,
+             tuple **     const terrain,
+             int          const cx,
+             int          const cy,
+             double       const radius) {
+/*----------------------------------------------------------------------------
+   Generate a regular (not tiny) crater.
+
+   Generate an impact feature of the correct size and shape.
+-----------------------------------------------------------------------------*/
+    int    const impactRadius = (int) MAX(2, (radius / 3));
+    int    const craterRadius = (int) radius;
+    double const rollmin      = 0.9;
+
+    int y;
+    unsigned int amptot, axelev;
+    unsigned int npatch;
+
+    /* Determine mean elevation around the impact area.
+       We assume the impact area is a fraction of the total crater size. */
+
+    for (y = cy - impactRadius, amptot = 0, npatch = 0;
+         y <= cy + impactRadius;
+         ++y) {
+        int x;
+        for (x = cx - impactRadius; x <= cx + impactRadius; ++x) {
+            amptot += terrainMod(pamP, terrain, x, y);
+            ++npatch;
+        }
+    }
+    assert(npatch > 0);
+    axelev = amptot / npatch;
+
+    for (y = cy - craterRadius; y <= cy + craterRadius; ++y) {
+        int const dysq = (cy - y) * (cy - y);
+
+        int x;
+
+        for (x = cx - craterRadius; x <= cx + craterRadius; ++x) {
+            int  const dxsq = (cx - x) * (cx - x);
+            double const cd = (dxsq + dysq) /
+                              (double) (craterRadius * craterRadius);
+            double const cd2 = cd * 2.25;
+            double const tcz = sqrt(DepthBias2) - sqrt(fabs(1 - cd2));
+            double cz;
+            double roll;
+
+            cz = MAX((cd2 > 1) ? 0.0 : -10, tcz);  /* Initial value */
+
+            cz *= pow(craterRadius, CdepthPower);
+            if (dysq == 0 && dxsq == 0 && ((int) cz) == 0) {
+                cz = cz < 0 ? -1 : 1;
+            }
+
+            roll = (((1 / (1 - MIN(rollmin, cd))) /
+                     (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
+
+            {
+                unsigned int av;
+                av = (axelev + cz) * (1 - roll) +
+                    (terrainMod(pamP, terrain, x, y) + cz) * roll;
+                av = MAX(1000, MIN(64000, av));
+                
+                *terrainModP(pamP, terrain, x, y) = av;
+            }
+        }
+    }
+}
+
+
+
+/* We should also have largeCrater() */
+
+
+
+static void
+plopCrater(struct pam * const pamP,
+           tuple **     const terrain,
+           int          const cx,
+           int          const cy,
+           double       const radius) {
+
+    if (radius < 3)
+        smallCrater (pamP, terrain, cx, cy, radius);
+    else
+        normalCrater(pamP, terrain, cx, cy, radius);
+}
+
+
+
+static void
+genCraters(struct CmdlineInfo const cmdline) {
+/*----------------------------------------------------------------------------
+   Generate cratered terrain
+-----------------------------------------------------------------------------*/
+    tuple ** terrain;    /* elevation array */
+    unsigned int row;
+    struct pam pam;
+
+    /* Allocate the elevation array and initialize it to mean surface
+       elevation.
+    */
+
+    pam.size   = sizeof(pam);
+    pam.len    = PAM_STRUCT_SIZE(tuple_type);
+    pam.file   = stdout;
+    pam.format = PAM_FORMAT;
+    pam.height = cmdline.height;
+    pam.width  = cmdline.width;
+    pam.depth  = 1;
+    pam.maxval = 65535;
+    pam.bytes_per_sample = 2;
+    STRSCPY(pam.tuple_type, "elevation");
+
+    terrain = pnm_allocpamarray(&pam);
+
+    for (row = 0; row < pam.height; ++row) {
+        unsigned int col;
+        for (col = 0; col < pam.width; ++col)
+            terrain[row][col][0] = pam.maxval / 2;
+    }
+    
+    if (cmdline.test)
+        plopCrater(&pam, terrain,
+                   pam.width/2, pam.height/2, (double) cmdline.radius);
+    else {
+        unsigned int const ncraters = cmdline.number; /* num of craters */
+        unsigned int l;
+
+        for (l = 0; l < ncraters; ++l) {
+            int const cx = cast((double) pam.width  - 1);
+            int const cy = cast((double) pam.height - 1);
+
+            /* Thanks, Rudy, for this equation that maps the uniformly
+               distributed numbers from cast() into an area-law distribution
+               as observed on cratered bodies.
+
+               Produces values within the interval:
+               0.56419 <= radius <= 56.419
+            */
+            double const radius = sqrt(1 / (M_PI * (1 - cast(0.9999))));
+
+            plopCrater(&pam, terrain, cx, cy, radius);
+
+            if (((l + 1) % 100000) == 0)
+                pm_message("%u craters generated of %u (%u%% done)",
+                           l + 1, ncraters, ((l + 1) * 100) / ncraters);
+        }
+    }
+
+    pnm_writepam(&pam, terrain);
+
+    pnm_freepamarray(terrain, &pam);
+
+    pm_close(stdout);
+}
+
+
+
+int
+main(int argc, const char ** argv) {
+
+    struct CmdlineInfo cmdline;
+
+    pm_proginit(&argc, argv);
+
+    parseCommandLine(argc, argv, &cmdline);
+
+    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
+
+    genCraters(cmdline);
+
+    return 0;
+}
+
+
+
diff --git a/generator/pamshadedrelief.c b/generator/pamshadedrelief.c
new file mode 100644
index 00000000..89996c83
--- /dev/null
+++ b/generator/pamshadedrelief.c
@@ -0,0 +1,250 @@
+/*=============================================================================
+                               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   /* 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;
+}
+
+
diff --git a/generator/pgmcrater b/generator/pgmcrater
new file mode 100755
index 00000000..1c22ed70
--- /dev/null
+++ b/generator/pgmcrater
@@ -0,0 +1,94 @@
+#!/bin/sh
+
+##############################################################################
+# This is essentially a Perl program.  We exec the Perl interpreter specifying
+# this same file as the Perl program and use the -x option to cause the Perl
+# interpreter to skip down to the Perl code.  The reason we do this instead of
+# just making /usr/bin/perl the script interpreter (instead of /bin/sh) is
+# that the user may have multiple Perl interpreters and the one he wants to
+# use is properly located in the PATH.  The user's choice of Perl interpreter
+# may be crucial, such as when the user also has a PERL5LIB environment
+# variable and it selects modules that work with only a certain main
+# interpreter program.
+#
+# An alternative some people use is to have /usr/bin/env as the script
+# interpreter.  We don't do that because we think the existence and
+# compatibility of /bin/sh is more reliable.
+#
+# Note that we aren't concerned about efficiency because the user who needs
+# high efficiency can use directly the programs that this program invokes.
+#
+##############################################################################
+
+exec perl -w -x -S -- "$0" "$@"
+
+#!/usr/bin/perl
+##############################################################################
+#  This is nothing but a compatibility interface for
+#  Pamcrater/Pamshadedrelief.  An old program coded to call Pgmcrater will
+#  continue working because this interface exists.  All new (or newly
+#  modified) programs should call Pamcrater and Pamshadedrelief instead.
+#
+#  In days past, Pamcrater and Pamshadedrelief did not exist.  Pgmcrater did
+#  both jobs together, with PGM output.
+##############################################################################
+
+use strict;
+
+use Getopt::Long;
+
+my @pgmcraterArgv = @ARGV;
+
+my $validOptions = GetOptions(
+    'number=i'     => \my $numberOpt,
+    'height=i'     => \my $heightOpt,
+    'ysize=i'      => \my $ysizeOpt,
+    'width=i'      => \my $widthOpt,
+    'xsize=i'      => \my $xsizeOpt,
+    'gamma=i'      => \my $gammaOpt,
+    'randomseed=i' => \my $randomseedOpt);
+
+if (!$validOptions) {
+    print STDERR "Invalid syntax\n";
+    exit(100);
+}
+
+my $pamcraterArgs;
+
+$pamcraterArgs = '';  # initial value
+
+if (defined($numberOpt)) {
+    $pamcraterArgs .= "'-number=$numberOpt' ";
+}
+
+if (defined($heightOpt)) {
+    $pamcraterArgs .= "'-height=$heightOpt' ";
+} elsif (defined($ysizeOpt)) {
+    $pamcraterArgs .= "'-height=$ysizeOpt' ";
+}
+
+if (defined($widthOpt)) {
+    $pamcraterArgs .= "'-width=$widthOpt' ";
+} elsif (defined($xsizeOpt)) {
+    $pamcraterArgs .= "'-width=$xsizeOpt' ";
+}
+
+if (defined($randomseedOpt)) {
+    $pamcraterArgs .= "'-randomseed=$randomseedOpt' ";
+}
+
+my $pamshadedreliefArgs;
+
+$pamshadedreliefArgs = '';  # initial value
+
+if (defined($gammaOpt)) {
+    $pamshadedreliefArgs .= "'-gamma=$gammaOpt' ";
+}
+
+my $termStatus =
+    system(
+        "pamcrater $pamcraterArgs | " .
+        "pamshadedrelief $pamshadedreliefArgs |" .
+        "pamtopnm");
+
+exit($termStatus == 0 ? 0 : 1);
diff --git a/generator/pgmcrater.c b/generator/pgmcrater.c
deleted file mode 100644
index 162e55bb..00000000
--- a/generator/pgmcrater.c
+++ /dev/null
@@ -1,479 +0,0 @@
-/*
-
-              Fractal cratering
-
-       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
-
-    The  algorithm  used  to  determine crater size is as described on
-    pages 31 and 32 of:
-
-    Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
-        Images, New York: Springer Verlag, 1988.
-
-    The  mathematical  technique  used  to calculate crater radii that
-    obey the proper area law distribution from a uniformly distributed
-    pseudorandom sequence was developed by Rudy Rucker.
-
-    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.
-
-                PLUGWARE!
-
-    If you like this kind of stuff, you may also enjoy "James  Gleick's
-    Chaos--The  Software"  for  MS-DOS,  available for $59.95 from your
-    local software store or directly from Autodesk, Inc., Attn: Science
-    Series,  2320  Marinship Way, Sausalito, CA 94965, USA.  Telephone:
-    (800) 688-2344 toll-free or, outside the  U.S. (415)  332-2344  Ext
-    4886.   Fax: (415) 289-4718.  "Chaos--The Software" includes a more
-    comprehensive   fractal    forgery    generator    which    creates
-    three-dimensional  landscapes  as  well as clouds and planets, plus
-    five more modules which explore other aspects of Chaos.   The  user
-    guide  of  more  than  200  pages includes an introduction by James
-    Gleick and detailed explanations by Rudy Rucker of the  mathematics
-    and algorithms used by each program.
-
-*/
-
-/* Modifications by Arjen Bax, 2001-06-21: Remove black vertical line at
-   right edge. Make craters wrap around the image (enables tiling of image).
- */
-
-#define _XOPEN_SOURCE   /* get M_PI in math.h */
-
-#include <assert.h>
-#include <math.h>
-
-#include "pm_c_util.h"
-#include "pgm.h"
-#include "mallocvar.h"
-#include "shhopt.h"
-
-
-struct CmdlineInfo {
-    /* All the information the user supplied in the command line,
-       in a form easy for the program to use.
-    */
-    unsigned int number;
-    unsigned int height;
-    unsigned int width;
-    float        gamma;
-    unsigned int randomseed;
-    unsigned int randomseedSpec;
-    unsigned int test;
-    unsigned int terrain;
-    unsigned int radius;
-
-};
-
-
-
-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 numberSpec, heightSpec, widthSpec, gammaSpec;
-
-    MALLOCARRAY_NOFAIL(option_def, 100);
-
-    option_def_index = 0;   /* incremented by OPTENT3 */
-    OPTENT3(0,   "number",   OPT_UINT,    &cmdlineP->number,
-            &numberSpec,      0);
-    OPTENT3(0,   "height",   OPT_UINT,    &cmdlineP->height,
-            &heightSpec,      0);
-    OPTENT3(0,   "ysize",    OPT_UINT,    &cmdlineP->height,
-            &heightSpec,      0);
-    OPTENT3(0,   "width",    OPT_UINT,    &cmdlineP->width,
-            &widthSpec,       0);
-    OPTENT3(0,   "xsize",    OPT_UINT,    &cmdlineP->width,
-            &widthSpec,       0);
-    OPTENT3(0,   "gamma",    OPT_FLOAT,   &cmdlineP->gamma,
-            &gammaSpec,       0);
-    OPTENT3(0,   "randomseed",  OPT_UINT, &cmdlineP->randomseed,
-            &cmdlineP->randomseedSpec,   0);
-    OPTENT3(0,   "test",     OPT_UINT,   &cmdlineP->radius,
-            &cmdlineP->test,      0);
-    OPTENT3(0,   "terrain",  OPT_FLAG,    NULL,
-            &cmdlineP->terrain,   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 (argc-1 > 0)
-        pm_error("There are no non-option arguments.  You specified %u",
-                 argc-1);
-
-    if (!heightSpec)
-        cmdlineP->height = 256;
-
-    if (cmdlineP->height == 0)
-        pm_error("-height must be positive");
-
-    if (!widthSpec)
-        cmdlineP->width = 256;
-
-    if (cmdlineP->width == 0)
-        pm_error("-width must be positive");
-
-    if (cmdlineP->test) {
-        if (numberSpec || cmdlineP->randomseedSpec)
-            pm_message("Test mode.  Only one fixed crater will be created.  "
-                       "-number and/or -randomseed ignored.");
-
-        if(MAX(cmdlineP->height, cmdlineP->width) * 2 < cmdlineP->radius)
-            pm_error("Radius (%u) too large", cmdlineP->radius);
-    } else {
-        if (!numberSpec)
-            cmdlineP->number = 50000;
-
-        if (cmdlineP->number == 0)
-            pm_error("-number must be positive");
-    }
-
-    if (cmdlineP->terrain) {
-        if (gammaSpec)
-            pm_message("Terrain elevation chart will be output.  "
-                       "-gamma argument (%f) ignored.", cmdlineP->gamma);
-    } else {
-        if (!gammaSpec)
-            cmdlineP->gamma = 1.0;
-
-        if (cmdlineP->gamma <= 0.0)
-            pm_error("gamma correction must be greater than 0");
-    }
-
-    free(option_def);
-}
-
-
-/* Definitions for obtaining random numbers. */
-
-/*  Display parameters  */
-
-static double const ImageGamma = 0.5;     /* Inherent gamma of mapped image */
-static double const arand = 32767.0;      /* Random number parameters */
-static double const CdepthPower = 1.5;    /* Crater depth power factor */
-static double DepthBias2 = 0.5;           /* Square of depth bias */
-static int const slopemin = -52;
-static int const slopemax = 52;
-
-static double const
-Cast(double const high) {
-
-    return high * ((rand() & 0x7FFF) / arand);
-}
-
-
-
-static unsigned int
-mod(int const t,
-    unsigned int const n) {
-
-    /* This is used to transform coordinates beyond bounds into ones
-       within: craters "wrap around" the edges.  This enables tiling
-       of the image.
-
-       Produces strange effects when crater radius is very large compared
-       to image size.
-    */
-
-    int m;
-    m = t % (int) n;
-    if (m < 0)
-        m += n;
-    return m;
-}
-
-
-static void
-generateSlopeMap(gray * const slopemap,
-                 double const dgamma) {
-
-    /* Prepare an array which maps the difference in altitude between two
-       adjacent points (slope) to shades of gray.  Used for output in
-       default (non-terrain) mode.   Uphill slopes are bright; downhill
-       slopes are dark.
-    */ 
-
-    int i;
-    double const gamma = dgamma * ImageGamma;
-
-    for (i = slopemin; i <= 0; i++) {   /* Negative, downhill, dark */
-        slopemap[i - slopemin] =
-            128 - 127.0 * pow(sin((M_PI / 2) * i / slopemin), gamma);
-    }
-    for (i = 0; i <= slopemax; i++) {   /* Positive, uphill, bright */
-        slopemap[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
-slopeToGrayval(int    const slope,
-               gray * const slopemap) {
-
-    return( slopemap[ MIN (MAX (slopemin, slope), slopemax) - slopemin] );
-
-}
-
-
-
-static void
-generateScreenImage(gray **       const terrain,
-                    unsigned int  const width,
-                    unsigned int  const height,
-                    double        const dgamma ) {
-
-    /* Convert a terrain elevation chart into a shaded image and output */
-
-    unsigned int row;
-    gray * const pixrow   = pgm_allocrow(width);  /* output row */
-    gray * const slopemap = pgm_allocrow(slopemax - slopemin +1);
-                            /* Slope to pixel grayval map */
-    gray const maxval = 255;
-
-    pgm_writepgminit(stdout, width, height, maxval, FALSE);
-
-    generateSlopeMap(slopemap, dgamma);
-
-    for (row = 0; row < height; ++row) {
-        unsigned int col;
-
-        for (col = 0; col < width -1; ++col) {
-            int const slope = terrain[row][col+1] - terrain[row][col];
-            pixrow[col] = slopeToGrayval(slope, slopemap);
-        }
-        /* Wrap around to determine shade of pixel on right edge */
-        pixrow[width -1] =
-            slopeToGrayval(terrain[row][0] - terrain[row][width-1], slopemap);
-
-        pgm_writepgmrow(stdout, pixrow, width, maxval, FALSE);
-    }
-
-    pgm_freerow(slopemap);
-    pgm_freerow(pixrow);
-
-}
-
-
-static void
-smallCrater(gray **      const terrain,
-            unsigned int const width,
-            unsigned int const height,
-            int          const cx,
-            int          const cy,
-            double       const g) {
-
-    /* If the crater is tiny, handle it specially. */
-
-    int x, y;
-    unsigned int amptot = 0, axelev;
-    unsigned int npatch = 0;
-
-    /* Set pixel to the average of its Moore neighbourhood. */
-
-    for (y = cy - 1; y <= cy + 1; y++) {
-        for (x = cx - 1; x <= cx + 1; x++) {
-            amptot += terrain[mod(y, height)][mod(x, width)];
-            npatch++;
-        }
-    }
-    axelev = amptot / npatch;
-
-    /* Perturb the mean elevation by a small random factor. */
-
-    x = (g >= 1) ? ((rand() >> 8) & 3) - 1 : 0;
-    terrain[mod(cy, height)][mod(cx, width)] = axelev + x;
-
-}
-
-
-
-static void
-normalCrater(gray **     const terrain,
-            unsigned int const width,
-            unsigned int const height,
-            int          const cx,
-            int          const cy,
-            double       const radius) {
-
-     /* Regular crater.  Generate an impact feature of the
-               correct size and shape. */
-
-    int x, y;
-    unsigned int amptot = 0, axelev;
-    unsigned int npatch = 0;
-    int const impactRadius = (int) MAX(2, (radius / 3));
-    int const craterRadius = (int) radius;
-    double const rollmin = 0.9;
-
-    /* Determine mean elevation around the impact area.
-       We assume the impact area is a fraction of the total crater size. */
-
-    for (y = cy - impactRadius; y <= cy + impactRadius; y++) {
-        for (x = cx - impactRadius; x <= cx + impactRadius; x++) {
-             amptot += terrain[mod(y, height)][mod(x, width)];
-             npatch++;
-        }
-    }
-    axelev = amptot / npatch;
-
-    for (y = cy - craterRadius; y <= cy + craterRadius; y++) {
-        int const dysq = (cy - y) * (cy - y);
-
-        for (x = cx - craterRadius; x <= cx + craterRadius; x++) {
-            int  const dxsq = (cx - x) * (cx - x);
-            double const cd = (dxsq + dysq) /
-                              (double) (craterRadius * craterRadius);
-            double const cd2 = cd * 2.25;
-            double const tcz = sqrt(DepthBias2) - sqrt(fabs(1 - cd2));
-            double cz = MAX((cd2 > 1) ? 0.0 : -10, tcz);  /* Initial value */
-            double roll;
-            unsigned int av;
-
-            cz *= pow(craterRadius, CdepthPower);
-            if (dysq == 0 && dxsq == 0 && ((int) cz) == 0) {
-                cz = cz < 0 ? -1 : 1;
-                }
-
-             roll = (((1 / (1 - MIN(rollmin, cd))) /
-                     (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
-
-             av = (axelev + cz) * (1 - roll) +
-                  (terrain[mod(y, height)][mod(x, width)] + cz) * roll;
-             av = MAX(1000, MIN(64000, av));
-
-             terrain[mod(y, height)][mod(x, width)] = av;
-        }
-    }
-}
-
-
-/* Todo:  We should also have largeCrater() */
-
-
-static void
-plopCrater(gray **         const terrain,
-              unsigned int const width,
-              unsigned int const height,
-              int          const cx,
-              int          const cy,
-              double       const radius) {
-
-        if (radius < 3)
-          smallCrater (terrain, width, height, cx, cy, radius);
-        else
-          normalCrater(terrain, width, height, cx, cy, radius);
-}
-
-
-
-static void
-genCraters(struct CmdlineInfo const cmdline) {
-/*----------------------------------------------------------------------------
-   Generate cratered terrain
------------------------------------------------------------------------------*/
-    unsigned int const width  = cmdline.width;  /* screen X size */
-    unsigned int const height = cmdline.height; /* screen Y size */
-    double       const dgamma = cmdline.gamma;  /* display gamma */
-    gray         const tmaxval = 65535; /* maxval of elevation array */
-
-    gray ** terrain;    /* elevation array */
-    unsigned int row, col;
-
-    /* Acquire the elevation array and initialize it to mean
-       surface elevation. */
-
-    terrain = pgm_allocarray(width, height);
-
-    for (row = 0; row < height; ++row) {
-        for (col = 0; col < width; ++col)
-            terrain[row][col] = tmaxval/2;
-    }
-
-    if ( cmdline.test )
-        plopCrater(terrain, width, height,
-                   width/2, height/2, (double) cmdline.radius);
-
-    else {
-        unsigned int const ncraters = cmdline.number; /* num of craters */
-        unsigned int l;
-
-        for (l = 0; l < ncraters; l++) {
-            int const cx = Cast((double) width  - 1);
-            int const cy = Cast((double) height - 1);
-
-            /* Thanks, Rudy, for this equation  that maps the uniformly
-               distributed  numbers  from   Cast   into   an   area-law
-               distribution as observed on cratered bodies.
-
-               Produces values within the interval:
-               0.56419 <= radius <= 56.419 */
-
-            double const radius = sqrt(1 / (M_PI * (1 - Cast(0.9999))));
-
-            plopCrater(terrain, width, height, cx, cy, radius);
-
-            if (((l + 1) % 5000) == 0)
-                pm_message("%u craters generated of %u (%u%% done)",
-                           l + 1, ncraters, ((l + 1) * 100) / ncraters);
-        }
-    }
-
-    if (cmdline.terrain)
-        pgm_writepgm(stdout, terrain, width, height, tmaxval, FALSE);
-    else
-        generateScreenImage(terrain, width, height, dgamma);
-
-    pgm_freearray(terrain, height);
-    pm_close(stdout);
-}
-
-
-
-int
-main(int argc, const char ** argv) {
-
-    struct CmdlineInfo cmdline;
-
-    pm_proginit(&argc, argv);
-
-    parseCommandLine(argc, argv, &cmdline);
-
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-    genCraters(cmdline);
-
-    exit(0);
-}