/* * This file is part of XPaint. * * Copyright 1996 by Torsten Martinsen * Copyright 1994 by Jer Johnson * * This is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. */ #include #include #include #include #include #include #include "xpaint.h" #include "misc.h" #include "image.h" #ifndef M_PI #define M_PI 3.14159265358979323846 #endif #if FEATURE_FRACTAL static void draw_grid(int x, int y, unsigned char col, unsigned char *grid); static void subdiv(int x1, int y1, int x2, int y2, unsigned char *grid); static void adjust(int x1, int y1, int x, int y, int x2, int y2, unsigned char *grid); static void fourn(float data[], int nn[], int ndim, int isign); static void spectralsynth(float **x, unsigned int n, double h); static int number, numcolors, width, height; /* * draw_plasma() and friends are stolen from Xtacy, written by Jer Johnson * (mpython@gnu.ai.mit.edu) who says that his work in turn is based on OOZE.C, * written by Jeff Clough. */ Image * draw_plasma(int w, int h) { int i, n; unsigned char *grid; Image *output; unsigned char *op; /* customizable parameters */ numcolors = 100; number = 10; width = w; n = numcolors / 5; numcolors = n * 5; output = ImageNewCmap(w, h, numcolors); grid = output->data; memset(grid, 0, w * h); /* * Make points at (0,0) (0, h-1) (w-1, 0) (w-1, h-1) of random colors */ draw_grid(0, 0, RANDOMI2(0, numcolors-1), grid); draw_grid(w - 1, 0, RANDOMI2(0, numcolors-1), grid); draw_grid(0, h - 1, RANDOMI2(0, numcolors-1), grid); draw_grid(w - 1, h - 1, RANDOMI2(0, numcolors-1), grid); subdiv(0, 0, w - 1, h - 1, grid); /* Build the colour map */ op = output->cmapData; for (i = 0; i < n; ++i) { *op++ = 255; *op++ = 255 * i / n; *op++ = 0; } for (i = 0; i < n; ++i) { *op++ = 255 - 255 * i / n; *op++ = 255; *op++ = 0; } for (i = 0; i < n; ++i) { *op++ = 0; *op++ = 255; *op++ = 255 * i / n; } for (i = 0; i < n; ++i) { *op++ = 255 * i / n; *op++ = 0; *op++ = 255; } for (i = 0; i < n; ++i) { *op++ = 255; *op++ = 0; *op++ = 255 - 255 * i / n; } #if 0 op = output->cmapData; for (i = 0; i < 5 * n; ++i) { printf("%d: %d %d %d\n", i, op[0], op[1], op[2]); op += 3; } #endif return output; } static void draw_grid(int x, int y, unsigned char col, unsigned char *grid) { grid[x + y * width] = col; } static void subdiv(int x1, int y1, int x2, int y2, unsigned char *grid) { int x, y; if ((x2 - x1 < 2) && (y2 - y1 < 2)) return; x = (x1 + x2) / 2; y = (y1 + y2) / 2; adjust(x1, y1, x, y1, x2, y1, grid); adjust(x2, y1, x2, y, x2, y2, grid); adjust(x1, y2, x, y2, x2, y2, grid); adjust(x1, y1, x1, y, x1, y2, grid); if (grid[x + y * width] == 0) { unsigned char spooge; spooge = ((int) grid[x1 + y1 * width] + (int) grid[x1 + y2 * width] + (int) grid[x2 + y1 * width] + (int) grid[x2 + y2 * width]) / 4; draw_grid(x, y, spooge, grid); } subdiv(x1, y1, x, y, grid); subdiv(x, y1, x2, y, grid); subdiv(x1, y, x, y2, grid); subdiv(x, y, x2, y2, grid); } static void adjust(int x1, int y1, int x, int y, int x2, int y2, unsigned char *grid) { int c, d, spoo; int clr1, clr2, clr3, clr4, horz, vert; if (grid[x + y * width] != 0) return; clr1 = grid[x1 + y1 * width]; clr2 = grid[x2 + y2 * width]; clr3 = grid[x1 + y2 * width]; clr4 = grid[x2 + y1 * width]; horz = x2 - x1; vert = y2 - y1; d = hypot(horz, vert); #if 1 /* 0 for band effect */ if ((spoo = RANDOMI2(0, 100)) > number) { if (spoo % 2 == 1) c = (((clr1 + clr2) / 2) + ((int) RANDOMI2(0, d))) % numcolors; else c = (((clr1 + clr2) / 2) - ((int) RANDOMI2(0, d))) % numcolors; } else #endif c = (clr1 + clr2 + clr3 + clr4) / 4; if (c < 0) c = -c; else if (c == 0) c = 1; draw_grid(x, y, c, grid); } /* * Fractal landscape generator from ppmforge.c, * written by John Walker (kelvin@Autodesk.com). */ #define Real(v, x, y) v[1 + (((x) * meshsize) + (y)) * 2] #define Imag(v, x, y) v[2 + (((x) * meshsize) + (y)) * 2] #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0])) /* customize: */ static double fracdim = 1.8; /* Fractal dimension - ]0,3] */ static double powscale = 1.0; /* Power law scaling exponent - ]0,3] */ static int meshsize = 256; /* FFT mesh size - generally >= max(w,h), but small numbers can give special effects */ /* colours */ static unsigned char pgnd[][3] = { {206, 205, 0}, {208, 207, 0}, {211, 208, 0}, {214, 208, 0}, {217, 208, 0}, {220, 208, 0}, {222, 207, 0}, {225, 205, 0}, {227, 204, 0}, {229, 202, 0}, {231, 199, 0}, {232, 197, 0}, {233, 194, 0}, {234, 191, 0}, {234, 188, 0}, {233, 185, 0}, {232, 183, 0}, {231, 180, 0}, {229, 178, 0}, {227, 176, 0}, {225, 174, 0}, {223, 172, 0}, {221, 170, 0}, {219, 168, 0}, {216, 166, 0}, {214, 164, 0}, {212, 162, 0}, {210, 161, 0}, {207, 159, 0}, {205, 157, 0}, {203, 156, 0}, {200, 154, 0}, {198, 152, 0}, {195, 151, 0}, {193, 149, 0}, {190, 148, 0}, {188, 147, 0}, {185, 145, 0}, {183, 144, 0}, {180, 143, 0}, {177, 141, 0}, {175, 140, 0}, {172, 139, 0}, {169, 138, 0}, {167, 137, 0}, {164, 136, 0}, {161, 135, 0}, {158, 134, 0}, {156, 133, 0}, {153, 132, 0}, {150, 132, 0}, {147, 131, 0}, {145, 130, 0}, {142, 130, 0}, {139, 129, 0}, {136, 128, 0}, {133, 128, 0}, {130, 127, 0}, {127, 127, 0}, {125, 127, 0}, {122, 127, 0}, {119, 127, 0}, {116, 127, 0}, {113, 127, 0}, {110, 128, 0}, {107, 128, 0}, {104, 128, 0}, {102, 127, 0}, {99, 126, 0}, {97, 124, 0}, {95, 122, 0}, {93, 120, 0}, {92, 117, 0}, {92, 114, 0}, {92, 111, 0}, {93, 108, 0}, {94, 106, 0}, {96, 104, 0}, {98, 102, 0}, {100, 100, 0}, {103, 99, 0}, {106, 99, 0}, {109, 99, 0}, {111, 100, 0}, {114, 101, 0}, {117, 102, 0}, {120, 103, 0}, {123, 102, 0}, {125, 102, 0}, {128, 100, 0}, {130, 98, 0}, {132, 96, 0}, {133, 94, 0}, {134, 91, 0}, {134, 88, 0}, {134, 85, 0}, {133, 82, 0}, {131, 80, 0}, {129, 78, 0} }; /* FOURN -- Multi-dimensional fast Fourier transform Called with arguments: data A one-dimensional array of floats (NOTE!!! NOT DOUBLES!!), indexed from one (NOTE!!! NOT ZERO!!), containing pairs of numbers representing the complex valued samples. The Fourier transformed results are returned in the same array. nn An array specifying the edge size in each dimension. THIS ARRAY IS INDEXED FROM ONE, AND ALL THE EDGE SIZES MUST BE POWERS OF TWO!!! ndim Number of dimensions of FFT to perform. Set to 2 for two dimensional FFT. isign If 1, a Fourier transform is done; if -1 the inverse transformation is performed. This function is essentially as given in Press et al., "Numerical Recipes In C", Section 12.11, pp. 467-470. */ static void fourn(float data[], int nn[], int ndim, int isign) { register int i1, i2, i3; int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2; int ibit, idim, k1, k2, n, nprev, nrem, ntot; float tempi, tempr; double theta, wi, wpi, wpr, wr, wtemp; #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr ntot = 1; for (idim = 1; idim <= ndim; idim++) ntot *= nn[idim]; nprev = 1; for (idim = ndim; idim >= 1; idim--) { n = nn[idim]; nrem = ntot / (n * nprev); ip1 = nprev << 1; ip2 = ip1 * n; ip3 = ip2 * nrem; i2rev = 1; for (i2 = 1; i2 <= ip2; i2 += ip1) { if (i2 < i2rev) { for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) { for (i3 = i1; i3 <= ip3; i3 += ip2) { i3rev = i2rev + i3 - i2; SWAP(data[i3], data[i3rev]); SWAP(data[i3 + 1], data[i3rev + 1]); } } } ibit = ip2 >> 1; while (ibit >= ip1 && i2rev > ibit) { i2rev -= ibit; ibit >>= 1; } i2rev += ibit; } ifp1 = ip1; while (ifp1 < ip2) { ifp2 = ifp1 << 1; theta = isign * (M_PI * 2) / (ifp2 / ip1); wtemp = sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = sin(theta); wr = 1.0; wi = 0.0; for (i3 = 1; i3 <= ifp1; i3 += ip1) { for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) { for (i2 = i1; i2 <= ip3; i2 += ifp2) { k1 = i2; k2 = k1 + ifp1; tempr = wr * data[k2] - wi * data[k2 + 1]; tempi = wr * data[k2 + 1] + wi * data[k2]; data[k2] = data[k1] - tempr; data[k2 + 1] = data[k1 + 1] - tempi; data[k1] += tempr; data[k1 + 1] += tempi; } } wr = (wtemp = wr) * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } ifp1 = ifp2; } nprev *= n; } } #undef SWAP /* * SPECTRALSYNTH -- Spectrally synthesised fractal motion in two dimensions. */ static void spectralsynth(float **x, unsigned int n, double h) { unsigned bl; int i, j, i0, j0, nsize[3]; double rad, phase, rcos, rsin; float *a; bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float); a = (float *) XtCalloc(bl, 1); *x = a; for (i = 0; i <= n / 2; i++) { for (j = 0; j <= n / 2; j++) { phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0; if (i != 0 || j != 0) { rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss(); } else { rad = 0; } rcos = rad * cos(phase); rsin = rad * sin(phase); Real(a, i, j) = rcos; Imag(a, i, j) = rsin; i0 = (i == 0) ? 0 : n - i; j0 = (j == 0) ? 0 : n - j; Real(a, i0, j0) = rcos; Imag(a, i0, j0) = -rsin; } } Imag(a, n / 2, 0) = 0; Imag(a, 0, n / 2) = 0; Imag(a, n / 2, n / 2) = 0; for (i = 1; i <= n / 2 - 1; i++) { for (j = 1; j <= n / 2 - 1; j++) { phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0; rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss(); rcos = rad * cos(phase); rsin = rad * sin(phase); Real(a, i, n - j) = rcos; Imag(a, i, n - j) = rsin; Real(a, n - i, j) = rcos; Imag(a, n - i, j) = -rsin; } } nsize[0] = 0; nsize[1] = nsize[2] = n; /* Dimension of frequency domain array */ fourn(a, nsize, 2, -1); /* Take inverse 2D Fourier transform */ } Image * draw_landscape(int w, int h, int clouds) { float *a; double rmin = HUGE_VAL, rmax = -HUGE_VAL, rmean, rrange; Image *output = ImageNew(w, h); unsigned char *op = output->data; int x, y, i, j; unsigned char *cp, *ap; double *u, *u1; unsigned int *bxf, *bxc; width = w; height = h; spectralsynth(&a, meshsize, 3.0 - fracdim); /* Apply power law scaling if non-unity scale is requested. */ if (powscale != 1.0) for (i = 0; i < meshsize; i++) for (j = 0; j < meshsize; j++) { double r = Real(a, i, j); if (r > 0) Real(a, i, j) = pow(r, powscale); } /* Compute extrema for autoscaling. */ for (i = 0; i < meshsize; i++) for (j = 0; j < meshsize; j++) { double r = Real(a, i, j); rmin = MIN(rmin, r); rmax = MAX(rmax, r); } rmean = (rmin + rmax) / 2; rrange = (rmax - rmin) / 2; for (i = 0; i < meshsize; i++) for (j = 0; j < meshsize; j++) Real(a, i, j) = (Real(a, i, j) - rmean) / rrange; u = (double *) XtMalloc((unsigned int) (width * sizeof(double))); u1 = (double *) XtMalloc((unsigned int) (width * sizeof(double))); bxf = (unsigned int *) XtMalloc((unsigned int) width * sizeof(unsigned int)); bxc = (unsigned int *) XtMalloc((unsigned int) width * sizeof(unsigned int)); /* Prescale the grid points into intensities. */ cp = (unsigned char *) XtMalloc(meshsize * meshsize + 1); ap = cp; for (i = 0; i < meshsize; i++) for (j = 0; j < meshsize; j++) *ap++ = (255.0 * (Real(a, i, j) + 1.0)) / 2.0; /* Hack: make an extra array entry */ *ap = ap[-1]; /* Fill the screen from the computed intensity grid by mapping screen points onto the grid, then calculating each pixel value by bilinear interpolation from the surrounding grid points. (N.b. the pictures would undoubtedly look better when generated with small grids if a more well-behaved interpolation were used.) Before we get started, precompute the line-level interpolation parameters and store them in an array so we don't have to do this every time around the inner loop. */ #define UPRJ(a,size) ((a)/((size)-1.0)) for (x = 0; x < width; x++) { double bx = (meshsize - 1) * UPRJ(x, width); bxf[x] = floor(bx); bxc[x] = bxf[x] + 1; u[x] = bx - bxf[x]; u1[x] = 1 - u[x]; } for (y = 0; y < height; y++) { double t, t1, by; int byf, byc; by = (meshsize - 1) * UPRJ(y, height); byf = floor(by) * meshsize; byc = byf + meshsize; t = by - floor(by); t1 = 1 - t; if (clouds) { /* Render the FFT output as clouds. */ for (x = 0; x < width; x++) { double r = t1 * u1[x] * cp[byf + bxf[x]] + t * u1[x] * cp[byc + bxf[x]] + t * u[x] * cp[byc + bxc[x]] + t1 * u[x] * cp[byf + bxc[x]]; unsigned char w = (r > 127.0) ? (255 * ((r - 127.0) / 128.0)) : 0; *op++ = w; /* blue trough white */ *op++ = w; *op++ = 255; } } else for (x = 0; x < width; x++) { double r = t1 * u1[x] * cp[byf + bxf[x]] + t * u1[x] * cp[byc + bxf[x]] + t * u[x] * cp[byc + bxc[x]] + t1 * u[x] * cp[byf + bxc[x]]; int ir, ig, ib; if (r >= 128) { int ix = ((r - 128) * (ELEMENTS(pgnd) - 1)) / 127; /* Land area. Look up colour based on elevation from precomputed colour map table. */ ir = pgnd[ix][0]; ig = pgnd[ix][1]; ib = pgnd[ix][2]; } else { /* Water. Generate clouds above water based on elevation. */ ir = ig = r > 64 ? (r - 64) * 4 : 0; ib = 255; } *op++ = ir; *op++ = ig; *op++ = ib; } } XtFree(cp); XtFree((char *) u); XtFree((char *) u1); XtFree((char *) bxf); XtFree((char *) bxc); XtFree((char *) a); return output; } #endif /* FEATURE_FRACTAL */