diff options
Diffstat (limited to 'generator/ppmforge.c')
| -rw-r--r-- | generator/ppmforge.c | 258 |
1 files changed, 134 insertions, 124 deletions
diff --git a/generator/ppmforge.c b/generator/ppmforge.c index 114f7f18..e32d2b59 100644 --- a/generator/ppmforge.c +++ b/generator/ppmforge.c @@ -34,12 +34,14 @@ #define _XOPEN_SOURCE 500 /* get M_PI in math.h */ #include <math.h> +#include <float.h> #include <assert.h> #include "pm_c_util.h" -#include "ppm.h" #include "mallocvar.h" +#include "rand.h" #include "shhopt.h" +#include "ppm.h" static double const hugeVal = 1e50; @@ -59,12 +61,8 @@ typedef struct { /* Definition for obtaining random numbers. */ -#define nrand 4 /* Gauss() sample count */ -#define Cast(low, high) ((low)+(((high)-(low)) * ((rand() & 0x7FFF) / arand))) - /* 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 */ @@ -92,6 +90,7 @@ struct CmdlineInfo { float ice; int saturation; unsigned int seed; + unsigned int seedSpec; int stars; unsigned int starsSpec; unsigned int width; @@ -172,6 +171,11 @@ parseCommandLine(int argc, const char **argv, if (cmdlineP->dimension <= 0.0) pm_error("-dimension must be greater than zero. " "You specified %f", cmdlineP->dimension); + else if (cmdlineP->dimension > 5.0 + FLT_EPSILON) + pm_error("-dimension must not be greater than 5. " + "Results are not interesting with higher numbers, so " + "we assume it is a mistake. " + "You specified %f", cmdlineP->dimension); } else cmdlineP->dimension = cmdlineP->clouds ? 2.15 : 2.4; @@ -238,7 +242,7 @@ parseCommandLine(int argc, const char **argv, static void -fourn(float * const data, +fourn(double * const data, const int * const nn, int const ndim, int const isign) { @@ -269,7 +273,7 @@ fourn(float * const data, 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 tempi, tempr; double theta, wi, wpi, wpr, wr, wtemp; #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr @@ -332,72 +336,64 @@ fourn(float * const data, nprev *= n; } } -#undef SWAP -static void -initgauss(unsigned int const seed) { +static double +cast(double const low, + double const high, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- - Initialize random number generators. - - As given in Peitgen & Saupe, page 77. + A random number in the range ['low', 'high']. -----------------------------------------------------------------------------*/ - /* Range of random generator */ - arand = pow(2.0, 15.0) - 1.0; - gaussadd = sqrt(3.0 * nrand); - gaussfac = 2 * gaussadd / (nrand * arand); - srand(seed); + return (low + (high - low) * pm_drand(randStP)); + } static double -gauss() { +randPhase(struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- - A Gaussian random number. - - As given in Peitgen & Saupe, page 77. + A random number in the range [0, 2 * PI). -----------------------------------------------------------------------------*/ - int i; - double sum; + return (2.0 * M_PI * pm_drand1(randStP)); - for (i = 1, sum = 0.0; i <= nrand; ++i) { - sum += (rand() & 0x7FFF); - } - return gaussfac * sum - gaussadd; } static void -spectralsynth(float ** const x, - unsigned int const n, - double const h) { +spectralsynth(double ** const aP, + unsigned int const n, + double const h, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Spectrally synthesized fractal motion in two dimensions. This algorithm is given under the name SpectralSynthesisFM2D on page 108 of Peitgen & Saupe. -----------------------------------------------------------------------------*/ - unsigned bl; + unsigned int const bl = ((((unsigned long) n) * n) + 1) * 2; + int i, j, i0, j0, nsize[3]; double rad, phase, rcos, rsin; - float *a; + double * a; + + MALLOCARRAY(a, bl); - 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).", + if (!a) { + pm_error("Cannot allocate %u x %u result array (%u doubles).", n, n, bl); } - *x = a; + for (i = 0; i < bl; ++i) + a[i] = 0.0; /* initial value */ for (i = 0; i <= n / 2; i++) { for (j = 0; j <= n / 2; j++) { - phase = 2 * M_PI * ((rand() & 0x7FFF) / arand); + phase = randPhase(randStP); if (i != 0 || j != 0) { - rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss(); + rad = pow( (double) (i*i + j*j), - (h + 1) / 2) * pm_gaussrand(randStP); } else { rad = 0; } @@ -416,8 +412,8 @@ spectralsynth(float ** const x, 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(); + phase = randPhase(randStP); + rad = pow((double) (i * i + j * j), -(h + 1) / 2) * pm_gaussrand(randStP); rcos = rad * cos(phase); rsin = rad * sin(phase); Real(a, i, n - j) = rcos; @@ -430,6 +426,8 @@ spectralsynth(float ** const x, nsize[0] = 0; nsize[1] = nsize[2] = n; /* Dimension of frequency domain array */ fourn(a, nsize, 2, -1); /* Take inverse 2D Fourier transform */ + + *aP = a; } @@ -469,22 +467,22 @@ temprgb(double const temp, static void -etoile(pixel * const pix) { +etoile(pixel * const pix, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Set a pixel in the starry sky. -----------------------------------------------------------------------------*/ - 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; - } + if ((pm_rand(randStP) % 1000) < starfraction) { + double const starQuality = 0.5; + /* Brightness distribution exponent */ + double const starIntensity = 8; + /* Brightness scale factor */ + double const starTintExp = 0.5; + /* Tint distribution exponent */ + double const v = + MIN(255.0, + starIntensity * pow(1 / (1 - cast(0, 0.9999, randStP)), + (double) starQuality)); /* We make a special case for star color of zero in order to prevent floating point roundoff which would otherwise @@ -493,13 +491,17 @@ etoile(pixel * const pix) { 256 shades in the image. */ if (starcolor == 0) { - int vi = v; + pixval const vi = v; PPM_ASSIGN(*pix, vi, vi, vi); } else { + double temp; + double r, g, b; + temp = 5500 + starcolor * - pow(1 / (1 - Cast(0, 0.9999)), StarTintExp) * - ((rand() & 7) ? -1 : 1); + pow(1 / (1 - cast(0, 0.9999, randStP)), starTintExp) * + ((pm_rand(randStP) & 7) ? -1 : 1); + /* Constrain temperature to a reasonable value: >= 2600K (S Cephei/R Andromedae), <= 28,000 (Spica). */ temp = MAX(2600, MIN(28000, temp)); @@ -535,9 +537,9 @@ atSat(double const x, static unsigned char * -makeCp(float * const a, - unsigned int const n, - pixval const maxval) { +makeCp(const double * const a, + unsigned int const n, + pixval const maxval) { /* Prescale the grid points into intensities. */ @@ -550,7 +552,7 @@ makeCp(float * const a, if (cp == NULL) pm_error("Unable to allocate %u bytes for cp array", n); - ap = cp; + ap = cp; /* initial value */ { unsigned int i; for (i = 0; i < n; i++) { @@ -565,17 +567,18 @@ makeCp(float * const a, static void -createPlanetStuff(bool const clouds, - 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) { +createPlanetStuff(bool const clouds, + const double * 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, + struct pm_randSt * const randStP) { double *u, *u1; unsigned int *bxf, *bxc; @@ -585,8 +588,8 @@ createPlanetStuff(bool const clouds, /* 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); + shang = hourspec ? hourangle : randPhase(randStP); + siang = inclspec ? inclangle : cast(-M_PI * 0.12, M_PI * 0.12, randStP); sunvecP->x = sin(shang) * cos(siang); sunvecP->y = sin(siang); @@ -594,7 +597,7 @@ createPlanetStuff(bool const clouds, /* Allow only 25% of random pictures to be crescents */ - if (!hourspec && ((rand() % 100) < 75)) { + if (!hourspec && ((pm_rand(randStP) % 100) < 75)) { flipped = (sunvecP->z < 0); sunvecP->z = fabs(sunvecP->z); } else @@ -634,7 +637,7 @@ createPlanetStuff(bool const clouds, pm_error("Cannot allocate %u element interpolation tables.", cols); { unsigned int j; - for (j = 0; j < cols; j++) { + for (j = 0; j < cols; ++j) { double const bx = (n - 1) * uprj(j, cols); bxf[j] = floor(bx); @@ -651,8 +654,9 @@ createPlanetStuff(bool const clouds, static void -generateStarrySkyRow(pixel * const pixels, - unsigned int const cols) { +generateStarrySkyRow(pixel * const pixels, + unsigned int const cols, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Generate a starry sky. Note that no FFT is performed; the output is generated directly from a power law @@ -660,8 +664,8 @@ generateStarrySkyRow(pixel * const pixels, -----------------------------------------------------------------------------*/ unsigned int j; - for (j = 0; j < cols; j++) - etoile(pixels + j); + for (j = 0; j < cols; ++j) + etoile(pixels + j, randStP); } @@ -867,21 +871,22 @@ limbDarken(int * const irP, 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) { +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, + struct pm_randSt * const randStP) { unsigned int const StarClose = 2; @@ -924,24 +929,25 @@ generatePlanetRow(pixel * const pixelrow, /* Left stars */ for (col = 0; (int)col < (int)(cols/2 - (lcos + StarClose)); ++col) - etoile(&pixelrow[col]); + etoile(&pixelrow[col], randStP); /* Right stars */ for (col = cols/2 + (lcos + StarClose); col < cols; ++col) - etoile(&pixelrow[col]); + etoile(&pixelrow[col], randStP); } 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) { +genplanet(bool const stars, + bool const clouds, + const double * const a, + unsigned int const cols, + unsigned int const rows, + unsigned int const n, + unsigned int const rseed, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Generate planet from elevation array. @@ -971,13 +977,13 @@ genplanet(bool const stars, clouds ? "clouds" : "planet", rseed, fracdim, powscale, meshsize); createPlanetStuff(clouds, a, n, &u, &u1, &bxf, &bxc, &cp, &sunvec, - cols, maxval); + cols, maxval, randStP); } pixelrow = ppm_allocrow(cols); for (row = 0; row < rows; ++row) { if (stars) - generateStarrySkyRow(pixelrow, cols); + generateStarrySkyRow(pixelrow, cols, randStP); else { double const by = (n - 1) * uprj(row, rows); int const byf = floor(by) * n; @@ -993,7 +999,8 @@ genplanet(bool const stars, generatePlanetRow(pixelrow, row, rows, cols, t, t1, u, u1, cp, bxc, bxf, byc, byf, sunvec, - maxval); + maxval, + randStP); } ppm_writeppmrow(stdout, pixelrow, cols, maxval, FALSE); } @@ -1010,9 +1017,9 @@ genplanet(bool const stars, static void -applyPowerLawScaling(float * const a, - int const meshsize, - double const powscale) { +applyPowerLawScaling(double * const a, + int const meshsize, + double const powscale) { /* Apply power law scaling if non-unity scale is requested. */ @@ -1023,7 +1030,7 @@ applyPowerLawScaling(float * const a, for (j = 0; j < meshsize; j++) { double const r = Real(a, i, j); if (r > 0) - Real(a, i, j) = pow(r, powscale); + Real(a, i, j) = MIN(hugeVal, pow(r, powscale)); } } } @@ -1032,10 +1039,10 @@ applyPowerLawScaling(float * const a, static void -computeExtremeReal(const float * const a, - int const meshsize, - double * const rminP, - double * const rmaxP) { +computeExtremeReal(const double * const a, + int const meshsize, + double * const rminP, + double * const rmaxP) { /* Compute extrema for autoscaling. */ @@ -1061,8 +1068,8 @@ computeExtremeReal(const float * const a, static void -replaceWithSpread(float * const a, - int const meshsize) { +replaceWithSpread(double * 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 @@ -1096,17 +1103,19 @@ planet(unsigned int const cols, /*---------------------------------------------------------------------------- Make a planet. -----------------------------------------------------------------------------*/ - float * a; + double * a; bool error; + struct pm_randSt randSt; - initgauss(rseed); + pm_randinit(&randSt); + pm_srand(&randSt, rseed); if (stars) { a = NULL; error = FALSE; } else { - spectralsynth(&a, meshsize, 3.0 - fracdim); - if (a == NULL) { + spectralsynth(&a, meshsize, 3.0 - fracdim, &randSt); + if (!a) { error = TRUE; } else { applyPowerLawScaling(a, meshsize, powscale); @@ -1117,7 +1126,7 @@ planet(unsigned int const cols, } } if (!error) - genplanet(stars, clouds, a, cols, rows, meshsize, rseed); + genplanet(stars, clouds, a, cols, rows, meshsize, rseed, &randSt); if (a != NULL) free(a); @@ -1159,7 +1168,8 @@ main(int argc, const char ** argv) { cols = (MAX(cmdline.height, cmdline.width) + 1) & (~1); rows = cmdline.height; - success = planet(cols, rows, cmdline.night, cmdline.clouds, cmdline.seed); + success = planet(cols, rows, cmdline.night, + cmdline.clouds, cmdline.seed); exit(success ? 0 : 1); } |