1 /*
2 * This file is part of XPaint.
3 *
4 * Copyright 1996 by Torsten Martinsen
5 * Copyright 1994 by Jer Johnson
6 *
7 * This is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 */
12
13 #include <X11/X.h>
14 #include <X11/Intrinsic.h>
15 #include <X11/Xlib.h>
16
17 #include <stdio.h>
18 #include <math.h>
19 #include <stdlib.h>
20
21 #ifdef __EMX__
22 #include <float.h>
23 #endif
24
25 #include "xpaint.h"
26 #include "misc.h"
27 #include "image.h"
28
29 #if FEATURE_FRACTAL
30
31 static void draw_grid(int x, int y, unsigned char col, unsigned char *grid);
32 static void subdiv(int x1, int y1, int x2, int y2, unsigned char *grid);
33 static void adjust(int x1, int y1, int x, int y, int x2, int y2,
34 unsigned char *grid);
35 static void fourn(float data[], int nn[], int ndim, int isign);
36 static void spectralsynth(float **x, unsigned int n, double h);
37
38 static int number, numcolors, width, height;
39 extern char fractalDensityStr[10];
40
41 /*
42 * draw_plasma() and friends are stolen from Xtacy, written by Jer Johnson
43 * (mpython@gnu.ai.mit.edu) who says that his work in turn is based on OOZE.C,
44 * written by Jeff Clough.
45 */
46 Image *
draw_plasma(int w,int h)47 draw_plasma(int w, int h)
48 {
49 int i, n;
50 unsigned char *grid;
51 Image *output;
52 unsigned char *op;
53
54
55 /* customizable parameters */
56 numcolors = 100;
57 number = 10;
58
59 width = w;
60
61 n = numcolors / 5;
62 numcolors = n * 5;
63 output = ImageNewCmap(w, h, numcolors);
64 grid = output->data;
65 memset(grid, 0, w * h);
66
67 /*
68 * Make points at (0,0) (0, h-1) (w-1, 0) (w-1, h-1) of random colors
69 */
70 draw_grid(0, 0, RANDOMI2(0, numcolors-1), grid);
71 draw_grid(w - 1, 0, RANDOMI2(0, numcolors-1), grid);
72 draw_grid(0, h - 1, RANDOMI2(0, numcolors-1), grid);
73 draw_grid(w - 1, h - 1, RANDOMI2(0, numcolors-1), grid);
74
75 subdiv(0, 0, w - 1, h - 1, grid);
76
77 /* Build the colour map */
78 op = output->cmapData;
79 for (i = 0; i < n; ++i) {
80 *op++ = 255;
81 *op++ = 255 * i / n;
82 *op++ = 0;
83 }
84 for (i = 0; i < n; ++i) {
85 *op++ = 255 - 255 * i / n;
86 *op++ = 255;
87 *op++ = 0;
88 }
89 for (i = 0; i < n; ++i) {
90 *op++ = 0;
91 *op++ = 255;
92 *op++ = 255 * i / n;
93 }
94 for (i = 0; i < n; ++i) {
95 *op++ = 255 * i / n;
96 *op++ = 0;
97 *op++ = 255;
98 }
99 for (i = 0; i < n; ++i) {
100 *op++ = 255;
101 *op++ = 0;
102 *op++ = 255 - 255 * i / n;
103 }
104 #if 0
105 op = output->cmapData;
106 for (i = 0; i < 5 * n; ++i) {
107 printf("%d: %d %d %d\n", i, op[0], op[1], op[2]);
108 op += 3;
109 }
110 #endif
111
112 return output;
113 }
114
115 static void
draw_grid(int x,int y,unsigned char col,unsigned char * grid)116 draw_grid(int x, int y, unsigned char col, unsigned char *grid)
117 {
118 grid[x + y * width] = col;
119 }
120
121 static void
subdiv(int x1,int y1,int x2,int y2,unsigned char * grid)122 subdiv(int x1, int y1, int x2, int y2, unsigned char *grid)
123 {
124 int x, y;
125
126 if ((x2 - x1 < 2) && (y2 - y1 < 2))
127 return;
128
129 x = (x1 + x2) / 2;
130 y = (y1 + y2) / 2;
131 adjust(x1, y1, x, y1, x2, y1, grid);
132 adjust(x2, y1, x2, y, x2, y2, grid);
133 adjust(x1, y2, x, y2, x2, y2, grid);
134 adjust(x1, y1, x1, y, x1, y2, grid);
135
136 if (grid[x + y * width] == 0) {
137 unsigned char spooge;
138
139 spooge = ((int) grid[x1 + y1 * width] +
140 (int) grid[x1 + y2 * width] +
141 (int) grid[x2 + y1 * width] +
142 (int) grid[x2 + y2 * width]) / 4;
143 draw_grid(x, y, spooge, grid);
144 }
145 subdiv(x1, y1, x, y, grid);
146 subdiv(x, y1, x2, y, grid);
147 subdiv(x1, y, x, y2, grid);
148 subdiv(x, y, x2, y2, grid);
149 }
150
151 static void
adjust(int x1,int y1,int x,int y,int x2,int y2,unsigned char * grid)152 adjust(int x1, int y1, int x, int y, int x2, int y2, unsigned char *grid)
153 {
154 int c, d, spoo;
155 int clr1, clr2, clr3, clr4, horz, vert;
156
157
158 if (grid[x + y * width] != 0)
159 return;
160
161 clr1 = grid[x1 + y1 * width];
162 clr2 = grid[x2 + y2 * width];
163 clr3 = grid[x1 + y2 * width];
164 clr4 = grid[x2 + y1 * width];
165 horz = x2 - x1;
166 vert = y2 - y1;
167 d = hypot(horz, vert);
168
169 #if 1 /* 0 for band effect */
170 if ((spoo = RANDOMI2(0, 100)) > number) {
171 if (spoo % 2 == 1)
172 c = (((clr1 + clr2) / 2) + ((int) RANDOMI2(0, d))) % numcolors;
173 else
174 c = (((clr1 + clr2) / 2) - ((int) RANDOMI2(0, d))) % numcolors;
175 } else
176 #endif
177 c = (clr1 + clr2 + clr3 + clr4) / 4;
178 if (c < 0)
179 c = -c;
180 else if (c == 0)
181 c = 1;
182 draw_grid(x, y, c, grid);
183 }
184
185 /*
186 * Fractal landscape generator from ppmforge.c,
187 * written by John Walker (kelvin@Autodesk.com).
188 */
189
190 #define Real(v, x, y) v[1 + (((x) * meshsize) + (y)) * 2]
191 #define Imag(v, x, y) v[2 + (((x) * meshsize) + (y)) * 2]
192
193 #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
194
195 /* customize: */
196 static double fracdim = 1.8; /* Fractal dimension - ]0,3] */
197 static double powscale = 1.0; /* Power law scaling exponent - ]0,3] */
198 static int meshsize = 256; /* FFT mesh size - generally >= max(w,h),
199 but small numbers can give special effects */
200 /* colours */
201
202 static unsigned char pgnd[][3] =
203 {
204 {206, 205, 0},
205 {208, 207, 0},
206 {211, 208, 0},
207 {214, 208, 0},
208 {217, 208, 0},
209 {220, 208, 0},
210 {222, 207, 0},
211 {225, 205, 0},
212 {227, 204, 0},
213 {229, 202, 0},
214 {231, 199, 0},
215 {232, 197, 0},
216 {233, 194, 0},
217 {234, 191, 0},
218 {234, 188, 0},
219 {233, 185, 0},
220 {232, 183, 0},
221 {231, 180, 0},
222 {229, 178, 0},
223 {227, 176, 0},
224 {225, 174, 0},
225 {223, 172, 0},
226 {221, 170, 0},
227 {219, 168, 0},
228 {216, 166, 0},
229 {214, 164, 0},
230 {212, 162, 0},
231 {210, 161, 0},
232 {207, 159, 0},
233 {205, 157, 0},
234 {203, 156, 0},
235 {200, 154, 0},
236 {198, 152, 0},
237 {195, 151, 0},
238 {193, 149, 0},
239 {190, 148, 0},
240 {188, 147, 0},
241 {185, 145, 0},
242 {183, 144, 0},
243 {180, 143, 0},
244 {177, 141, 0},
245 {175, 140, 0},
246 {172, 139, 0},
247 {169, 138, 0},
248 {167, 137, 0},
249 {164, 136, 0},
250 {161, 135, 0},
251 {158, 134, 0},
252 {156, 133, 0},
253 {153, 132, 0},
254 {150, 132, 0},
255 {147, 131, 0},
256 {145, 130, 0},
257 {142, 130, 0},
258 {139, 129, 0},
259 {136, 128, 0},
260 {133, 128, 0},
261 {130, 127, 0},
262 {127, 127, 0},
263 {125, 127, 0},
264 {122, 127, 0},
265 {119, 127, 0},
266 {116, 127, 0},
267 {113, 127, 0},
268 {110, 128, 0},
269 {107, 128, 0},
270 {104, 128, 0},
271 {102, 127, 0},
272 {99, 126, 0},
273 {97, 124, 0},
274 {95, 122, 0},
275 {93, 120, 0},
276 {92, 117, 0},
277 {92, 114, 0},
278 {92, 111, 0},
279 {93, 108, 0},
280 {94, 106, 0},
281 {96, 104, 0},
282 {98, 102, 0},
283 {100, 100, 0},
284 {103, 99, 0},
285 {106, 99, 0},
286 {109, 99, 0},
287 {111, 100, 0},
288 {114, 101, 0},
289 {117, 102, 0},
290 {120, 103, 0},
291 {123, 102, 0},
292 {125, 102, 0},
293 {128, 100, 0},
294 {130, 98, 0},
295 {132, 96, 0},
296 {133, 94, 0},
297 {134, 91, 0},
298 {134, 88, 0},
299 {134, 85, 0},
300 {133, 82, 0},
301 {131, 80, 0},
302 {129, 78, 0}
303 };
304
305
306 /* FOURN -- Multi-dimensional fast Fourier transform
307
308 Called with arguments:
309
310 data A one-dimensional array of floats (NOTE!!! NOT
311 DOUBLES!!), indexed from one (NOTE!!! NOT ZERO!!),
312 containing pairs of numbers representing the complex
313 valued samples. The Fourier transformed results are
314 returned in the same array.
315
316 nn An array specifying the edge size in each dimension.
317 THIS ARRAY IS INDEXED FROM ONE, AND ALL THE EDGE
318 SIZES MUST BE POWERS OF TWO!!!
319
320 ndim Number of dimensions of FFT to perform. Set to 2 for
321 two dimensional FFT.
322
323 isign If 1, a Fourier transform is done; if -1 the inverse
324 transformation is performed.
325
326 This function is essentially as given in Press et al., "Numerical
327 Recipes In C", Section 12.11, pp. 467-470.
328 */
329
330 static void
fourn(float data[],int nn[],int ndim,int isign)331 fourn(float data[], int nn[], int ndim, int isign)
332 {
333 register int i1, i2, i3;
334 int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
335 int ibit, idim, k1, k2, n, nprev, nrem, ntot;
336 float tempi, tempr;
337 double theta, wi, wpi, wpr, wr, wtemp;
338
339 #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr
340
341 ntot = 1;
342 for (idim = 1; idim <= ndim; idim++)
343 ntot *= nn[idim];
344 nprev = 1;
345 for (idim = ndim; idim >= 1; idim--) {
346 n = nn[idim];
347 nrem = ntot / (n * nprev);
348 ip1 = nprev << 1;
349 ip2 = ip1 * n;
350 ip3 = ip2 * nrem;
351 i2rev = 1;
352 for (i2 = 1; i2 <= ip2; i2 += ip1) {
353 if (i2 < i2rev) {
354 for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
355 for (i3 = i1; i3 <= ip3; i3 += ip2) {
356 i3rev = i2rev + i3 - i2;
357 SWAP(data[i3], data[i3rev]);
358 SWAP(data[i3 + 1], data[i3rev + 1]);
359 }
360 }
361 }
362 ibit = ip2 >> 1;
363 while (ibit >= ip1 && i2rev > ibit) {
364 i2rev -= ibit;
365 ibit >>= 1;
366 }
367 i2rev += ibit;
368 }
369 ifp1 = ip1;
370 while (ifp1 < ip2) {
371 ifp2 = ifp1 << 1;
372 theta = isign * (M_PI * 2) / (ifp2 / ip1);
373 wtemp = sin(0.5 * theta);
374 wpr = -2.0 * wtemp * wtemp;
375 wpi = sin(theta);
376 wr = 1.0;
377 wi = 0.0;
378 for (i3 = 1; i3 <= ifp1; i3 += ip1) {
379 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
380 for (i2 = i1; i2 <= ip3; i2 += ifp2) {
381 k1 = i2;
382 k2 = k1 + ifp1;
383 tempr = wr * data[k2] - wi * data[k2 + 1];
384 tempi = wr * data[k2 + 1] + wi * data[k2];
385 data[k2] = data[k1] - tempr;
386 data[k2 + 1] = data[k1 + 1] - tempi;
387 data[k1] += tempr;
388 data[k1 + 1] += tempi;
389 }
390 }
391 wr = (wtemp = wr) * wpr - wi * wpi + wr;
392 wi = wi * wpr + wtemp * wpi + wi;
393 }
394 ifp1 = ifp2;
395 }
396 nprev *= n;
397 }
398 }
399 #undef SWAP
400
401 /*
402 * SPECTRALSYNTH -- Spectrally synthesised fractal motion in two dimensions.
403 */
404 static void
spectralsynth(float ** x,unsigned int n,double h)405 spectralsynth(float **x, unsigned int n, double h)
406 {
407 unsigned bl;
408 int i, j, i0, j0, nsize[3];
409 double rad, phase, rcos, rsin;
410 float *a;
411
412 bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float);
413 a = (float *) XtCalloc(bl, 1);
414 *x = a;
415
416 for (i = 0; i <= n / 2; i++) {
417 for (j = 0; j <= n / 2; j++) {
418 phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0;
419 if (i != 0 || j != 0) {
420 rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
421 } else {
422 rad = 0;
423 }
424 rcos = rad * cos(phase);
425 rsin = rad * sin(phase);
426 Real(a, i, j) = rcos;
427 Imag(a, i, j) = rsin;
428 i0 = (i == 0) ? 0 : n - i;
429 j0 = (j == 0) ? 0 : n - j;
430 Real(a, i0, j0) = rcos;
431 Imag(a, i0, j0) = -rsin;
432 }
433 }
434 Imag(a, n / 2, 0) = 0;
435 Imag(a, 0, n / 2) = 0;
436 Imag(a, n / 2, n / 2) = 0;
437 for (i = 1; i <= n / 2 - 1; i++) {
438 for (j = 1; j <= n / 2 - 1; j++) {
439 phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0;
440 rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
441 rcos = rad * cos(phase);
442 rsin = rad * sin(phase);
443 Real(a, i, n - j) = rcos;
444 Imag(a, i, n - j) = rsin;
445 Real(a, n - i, j) = rcos;
446 Imag(a, n - i, j) = -rsin;
447 }
448 }
449
450 nsize[0] = 0;
451 nsize[1] = nsize[2] = n; /* Dimension of frequency domain array */
452 fourn(a, nsize, 2, -1); /* Take inverse 2D Fourier transform */
453 }
454
455 Image *
draw_landscape(int w,int h,int clouds)456 draw_landscape(int w, int h, int clouds)
457 {
458 float *a;
459 int density = atoi(fractalDensityStr);
460 double rmin = HUGE_VAL, rmax = -HUGE_VAL, rmean, rrange;
461 Image *output = ImageNew(w, h);
462 unsigned char *op = output->data;
463 int x, y, i, j;
464 unsigned char *cp, *ap;
465 double *u, *u1;
466 unsigned int *bxf, *bxc;
467
468 width = w;
469 height = h;
470 spectralsynth(&a, meshsize, 3.0 - fracdim);
471
472 /* Apply power law scaling if non-unity scale is requested. */
473
474 if (powscale != 1.0)
475 for (i = 0; i < meshsize; i++)
476 for (j = 0; j < meshsize; j++) {
477 double r = Real(a, i, j);
478
479 if (r > 0)
480 Real(a, i, j) = pow(r, powscale);
481 }
482 /* Compute extrema for autoscaling. */
483 for (i = 0; i < meshsize; i++)
484 for (j = 0; j < meshsize; j++) {
485 double r = Real(a, i, j);
486
487 rmin = MIN(rmin, r);
488 rmax = MAX(rmax, r);
489 }
490
491 rmean = (rmin + rmax) / 2;
492 rrange = (rmax - rmin) / 2;
493 for (i = 0; i < meshsize; i++)
494 for (j = 0; j < meshsize; j++)
495 Real(a, i, j) = (Real(a, i, j) - rmean) / rrange;
496
497 u = (double *) XtMalloc((unsigned int) (width * sizeof(double)));
498 u1 = (double *) XtMalloc((unsigned int) (width * sizeof(double)));
499 bxf = (unsigned int *) XtMalloc((unsigned int) width *
500 sizeof(unsigned int));
501 bxc = (unsigned int *) XtMalloc((unsigned int) width *
502 sizeof(unsigned int));
503
504 /* Prescale the grid points into intensities. */
505
506 cp = (unsigned char *) XtMalloc(meshsize * meshsize + 1);
507 ap = cp;
508 for (i = 0; i < meshsize; i++)
509 for (j = 0; j < meshsize; j++)
510 *ap++ = (255.0 * (Real(a, i, j) + 1.0)) / 2.0;
511 /* Hack: make an extra array entry */
512 *ap = ap[-1];
513
514 /* Fill the screen from the computed intensity grid by mapping screen points
515 onto the grid, then calculating each pixel value by bilinear
516 interpolation from the surrounding grid points. (N.b. the pictures would
517 undoubtedly look better when generated with small grids if a more
518 well-behaved interpolation were used.)
519
520 Before we get started, precompute the line-level interpolation
521 parameters and store them in an array so we don't have to do
522 this every time around the inner loop. */
523
524 #define UPRJ(a,size) (((a))? (a)/((size)-1.0) : 0)
525
526 for (x = 0; x < width; x++) {
527 double bx = (meshsize - 1) * UPRJ(x, width);
528
529 bxf[x] = floor(bx);
530 bxc[x] = bxf[x] + 1;
531 u[x] = bx - bxf[x];
532 u1[x] = 1 - u[x];
533 }
534
535 for (y = 0; y < height; y++) {
536 double t, t1, by;
537 int byf, byc;
538
539 by = (meshsize - 1) * UPRJ(y, height);
540 byf = floor(by) * meshsize;
541 byc = byf + meshsize;
542 t = by - floor(by);
543 t1 = 1 - t;
544
545 if (clouds) {
546 /* Render the FFT output as clouds. */
547
548 for (x = 0; x < width; x++) {
549 double r = t1 * u1[x] * cp[byf + bxf[x]] +
550 t * u1[x] * cp[byc + bxf[x]] +
551 t * u[x] * cp[byc + bxc[x]] +
552 t1 * u[x] * cp[byf + bxc[x]];
553 unsigned char w = (r > 127.0) ? (255 * ((r - 127.0) / 128.0)) : 0;
554
555 w = (255*(100-density) + density * (unsigned int)w) / 100;
556 /* blue of some density through white */
557 *op++ = w;
558 *op++ = w;
559 *op++ = 255;
560 }
561 } else
562 for (x = 0; x < width; x++) {
563 double r = t1 * u1[x] * cp[byf + bxf[x]] +
564 t * u1[x] * cp[byc + bxf[x]] +
565 t * u[x] * cp[byc + bxc[x]] +
566 t1 * u[x] * cp[byf + bxc[x]];
567 int ir, ig, ib;
568
569 if (r >= 128) {
570 int ix = ((r - 128) * (ELEMENTS(pgnd) - 1)) / 127;
571
572 /* Land area. Look up colour based on elevation from
573 precomputed colour map table. */
574
575 ir = pgnd[ix][0];
576 ig = pgnd[ix][1];
577 ib = pgnd[ix][2];
578 } else {
579 /* Water. Generate clouds above water based on elevation. */
580 ir = ig = density * ((r > 64) ? (r - 64)*4 : 0) / 100;
581 ib = 255;
582 }
583
584 *op++ = ir;
585 *op++ = ig;
586 *op++ = ib;
587 }
588 }
589
590 XtFree((char *) cp);
591 XtFree((char *) u);
592 XtFree((char *) u1);
593 XtFree((char *) bxf);
594 XtFree((char *) bxc);
595 XtFree((char *) a);
596
597 return output;
598 }
599
600 #endif /* FEATURE_FRACTAL */
601