Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!tut.cis.ohio-state.edu!rutgers!ucsd!nprdc!malloy
From: malloy@nprdc.arpa (Sean Malloy)
Newsgroups: comp.graphics
Subject: Re: Fractals Demo
Message-ID: <1858@skinner.nprdc.arpa>
Date: 3 May 89 13:53:35 GMT
References: <12725@pasteur.Berkeley.EDU> <20@steven.COM> <11@amcom.UUCP> <1472@cunixc.cc.columbia.edu>
Reply-To: malloy@nprdc.arpa (Sean Malloy)
Distribution: na
Organization: Navy Personnel R&D Center, San Diego
Lines: 440

In article <12725@pasteur.Berkeley.EDU>, dunross@petrus.berkeley.edu (Ian Dunross) writes:
>    I'm looking for a sample C program that will generate 3-D fractal
>    surfaces.  Ideally, it should run under X (X10 or X11), but any
>    sort of implementation will do.

Finally managing to find where I hid it, here is a program that
generates the stock triangular subdivision fractal terrain, with
command-line options to specify the number of subdivisions, the seed
to use to initialize the random number generator, and whether or not
to clip at 'sea level'. The program outputs PostScript code for
printing the terrain.

The syntax of the command is

	fractal [-c] [-n#] [-s#]

The -c option turns on clipping at sea level.

The -n option, with a number from 1 to 7, specifies the number of
subdivisions that are performed on the terrain. The more subdivisions,
the less coarse the generated terrain, but the more time the program
takes. If this option is not specified, the program assumes 7
subdivisions.

The -s option initializes the random number generator to the number
following it. When printed, the output from the program prints the
seed on the bottom of the page, so you can duplicate runs.

The program sends information on its progress to stderr, and the
PostScript code to stdout, so you can capture the output for later
printing by redirecting the output to a file, or simply pipe it to
your print spooler.

If you are using a random number generator that returns 16-bit random
values, you will have to change the value of RAND_MAX.

The code is derived from a BASIC program published in a PC-oriented
magazine a number of years ago; unfortunately, I don't recall the name
of the magazine.


 Sean Malloy					| "The proton absorbs a photon
 Navy Personnel Research & Development Center	| and emits two morons, a
 San Diego, CA 92152-6800			| lepton, a boson, and a
 malloy@nprdc.navy.mil				| boson's mate. Why did I ever
						| take high-energy physics?"

----8<--------8<--------8<--------8<--------8<--------8<--------8<--------8<----
#include <math.h>
#include <stdio.h>
#include <sys/types.h>
#include <sys/timeb.h>

#define pi 3.14159263589793238462
#define bool int
#define TRUE 1
#define FALSE 0
#define VERSION "1.0"
#define RAND_MAX (2147483647.0)


int mx, my;
int fractures = 7;
int mode = 0;
int seed = 0;
bool doflood = FALSE;
double l,ds;
double terrain[257][129];
double xx,yy,zz,x0;
double x2,y2,z2,rh,vt;
float lastx,lasty;

init()
{
	struct timeb grap;

	if (seed == 0)
	{
		ftime(&grap);
		seed = (int)((grap.time%100000)+grap.millitm);
	}
	(void)srand(seed);
	printf("%%!\n");
	printf("%% PS_Fractal - V%s; seed=%d\n\n",VERSION,seed);
	printf("90 rotate 0 -612 translate\n");
	printf("100 100 translate\n0 setgray\n");
	printf("\n/Helvetica-Bold findfont 10 scalefont setfont\n");
	printf("0 0 moveto\n");
	printf("(PS_Fractal - V%s; seed=%d) show\n\n",VERSION,seed);
}

main(argc,argv)
int argc;
char *argv[];
{
	char *c, *prognam;
	bool skip = FALSE;

	prognam = argv[0];
	while (--argc > 0 && (*++argv)[0] == '-')
	{
		for (c = argv[0]+1, skip=FALSE; (*c != '\0') && (!(skip)); c++)
			switch (*c)
			{
			case 'c':       /* set lisp output */
				doflood = TRUE;
				break;
			case 's':       /* set random seed */
				seed = atoi(&(*++c));
				skip = TRUE;
				break;
			case 'n':       /* set random seed */
				fractures = atoi(&(*++c));
				skip = TRUE;
				break;
			default:
			case '?':
				fprintf(stderr, "%s: Usage: %s [-c] [-n#] [-s#] \n",
					prognam,prognam);
				return (1);
			}
	}
	init();
    fracture();
    plot();
}

fracture()
{
    int i,n,sk,ib;
    double power();

    for (i = 0; i <= 257; i++)
	for (n = 0; n <= 129; n++)
	    terrain[i][n] = 0.0;
    ds = 2.0;
    for (n = 1; n <= fractures; n++)
	ds = ds + power(2.0, n-1);
    mx = ds - 1;
    my = mx / 2;
    rh = pi * 30.0/180.0;
    vt = rh * 1.2;
    for (i = 1; i<= fractures; i++)
    {
		fprintf(stderr,"Computing fracture level %d\n",i);
		l = 2000.0 / power(1.8,i);
		ib = mx / power(2.0,i);
		sk = ib * 2;
		assignx(sk,ib);
		assigny(sk,ib);
		assigndiag(sk,ib);
    }
}

assignx(sk,ib)
int sk,ib;
{
    int xe,ye;
    double d1,d2,newval;
    double buckle(), getarray();

    for (ye = 0; ye <= mx-1; ye += sk)
    {
		for (xe = ye+ib; xe <= mx; xe += sk)
		{
			d1 = getarray(xe-ib,ye);
			d2 = getarray(xe+ib,ye);
			newval = (d1+d2)/2.0 + buckle(l/2);
			putarray(xe,ye,newval);
		}
    }
}

assigny(sk,ib)
int sk,ib;
{
    int xe,ye;
    double d1,d2,newval;
    double buckle(), getarray();

    for (xe = mx; xe >= 1; xe -= sk)
    {
		for (ye = ib; ye <= xe; ye += sk)
		{
			d1 = getarray(xe,ye+ib);
			d2 = getarray(xe,ye-ib);
			newval = (d1+d2)/2.0 + buckle(l/2);
			putarray(xe,ye,newval);
		}
    }
}

assigndiag(sk,ib)
int sk,ib;
{
    int xe,ye;
    double d1,d2,newval;
    double buckle(), getarray();

    for (xe = 0; xe <= mx-1; xe += sk)
    {
		for (ye = ib; ye <= mx-xe; ye += sk)
		{
			d1 = getarray(xe+ye-ib,ye-ib);
			d2 = getarray(xe+ye+ib,ye+ib);
			newval = (d1+d2)/2.0 + buckle(l/2);
			putarray(xe+ye,ye,newval);
		}
    }
}

double buckle(max)
double max;
{
    int temp;

    return((max*(((double) rand())/(double)(RAND_MAX)))-(max/2.0));
}

double getarray(x,y)
int x,y;
{

    if (y>my)
		return(terrain[mx-x][mx+1-y]);
    else
		return(terrain[x][y]);
}

putarray(x,y,val)
int x,y;
double val;
{

    if (y>my)
		terrain[mx-x][mx+1-y] = val;
    else
		terrain[x][y] = val;
}

double power(x,n)
int n;
double x;
{
    double p;

    for (p = 1; n > 0; --n)
		p *= x;
    return(p);
}

plot()
{
    int ax,ay,ex,ey,c;
    short len,control,ierr;

/* plot initialization */
	fprintf(stderr,"Generating output\n");
	printf("\n0 setgray\n");
	mode = 1;
	x2 = 0.0;
	y2 = 0.0;
	z2 = 0.0;
	fprintf(stderr,"Writing first axis\n");
	for (ax = 0; ax <= mx; ax++)
	{
		x0 = -999.0;
		for (ay = 0; ay <= ax; ay++)
		{
			zz = getarray(ax,ay);
			yy = ((double) ay)/((double) mx) * 1000.0;
			xx = ((double) ax)/((double) mx) * 1000.0 - yy/2.0;
			putpoint();
		}
		printf("stroke\n");
	}
	printf("\n0 setgray\n");
	mode = 1;
	fprintf(stderr,"Writing second axis\n");
	for (ay = 0; ay <= mx; ay++)
	{
		x0 = -999.0;
		for (ax = ay; ax <= mx; ax++)
		{
			zz = getarray(ax,ay);
			yy = ((double) ay)/((double) mx) * 1000.0;
			xx = ((double) ax)/((double) mx) * 1000.0 - yy/2.0;
			putpoint();
		}
		printf("stroke\n");
	}
	printf("\n0 setgray\n");
	mode = 1;
	fprintf(stderr,"Writing third axis\n");
	for (ex = 0; ex <= mx; ex++)
	{
		x0 = -999.0;
			for (ey = 0; ey <= mx-ex; ey++)
		{
			zz = getarray(ex+ey,ey);
			yy = ((double) ey)/((double) mx) * 1000.0;
			xx = ((double) ex+ey)/((double) mx) * 1000.0 - yy/2.0;
			putpoint();
		}
		printf("stroke\n");
	}
	printf("\nshowpage\n");
	fprintf(stderr,"Done.\n");
}

putpoint()
{

    if (doflood == TRUE) flood();
    viewpoint();
}

flood()
{
    double w3,xt,yt,zt;
    double x3,y3,z3;

    if (x0 == -999.0)
    {
/* start of line in plot */
	x2 = xx;
	y2 = yy;
	z2 = zz;
	if (zz <= 0.0)
	{
/* point is below 'sea level' */
	    zz = 0.0;
		if (mode != 0)
			printf("stroke 0.5 setgray \n%.2f %.2f moveto\n",lastx,lasty);
		mode = 0;
	}
	else
	{
/* point is above 'sea level' */
		if (mode != 1)
			printf("stroke 0 setgray \n%.2f %.2f moveto\n",lastx,lasty);
		mode = 1;
	}
    }
    else if ((zz > 0.0) && (z2 > 0.0))
    {
/* both the new and previous points are above 'sea level' */
		if (mode != 1)
			printf("stroke 0 setgray \n%.2f %.2f moveto\n",lastx,lasty);
		mode = 1;
		x2 = xx;
		y2 = yy;
		z2 = zz;
    }
    else if ((zz <= 0) && (z2 <= 0))
    {
/* both the new and previous points are below 'sea level' */
		x2 = xx;
		y2 = yy;
		z2 = zz;
		zz = 0.0;
		if (mode != 0)
			printf("stroke 0.5 setgray \n%.2f %.2f moveto\n",lastx,lasty);
		mode = 0;
    }
    else
    {
/* a line between the points crosses 'sea level' */
		w3 = zz/(zz-z2);
		x3 = (x2-xx)*w3+xx;
		y3 = (y2-yy)*w3+yy;
		z3 = 0.0;
		xt = xx;
		yt = yy;
		zt = zz;
		if (zt < 0.0)
		{
/* going down into the water */
			xx = x3;
			yy = y3;
			zz = z3;
			viewpoint();
			if (mode != 0)
				printf("stroke 0.5 setgray \n%.2f %.2f moveto\n",lastx,lasty);
			mode = 0;
			xx = xt;
			yy = yt;
			zz = 0.0;
			x2 = xt;
			y2 = yt;
			z2 = zt;
		}
		else
		{
/* coming up out of the water */
			xx = x3;
			yy = y3;
			zz = z3;
			viewpoint();
			if (mode != 1)
				printf("stroke 0 setgray \n%.2f %.2f moveto\n",lastx,lasty);
			mode = 1;
			xx = xt;
			yy = yt;
			zz = 0.0;
			x2 = xt;
			y2 = yt;
			z2 = zt;
		}
    }
}

viewpoint()
{
    double r1,r2,rd,ra;
    double xp,yp,zp;
    float xplot,yplot;

/* rotate display for viewing */
/* tilt display for relief */
/* combined rotate/tilt algorithm */
    r1 = pi/3.0;	/* rotate angle */
    r2 = pi/3.5;	/* tilt angle */
    xp = xx;
    yp = yy;
    zp = zz / 2.0;
    xx = xp*cos(r1) - yp*sin(r1);
    yy = (xp*sin(r1) + yp*cos(r1))*cos(r2) + zp*sin(r2);
    zz = zp*cos(r2) - (xp*sin(r1) + yp*cos(r1))*sin(r2);
    xplot = ((float) xx/2.5 + 300.0);
    yplot = ((float) yy/2.5);
/* plot the point */
    if (x0 == -999.0)
		printf("%.2f %.2f moveto\n",xplot,yplot);
    else
		printf("%.2f %.2f lineto\n",xplot,yplot);
	lastx = xplot;
	lasty = yplot;
    x0 = yy;
}