home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsm / netpbmsca / pgm / c / pgmcrater < prev    next >
Encoding:
Text File  |  1993-10-04  |  10.2 KB  |  372 lines
  1. /*
  2.  
  3.               Fractal cratering
  4.  
  5.        Designed and implemented in November of 1989 by:
  6.  
  7.         John Walker
  8.         Autodesk SA
  9.         Avenue des Champs-Montants 14b
  10.         CH-2074 MARIN
  11.         Switzerland
  12.         Usenet: kelvin@Autodesk.com
  13.         Fax:    038/33 88 15
  14.         Voice:  038/33 76 33
  15.  
  16.     The  algorithm  used  to  determine crater size is as described on
  17.     pages 31 and 32 of:
  18.  
  19.     Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
  20.         Images, New York: Springer Verlag, 1988.
  21.  
  22.     The  mathematical  technique  used    to calculate crater radii that
  23.     obey the proper area law distribution from a uniformly distributed
  24.     pseudorandom sequence was developed by Rudy Rucker.
  25.  
  26.     Permission    to  use, copy, modify, and distribute this software and
  27.     its documentation  for  any  purpose  and  without    fee  is  hereby
  28.     granted,  without any conditions or restrictions.  This software is
  29.     provided "as is" without express or implied warranty.
  30.  
  31.                 PLUGWARE!
  32.  
  33.     If you like this kind of stuff, you may also enjoy "James  Gleick's
  34.     Chaos--The  Software"  for  MS-DOS,  available for $59.95 from your
  35.     local software store or directly from Autodesk, Inc., Attn: Science
  36.     Series,  2320  Marinship Way, Sausalito, CA 94965, USA.  Telephone:
  37.     (800) 688-2344 toll-free or, outside the  U.S. (415)  332-2344  Ext
  38.     4886.   Fax: (415) 289-4718.  "Chaos--The Software" includes a more
  39.     comprehensive   fractal    forgery      generator    which    creates
  40.     three-dimensional  landscapes  as  well as clouds and planets, plus
  41.     five more modules which explore other aspects of Chaos.   The  user
  42.     guide  of  more  than  200    pages includes an introduction by James
  43.     Gleick and detailed explanations by Rudy Rucker of the  mathematics
  44.     and algorithms used by each program.
  45.  
  46. */
  47.  
  48. #include <math.h>
  49. #include "pgm.h"
  50.  
  51. static void gencraters ARGS((void));
  52. static void initseed ARGS((void));
  53.  
  54. #ifndef M_PI
  55. #define M_PI    3.14159265358979323846
  56. #endif
  57.  
  58. /* Definitions for obtaining random numbers. */
  59.  
  60. #define Cast(low, high) ((low)+((high)-(low)) * ((random() & 0x7FFF) / arand))
  61.  
  62. /*  Data types    */
  63.  
  64. typedef int Boolean;
  65. #define FALSE 0
  66. #define TRUE 1
  67.  
  68. #define V   (void)
  69.  
  70. /*  Display parameters    */
  71.  
  72. #define SCRX    screenxsize          /* Screen width */
  73. #define SCRY    screenysize          /* Screen height */
  74. #define SCRGAMMA 1.0              /* Display gamma */
  75.  
  76. /*  Local variables  */
  77.  
  78. #define ImageGamma  0.5           /* Inherent gamma of mapped image */
  79.  
  80. static int screenxsize = 256;          /* Screen X size */
  81. static int screenysize = 256;          /* Screen Y size */
  82. static double dgamma = SCRGAMMA;      /* Display gamma */
  83. static double arand = 32767.0;          /* Random number parameters */
  84. static long ncraters = 50000L;          /* Number of craters to generate */
  85. static double CdepthPower = 1.5;      /* Crater depth power factor */
  86. static double DepthBias = 0.707107;   /* Depth bias */
  87.  
  88. /*  INITSEED  --  Generate initial random seed, if needed.  */
  89.  
  90. static void initseed()
  91. {
  92.     int i;
  93.  
  94.     i = time((long *) 0) * 0xF37C;
  95.     srandom(i);
  96.     for (i = 0; i < 7; i++) {
  97.     V random();
  98.     }
  99. }
  100.  
  101. /*  GENCRATERS    --  Generate cratered terrain.    */
  102.  
  103. static void gencraters()
  104. {
  105.     int i, j, x, y;
  106.     long l;
  107.     unsigned short *aux;
  108.     int slopemin = -52, slopemax = 52;
  109. #define RGBQuant    255
  110.     unsigned char *slopemap;   /* Slope to pixel map */
  111.     gray *pixels;           /* Pixel vector */
  112.  
  113. #define Auxadr(x, y)  ((unsigned short *) (aux + ((((y)) * SCRX) + (x))))
  114.  
  115.     /* Acquire the elevation array and initialise it to mean
  116.        surface elevation. */
  117.  
  118.     aux = (unsigned short *) malloc(SCRX * SCRY * sizeof(short));
  119.     if (aux == (unsigned short *) 0) {
  120.         pm_error("out of memory allocating elevation array");
  121.     }
  122.  
  123.     /* Acquire the elevation buffer and initialise to mean
  124.        initial elevation. */
  125.  
  126.     for (i = 0; i < SCRY; i++) {
  127.     unsigned short *zax = aux + (((long) SCRX) * i);
  128.  
  129.     for (j = 0; j < SCRX; j++) {
  130.         *zax++ = 32767;
  131.     }
  132.     }
  133.  
  134.     /* Every time we go around this loop we plop another crater
  135.        on the surface.    */
  136.  
  137.     for (l = 0; l < ncraters; l++) {
  138.     double g;
  139.     int cx = Cast(0.0, ((double) SCRX - 1)),
  140.         cy = Cast(0.0, ((double) SCRY - 1)),
  141.         gx, gy, x, y;
  142.     unsigned long amptot = 0, axelev;
  143.     unsigned int npatch = 0;
  144.  
  145.  
  146.     /* Phase 1.  Compute the mean elevation of the impact
  147.              area.  We assume the impact area is a
  148.              fraction of the total crater size. */
  149.  
  150.     /* Thanks, Rudy, for this equation  that maps the uniformly
  151.        distributed    numbers  from    Cast   into   an   area-law
  152.        distribution as observed on cratered bodies. */
  153.  
  154.     g = sqrt(1 / (M_PI * (1 - Cast(0, 0.9999))));
  155.  
  156.     /* If the crater is tiny, handle it specially. */
  157.  
  158.     if (g < 3) {
  159.  
  160.        /* Set pixel to the average of its Moore neighbourhood. */
  161.  
  162.         for (y = max(0, cy - 1); y <= min(SCRY - 1, cy + 1); y++) {
  163.         int sx = max(0, cx - 1);
  164.         unsigned short *a = Auxadr(sx, y);
  165.  
  166.         for (x = sx; x <= min(SCRX - 1, cx + 1); x++) {
  167.             amptot += *a++;
  168.             npatch++;
  169.         }
  170.         }
  171.         axelev = amptot / npatch;
  172.  
  173.         /* Perturb the mean elevation by a small random factor. */
  174.  
  175.         x = (g >= 1) ? ((random() >> 8) & 3) - 1 : 0;
  176.         *Auxadr(cx, cy) = axelev + x;
  177.  
  178.         /* Jam repaint sizes to correct patch. */
  179.  
  180.         gx = 1;
  181.         gy = 0;
  182.  
  183.     } else {
  184.  
  185.         /* Regular crater.    Generate an impact feature of the
  186.            correct size and shape. */
  187.  
  188.         /* Determine mean elevation around the impact area. */
  189.  
  190.         gx = max(2, (g / 3));
  191.         gy = max(2, g / 3);
  192.  
  193.         for (y = max(0, cy - gy); y <= min(SCRY - 1, cy + gy); y++) {
  194.         int sx = max(0, cx - gx);
  195.         unsigned short *a = Auxadr(sx, y);
  196.  
  197.         for (x = sx; x <= min(SCRX - 1, cx + gx); x++) {
  198.             amptot += *a++;
  199.             npatch++;
  200.         }
  201.         }
  202.         axelev = amptot / npatch;
  203.  
  204.         gy = max(2, g);
  205.         g = gy;
  206.         gx = max(2, g);
  207.  
  208.         for (y = max(0, cy - gy); y <= min(SCRY - 1, cy + gy); y++) {
  209.         int sx = max(0, cx - gx);
  210.         unsigned short *ax = Auxadr(sx, y);
  211.         double dy = (cy - y) / (double) gy,
  212.                dysq = dy * dy;
  213.  
  214.         for (x = sx; x <= min(SCRX - 1, cx + gx); x++) {
  215.             double dx = ((cx - x) / (double) gx),
  216.                cd = (dx * dx) + dysq,
  217.                cd2 = cd * 2.25,
  218.                tcz = DepthBias - sqrt(fabs(1 - cd2)),
  219.                cz = max((cd2 > 1) ? 0.0 : -10, tcz),
  220.                roll, iroll;
  221.             unsigned short av;
  222.  
  223.             cz *= pow(g, CdepthPower);
  224.             if (dy == 0 && dx == 0 && ((int) cz) == 0) {
  225.                cz = cz < 0 ? -1 : 1;
  226.             }
  227.  
  228. #define         rollmin 0.9
  229.             roll = (((1 / (1 - min(rollmin, cd))) /
  230.                  (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
  231.             iroll = 1 - roll;
  232.  
  233.             av = (axelev + cz) * iroll + (*ax + cz) * roll;
  234.             av = max(1000, min(64000, av));
  235.             *ax++ = av;
  236.         }
  237.         }
  238.      }
  239.     if ((l % 5000) == 4999) {
  240.         pm_message( "%ld craters generated of %ld (%ld%% done)",
  241.         l + 1, ncraters, ((l + 1) * 100) / ncraters);
  242.     }
  243.     }
  244.  
  245.     i = max((slopemax - slopemin) + 1, 1);
  246.     slopemap = (unsigned char *) malloc(i * sizeof(unsigned char));
  247.     if (slopemap == (unsigned char *) 0) {
  248.         pm_error("out of memory allocating slope map");
  249.     }
  250.     for (i = slopemin; i <= slopemax; i++) {
  251.  
  252.         /* Confused?   OK,  we're using the  left-to-right slope to
  253.        calculate a shade based on the sine of  the    angle  with
  254.        respect  to the vertical (light incident from the left).
  255.        Then, with one exponentiation, we account for  both    the
  256.        inherent   gamma   of   the     image    (ad-hoc),  and    the
  257.        user-specified display gamma, using the identity:
  258.  
  259.          (x^y)^z = (x^(y*z))             */
  260.  
  261.     slopemap[i - slopemin] = i > 0 ?
  262.         (128 + 127.0 *
  263.         pow(sin((M_PI / 2) * i / slopemax),
  264.                dgamma * ImageGamma)) :
  265.         (128 - 127.0 *
  266.         pow(sin((M_PI / 2) * i / slopemin),
  267.                dgamma * ImageGamma));
  268.     }
  269.  
  270.     /* Generate the screen image. */
  271.  
  272.     pgm_writepgminit(stdout, SCRX, SCRY, RGBQuant, FALSE);
  273.     pixels = pgm_allocrow(SCRX);
  274.  
  275.     for (y = 0; y < SCRY; y++) {
  276.     unsigned short *ax = Auxadr(0, y);
  277.     gray *pix = pixels;
  278.  
  279.     for (x = 0; x < SCRX - 1; x++) {
  280.         int j = ax[1] - ax[0];
  281.         j = min(max(slopemin, j), slopemax);
  282.         *pix++ = slopemap[j - slopemin];
  283.         ax++;
  284.     }
  285.     pgm_writepgmrow(stdout, pixels, SCRX, RGBQuant, FALSE);
  286.     }
  287.     pm_close(stdout);
  288.     pgm_freerow(pixels);
  289.  
  290. #undef Auxadr
  291. #undef Scradr
  292.     free((char *) slopemap);
  293.     free((char *) aux);
  294. }
  295.  
  296. /*  MAIN  --  Main program.  */
  297.  
  298. int main(argc, argv)
  299.   int argc;
  300.   char *argv[];
  301. {
  302.     int i;
  303.     Boolean gammaspec = FALSE, numspec = FALSE,
  304.         widspec = FALSE, hgtspec = FALSE;
  305.     char *usage = "[-number <n>] [-width|-xsize <w>]\n\
  306.                   [-height|-ysize <h>] [-gamma <f>]";
  307.  
  308.     DepthBias = sqrt(0.5);          /* Get exact value for depth bias */
  309.  
  310.  
  311.     pgm_init(&argc, argv);
  312.  
  313.     i = 1;
  314.     while ((i < argc) && (argv[i][0] == '-') && (argv[i][1] != '\0')) {
  315.         if (pm_keymatch(argv[i], "-gamma", 2)) {
  316.         if (gammaspec) {
  317.                 pm_error("already specified gamma correction");
  318.         }
  319.         i++;
  320.             if ((i == argc) || (sscanf(argv[i], "%lf", &dgamma)  != 1))
  321.         pm_usage(usage);
  322.         if (dgamma <= 0.0) {
  323.                 pm_error("gamma correction must be greater than 0");
  324.         }
  325.         gammaspec = TRUE;
  326.         } else if (pm_keymatch(argv[i], "-number", 2)) {
  327.         if (numspec) {
  328.                 pm_error("already specified number of craters");
  329.         }
  330.         i++;
  331.             if ((i == argc) || (sscanf(argv[i], "%ld", &ncraters) != 1))
  332.         pm_usage(usage);
  333.         if (ncraters <= 0) {
  334.                 pm_error("number of craters must be greater than 0!");
  335.         }
  336.         numspec = TRUE;
  337.         } else if (pm_keymatch(argv[i], "-xsize", 2) ||
  338.                    pm_keymatch(argv[i], "-width", 2)) {
  339.         if (widspec) {
  340.                 pm_error("already specified a width/xsize");
  341.         }
  342.         i++;
  343.             if ((i == argc) || (sscanf(argv[i], "%d", &screenxsize) != 1))
  344.         pm_usage(usage);
  345.         if (screenxsize <= 0) {
  346.                 pm_error("screen width must be greater than 0");
  347.         }
  348.         widspec = TRUE;
  349.         } else if (pm_keymatch(argv[i], "-ysize", 2) ||
  350.                    pm_keymatch(argv[i], "-height", 2)) {
  351.         if (hgtspec) {
  352.                 pm_error("already specified a height/ysize");
  353.         }
  354.         i++;
  355.             if ((i == argc) || (sscanf(argv[i], "%d", &screenysize) != 1))
  356.         pm_usage(usage);
  357.         if (screenxsize <= 0) {
  358.                 pm_error("screen height must be greater than 0");
  359.         }
  360.         hgtspec = TRUE;
  361.     } else {
  362.         pm_usage(usage);
  363.     }
  364.     i++;
  365.     }
  366.  
  367.     initseed();
  368.     gencraters();
  369.  
  370.     exit(0);
  371. }
  372.