Add -terrain, -test

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

1 files changed, 252 insertions, 196 deletions

diff --git a/generator/pgmcrater.c b/generator/pgmcrater.c
index 25b7b817..162e55bb 100644
--- a/generator/pgmcrater.c
+++ b/generator/pgmcrater.c
@@ -45,8 +45,8 @@
 
 */
 
-/* Modifications by Arjen Bax, 2001-06-21: Remove black vertical line at right
- * edge. Make craters wrap around the image (enables tiling of image).
+/* 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 */
@@ -70,6 +70,10 @@ struct CmdlineInfo {
     float        gamma;
     unsigned int randomseed;
     unsigned int randomseedSpec;
+    unsigned int test;
+    unsigned int terrain;
+    unsigned int radius;
+
 };
 
 
@@ -102,10 +106,14 @@ parseCommandLine(int argc, const char ** const argv,
             &widthSpec,       0);
     OPTENT3(0,   "xsize",    OPT_UINT,    &cmdlineP->width,
             &widthSpec,       0);
-    OPTENT3(0,   "gamma",    OPT_FLOAT,    &cmdlineP->gamma,
+    OPTENT3(0,   "gamma",    OPT_FLOAT,   &cmdlineP->gamma,
             &gammaSpec,       0);
-    OPTENT3(0,   "randomseed",   OPT_UINT,    &cmdlineP->randomseed,
-            &cmdlineP->randomseedSpec,      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 */
@@ -118,12 +126,6 @@ parseCommandLine(int argc, const char ** const argv,
         pm_error("There are no non-option arguments.  You specified %u",
                  argc-1);
 
-    if (!numberSpec)
-        cmdlineP->number = 50000;
-
-    if (cmdlineP->number == 0)
-        pm_error("-number must be positive");
-
     if (!heightSpec)
         cmdlineP->height = 256;
 
@@ -136,267 +138,327 @@ parseCommandLine(int argc, const char ** const argv,
     if (cmdlineP->width == 0)
         pm_error("-width must be positive");
 
-    if (!gammaSpec)
-        cmdlineP->gamma = 1.0;
+    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->gamma <= 0.0)
-        pm_error("gamma correction must be greater than 0");
+    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  */
 
-#define SCRX    screenxsize       /* Screen width */
-#define SCRY    screenysize       /* Screen height */
-#define SCRGAMMA 1.0              /* Display gamma */
-
-#define RGBQuant    255
-
-
 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 DepthBias;  /* sqrt(.5) */
-
+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 low,
-     double const high) {
+Cast(double const high) {
 
-    return low + (high - low) * ((rand() & 0x7FFF) / arand);
+    return high * ((rand() & 0x7FFF) / arand);
 }
 
 
 
-static int
-modulo(int const t,
-       int const n) {
+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;
+}
 
-    assert(n > 0);
 
-    m = t % n;
+static void
+generateSlopeMap(gray * const slopemap,
+                 double const dgamma) {
 
-    while (m < 0) {
-        m += n;
+    /* 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);
     }
 
-    return m;
+    /* 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) {
 
-#define Auxadr(x, y) &aux[modulo(y, screenysize)*screenxsize+modulo(x, screenxsize)]
+    return( slopemap[ MIN (MAX (slopemin, slope), slopemax) - slopemin] );
+
+}
 
 
 
 static void
-generateScreenImage(const unsigned short * const aux,
-                    unsigned int           const screenxsize,
-                    unsigned int           const screenysize,
-                    unsigned char *        const slopemap) {
+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 * pixrow;
+    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);
 
-    pgm_writepgminit(stdout, screenxsize, screenysize, RGBQuant, FALSE);
-    pixrow = pgm_allocrow(screenxsize);
+    generateSlopeMap(slopemap, dgamma);
 
-    for (row = 0; row < screenysize; ++row) {
+    for (row = 0; row < height; ++row) {
         unsigned int col;
 
-        for (col = 0; col < screenxsize; ++col) {
-            int j;
-            j = *Auxadr(col+1, row) - *Auxadr(col, row);
-            j = MIN(MAX(slopemin, j), slopemax);
-            pixrow[col] = slopemap[j - slopemin];
+        for (col = 0; col < width -1; ++col) {
+            int const slope = terrain[row][col+1] - terrain[row][col];
+            pixrow[col] = slopeToGrayval(slope, slopemap);
         }
-        pgm_writepgmrow(stdout, pixrow, screenxsize, RGBQuant, FALSE);
+        /* 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);
     }
-    pm_close(stdout);
+
+    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
-gencraters(struct CmdlineInfo const cmdline) {
-/*----------------------------------------------------------------------------
-   Generate cratered terrain
------------------------------------------------------------------------------*/
-    unsigned int const screenxsize = cmdline.width;  /* screen X size */
-    unsigned int const screenysize = cmdline.height; /* screen Y size */
-    double       const dgamma =      cmdline.gamma;  /* display gamma */
-    unsigned int const ncraters =    cmdline.number; /* num craters to gen */
+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 i, j;
-    unsigned int l;
-    unsigned short * aux;
-    unsigned char * slopemap;   /* Slope to pixel map */
+    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;
 
-    /* Acquire the elevation array and initialize it to mean
-       surface elevation. */
+    /* Determine mean elevation around the impact area.
+       We assume the impact area is a fraction of the total crater size. */
 
-    MALLOCARRAY(aux, SCRX * SCRY);
-    if (aux == NULL) 
-        pm_error("out of memory allocating elevation array");
+    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;
+                }
 
-    /* Acquire the elevation buffer and initialize to mean
-       initial elevation. */
+             roll = (((1 / (1 - MIN(rollmin, cd))) /
+                     (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
 
-    for (i = 0; i < SCRY; i++) {
-        unsigned short *zax = aux + (((long) SCRX) * i);
+             av = (axelev + cz) * (1 - roll) +
+                  (terrain[mod(y, height)][mod(x, width)] + cz) * roll;
+             av = MAX(1000, MIN(64000, av));
 
-        for (j = 0; j < SCRX; j++) {
-            *zax++ = 32767;
+             terrain[mod(y, height)][mod(x, width)] = av;
         }
     }
+}
 
-    /* Every time we go around this loop we plop another crater
-       on the surface.  */
 
-    for (l = 0; l < ncraters; l++) {
-        double g;
-        int cx = Cast(0.0, ((double) SCRX - 1)),
-            cy = Cast(0.0, ((double) SCRY - 1)),
-            gx, gy, x, y;
-        unsigned int amptot = 0, axelev;
-        unsigned int npatch = 0;
+/* Todo:  We should also have largeCrater() */
 
 
-        /* Phase 1.  Compute the mean elevation of the impact
-           area.  We assume the impact area is a
-           fraction of the total crater size. */
+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);
+}
 
-        /* Thanks, Rudy, for this equation  that maps the uniformly
-           distributed  numbers  from   Cast   into   an   area-law
-           distribution as observed on cratered bodies. */
 
-        g = sqrt(1 / (M_PI * (1 - Cast(0, 0.9999))));
 
-        /* If the crater is tiny, handle it specially. */
+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 */
 
-        if (g < 3) {
+    gray ** terrain;    /* elevation array */
+    unsigned int row, col;
 
-            /* Set pixel to the average of its Moore neighbourhood. */
+    /* Acquire the elevation array and initialize it to mean
+       surface elevation. */
 
-            for (y = cy - 1; y <= cy + 1; y++) {
-                for (x = cx - 1; x <= cx + 1; x++) {
-                    amptot += *Auxadr(x, y);
-                    npatch++;
-                }
-            }
-            axelev = amptot / npatch;
+    terrain = pgm_allocarray(width, height);
 
-            /* Perturb the mean elevation by a small random factor. */
+    for (row = 0; row < height; ++row) {
+        for (col = 0; col < width; ++col)
+            terrain[row][col] = tmaxval/2;
+    }
 
-            x = (g >= 1) ? ((rand() >> 8) & 3) - 1 : 0;
-            *Auxadr(cx, cy) = axelev + x;
+    if ( cmdline.test )
+        plopCrater(terrain, width, height,
+                   width/2, height/2, (double) cmdline.radius);
 
-            /* Jam repaint sizes to correct patch. */
+    else {
+        unsigned int const ncraters = cmdline.number; /* num of craters */
+        unsigned int l;
 
-            gx = 1;
-            gy = 0;
+        for (l = 0; l < ncraters; l++) {
+            int const cx = Cast((double) width  - 1);
+            int const cy = Cast((double) height - 1);
 
-        } else {
+            /* Thanks, Rudy, for this equation  that maps the uniformly
+               distributed  numbers  from   Cast   into   an   area-law
+               distribution as observed on cratered bodies.
 
-            /* Regular crater.  Generate an impact feature of the
-               correct size and shape. */
+               Produces values within the interval:
+               0.56419 <= radius <= 56.419 */
 
-            /* Determine mean elevation around the impact area. */
+            double const radius = sqrt(1 / (M_PI * (1 - Cast(0.9999))));
 
-            gx = MAX(2, (g / 3));
-            gy = MAX(2, g / 3);
+            plopCrater(terrain, width, height, cx, cy, radius);
 
-            for (y = cy - gy; y <= cy + gy; y++) {
-                for (x = cx-gx; x <= cx + gx; x++) {
-                    amptot += *Auxadr(x,y);
-                    npatch++;
-                }
-            }
-            axelev = amptot / npatch;
-
-            gy = MAX(2, g);
-            g = gy;
-            gx = MAX(2, g);
-
-            for (y = cy - gy; y <= cy + gy; y++) {
-                double dy = (cy - y) / (double) gy,
-                    dysq = dy * dy;
-
-                for (x = cx - gx; x <= cx + gx; x++) {
-                    double dx = ((cx - x) / (double) gx),
-                        cd = (dx * dx) + dysq,
-                        cd2 = cd * 2.25,
-                        tcz = DepthBias - sqrt(fabs(1 - cd2)),
-                        cz = MAX((cd2 > 1) ? 0.0 : -10, tcz),
-                        roll, iroll;
-                    unsigned short av;
-
-                    cz *= pow(g, CdepthPower);
-                    if (dy == 0 && dx == 0 && ((int) cz) == 0) {
-                        cz = cz < 0 ? -1 : 1;
-                    }
-
-#define         rollmin 0.9
-                    roll = (((1 / (1 - MIN(rollmin, cd))) /
-                             (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
-                    iroll = 1 - roll;
-
-                    av = (axelev + cz) * iroll + (*Auxadr(x,y) + cz) * roll;
-                    av = MAX(1000, MIN(64000, av));
-                    *Auxadr(x,y) = av;
-                }
-            }
-        }
-        if ((l % 5000) == 4999) {
-            pm_message( "%u craters generated of %u (%u%% done)",
-                        l + 1, ncraters, ((l + 1) * 100) / ncraters);
+            if (((l + 1) % 5000) == 0)
+                pm_message("%u craters generated of %u (%u%% done)",
+                           l + 1, ncraters, ((l + 1) * 100) / ncraters);
         }
     }
 
-    i = MAX((slopemax - slopemin) + 1, 1);
-    MALLOCARRAY(slopemap, i);
-    if (slopemap == NULL)
-        pm_error("out of memory allocating slope map");
-
-    for (i = slopemin; i <= slopemax; i++) {
-
-        /* 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))             */
-
-        slopemap[i - slopemin] = i > 0 ?
-            (128 + 127.0 *
-             pow(sin((M_PI / 2) * i / slopemax),
-                 dgamma * ImageGamma)) :
-            (128 - 127.0 *
-             pow(sin((M_PI / 2) * i / slopemin),
-                 dgamma * ImageGamma));
-    }
+    if (cmdline.terrain)
+        pgm_writepgm(stdout, terrain, width, height, tmaxval, FALSE);
+    else
+        generateScreenImage(terrain, width, height, dgamma);
 
-    generateScreenImage(aux, screenxsize, screenysize, slopemap);
-
-    free((char *) slopemap);
-    free((char *) aux);
+    pgm_freearray(terrain, height);
+    pm_close(stdout);
 }
 
 
@@ -411,13 +473,7 @@ main(int argc, const char ** argv) {
     parseCommandLine(argc, argv, &cmdline);
 
     srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
-    DepthBias = sqrt(0.5);
-
-    gencraters(cmdline);
+    genCraters(cmdline);
 
     exit(0);
 }
-
-
-