/* 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 #include #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; }; 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); 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 (!numberSpec) cmdlineP->number = 50000; if (cmdlineP->number == 0) pm_error("-number must be positive"); 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 (!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 int const slopemin = -52; static int const slopemax = 52; static double const Cast(double const low, double const high) { return low + (high - low) * ((rand() & 0x7FFF) / arand); } static int modulo(int const t, int const n) { int m; assert(n > 0); m = t % n; while (m < 0) { m += n; } return m; } #define Auxadr(x, y) &aux[modulo(y, screenysize)*screenxsize+modulo(x, screenxsize)] static void generateScreenImage(const unsigned short * const aux, unsigned int const screenxsize, unsigned int const screenysize, unsigned char * const slopemap) { unsigned int row; gray * pixrow; pgm_writepgminit(stdout, screenxsize, screenysize, RGBQuant, FALSE); pixrow = pgm_allocrow(screenxsize); for (row = 0; row < screenysize; ++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]; } pgm_writepgmrow(stdout, pixrow, screenxsize, RGBQuant, FALSE); } pm_close(stdout); pgm_freerow(pixrow); } 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 */ int i, j; unsigned int l; unsigned short * aux; unsigned char * slopemap; /* Slope to pixel map */ /* Acquire the elevation array and initialize it to mean surface elevation. */ MALLOCARRAY(aux, SCRX * SCRY); if (aux == NULL) pm_error("out of memory allocating elevation array"); /* Acquire the elevation buffer and initialize to mean initial elevation. */ for (i = 0; i < SCRY; i++) { unsigned short *zax = aux + (((long) SCRX) * i); for (j = 0; j < SCRX; j++) { *zax++ = 32767; } } /* 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; /* Phase 1. Compute the mean elevation of the impact area. We assume the impact area is a fraction of the total crater size. */ /* 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. */ if (g < 3) { /* 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 += *Auxadr(x, y); npatch++; } } axelev = amptot / npatch; /* Perturb the mean elevation by a small random factor. */ x = (g >= 1) ? ((rand() >> 8) & 3) - 1 : 0; *Auxadr(cx, cy) = axelev + x; /* Jam repaint sizes to correct patch. */ gx = 1; gy = 0; } else { /* Regular crater. Generate an impact feature of the correct size and shape. */ /* Determine mean elevation around the impact area. */ gx = MAX(2, (g / 3)); gy = MAX(2, g / 3); 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); } } 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)); } generateScreenImage(aux, screenxsize, screenysize, slopemap); free((char *) slopemap); free((char *) aux); } 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()); DepthBias = sqrt(0.5); gencraters(cmdline); exit(0); }