about summary refs log tree commit diff
path: root/generator/pgmcrater.c
blob: 1833e60478577a00ffff51058aaa36707620b7a9 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
/*

			  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"

static void gencraters ARGS((void));
static void initseed ARGS((void));

/* Definitions for obtaining random numbers. */

#define Cast(low, high) ((low)+((high)-(low)) * ((rand() & 0x7FFF) / arand))

/*  Data types	*/

typedef int Boolean;
#define FALSE 0
#define TRUE 1

#define V   (void)

/*  Display parameters	*/

#define SCRX	screenxsize	      /* Screen width */
#define SCRY	screenysize	      /* Screen height */
#define SCRGAMMA 1.0		      /* Display gamma */

/*  Local variables  */

#define ImageGamma  0.5 	      /* Inherent gamma of mapped image */

static int screenxsize = 256;	      /* Screen X size */
static int screenysize = 256;	      /* Screen Y size */
static double dgamma = SCRGAMMA;      /* Display gamma */
static double arand = 32767.0;	      /* Random number parameters */
static long ncraters = 50000L;	      /* Number of craters to generate */
static double CdepthPower = 1.5;      /* Crater depth power factor */
static double DepthBias = 0.707107;   /* Depth bias */

static int modulo(int t, int n)
{
    int m;
    assert(n>0);
    m = t % n;
    while (m<0) {
	m+=n;
    }
    return m;
}

/*  INITSEED  --  Generate initial random seed, if needed.  */

static void initseed()
{
    int i;

    i = time(NULL) * 0xF37C;
    srand(i);
    for (i = 0; i < 7; i++) 
        V rand();
}

/*  GENCRATERS	--  Generate cratered terrain.	*/

static void gencraters()
{
    int i, j, x, y;
    long l;
    unsigned short *aux;
    int slopemin = -52, slopemax = 52;
#define RGBQuant    255
    unsigned char *slopemap;   /* Slope to pixel map */
    gray *pixels;	       /* Pixel vector */

#define Auxadr(x, y)  &aux[modulo(y, SCRY)*SCRX+modulo(x, SCRX)]

    /* 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 long 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 0
	fprintf(stderr, "crater=%lu ", (unsigned long)l);
	fprintf(stderr, "cx=%d ", cx);
	fprintf(stderr, "cy=%d ", cy);
	fprintf(stderr, "size=%g\n", g);
#endif

	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( "%ld craters generated of %ld (%ld%% 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));
    }

    /* Generate the screen image. */

    pgm_writepgminit(stdout, SCRX, SCRY, RGBQuant, FALSE);
    pixels = pgm_allocrow(SCRX);

    for (y = 0; y < SCRY; y++) {
	gray *pix = pixels;

	for (x = 0; x < SCRX; x++) {
	    int j = *Auxadr(x+1, y) - *Auxadr(x, y);
	    j = MIN(MAX(slopemin, j), slopemax);
	    *pix++ = slopemap[j - slopemin];
	}
	pgm_writepgmrow(stdout, pixels, SCRX, RGBQuant, FALSE);
    }
    pm_close(stdout);
    pgm_freerow(pixels);

#undef Auxadr
#undef Scradr
    free((char *) slopemap);
    free((char *) aux);
}

/*  MAIN  --  Main program.  */

int main(argc, argv)
  int argc;
  char *argv[];
{
    int i;
    Boolean gammaspec = FALSE, numspec = FALSE,
	    widspec = FALSE, hgtspec = FALSE;
    const char * const usage = "[-number <n>] [-width|-xsize <w>]\n\
                  [-height|-ysize <h>] [-gamma <f>]";

    DepthBias = sqrt(0.5);	      /* Get exact value for depth bias */


    pgm_init(&argc, argv);

    i = 1;
    while ((i < argc) && (argv[i][0] == '-') && (argv[i][1] != '\0')) {
        if (pm_keymatch(argv[i], "-gamma", 2)) {
	    if (gammaspec) {
                pm_error("already specified gamma correction");
	    }
	    i++;
            if ((i == argc) || (sscanf(argv[i], "%lf", &dgamma)  != 1))
		pm_usage(usage);
	    if (dgamma <= 0.0) {
                pm_error("gamma correction must be greater than 0");
	    }
	    gammaspec = TRUE;
        } else if (pm_keymatch(argv[i], "-number", 2)) {
	    if (numspec) {
                pm_error("already specified number of craters");
	    }
	    i++;
            if ((i == argc) || (sscanf(argv[i], "%ld", &ncraters) != 1))
		pm_usage(usage);
	    if (ncraters <= 0) {
                pm_error("number of craters must be greater than 0!");
	    }
	    numspec = TRUE;
        } else if (pm_keymatch(argv[i], "-xsize", 2) ||
                   pm_keymatch(argv[i], "-width", 2)) {
	    if (widspec) {
                pm_error("already specified a width/xsize");
	    }
	    i++;
            if ((i == argc) || (sscanf(argv[i], "%d", &screenxsize) != 1))
		pm_usage(usage);
	    if (screenxsize <= 0) {
                pm_error("screen width must be greater than 0");
	    }
	    widspec = TRUE;
        } else if (pm_keymatch(argv[i], "-ysize", 2) ||
                   pm_keymatch(argv[i], "-height", 2)) {
	    if (hgtspec) {
                pm_error("already specified a height/ysize");
	    }
	    i++;
            if ((i == argc) || (sscanf(argv[i], "%d", &screenysize) != 1))
		pm_usage(usage);
	    if (screenxsize <= 0) {
                pm_error("screen height must be greater than 0");
	    }
	    hgtspec = TRUE;
	} else {
	    pm_usage(usage);
	}
	i++;
    }

    initseed();
    gencraters();

    exit(0);
}