/* -------------------------------------------------------------- Gforge as a whole, and all routines not commented otherwise, is Copyright 1995 by John Beale The crater code is Copyright 1995 by Heiko Eissfeldt (see crater.c). This program 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. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. ------------------------------------------------------------- Much of the base gforge code was taken from "ppmforge.c" (a "fractal forgery generator") for the PBMPLUS toolkit. From those comments: ppmforge.c was 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. It was adapted in September of 1991 for use with Jef Poskanzer's raster toolkit. Author of ppmforge.c: John Walker kelvin@fourmilab.ch (current as of May 26 1995) References cited in the ppmforge 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. */ #include #include #include #include /* #include */ #include "pp.h" #include "gforge.h" #include "crater.h" /* from Heiko Eissfeldt */ #include "fftn.h" #define GVERSION "gforge v1.3a 18MAY96" #define GAUTHOR "John Beale " void f_filter(float *a, int n, double center, double Q, int sign); void rescale(float **a, int meshsize, float min, float max); void findpeak(float *a, int meshsize, int *imax, int *jmax); void fac_warn(int mesh); /* Local variables */ char output_filename[80]; /* output filename string */ char input_filename[80]; char *string; /* temp string */ static double arand, gaussadd, gaussfac; /* Gaussian random parameters */ static double fracdim; /* Fractal dimension */ static double powscale; /* Power law scaling exponent */ static int meshsize; /* FFT mesh size */ static int screenxsize, screenysize; /* output image size in pixels */ static int xstart, ystart; /* offset in data array for image */ static float xfrac, yfrac; /* peak location in image */ static double dim[10], dscale[10]; /* 1/f powers and relative scales */ static int d_factors; /* number of 1/f powers defined */ static unsigned int rseed = 3; /* Current random seed */ static image_data img; /* struct for info to file open */ static Boolean seedspec = FALSE; /* Did the user specify a seed ? */ static Boolean peakspec = FALSE; /* Specified fixed peak location ? */ static Boolean include_craters = FALSE; /* add craters to surface */ static Boolean wrap = TRUE; /* wrap craters around at edges */ static double c_density = 1.0; /* surface crater density factor */ static double ch_scale = 1.0; /* crater amplitude scaling factor */ static Boolean limitspec = FALSE; /* set different limits? */ static double limit_hi = 1.0; static double limit_lo = 0.0; /* limits for scaling before pow() */ static Boolean bpspec = FALSE; /* specify bandpass filter ? */ static double bpcenter, bpQ; /* bandpass parameters */ static Boolean brspec = FALSE; /* specify bandreject fileter ? */ static double brcenter, brQ; /* bandreject parameters */ static Boolean hpspec = FALSE; /* highpass filter ? */ static double hpcut, hporder; static Boolean lpspec = FALSE; /* lowpass filter ? */ static double lpcut, lporder; /* low-pass cutoff frequency, order */ /* INITGAUSS -- Initialise random number generators. As given in Peitgen & Saupe, page 77. */ static void initgauss(unsigned int seed) { int stemp; /* Range of random generator */ arand = pow(2.0, 15.0) - 1.0; gaussadd = sqrt(3.0 * Nrand); gaussfac = 2 * gaussadd / (Nrand * arand); /* srand(seed); */ stemp = seed; /* observe that 'seed' is unsigned, ie nonnegative */ seed_ran1(stemp); /* seed the generator */ } /* 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 += (ran1() * 0x7FFF); } return gaussfac * sum - gaussadd; } /* fill array with 1/f gaussian noise */ static void fillarray(float *a, int n, double h, double scale) { int x,y, k, i0, j0, rank, nx, ny, xcent, ycent; double rad, phase, rcos, rsin; /* printf("h = %f scale = %f\n",h, scale); */ nx = n; ny = n; /* horizontal and vertical dimension of array */ xcent = (int)(nx / 2.0 - 0.5); /* center dimensions of array */ ycent = (int)(ny / 2.0 - 0.5); printf("filling random array..."); fflush(stdout); /* fill in mx. in order of radius, so we can generate higher resolutions with the same overall aspect (if seed is held fixed) */ for (rank = 0; rank <= xcent; rank++) { /* fill quadrants 2 and 4 */ for (k=0;k<=rank;k++) { x = k; y = rank; phase = 2 * M_PI * ((ran1() * 0x7FFF) / arand); if ((x == 0) && (y == 0)) rad = 0; else rad = pow((double) (x*x + y*y), -(h+1) / 2) * gauss(); rcos = rad * cos(phase)*scale; rsin = rad * sin(phase)*scale; Real(a, x, y) += rcos; Imag(a, x, y) += rsin; if (!((x == 0) && (y == 0))) { i0 = nx-x-1; j0 = ny-y-1; Real(a, i0,j0) += rcos; Imag(a, i0,j0) += rsin; } x = rank; y = k; phase = 2 * M_PI * ((ran1() * 0x7FFF) / arand); if ((x == 0) && (y == 0)) rad = 0; else rad = pow((double) (x*x + y*y), -(h+1) / 2) * gauss(); rcos = rad * cos(phase)*scale; rsin = rad * sin(phase)*scale; Real(a, x, y) += rcos; Imag(a, x, y) += rsin; if (!((x == 0) && (y == 0))) { i0 = nx-x-1; j0 = ny-y-1; Real(a, i0,j0) += rcos; Imag(a, i0,j0) += rsin; } } /* end for k */ /* now handle quadrants 1 and 3 */ for (k=0;k<=rank;k++) { x = k; y = rank; phase = 2 * M_PI * ((ran1() * 0x7FFF) / arand); if ((x == 0) && (y == 0)) rad = 0; else rad = pow((double) (x*x + y*y), -(h+1) / 2) * gauss(); rcos = rad * cos(phase)*scale; rsin = rad * sin(phase)*scale; Real(a, x, ny-y-1) += rcos; Imag(a, x, ny-y-1) += rsin; Real(a, nx-x-1, y) += rcos; Imag(a, nx-x-1, y) += rsin; x = rank; y = k; phase = 2 * M_PI * ((ran1() * 0x7FFF) / arand); if ((x == 0) && (y == 0)) rad = 0; else rad = pow((double) (x*x + y*y), -(h+1) / 2) * gauss(); rcos = rad * cos(phase)*scale; rsin = rad * sin(phase)*scale; Real(a, x, ny-y-1) += rcos; Imag(a, x, ny-y-1) += rsin; Real(a, nx-x-1, y) += rcos; Imag(a, nx-x-1, y) += rsin; } /* end for k */ } /* end for rank */ Imag(a, nx / 2, 0) = 0; Imag(a, 0, ny / 2) = 0; Imag(a, nx / 2, ny / 2) = 0; } /* end fillarray() */ /* f_filter() -- filter array of noise with peak or notch filter */ /* since this happens before the inverse FFT, it is */ /* a frequency-domain filter. */ /* f_type 1 bpass, -1 breject, 2 lopass, -2 hipass */ void f_filter(float *a, int n, double center, double Q, int f_type) { int i, j, i0, j0; double rad, fac, p, sfac; int xsize, ysize; xsize = n; ysize = n; if (center == 0.0) center = -0.00001; /* avoid a singularity */ /* sfac is radius scaling factor */ sfac = 1.0/sqrt((double)(xsize*xsize/4 + ysize*ysize/4)); printf("filtering array.."); fflush(stdout); for (i = 0; i <= n / 2; i++) { /* do quadrants 2 and 4 */ for (j = 0; j <= n / 2; j++) { if (i != 0 || j != 0) { rad = sqrt((double) (i * i + j * j)) * sfac; } else { rad = 0; } p = 1.0 / pow(Q * center, 2); if (abs(f_type)==1) fac = p / (p + pow((1.0-rad/center),2)) ; /* bandpass/rej. */ else fac = 1.0 / (1.0 + pow((rad/center),Q) ); /* lo/hi-pass */ if (f_type < 0) fac = (1.0 - fac); /* invert filter */ Real(a, i, j) *= fac; Imag(a, i, j) *= fac; i0 = (i == 0) ? 0 : n - i; j0 = (j == 0) ? 0 : n - j; Real(a, i0, j0) *= fac; Imag(a, i0, j0) *= fac; } } printf("."); fflush(stdout); 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++) { /* do quadrants 1 and 3 */ for (j = 1; j <= n / 2 - 1; j++) { rad = sfac * sqrt((double) (i * i + j * j)); p = 1.0 / pow(Q * center, 2); if (abs(f_type)==1) fac = p / (p + pow((1.0-rad/center),2)) ; /* bandpass/rej. */ else fac = 1.0 / (1.0 + pow((rad/center),Q) ); /* lo/hi-pass */ if (f_type < 0) fac = (1.0 - fac); /* invert filter */ Real(a, i, n - j) *= fac; Imag(a, i, n - j) *= fac; Real(a, n - i, j) *= fac; Imag(a, n - i, j) *= fac; } } } /* end f_filter() */ /* SPECTRALSYNTH -- Spectrally synthesised fractal motion in two dimensions. This algorithm is given under the name SpectralSynthesisFM2D on page 108 of Peitgen & Saupe. */ static void spectralsynth(float **x, unsigned int n) { int i,j,nsize[2]; unsigned long bl; float *a; bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float); a = (float *) malloc( (size_t) (n+4)*(n+4)*2 * sizeof(float)); if (a == (float *) 0) { fprintf(stderr,"Cannot allocate %d x %d result array (%lu bytes).\n", n, n, bl); exit(1); } else printf("Success allocating %d x %d result array (%lu bytes).\n", n, n, bl); *x = a; for (j=0;j 0) return(1.0); else if (x < 0) return(-1.0); else return(0); } void rescale(float **a, int meshsize, float min, float max) { int i,j; float r; float rmin, rmax, scale, offset; rmin = Real(a[0],0,0); rmax = rmin; for (i = 0; i < meshsize; i++) { /* find current rmin, rmax */ for (j = 0; j < meshsize; j++) { r = Real(a[0], i, j); rmin = min(rmin, r); rmax = max(rmax, r); } } scale = (max-min) / (rmax - rmin); offset = min; /* rescale all data to lie in the range [min..max] */ for (i = 0; i < meshsize; i++) { for (j = 0; j < meshsize; j++) { Real(a[0], i, j) = (scale * (Real(a[0], i, j)-rmin))+ offset; } } /* printf("\nrescale: min=%f max = %f\n",rmin,rmax); */ } /* end rescale() */ /* findpeak() -- find peak location */ void findpeak(float *a, int meshsize, int *imax, int *jmax) { int i,j; float r; float rmax; rmax = Real(a,0,0); imax[0] = 0; jmax[0] = 0; for (i = 0; i < meshsize; i++) { /* find current rmin, rmax */ for (j = 0; j < meshsize; j++) { r = Real(a, i, j); if (r > rmax) { rmax = r; imax[0] = i; jmax[0] = j; } } } } /* end findpeak() */ /* Warn if meshsize is prime or has large prime factor. */ void fac_warn(int mesh) { int fac, flimit,num, maxf; #define div(a,b) ( (b)*floor((a)/(b)) == (a)) fac = 2; num = mesh; maxf = 0; flimit = (1+num/2); do { if (div(num,fac)) { num /= fac; if (fac > maxf) maxf = fac; } else { fac++; } } while (fac < flimit); if ((num > 1) && (mesh > 150)) { printf("\nAdvisory: your meshsize (%d) is a prime number.\n",mesh); printf("If the IFFT is too slow, choose a different meshsize.\n\n"); } else { if (maxf > 150) { printf("\nAdvisory: meshsize %d has a large prime factor (%d).\n",mesh,maxf); printf("If the IFFT is too slow, choose a different meshsize.\n\n"); } } /* end else */ } /* end fac_warn() */ /* MAIN -- Main gforge program: read command line, set vars. */ int main(int argc, char **argv) { int i; char *usage = " Usage: \n\ gforge [-mesh ] [-dimension [-adim ]] [-power ] \n\ [-peak ] [-seed ] [-name ] [-craters [ ]] \n\ [-limit ] [-wrapoff] [-type tga|pgm|pg8|oct|mat|png] \n\ [-bpf ] [-brf ] [-lpf ] [-hpf ] [-version]\n\ -------\n\ Default behavior is to generate a 128x128 TGA-format heightfield.\n\ See the Readme file and the man page for more details.\n"; Boolean dimspec = FALSE, meshspec = FALSE, powerspec = FALSE, namespec = FALSE, typespec = FALSE, adimspec = FALSE, in_namespec = FALSE; float ftmp; /* temporary float value */ /* -------------------------------------------------- */ /* scan command line args */ i=1; d_factors = 1; /* while ((i < argc) && (argv[i][1] != '\0')) { */ while ((i < argc) && (argv[i][1] != '\0')) { if (pm_keymatch(argv[i], "-dimension", 2)) { i++; if ((i == argc) || (sscanf(argv[i], "%lf", &dim[0]) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} if (dim[0] <= 0.0) { pperror("fractal dimension must be greater than 0\n"); } dscale[0] = 1; dimspec = TRUE; } else if (pm_keymatch(argv[i], "-adim", 2)) { if (d_factors > 9) pperror("sheesh! maximum of 10 dimensions\n"); i++; if ((i == argc) || (sscanf(argv[i], "%lf", &dim[d_factors]) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} if (dim[d_factors] < 0.0) pperror("require dimension > 0\n"); i++; if ((i==argc) || (sscanf(argv[i],"%lf",&dscale[d_factors]) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} adimspec = TRUE; /* may be set more than once */ d_factors++; } else if (pm_keymatch(argv[i], "-bpfilter", 3)) { i++; if ((i == argc) || (sscanf(argv[i],"%lf", &bpcenter) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} i++; if ((i == argc) || (sscanf(argv[i],"%lf", &bpQ) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} bpspec = TRUE; } else if (pm_keymatch(argv[i], "-brfilter", 3)) { i++; if ((i == argc) || (sscanf(argv[i],"%lf", &brcenter) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} i++; if ((i == argc) || (sscanf(argv[i],"%lf", &brQ) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} brspec = TRUE; } else if (pm_keymatch(argv[i], "-hpfilter", 3)) { i++; if ((i == argc) || (sscanf(argv[i],"%lf", &hpcut) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} i++; if ((i == argc) || (sscanf(argv[i],"%lf", &hporder) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} hpspec = TRUE; } else if (pm_keymatch(argv[i], "-lpfilter", 3)) { i++; if ((i == argc) || (sscanf(argv[i],"%lf", &lpcut) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} i++; if ((i == argc) || (sscanf(argv[i],"%lf", &lporder) != 1)) fprintf(stderr,GVERSION"%s\n",usage); lpspec = TRUE; } else if (pm_keymatch(argv[i], "-mesh", 2)) { i++; if ((i == argc) || (sscanf(argv[i], "%d", &meshsize) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} if (meshsize <= 0) pperror("A mesh size greater than zero is suggested!\n"); meshspec = TRUE; } else if (pm_keymatch(argv[i], "-power", 3)) { i++; if ((i == argc) || (sscanf(argv[i], "%lf", &powscale) != 1)) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} if (powscale <= 0.0) { pperror("power factor must be greater than 0\n"); } powerspec = TRUE; } else if (pm_keymatch(argv[i], "-limit", 2)) { i++; if ((i == argc) || (sscanf(argv[i], "%lf", &limit_lo) != 1)) {fprintf(stderr,GVERSION"%s\n",usage); exit(1);} i++; if ((i == argc) || (sscanf(argv[i], "%lf", &limit_hi) != 1)) {fprintf(stderr,GVERSION"%s\n",usage); exit(1);} if (limit_hi <= limit_lo) { fprintf(stderr,"Got limit low: %1.1e and high: %1.1e\n", limit_lo, limit_hi); pperror("High limit must be greater than low limit.\n"); } limitspec = TRUE; } else if (pm_keymatch(argv[i], "-seed", 2)) { if (seedspec) { pperror("already specified a random seed\n"); } i++; if ((i == argc) || (sscanf(argv[i], "%d", &rseed) != 1)) {fprintf(stderr,GVERSION"%s\n",usage); exit(1); } seedspec = TRUE; } else if (pm_keymatch(argv[i], "-peak", 3)) { if (peakspec) pperror("already specified peak location\n"); i++; if ((i == argc) || (sscanf(argv[i], "%f", &xfrac) != 1)) {fprintf(stderr,GVERSION"%s\n",usage); exit(1);} i++; if ((i == argc) || (sscanf(argv[i], "%f", &yfrac) != 1)) pperror("-peak requires two arguments, xfrac and yfrac.\n"); if ((xfrac<0)||(xfrac>1)||(yfrac<0)||(yfrac>1)) pperror("-peak arguments lie in the range [0..1]\n"); peakspec = TRUE; } else if (pm_keymatch(argv[i], "-name", 2)) { i++; if (i == argc) { fprintf(stderr,GVERSION"%s\n",usage); exit(1); } strcpy(output_filename,argv[i]); namespec = TRUE; } else if (pm_keymatch(argv[i], "-input", 2)) { i++; if (i == argc) {fprintf(stderr,GVERSION"%s\n",usage); exit(1);} strcpy(input_filename,argv[i]); in_namespec = TRUE; } else if (pm_keymatch(argv[i], "-craters", 1)) { include_craters = TRUE; i++; if ((i==argc) || (sscanf(argv[i],"%f",&ftmp) != 1) || (ftmp==0)) i--; /* oops, didn't have a numerical argument */ else { c_density = ftmp; /* crater density argument */ i++; /* check for 2nd arg. */ if ((i==argc) || (sscanf(argv[i],"%f",&ftmp) != 1) || (ftmp==0)) i--; /* oops, didn't have another numerical argument */ else ch_scale = ftmp; /* crater amplitude argument */ } } else if (pm_keymatch(argv[i], "-wrapoff",2)) { wrap = FALSE; } else if (pm_keymatch(argv[i], "-version",2)) { fprintf(stderr,GVERSION" - "GAUTHOR"\n"); exit(1); } else if (pm_keymatch(argv[i], "-type", 2)) { if (typespec) { pperror("already specified a filetype\n"); } i++; if (i == argc) {fprintf(stderr,GVERSION"%s\n",usage);exit(1);} if (pm_keymatch(argv[i], "TGA", 3)) pp_filetype = TGA; else if (pm_keymatch(argv[i], "PGM", 3)) pp_filetype = PGM; else if (pm_keymatch(argv[i], "PG8", 3)) pp_filetype = PG8; else if (pm_keymatch(argv[i], "OCT", 3)) pp_filetype = OCT; else if (pm_keymatch(argv[i], "MAT", 3)) pp_filetype = MAT; else if (pm_keymatch(argv[i], "PNG", 3)) pp_filetype = PNG; else pperror("Don't understand that file type.\n"); typespec = TRUE; } else { fprintf(stderr,"%s not understood.\n",argv[i]); fprintf(stderr,GVERSION"%s\n",usage); exit(1); } i++; } /* end while */ /* -------------------------------------------------- */ fac_warn(meshsize); /* warn if meshsize has large prime factor */ /* Set defaults when explicit specifications were not given. */ if (!dimspec) { if (adimspec) pperror("use of -adim requires also the -dim option\n"); dim[0] = 2.15; dscale[0] = 1.0; d_factors = 1; /* number of powers defined */ fracdim = 2.15; } if (!powerspec) { powscale = 1.2; } if (!meshspec) { meshsize = 128; } if (!typespec) { pp_filetype = TGA; } if (!namespec) { if ((pp_filetype == PGM) || (pp_filetype == PG8) ) strcpy(output_filename,"output.pgm"); else if (pp_filetype == TGA) strcpy(output_filename,"output.tga"); else if (pp_filetype == PNG) strcpy(output_filename,"output.png"); else if (pp_filetype == OCT) strcpy(output_filename,"output.oct"); else if (pp_filetype == MAT) strcpy(output_filename,"output.mat"); else strcpy(output_filename,"output.dat"); /* unknown type (!?) */ } if (!wrap && peakspec) { pperror("-wrapoff and -peak options are mutually exclusive!\n"); } screenxsize = meshsize; screenysize = meshsize; /* print out parameters so we know how we generated this file */ printf("Parameters: -dim %.2f -power %.2f -mesh %d -name %s\n", dim[0], powscale, meshsize, output_filename); if (adimspec) { for (i=1;i