From 1fd361a1ea06e44286c213ca1f814f49306fdc43 Mon Sep 17 00:00:00 2001 From: giraffedata Date: Sat, 19 Aug 2006 03:12:28 +0000 Subject: Create Subversion repository git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@1 9d0c8265-081b-0410-96cb-a4ca84ce46f8 --- generator/ppmforge.c | 1182 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1182 insertions(+) create mode 100644 generator/ppmforge.c (limited to 'generator/ppmforge.c') diff --git a/generator/ppmforge.c b/generator/ppmforge.c new file mode 100644 index 00000000..64b1ad79 --- /dev/null +++ b/generator/ppmforge.c @@ -0,0 +1,1182 @@ +/* + + Fractal forgery generator for the PPM toolkit + + Originally designed and implemented in December of 1989 by + John Walker as a stand-alone program for the Sun and MS-DOS + under Turbo C. Adapted in September of 1991 for use with Jef + Poskanzer's raster toolkit. + + References cited in the comments are: + + Foley, J. D., and Van Dam, A., Fundamentals of Interactive + Computer Graphics, Reading, Massachusetts: Addison + Wesley, 1984. + + Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal + Images, New York: Springer Verlag, 1988. + + Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, + W. T., Numerical Recipes In C, New Rochelle: Cambridge + University Press, 1988. + + Author: + John Walker + http://www.fourmilab.ch/ + + 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. + +*/ + +#define _XOPEN_SOURCE /* get M_PI in math.h */ + +#include +#include + +#include "pm_c_util.h" +#include "ppm.h" +#include "mallocvar.h" + +#ifdef VMS +static double const hugeVal = HUGE_VAL; +#else +static double const hugeVal = 1e50; +#endif + +/* Definitions used to address real and imaginary parts in a two-dimensional + array of complex numbers as stored by fourn(). */ + +#define Real(v, x, y) v[1 + (((x) * meshsize) + (y)) * 2] +#define Imag(v, x, y) v[2 + (((x) * meshsize) + (y)) * 2] + +/* Co-ordinate indices within arrays. */ + +typedef struct { + double x; + double y; + double z; +} vector; + +/* Definition for obtaining random numbers. */ + +#define nrand 4 /* Gauss() sample count */ +#define Cast(low, high) ((low)+(((high)-(low)) * ((rand() & 0x7FFF) / arand))) + +/* prototypes */ +static void fourn ARGS((float data[], int nn[], int ndim, int isign)); +static void initgauss ARGS((unsigned int seed)); +static double gauss ARGS((void)); +static void spectralsynth ARGS((float **x, unsigned int n, double h)); +static void temprgb ARGS((double temp, double *r, double *g, double *b)); +static void etoile ARGS((pixel *pix)); +/* Local variables */ + +static double arand, gaussadd, gaussfac; /* Gaussian random parameters */ +static double fracdim; /* Fractal dimension */ +static double powscale; /* Power law scaling exponent */ +static int meshsize = 256; /* FFT mesh size */ +static unsigned int seedarg; /* Seed specified by user */ +static bool seedspec = FALSE; /* Did the user specify a seed ? */ +static bool clouds = FALSE; /* Just generate clouds */ +static bool stars = FALSE; /* Just generate stars */ +static int screenxsize = 256; /* Screen X size */ +static int screenysize = 256; /* Screen Y size */ +static double inclangle, hourangle; /* Star position relative to planet */ +static bool inclspec = FALSE; /* No inclination specified yet */ +static bool hourspec = FALSE; /* No hour specified yet */ +static double icelevel; /* Ice cap theshold */ +static double glaciers; /* Glacier level */ +static int starfraction; /* Star fraction */ +static int starcolor; /* Star color saturation */ + +/* 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(data, nn, ndim, isign) + float data[]; + int nn[], ndim, 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 + +/* INITGAUSS -- Initialize random number generators. As given in + Peitgen & Saupe, page 77. */ + +static void initgauss(seed) + unsigned int seed; +{ + /* Range of random generator */ + arand = pow(2.0, 15.0) - 1.0; + gaussadd = sqrt(3.0 * nrand); + gaussfac = 2 * gaussadd / (nrand * arand); + srand(seed); +} + +/* GAUSS -- Return a Gaussian random number. As given in Peitgen + & Saupe, page 77. */ + +static double gauss() +{ + int i; + double sum = 0.0; + + for (i = 1; i <= nrand; i++) { + sum += (rand() & 0x7FFF); + } + return gaussfac * sum - gaussadd; +} + +/* SPECTRALSYNTH -- Spectrally synthesized fractal motion in two + dimensions. This algorithm is given under the + name SpectralSynthesisFM2D on page 108 of + Peitgen & Saupe. */ + +static void spectralsynth(x, n, h) + 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 *) calloc(bl, 1); + if (a == (float *) 0) { + pm_error("Cannot allocate %d x %d result array (% d bytes).", + n, n, bl); + } + *x = a; + + for (i = 0; i <= n / 2; i++) { + for (j = 0; j <= n / 2; j++) { + phase = 2 * M_PI * ((rand() & 0x7FFF) / arand); + 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 * ((rand() & 0x7FFF) / arand); + 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 */ +} + + +static unsigned int +initseed(void) { + /* Generate initial random seed. */ + + int i; + + i = time(NULL) ^ 0xF37C; + srand(i); + for (i = 0; i < 7; ++i) + rand(); + return rand(); +} + + + +/* TEMPRGB -- Calculate the relative R, G, and B components for a + black body emitting light at a given temperature. + The Planck radiation equation is solved directly for + the R, G, and B wavelengths defined for the CIE 1931 + Standard Colorimetric Observer. The color + temperature is specified in degrees Kelvin. */ + +static void temprgb(temp, r, g, b) + double temp; + double *r, *g, *b; +{ + double c1 = 3.7403e10, + c2 = 14384.0, + er, eg, eb, es; + +/* Lambda is the wavelength in microns: 5500 angstroms is 0.55 microns. */ + +#define Planck(lambda) ((c1 * pow((double) lambda, -5.0)) / \ + (pow(M_E, c2 / (lambda * temp)) - 1)) + + er = Planck(0.7000); + eg = Planck(0.5461); + eb = Planck(0.4358); +#undef Planck + + es = 1.0 / MAX(er, MAX(eg, eb)); + + *r = er * es; + *g = eg * es; + *b = eb * es; +} + +/* ETOILE -- Set a pixel in the starry sky. */ + +static void etoile(pix) + pixel *pix; +{ + if ((rand() % 1000) < starfraction) { +#define StarQuality 0.5 /* Brightness distribution exponent */ +#define StarIntensity 8 /* Brightness scale factor */ +#define StarTintExp 0.5 /* Tint distribution exponent */ + double v = StarIntensity * pow(1 / (1 - Cast(0, 0.9999)), + (double) StarQuality), + temp, + r, g, b; + + if (v > 255) { + v = 255; + } + + /* We make a special case for star color of zero in order to + prevent floating point roundoff which would otherwise + result in more than 256 star colors. We can guarantee + that if you specify no star color, you never get more than + 256 shades in the image. */ + + if (starcolor == 0) { + int vi = v; + + PPM_ASSIGN(*pix, vi, vi, vi); + } else { + temp = 5500 + starcolor * + pow(1 / (1 - Cast(0, 0.9999)), StarTintExp) * + ((rand() & 7) ? -1 : 1); + /* Constrain temperature to a reasonable value: >= 2600K + (S Cephei/R Andromedae), <= 28,000 (Spica). */ + temp = MAX(2600, MIN(28000, temp)); + temprgb(temp, &r, &g, &b); + PPM_ASSIGN(*pix, (int) (r * v + 0.499), + (int) (g * v + 0.499), + (int) (b * v + 0.499)); + } + } else { + PPM_ASSIGN(*pix, 0, 0, 0); + } +} + + + +static double +uprj(unsigned int const a, + unsigned int const size) { + + return (double)a/(size-1); +} + + + +static double +atSat(double const x, + double const y, + double const dsat) { + + return x*(1.0-dsat) + y*dsat; +} + + + +static unsigned char * +makeCp(float * const a, + unsigned int const n, + pixval const maxval) { + + /* Prescale the grid points into intensities. */ + + unsigned char * cp; + unsigned char * ap; + + if (UINT_MAX / n < n) + pm_error("arithmetic overflow squaring %u", n); + cp = malloc(n * n); + if (cp == NULL) + pm_error("Unable to allocate %u bytes for cp array", n); + + ap = cp; + { + unsigned int i; + for (i = 0; i < n; i++) { + unsigned int j; + for (j = 0; j < n; j++) + *ap++ = ((double)maxval * (Real(a, i, j) + 1.0)) / 2.0; + } + } + return cp; +} + + + +static void +createPlanetStuff(float * const a, + unsigned int const n, + double ** const uP, + double ** const u1P, + unsigned int ** const bxfP, + unsigned int ** const bxcP, + unsigned char ** const cpP, + vector * const sunvecP, + unsigned int const cols, + pixval const maxval) { + + double *u, *u1; + unsigned int *bxf, *bxc; + unsigned char * cp; + double shang, siang; + bool flipped; + + /* Compute incident light direction vector. */ + + shang = hourspec ? hourangle : Cast(0, 2 * M_PI); + siang = inclspec ? inclangle : Cast(-M_PI * 0.12, M_PI * 0.12); + + sunvecP->x = sin(shang) * cos(siang); + sunvecP->y = sin(siang); + sunvecP->z = cos(shang) * cos(siang); /* initial value */ + + /* Allow only 25% of random pictures to be crescents */ + + if (!hourspec && ((rand() % 100) < 75)) { + flipped = (sunvecP->z < 0); + sunvecP->z = fabs(sunvecP->z); + } else + flipped = FALSE; + + if (!clouds) { + pm_message( + " -inclination %.0f -hour %d -ice %.2f -glaciers %.2f", + (siang * (180.0 / M_PI)), + (int) (((shang * (12.0 / M_PI)) + 12 + + (flipped ? 12 : 0)) + 0.5) % 24, icelevel, glaciers); + pm_message(" -stars %d -saturation %d.", + starfraction, starcolor); + } + + cp = makeCp(a, n, maxval); + + /* 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.) + + Also compute the line-level interpolation parameters that + caller will need every time around his inner loop. + */ + + MALLOCARRAY(u, cols); + MALLOCARRAY(u1, cols); + MALLOCARRAY(bxf, cols); + MALLOCARRAY(bxc, cols); + + if (u == NULL || u1 == NULL || bxf == NULL || bxc == NULL) + pm_error("Cannot allocate %d element interpolation tables.", cols); + { + unsigned int j; + for (j = 0; j < cols; j++) { + double const bx = (n - 1) * uprj(j, cols); + + bxf[j] = floor(bx); + bxc[j] = bxf[j] + 1; + u[j] = bx - bxf[j]; + u1[j] = 1 - u[j]; + } + } + *uP = u; *u1P = u1; + *bxfP = bxf; *bxcP = bxc; + *cpP = cp; +} + + + +static void +generateStarrySkyRow(pixel * const pixels, + unsigned int const cols) { +/*---------------------------------------------------------------------------- + Generate a starry sky. Note that no FFT is performed; + the output is generated directly from a power law + mapping of a pseudorandom sequence into intensities. +-----------------------------------------------------------------------------*/ + unsigned int j; + + for (j = 0; j < cols; j++) + etoile(pixels + j); +} + + + +static void +generateCloudRow(pixel * const pixels, + unsigned int const cols, + double const t, + double const t1, + double * const u, + double * const u1, + unsigned char * const cp, + unsigned int * const bxc, + unsigned int * const bxf, + int const byc, + int const byf, + pixval const maxval) { + + /* Render the FFT output as clouds. */ + + unsigned int j; + + for (j = 0; j < cols; j++) { + double r; + pixval w; + + r = 0.0; /* initial value */ + /* Note that where t1 and t are zero, the cp[] element + referenced below does not exist. + */ + if (t1 > 0.0) + r += t1 * u1[j] * cp[byf + bxf[j]] + + t1 * u[j] * cp[byf + bxc[j]]; + if (t > 0.0) + r += t * u1[j] * cp[byc + bxf[j]] + + t * u[j] * cp[byc + bxc[j]]; + + w = (r > 127.0) ? (maxval * ((r - 127.0) / 128.0)) : 0; + + PPM_ASSIGN(*(pixels + j), w, w, maxval); + } +} + + + +static void +makeLand(int * const irP, + int * const igP, + int * const ibP, + double const r) { +/*---------------------------------------------------------------------------- + Land area. Look up color based on elevation from precomputed + color map table. +-----------------------------------------------------------------------------*/ + 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} + }; + + unsigned int const ix = ((r - 128) * (ARRAY_SIZE(pgnd) - 1)) / 127; + + *irP = pgnd[ix][0]; + *igP = pgnd[ix][1]; + *ibP = pgnd[ix][2]; +} + + + +static void +makeWater(int * const irP, + int * const igP, + int * const ibP, + double const r, + pixval const maxval) { + + /* Water. Generate clouds above water based on elevation. */ + + *irP = *igP = r > 64 ? (r - 64) * 4 : 0; + *ibP = maxval; +} + + + +static void +addIce(int * const irP, + int * const igP, + int * const ibP, + double const r, + double const azimuth, + double const icelevel, + double const glaciers, + pixval const maxval) { + + /* Generate polar ice caps. */ + + double const icet = pow(fabs(sin(azimuth)), 1.0 / icelevel) - 0.5; + double const ice = MAX(0.0, + (icet + glaciers * MAX(-0.5, (r - 128) / 128.0))); + if (ice > 0.125) { + *irP = maxval; + *igP = maxval; + *ibP = maxval; + } +} + + + +static void +limbDarken(int * const irP, + int * const igP, + int * const ibP, + unsigned int const col, + unsigned int const row, + unsigned int const cols, + unsigned int const rows, + vector const sunvec, + pixval const maxval) { + + /* With Gcc 2.95.3 compiler optimization level > 1, I have seen this + function confuse all the variables and ultimately generate a + completely black image. Adding an extra reference to 'rows' seems + to put things back in order, and the assert() below does that. + Take it out, and the problem comes back! 04.02.21. + */ + + /* Apply limb darkening by cosine rule. */ + + double const atthick = 1.03; + double const atSatFac = 1.0; + double const athfac = sqrt(atthick * atthick - 1.0); + /* Atmosphere thickness as a percentage of planet's diameter */ + + double const dy = 2 * ((double)rows/2 - row) / rows; + double const dysq = dy * dy; + /* Note: we are in fact normalizing this horizontal position by the + vertical size of the picture. And we know rows >= cols. + */ + double const dx = 2 * ((double)cols/2 - col) / rows; + double const dxsq = dx * dx; + + double const ds = MIN(1.0, sqrt(dxsq + dysq)); + + /* Calculate atmospheric absorption based on the thickness of + atmosphere traversed by light on its way to the surface. + */ + double const dsq = ds * ds; + double const dsat = atSatFac * ((sqrt(atthick * atthick - dsq) - + sqrt(1.0 * 1.0 - dsq)) / athfac); + + assert(rows >= cols); /* An input requirement */ + + *irP = atSat(*irP, maxval/2, dsat); + *igP = atSat(*igP, maxval/2, dsat); + *ibP = atSat(*ibP, maxval, dsat); + { + double const PlanetAmbient = 0.05; + + double const sqomdysq = sqrt(1.0 - dysq); + double const svx = sunvec.x; + double const svy = sunvec.y * dy; + double const svz = sunvec.z * sqomdysq; + double const di = + MAX(0, MIN(1.0, svx * dx + svy + svz * sqrt(1.0 - dxsq))); + double const inx = PlanetAmbient * 1.0 + (1.0 - PlanetAmbient) * di; + + *irP *= inx; + *igP *= inx; + *ibP *= inx; + } +} + + + +static void +generatePlanetRow(pixel * const pixelrow, + unsigned int const row, + unsigned int const rows, + unsigned int const cols, + double const t, + double const t1, + double * const u, + double * const u1, + unsigned char * const cp, + unsigned int * const bxc, + unsigned int * const bxf, + int const byc, + int const byf, + vector const sunvec, + pixval const maxval) { + + unsigned int const StarClose = 2; + + double const azimuth = asin(((((double) row) / (rows - 1)) * 2) - 1); + unsigned int const lcos = (rows / 2) * fabs(cos(azimuth)); + + unsigned int col; + + for (col = 0; col < cols; ++col) + PPM_ASSIGN(pixelrow[col], 0, 0, 0); + + for (col = cols/2 - lcos; col <= cols/2 + lcos; ++col) { + double r; + int ir, ig, ib; + + r = 0.0; /* initial value */ + + /* Note that where t1 and t are zero, the cp[] element + referenced below does not exist. + */ + if (t1 > 0.0) + r += t1 * u1[col] * cp[byf + bxf[col]] + + t1 * u[col] * cp[byf + bxc[col]]; + if (t > 0.0) + r += t * u1[col] * cp[byc + bxf[col]] + + t * u[col] * cp[byc + bxc[col]]; + + if (r >= 128) + makeLand(&ir, &ig, &ib, r); + else + makeWater(&ir, &ig, &ib, r, maxval); + + addIce(&ir, &ig, &ib, r, azimuth, icelevel, glaciers, maxval); + + limbDarken(&ir, &ig, &ib, col, row, cols, rows, sunvec, maxval); + + PPM_ASSIGN(pixelrow[col], ir, ig, ib); + } + + /* Left stars */ + + for (col = 0; (int)col < (int)(cols/2 - (lcos + StarClose)); ++col) + etoile(&pixelrow[col]); + + /* Right stars */ + + for (col = cols/2 + (lcos + StarClose); col < cols; ++col) + etoile(&pixelrow[col]); +} + + + +static void +genplanet(bool const stars, + bool const clouds, + float * const a, + unsigned int const cols, + unsigned int const rows, + unsigned int const n, + unsigned int const rseed) { +/*---------------------------------------------------------------------------- + Generate planet from elevation array. + + If 'stars' is true, a is undefined. Otherwise, it is defined. +-----------------------------------------------------------------------------*/ + pixval const maxval = PPM_MAXMAXVAL; + + unsigned char *cp; + double *u, *u1; + unsigned int *bxf, *bxc; + + pixel *pixelrow; + unsigned int row; + + vector sunvec; + + ppm_writeppminit(stdout, cols, rows, maxval, FALSE); + + if (stars) { + pm_message("night: -seed %d -stars %d -saturation %d.", + rseed, starfraction, starcolor); + cp = NULL; + u = NULL; u1 = NULL; + bxf = NULL; bxc = NULL; + } else { + pm_message("%s: -seed %d -dimension %.2f -power %.2f -mesh %d", + clouds ? "clouds" : "planet", + rseed, fracdim, powscale, meshsize); + createPlanetStuff(a, n, &u, &u1, &bxf, &bxc, &cp, &sunvec, + cols, maxval); + } + + pixelrow = ppm_allocrow(cols); + for (row = 0; row < rows; ++row) { + if (stars) + generateStarrySkyRow(pixelrow, cols); + else { + double const by = (n - 1) * uprj(row, rows); + int const byf = floor(by) * n; + int const byc = byf + n; + double const t = by - floor(by); + double const t1 = 1 - t; + + if (clouds) + generateCloudRow(pixelrow, cols, + t, t1, u, u1, cp, bxc, bxf, byc, byf, + maxval); + else + generatePlanetRow(pixelrow, row, rows, cols, + t, t1, u, u1, cp, bxc, bxf, byc, byf, + sunvec, + maxval); + } + ppm_writeppmrow(stdout, pixelrow, cols, maxval, FALSE); + } + pm_close(stdout); + + ppm_freerow(pixelrow); + if (cp) free(cp); + if (u) free(u); + if (u1) free(u1); + if (bxf) free(bxf); + if (bxc) free(bxc); +} + + + +static void +applyPowerLawScaling(float * const a, + int const meshsize, + double const powscale) { + + /* Apply power law scaling if non-unity scale is requested. */ + + if (powscale != 1.0) { + unsigned int i; + for (i = 0; i < meshsize; i++) { + unsigned int j; + for (j = 0; j < meshsize; j++) { + double const r = Real(a, i, j); + if (r > 0) + Real(a, i, j) = pow(r, powscale); + } + } + } +} + + + +static void +computeExtremeReal(const float * const a, + int const meshsize, + double * const rminP, + double * const rmaxP) { + + /* Compute extrema for autoscaling. */ + + double rmin, rmax; + unsigned int i; + + rmin = hugeVal; + rmax = -hugeVal; + + for (i = 0; i < meshsize; i++) { + unsigned int j; + for (j = 0; j < meshsize; j++) { + double r = Real(a, i, j); + + rmin = MIN(rmin, r); + rmax = MAX(rmax, r); + } + } + *rminP = rmin; + *rmaxP = rmax; +} + + + +static void +replaceWithSpread(float * const a, + int const meshsize) { +/*---------------------------------------------------------------------------- + Replace the real part of each element of the 'a' array with a + measure of how far the real is from the middle; sort of a standard + deviation. +-----------------------------------------------------------------------------*/ + double rmin, rmax; + double rmean, rrange; + unsigned int i; + + computeExtremeReal(a, meshsize, &rmin, &rmax); + + rmean = (rmin + rmax) / 2; + rrange = (rmax - rmin) / 2; + + for (i = 0; i < meshsize; i++) { + unsigned int j; + for (j = 0; j < meshsize; j++) { + Real(a, i, j) = (Real(a, i, j) - rmean) / rrange; + } + } +} + + + +static bool +planet(unsigned int const cols, + unsigned int const rows, + bool const stars, + bool const clouds) { +/*---------------------------------------------------------------------------- + Make a planet. +-----------------------------------------------------------------------------*/ + float * a; + bool error; + unsigned int rseed; /* Current random seed */ + + if (seedspec) + rseed = seedarg; + else + rseed = initseed(); + + initgauss(rseed); + + if (stars) { + a = NULL; + error = FALSE; + } else { + spectralsynth(&a, meshsize, 3.0 - fracdim); + if (a == NULL) { + error = TRUE; + } else { + applyPowerLawScaling(a, meshsize, powscale); + + replaceWithSpread(a, meshsize); + + error = FALSE; + } + } + if (!error) + genplanet(stars, clouds, a, cols, rows, meshsize, rseed); + + if (a != NULL) + free(a); + + return !error; +} + + + + +int +main(int argc, char ** argv) { + + bool success; + int i; + const char * const usage = "\n\ +[-width|-xsize ] [-height|-ysize ] [-mesh ]\n\ +[-clouds] [-dimension ] [-power ] [-seed ]\n\ +[-hour ] [-inclination|-tilt ] [-ice ] [-glaciers ]\n\ +[-night] [-stars ] [-saturation ]"; + bool dimspec = FALSE, meshspec = FALSE, powerspec = FALSE, + widspec = FALSE, hgtspec = FALSE, icespec = FALSE, + glacspec = FALSE, starspec = FALSE, starcspec = FALSE; + + int cols, rows; /* Dimensions of our output image */ + + ppm_init(&argc, argv); + i = 1; + + while ((i < argc) && (argv[i][0] == '-') && (argv[i][1] != '\0')) { + + if (pm_keymatch(argv[i], "-clouds", 2)) { + clouds = TRUE; + } else if (pm_keymatch(argv[i], "-night", 2)) { + stars = TRUE; + } else if (pm_keymatch(argv[i], "-dimension", 2)) { + if (dimspec) { + pm_error("already specified a dimension"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &fracdim) != 1)) + pm_usage(usage); + if (fracdim <= 0.0) { + pm_error("fractal dimension must be greater than 0"); + } + dimspec = TRUE; + } else if (pm_keymatch(argv[i], "-hour", 3)) { + if (hourspec) { + pm_error("already specified an hour"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &hourangle) != 1)) + pm_usage(usage); + hourangle = (M_PI / 12.0) * (hourangle + 12.0); + hourspec = TRUE; + } else if (pm_keymatch(argv[i], "-inclination", 3) || + pm_keymatch(argv[i], "-tilt", 2)) { + if (inclspec) { + pm_error("already specified an inclination/tilt"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &inclangle) != 1)) + pm_usage(usage); + inclangle = (M_PI / 180.0) * inclangle; + inclspec = TRUE; + } else if (pm_keymatch(argv[i], "-mesh", 2)) { + unsigned int j; + + if (meshspec) { + pm_error("already specified a mesh size"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%d", &meshsize) != 1)) + pm_usage(usage); + + /* Force FFT mesh to the next larger power of 2. */ + + for (j = meshsize; (j & 1) == 0; j >>= 1) ; + + if (j != 1) { + for (j = 2; j < meshsize; j <<= 1) ; + meshsize = j; + } + meshspec = TRUE; + } else if (pm_keymatch(argv[i], "-power", 2)) { + if (powerspec) { + pm_error("already specified a power factor"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &powscale) != 1)) + pm_usage(usage); + if (powscale <= 0.0) { + pm_error("power factor must be greater than 0"); + } + powerspec = TRUE; + } else if (pm_keymatch(argv[i], "-ice", 3)) { + if (icespec) { + pm_error("already specified ice cap level"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &icelevel) != 1)) + pm_usage(usage); + if (icelevel <= 0.0) { + pm_error("ice cap level must be greater than 0"); + } + icespec = TRUE; + } else if (pm_keymatch(argv[i], "-glaciers", 2)) { + if (glacspec) { + pm_error("already specified glacier level"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%lf", &glaciers) != 1)) + pm_usage(usage); + if (glaciers <= 0.0) { + pm_error("glacier level must be greater than 0"); + } + glacspec = TRUE; + } else if (pm_keymatch(argv[i], "-stars", 3)) { + if (starspec) { + pm_error("already specified a star fraction"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%d", &starfraction) != 1)) + pm_usage(usage); + starspec = TRUE; + } else if (pm_keymatch(argv[i], "-saturation", 3)) { + if (starcspec) { + pm_error("already specified a star color saturation"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%d", &starcolor) != 1)) + pm_usage(usage); + starcspec = TRUE; + } else if (pm_keymatch(argv[i], "-seed", 3)) { + if (seedspec) { + pm_error("already specified a random seed"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%u", &seedarg) != 1)) + pm_usage(usage); + seedspec = 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); + widspec = TRUE; + } else if (pm_keymatch(argv[i], "-ysize", 2) || + pm_keymatch(argv[i], "-height", 3)) { + if (hgtspec) { + pm_error("already specified a height/ysize"); + } + i++; + if ((i == argc) || (sscanf(argv[i], "%d", &screenysize) != 1)) + pm_usage(usage); + hgtspec = TRUE; + } else { + pm_usage(usage); + } + i++; + } + + + /* Set defaults when explicit specifications were not given. + + The default fractal dimension and power scale depend upon + whether we are generating a planet or clouds. + */ + + if (!dimspec) { + fracdim = clouds ? 2.15 : 2.4; + } + if (!powerspec) { + powscale = clouds ? 0.75 : 1.2; + } + if (!icespec) { + icelevel = 0.4; + } + if (!glacspec) { + glaciers = 0.75; + } + if (!starspec) { + starfraction = 100; + } + if (!starcspec) { + starcolor = 125; + } + + /* Force screen to be at least as wide as it is high. Long, + skinny screens cause crashes because picture width is + calculated based on height. + */ + + cols = (MAX(screenysize, screenxsize) + 1) & (~1); + rows = screenysize; + + success = planet(cols, rows, stars, clouds); + + exit(success ? 0 : 1); +} -- cgit 1.4.1