From e8c20a06e194739849fe7fcd8413f49c87cda203 Mon Sep 17 00:00:00 2001 From: giraffedata Date: Thu, 24 Jul 2014 01:45:08 +0000 Subject: Add -terrain, -test git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@2227 9d0c8265-081b-0410-96cb-a4ca84ce46f8 --- generator/pgmcrater.c | 448 ++++++++++++++++++++++++++++---------------------- 1 file 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); } - - - -- cgit 1.4.1