chavez@harvard.ARPA (R. Martin Chavez) (07/02/85)
I've received a fair number of inquiries regarding
my fractal generator. I should emphasize that
(a) the code requires oodles of CPU cycles if
you want to compute at high resolution;
(b) your frame buffer should have at least
six bits of color;
(c) I don't have time to maintain and debug
the code, but I will be happy to answer questions; and
(d) you should probably get a copy of Mandelbrot's
"The Fractal Geometry of Nature" if you want to generate
some serious fractals (though I have shell scripts to get
you started.) The choice of polynomial iterates is critical;
Mandelbrot describes the theory in some detail. There
was also in article in an ACM SIGGRAPH (two [three?] years
ago) that outlines the theory of 3D fractals. (One
of the authors was Mandelbrot's co-worker at IBM; to
him is due the inspiration for this code, even though I
can't remember his name at the moment.)
I do not take responsibility for any bugs you
may find, nor will I sympathize with you when your
system manager asks you why you created a 50Mb file.
The code is here presented exactly as it was
one year ago. The file fractal.ms contains some
sketchy documentation; you'll have to compile the files
yourself, as I haven't included a makefile. The
shell script in 'quaternion.sh' once produced a
memorable picture that I immediately christened 'alien';
I will gladly send the pixel map to any disbelievers
out there.
The program 'display.c' produces two files,
picture.{col,img}. The .col file is a color table,
and the .img file is a pixel map (one byte per pixel)
scanning from left to right and from bottom to top.
The constants in graphics.h (e.g., SCANLINES and
SCANWIDTH) will need to be changed. You may also
want to experiment with colors other than blue.
I attach only two stipulations to the distribution
of this code: (1) keep my name associated with it, and
(2) send me any interesting modifications that you devise.
I will be happy to post a revised (and cleaned up) shar
file once people have hacked away a bit. This code is
a real MESS, but it wasn't too bad for a two-week academic
term project. (After 150+ hours of hacking last spring,
I haven't been able to bear the sight of this code until
now.)
Again, my apologies for the crufty code, but
it *did* work in the not-too-distant past. Enjoy!
R. Martin Chavez
Q-4 Nuclear Safeguards
Los Alamos National Laboratory (rmc@lanl.ARPA)
# -------------------- TEAR HERE ----------------------------
echo x - belt.in
sed s/.// >belt.in <<\End.belt.in
-400
-500000
--2. 2. -2. 2. 0. 6.28
-0. 0. 0. 100. 100. 100.
-0
End.belt.in
echo x - belt.sh
sed s/.// >belt.sh <<\End.belt.sh
-#! /bin/csh
-cd ~/fractal
-wind belt.f <belt.in >&belt.out
-display belt.f <dbelt
End.belt.sh
echo x - circle.c
sed s/.// >circle.c <<\End.circle.c
-#include "complex.h"
-#include <stdio.h>
-#include <ctype.h>
-
-#define MAX_ITERATION 100
-
-#define bool char
-#define TRUE (1)
-#define FALSE (0)
-#define CORNER 0
-#define EDGE 1
-#define FACE 2
-
-#define TESTED 0
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-#define NBUF 26
-
-#define round(x) ((int) ((x)+0.5))
-#define EMPTY (CUBE *) 0
-#define END (struct fifo *) 0
-
-typedef struct point {
- int x, y, z;
-} CUBE;
-
-struct fifo {
- CUBE *p;
- struct fifo *next;
-} *top, *bottom;
-
-FILE *grid, *fopen();
-int xres, yres, zres, stringency;
-float xlow, xmax, ylow, ymax, zlow, zmax, xinc, yinc, zinc;
-
-char *malloc();
-
-char *buffer_ring[NBUF];
-bool altered[NBUF];
-int block_ptr[NBUF];
-int current_buffer;
-
-int memrefs = 0, faults = 0, bpts = 0, ctests = 0,
- stack_size = 0, max_stack_size = 0;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- FILE *clear();
-
- CUBE *locate(), p;
- bool toggled();
-
- float x0, y0, z0, xf, yf, zf;
- float size;
- int i;
-
- printf( "Enter the x, y, and z resolutions.\n" );
- scanf( "%d %d %d", &xres, &yres, &zres );
- printf( "Enter the x, y, and z limits.\n" );
- scanf( "%f %f %f %f %f %f", &xlow, &xmax, &ylow, &ymax, &zlow, &zmax );
- printf( "Enter the coordinates of a segment that crosses the boundary.\n");
- scanf( "%f %f %f %f %f %f", &x0, &y0, &z0, &xf, &yf, &zf );
- printf( "Stringency? (0 -- corner, 1 -- edge, 2 -- face)\n" );
- scanf( "%d", &stringency );
-
- xinc = (xmax-xlow)/xres;
- yinc = (ymax-ylow)/yres;
- zinc = (zmax-zlow)/zres;
- size = xres*yres*zres/2 + 1;
- grid = clear( argv[1], size );
- top = (bottom = END );
-
- assert( locate( x0, y0, z0, xf, yf, zf) );
- printf( "Boundary points:\n\n" );
-
- for ( i = 0; i < NBUF; i++ )
- if (buffer_ring[i] != NULL && altered[i])
- {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
-
- if ( xres < 80 )
- for ( p.z = 0; p.z < zres; p.z++ )
- {
- printf( "\nz = %d\n", p.z );
- for ( p.y = yres-1; p.y > 0; p.y-- )
- {
- for ( p.x = 0; p.x < xres; p.x++ )
- {
- if (toggled( BOUNDARY, &p )) putchar('1');
- else putchar('0');
- }
- putchar('\n');
- }
- }
-
- printf( "\nStatistics:\n----------\n\n" );
- printf( "\tGrid size: %dx%dx%d = %d\n",
- xres, yres, zres, xres*yres*zres );
- printf( "\tMemory references:\t%d\n\tPage faults:\t\t%d\n",
- memrefs, faults );
- printf( "\tBoundary points:\t%d\n\tConvergence tests:\t%d\n",
- bpts, ctests );
- printf( "\tMaximum stack size:\t%d\n\tBitmap size (bytes):\t%d\n",
- max_stack_size, (int) size );
- fclose( grid );
-}
-
-FILE *clear( fname, size )
-char *fname;
-float size;
-{
- FILE *buffer;
- char *empty;
- int i;
-
- empty = malloc( 1024 );
- buffer = fopen( fname, "w+" );
- if ( buffer == NULL ) {
- printf( "fractal: can't open %s\n", fname );
- exit( 1 );
- }
-
- for ( i = 1; i <= size/1024 + 1; i++ )
- fwrite( empty, 1024, 1, buffer );
-
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ )
- {
- buffer_ring[i] = NULL;
- altered[i] = FALSE;
- block_ptr[i] = -1;
- }
-
- cfree( empty );
- return( buffer );
-}
-
-push( p )
-CUBE *p;
-{
- struct fifo *new;
-
- new = (struct fifo *) malloc( sizeof( struct fifo ) );
- new->next = END;
- new->p = p;
-
- set( INSPECTED, p );
- set( BOUNDARY, p );
-
- if (top == END) top = new;
- else bottom->next = new;
-
- bottom = new;
-
- ++stack_size;
- if ( stack_size > max_stack_size ) max_stack_size = stack_size;
-
- return;
-}
-
-CUBE *pop( )
-{
- CUBE *value;
- struct fifo *release;
-
- if ( top == END )
- return( EMPTY );
-
- value = top->p;
- release = top;
- top = top->next;
- cfree( release->p );
- cfree( release );
-
- --stack_size;
- return( value );
-}
-
-cartesian( p, x, y, z )
-CUBE *p;
-float *x, *y, *z;
-{
- *x = xlow + p->x*xinc;
- *y = ylow + p->y*yinc;
- *z = zlow + p->z*zinc;
-}
-
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
-
- nibble = xres*(yres*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
-
-set( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
- if ( item == BOUNDARY ) ++bpts;
- if ( (bpts > 0) && (bpts-(bpts/1000)*1000 == 0) )
- printf( "%d boundary points processed." );
-
- nibble = xres*(yres*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- {
- buffer_ring[i][byte] |= mask;
- altered[i] = TRUE;
- return;
- }
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- buffer_ring[i][byte] |= mask;
- block_ptr[i] = block;
- altered[i] = TRUE;
-
- return;
-}
-
-bool interior( p )
-CUBE *p;
-{
- bool hit, toggled(), out_of_range();
- float x, y, z;
-
- if ( out_of_range( p ) ) return( FALSE );
- cartesian( p, &x, &y, &z );
- if ( toggled( TESTED, p ) ) return( toggled( STATUS, p ) );
- else {
- ++ctests;
- set( TESTED, p );
- if ( hit=( sqrt( x*x + y*y + z*z ) <= 1.0) ) set( STATUS, p );
- }
-
- return( hit );
-}
-
-assert( p )
-CUBE *p;
-{
- CUBE *n, *neighbor(), *pop();
- bool interior();
- int index;
-
- push( p );
-
- while ( (p = pop( )) != EMPTY )
- {
- index = 1;
- while( n = neighbor( p, &index ) )
- if ( interior( n ) && !toggled( INSPECTED, n ) ) examine( n );
- else cfree( n );
- }
-}
-
-examine( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
- int index;
-
- index = 1;
- while( n = neighbor( p, &index ) ) {
- if( !interior( n ) ) { push( p ); cfree( n ); break; }
- cfree( n );
- }
- set( INSPECTED, p );
-}
-
-CUBE *neighbor( p, index )
-CUBE *p;
-int *index;
-{
- bool out_of_range();
- CUBE *n;
-
- n = (CUBE *) malloc( sizeof( CUBE ) );
- n->x= p->x; n->y = p->y; n->z = p->z;
-
- if ( (stringency == FACE && *index > 6) ||
- (stringency == EDGE && *index > 18) ||
- (stringency == CORNER && *index > 26) )
- return( EMPTY );
-
- switch( *index ) {
- case 1: { --n->x; break; }
- case 2: { ++n->x; break; }
- case 3: { --n->y; break; }
- case 4: { ++n->y; break; }
- case 5: { --n->z; break; }
- case 6: { ++n->z; break; }
-
- case 7: { --n->z; --n->y; break; }
- case 8: { ++n->z; --n->y; break; }
- case 9: { --n->x; --n->y; break; }
- case 10: { ++n->x; --n->y; break; }
- case 11: { --n->z; --n->x; break; }
- case 12: { ++n->z; --n->x; break; }
- case 13: { --n->z; ++n->x; break; }
- case 14: { ++n->z; ++n->x; break; }
- case 15: { --n->z; ++n->y; break; }
- case 16: { ++n->z; ++n->y; break; }
- case 17: { --n->x; ++n->y; break; }
- case 18: { ++n->x; ++n->y; break; }
-
- case 19: { ++n->x; ++n->y; ++n->z; break; }
- case 20: { ++n->x; ++n->y; --n->z; break; }
- case 21: { ++n->x; --n->y; ++n->z; break; }
- case 22: { ++n->x; --n->y; --n->z; break; }
- case 23: { --n->x; ++n->y; ++n->z; break; }
- case 24: { --n->x; ++n->y; --n->z; break; }
- case 25: { --n->x; --n->y; ++n->z; break; }
- case 26: { --n->x; --n->y; --n->z; break; }
- }
-
- *index += 1;
- if ( out_of_range( n ) )
- {
- cfree( n );
- return( neighbor( p, index ) );
- }
-
- return( n );
-}
-
-bool out_of_range( p )
-CUBE *p;
-{
- return ( (p->x < 0) || (p->x > xres) ||
- (p->y < 0) || (p->y > yres) ||
- (p->z < 0) || (p->z > zres) );
-}
-
-bool converges( p )
-CUBE *p;
-{
- complex *iterate, *hold;
- float x, y, z;
- int index;
-
- index = 1;
- iterate = zero();
- hold = zero();
- cartesian( p, &x, &y, &z );
- iterate->r = x;
- iterate->i = y;
-
- while( magnitude( iterate ) < BIG && index++ <= MAX_ITERATION )
- {
- mult( iterate, iterate, hold );
- subtract( iterate, hold, iterate );
- scale( iterate, z );
- }
-
- cfree( hold );
- cfree( iterate );
-
- return( index > MAX_ITERATION );
-}
-
-CUBE *bitmap( x, y, z )
-float x, y, z;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = round( (x-xlow)/xinc );
- new->y = round( (y-ylow)/yinc );
- new->z = round( (z-zlow)/zinc );
-
- return( new );
-}
-
-CUBE *bisect( left, right )
-CUBE *left, *right;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = floor( (left->x + right->x)/2. );
- new->y = floor( (left->y + right->y)/2. );
- new->z = floor( (left->z + right->z)/2. );
-
- return( new );
-}
-
-CUBE *locate( x0, y0, z0, xf, yf, zf )
-float x0, y0, z0, xf, yf, zf;
-{
- bool interior();
- CUBE *in, *out, *middle, *bitmap(), *bisect();
-
- in = bitmap( x0, y0, z0 );
- out = bitmap( xf, yf, zf );
-
- if ( ! (interior( in ) && !interior( out ) ) )
- {
- printf( "The given segment does not span the boundary. ");
- exit( 1 );
- }
- do {
- middle = bisect( in, out );
- if (interior( middle ) )
- {
- cfree( in );
- in = middle;
- }
- else {
- cfree( out );
- out = middle;
- }
- } while ( abs( in->x - out->x ) > 1 ||
- abs( in->y - out->y ) > 1 ||
- abs( in->z - out->z ) > 1 );
- return( in );
-}
End.circle.c
echo x - complex.c
sed s/.// >complex.c <<\End.complex.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The algorithm is by no means complete, and there
- are many known bugs in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- in CS175 (Computer Graphics) in the spring of 1984.
-
- For complete documentation, please refer to the companion file
- ~chavez/fractal/fractal.doc on hu-midgard on harvard.
-
- The routines in "complex.c" provide a general interface for
- algebraic manipulation of complex numbers. In order to use
- the package, "#include complex.h" and reference the object
- file "complex.o" in the C compilation.
-
-*/
-
-#include <math.h>
-
-#define BIG 1.e18
-
-typedef struct complex_number {
- float r, i;
-} complex;
-
-float magnitude( z )
-complex *z;
-{
- if ( fabs( z->r ) < BIG && fabs( z->i ) < BIG )
- return( sqrt( z->r*z->r + z->i*z->i ) );
- else return( BIG );
-}
-
-complex *zero( )
-{
- complex *value;
-
- value = (complex *) malloc( sizeof( complex ) );
-
- value->r = 0.0;
- value->i = 0.0;
-
- return( value );
-}
-
-mult( a, b, c )
-complex *a, *b, *c;
-{
- float x, y;
-
- x = a->r*b->r - a->i*b->i;
- y = a->i*b->r + a->r*b->i;
-
- c->r = x;
- c->i = y;
-}
-
-scale( z, n )
-complex *z;
-float n;
-{
- z->r = n*z->r;
- z->i = n*z->i;
-}
-
-add( a, b, c )
-complex *a, *b, *c;
-{
- c->r = a->r + b->r;
- c->i = a->i + b->i;
-}
-
-subtract( a, b, c )
-complex *a, *b, *c;
-{
- c->r = a->r - b->r;
- c->i = a->i - b->i;
-}
-
-move( a, b )
-complex *a, *b;
-{
- b->r = a->r;
- b->i = a->i;
-}
End.complex.c
echo x - complex.h
sed s/.// >complex.h <<\End.complex.h
-#include <math.h>
-
-#define BIG 1.0e18
-
-typedef struct complex_number {
- float r, i;
-} complex;
-
-float magnitude();
-complex *zero();
End.complex.h
echo x - display.c
sed s/.// >display.c <<\End.display.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The algorithm is by no means complete, and there
- are many known bugs in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- in CS175 (Computer Graphics) in the spring of 1984.
-
- For complete documentation, please refer to the companion file
- ~chavez/fractal/fractal.doc on hu-midgard on harvard.
-
- This program, "display.c", is the image processor for the fractal
- display system. The command-line invocation takes one argument,
- which is the name of the fractal grid file (as produced by
- fractal.c and companions.) Compile the display processor
- with the csh command
-
- cc display.c -lm -o display
-
- The input parameters are self-explanatory. Display.c generates
- pixel and color files suitable for viewing on the JUPITER system.
-
-*/
-
-#include "graphics.h"
-#include <math.h>
-#include <stdio.h>
-
-#define MAX_COLORS 256 /* size of the color table */
-#define bool char /* my favorite data type */
-#define TESTED 0 /* bits for the fractal grid */
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-
-#define NBUF 26 /* number of buffers in the virtual scheme */
-#define KI 2. /* constant coefficients for Lambert shading */
-#define KR 1.
-#define legal(x) (((x)>=0) && ((x)<size)) /* legal grid index? */
-#define round(x) ((int) (x+0.5))
-
-typedef short unsigned int twobytes; /* data type for z buffer */
-
-typedef struct scaled_color { /* data type for the color table */
- unsigned int r, g, b;
-} SCOLOR;
-
-typedef struct fractal_surface { /* display primitives */
- unsigned int x, y, z;
-} CUBE;
-
-FILE *grid; /* fractal data file */
-char *malloc();
-
-char *buffer_ring[NBUF]; /* various storage areas for the */
-bool altered[NBUF]; /* virtual memory scheme */
-int block_ptr[NBUF];
-int current_buffer;
-
-int res, yon, size; /* resolution, yon plane, res^3 */
-char *pixel; /* raster data area */
-twobytes *z_buffer; /* z buffer */
-int color_index = 0; /* index into color table */
-SCOLOR *color_table[MAX_COLORS]; /* color table */
-SCOLOR *shade; /* the fractal's base color */
-
-main( argc, argv )
-int argc;
-char *argv[];
-{
- char *malloc(), *calloc();
- SCOLOR *tweak(), *make_color();
- CUBE *project();
- bool toggled();
- char lighting();
-
- int view; /* see fractal.doc for more */
- int i, j, bit, counter, x_offset, y_offset;
- CUBE p, *image;
-
- FILE *pixel_file, *color_file; /* output files */
- char fname[64], pixel_fname[64], color_fname[64];
-
- /* read the command line and open the input data file */
- grid = fopen( argv[1], "r" );
- if ( grid == NULL ) {
- printf( "display: can't open %s\n", argv[1] );
- exit( 1 );
- }
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ ) /* initialize the virtual buffers */
- {
- altered[i] = FALSE;
- block_ptr[i] = -1;
- buffer_ring[i] = NULL;
- }
-
- /* setup dialogue */
- printf( "Viewing direction?\n" );
- scanf( "%d", &view );
- printf( "Grid size?\n" );
- scanf( "%d", &res );
- printf( "File name for fractal image?\n" );
- scanf( "%s", fname );
- sprintf( pixel_fname, "%s.img", fname );
- sprintf( color_fname, "%s.col", fname );
- pixel_file = fopen( pixel_fname, "w" );
- color_file = fopen( color_fname, "w" );
- if ( pixel_file == NULL || color_file == NULL ) {
- printf( "display: Couldn't open %s files\n", fname );
- exit( 1 );
- }
-
- /* we don't use the whole screen -- stay in the middle */
- x_offset = (SCANWIDTH-res)/2;
- y_offset = (SCANLINES-res)/2;
- size = res*res;
- yon = res+1;
-
- /* allocate the pixel array and z buffer */
- pixel = malloc( SCANLINES*SCANWIDTH );
- z_buffer = (twobytes *)
- calloc( size, sizeof( twobytes ) );
- color_table[0] = make_color( 20, 20, 20 );
-
- /* specify the base shading for the fractal */
- shade = make_color( 0, 0, 100 );
-
- /* initialize the z buffer to the yon clipping plane */
- for ( i = 0; i < size; i++ )
- z_buffer[i] = yon;
-
- /* loop through the grid and toss into the z buffer */
- for ( p.z = 0; p.z < res; ++p.z )
- for ( p.y = 0; p.y < res; ++p.y )
- for ( p.x = 0; p.x < res; ++p.x )
- if ( toggled( BOUNDARY, &p ) )
- {
- image = project( &p, view );
- bit = image->y*res + image->x;
- if ( image->z < z_buffer[bit] )
- z_buffer[bit] = image->z;
- cfree( image );
- }
-
- /* shade the pixels that came out on top */
- for ( p.x = 0; p.x < res; ++p.x )
- for ( p.y = 0; p.y < res; ++p.y )
- if ( z_buffer[bit = p.y*res + p.x] != yon )
- pixel[(y_offset+p.y)*SCANWIDTH+x_offset+p.x] = lighting( bit );
-
- /* dump the pixel map and color table */
- printf( "Write out the pixel map...\n" );
- for ( i = 0; i <= color_index; i++ )
- fprintf( color_file, "%d %d %d\n",
- color_table[i]->r, color_table[i]->g, color_table[i]->b );
- fclose( color_file );
-
- fwrite( pixel, sizeof( char ), SCANLINES*SCANWIDTH, pixel_file );
- fclose( pixel_file );
-
- printf( "All done.\n" );
-}
-
-/* perform one of six orthogonal projections on a grid point */
-/* Standard view: center of interest at center of grid,
- view reference point at z = -infinity */
-CUBE *project( p, view )
-CUBE *p;
-int view;
-{
- CUBE *image;
-
- image = (CUBE *) malloc( sizeof( CUBE ) );
-
- image->x = p->x; image->y = p->y; image->z = p->z;
-
- switch ( view ) {
-
- /* x = +infinity */
- case 2: {
- image->x = p->z; image->z = res-p->x; break; }
-
- /* z = +infinity */
- case 3: {
- image->x = res-p->x; image->z = res-p->z; break; }
-
- /* x = -infinity */
- case 4: {
- image->x = res-p->z; image->z = p->x; break; }
-
- /* y = +infinity */
- case 5: {
- image->y = p->z; image->z = res-p->y; break; }
-
- /* y = -infinity */
- case 6: {
- image->y = res-p->z; image->z = p->y; break; }
-
- }
-
- return( image );
-}
-
-/* examine a pixel in the z buffer and compute the shading */
-char lighting( bit )
-int bit;
-{
- char lookup();
- SCOLOR *intensity();
-
- float xnorm, ynorm, gamma;
- int y, ybit, x, xbit;
-
- /* gamma is the z-component of the surface normal and therefore
- the dot product of the incident light and the normal */
- gamma = 0.0;
- /* look at each neighboring pixel */
- for ( y = -1; y <= 1; y += 2 )
- {
- ybit = bit + y*res;
- if ( !legal( ybit ) || z_buffer[ybit] == yon )
- break; /* don't shade border pixels */
- for ( x = -1; x <= 1; x += 2 )
- {
- xbit = bit + x;
- if ( !legal( xbit ) || z_buffer[xbit] == yon )
- break;
- /* set up the cross product */
- xnorm = (z_buffer[bit]-z_buffer[xbit])*y;
- ynorm = (z_buffer[bit]-z_buffer[ybit])*x;
- gamma += 1./sqrt( xnorm*xnorm + ynorm*ynorm + 1. );
- }
- if ( gamma == 0.0 ) break;
- }
-
- /* go compute the intensity based on the average normal
- z-component and check out the color table */
- return( lookup( intensity( gamma/4., z_buffer[bit] ) ) );
-}
-
-/* turn gamma and depth into a scaled color */
-SCOLOR *intensity( gamma, depth )
-float gamma;
-twobytes depth;
-{
- SCOLOR *illumination;
- float scaled_depth, it, red, green, blue;
-
- illumination = (SCOLOR *) malloc( sizeof( SCOLOR ) );
-
- scaled_depth = ( (float) depth ) / res;
- it = (KI/(KR+scaled_depth)) * gamma;
- if (it < 0.0) it = 0.0;
-
- red = it * shade->r;
- green = it * shade->g;
- blue = it * shade->b;
-
- /* clamp at maximum saturation */
- illumination->r = 255;
- illumination->b = 255;
- illumination->g = 255;
-
- if ( red < 255 ) illumination->r = red;
- if ( green < 255 ) illumination->g = green;
- if ( blue < 255 ) illumination->b = blue;
-
- return( illumination );
-}
-
-/* look up a scaled color in the color table;
- make up a new entry if it's not already there */
-char lookup( illumination )
-SCOLOR *illumination;
-{
- int i;
-
- for ( i = 0; i <= color_index; i++ )
- if ( (illumination->r == color_table[i]->r) &&
- (illumination->g == color_table[i]->g) &&
- (illumination->b == color_table[i]->b) ) return( (char) i );
-
- color_table[i] = illumination;
-
- color_index += 1;
- return( (char) i );
-}
-
-/* turn a COLOR into an SCOLOR */
-SCOLOR *tweak( c )
-COLOR *c;
-{
- SCOLOR *new;
-
- new = (SCOLOR *) malloc( sizeof( SCOLOR ) );
- new->r = 255;
- new->g = 255;
- new->b = 255;
-
- if ( c->red <= 1.0 ) new->r = c->red*255;
- if ( c->green <= 1.0 ) new->g = c->green*255;
- if ( c->blue <= 1.0 ) new->b = c->blue*255;
-
- return( new );
-}
-
-/* take r, g, b and make a new SCOLOR record */
-SCOLOR *make_color( r, g, b )
-unsigned int r, g, b;
-{
- SCOLOR *new;
-
- new = (SCOLOR *) malloc( sizeof( SCOLOR ) );
-
- new->r = r;
- new->g = g;
- new->b = b;
-
- return( new );
-}
-
-/* check to see whether a boundary point lives here */
-/* dfractal.c has fuller documentation of toggled */
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
End.display.c
echo x - dragon.c
sed s/.// >dragon.c <<\End.dragon.c
-#include "complex.h"
-
-complex *mu;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- int n = 0;
- complex *seed;
-
- seed = zero();
- mu = zero();
-
- printf( "\nMu? ");
- scanf( "%f %f", &mu->r, &mu->i );
-
- printf( "\nSeed? " );
- scanf( "%f %f", &seed->r, &seed->i );
- while (1) {
- printf( "\nIteration? " );
- scanf( "%d", &n );
- iterate( n, seed );
- printf( "\nResult: %.10f + %.10fi\n", seed->r, seed->i );
- }
-}
-
-iterate( n, seed )
-int n;
-complex *seed;
-{
- int i;
- double pi, theta;
- complex *one;
-
- one = zero();
- for ( i=1; i<=n; i++ ) {
- mult( seed, seed, seed );
- subtract( seed, mu, seed );
- }
-}
End.dragon.c
echo x - fractal.c
sed s/.// >fractal.c <<\End.fractal.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The information is by no means complete, and there
- are some imperfections in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- in CS175 (Computer Graphics) in the spring of 1984.
-
- The file "fractal.c" includes a boolean convergence evaluator
- that generates a parametric family of fractal curves, where
- lambda (calculated as the z coordinate) takes on real values
- from one to four.
-
- Compile "fractal.c" with the csh command
-
- cc fractal.c complex.o -lm -o fractal
-
- and invoke it with one argument, the name of the fractal grid
- file. Fractal.c generates data suitable for use with the
- "display.c" image processor.
-
-*/
-
-#include "complex.h"
-#include <stdio.h>
-#include <ctype.h>
-
-#define MAX_ITERATION 100
-#define EPSILON 1.e-4
-
-#define bool char
-#define TRUE (1)
-#define FALSE (0)
-#define CORNER 0
-#define EDGE 1
-#define FACE 2
-
-#define TESTED 0
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-#define NBUF 26
-
-#define round(x) ((int) ((x)+0.5))
-#define pack(p) ((long) (res*(res*p->z+p->y)+p->x))
-#define mod(x,y) ((x)-((x)/(y)*(y)))
-#define EMPTY (CUBE *) 0
-
-typedef struct point {
- int x, y, z;
-} CUBE;
-
-long *stack;
-
-FILE *grid, *fopen();
-int res, res2, stringency;
-long size;
-float xlow, xmax, ylow, ymax, zlow, zmax, xinc, yinc, zinc;
-
-char *malloc();
-
-char *buffer_ring[NBUF];
-bool altered[NBUF];
-int block_ptr[NBUF];
-int current_buffer;
-
-complex *hold, *iterate;
-
-int memrefs = 0, faults = 0, bpts = 0, ctests = 0,
- stack_ptr = 0, max_stack_size = 0;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- FILE *clear();
-
- CUBE *locate(), p;
- bool toggled();
-
- float x0, y0, z0, xf, yf, zf;
- int i;
-
- printf( "Resolution?\n" );
- scanf( "%d", &res );
- printf( "Maximum stack size?\n" );
- scanf( "%d", &max_stack_size );
- printf( "Enter the x, y, and z limits.\n" );
- scanf( "%f %f %f %f %f %f", &xlow, &xmax, &ylow, &ymax, &zlow, &zmax );
- printf( "Enter the coordinates of a segment that crosses the boundary.\n");
- scanf( "%f %f %f %f %f %f", &x0, &y0, &z0, &xf, &yf, &zf );
- printf( "Stringency? (0 -- corner, 1 -- edge, 2 -- face)\n" );
- scanf( "%d", &stringency );
-
- xinc = (xmax-xlow)/res;
- yinc = (ymax-ylow)/res;
- zinc = (zmax-zlow)/res;
- res2 = res*res;
- size = res*res2/2 + 1;
- grid = clear( argv[1], size );
- stack = (long *) calloc( max_stack_size, sizeof( long ) );
- if ( stack == NULL ) {
- printf( "fractal: couldn't create stack\n" );
- exit( 1 );
- }
-
- hold = zero();
- iterate = zero();
-
- assert( locate( x0, y0, z0, xf, yf, zf) );
- printf( "Boundary points:\n\n" );
-
- for ( i = 0; i < NBUF; i++ )
- if (buffer_ring[i] != NULL && altered[i])
- {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
-
- if ( res < 80 )
- for ( p.z = 0; p.z < res; p.z++ )
- {
- printf( "\nz = %d\n", p.z );
- for ( p.y = res-1; p.y > 0; p.y-- )
- {
- for ( p.x = 0; p.x < res; p.x++ )
- {
- if (toggled( BOUNDARY, &p )) putchar('1');
- else putchar('0');
- }
- putchar('\n');
- }
- }
-
- printf( "\nStatistics:\n----------\n\n" );
- printf( "\tGrid size: %dx%dx%d = %d\n",
- res, res, res, res*res*res );
- printf( "\tMemory references:\t%d\n\tPage faults:\t\t%d\n",
- memrefs, faults );
- printf( "\tBoundary points:\t%d\n\tConvergence tests:\t%d\n",
- bpts, ctests );
- printf( "\tBitmap size (bytes):\t%d\n",
- size );
- fclose( grid );
-}
-
-FILE *clear( fname, size )
-char *fname;
-long size;
-{
- FILE *buffer;
- char *empty;
- int i;
-
- empty = malloc( 1024 );
- buffer = fopen( fname, "w+" );
- if ( buffer == NULL ) {
- printf( "fractal: can't open %s\n", fname );
- exit( 1 );
- }
-
- for ( i = 1; i <= size/1024 + 1; i++ )
- fwrite( empty, 1024, 1, buffer );
-
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ )
- {
- buffer_ring[i] = NULL;
- altered[i] = FALSE;
- block_ptr[i] = -1;
- }
-
- cfree( empty );
- return( buffer );
-}
-
-cartesian( p, x, y, z )
-CUBE *p;
-float *x, *y, *z;
-{
- *x = xlow + p->x*xinc;
- *y = ylow + p->y*yinc;
- *z = zlow + p->z*zinc;
-}
-
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
-
-set( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
- if ( item == BOUNDARY ) {
- ++bpts;
- if ( (bpts > 0) && (bpts-(bpts/1000)*1000 == 0) )
- printf( "\n%d boundary points processed.\n", bpts );
- }
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- {
- buffer_ring[i][byte] |= mask;
- altered[i] = TRUE;
- return;
- }
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- buffer_ring[i][byte] |= mask;
- block_ptr[i] = block;
- altered[i] = TRUE;
-
- return;
-}
-
-bool interior( p )
-CUBE *p;
-{
- bool hit, converges(), toggled(), out_of_range();
-
- if ( out_of_range( p ) ) return( FALSE );
- if ( toggled( TESTED, p ) ) return( toggled( STATUS, p ) );
- else {
- ++ctests;
- set( TESTED, p );
- if ( (hit=converges( p )) ) set( STATUS, p );
- }
-
- return( hit );
-}
-
-assert( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
-
- int index;
- long nibble;
-
- while ( stack_ptr >= 0 )
- {
- index = 1;
- while( n = neighbor( p, &index ) )
- if ( interior( n ) && !toggled( INSPECTED, n ) ) examine( n );
- else cfree( n );
-
- if ( stack_ptr > 0 ) {
- nibble = stack[--stack_ptr];
- p->z = nibble/res2;
- p->y = nibble/res - res*p->z;
- p->x = nibble - (res*(res*p->z + p->y));
- }
- else --stack_ptr;
- }
-
- cfree( p );
-}
-
-examine( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
- int index;
-
- index = 1;
- while( n = neighbor( p, &index ) ) {
- if( !interior( n ) ) {
- stack[stack_ptr++] = pack( p );
- set( INSPECTED, p );
- set( BOUNDARY, p );
- cfree( n );
- break;
- }
- cfree( n );
- }
-
- set( INSPECTED, p );
-
- if ( n != EMPTY ) cfree( n );
- cfree( p );
-}
-
-CUBE *neighbor( p, index )
-CUBE *p;
-int *index;
-{
- bool out_of_range();
- CUBE *n;
-
- n = (CUBE *) malloc( sizeof( CUBE ) );
- n->x= p->x; n->y = p->y; n->z = p->z;
-
- if ( (stringency == FACE && *index > 6) ||
- (stringency == EDGE && *index > 18) ||
- (stringency == CORNER && *index > 26) )
- return( EMPTY );
-
- switch( *index ) {
- case 1: { --n->x; break; }
- case 2: { ++n->x; break; }
- case 3: { --n->y; break; }
- case 4: { ++n->y; break; }
- case 5: { --n->z; break; }
- case 6: { ++n->z; break; }
-
- case 7: { --n->z; --n->y; break; }
- case 8: { ++n->z; --n->y; break; }
- case 9: { --n->x; --n->y; break; }
- case 10: { ++n->x; --n->y; break; }
- case 11: { --n->z; --n->x; break; }
- case 12: { ++n->z; --n->x; break; }
- case 13: { --n->z; ++n->x; break; }
- case 14: { ++n->z; ++n->x; break; }
- case 15: { --n->z; ++n->y; break; }
- case 16: { ++n->z; ++n->y; break; }
- case 17: { --n->x; ++n->y; break; }
- case 18: { ++n->x; ++n->y; break; }
-
- case 19: { ++n->x; ++n->y; ++n->z; break; }
- case 20: { ++n->x; ++n->y; --n->z; break; }
- case 21: { ++n->x; --n->y; ++n->z; break; }
- case 22: { ++n->x; --n->y; --n->z; break; }
- case 23: { --n->x; ++n->y; ++n->z; break; }
- case 24: { --n->x; ++n->y; --n->z; break; }
- case 25: { --n->x; --n->y; ++n->z; break; }
- case 26: { --n->x; --n->y; --n->z; break; }
- }
-
- *index += 1;
- if ( out_of_range( n ) )
- {
- cfree( n );
- return( neighbor( p, index ) );
- }
-
- return( n );
-}
-
-bool out_of_range( p )
-CUBE *p;
-{
- return ( (p->x < 0) || (p->x > res) ||
- (p->y < 0) || (p->y > res) ||
- (p->z < 0) || (p->z > res) );
-}
-
-CUBE *bitmap( x, y, z )
-float x, y, z;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = round( (x-xlow)/xinc );
- new->y = round( (y-ylow)/yinc );
- new->z = round( (z-zlow)/zinc );
-
- return( new );
-}
-
-CUBE *bisect( left, right )
-CUBE *left, *right;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = floor( (left->x + right->x)/2. );
- new->y = floor( (left->y + right->y)/2. );
- new->z = floor( (left->z + right->z)/2. );
-
- return( new );
-}
-
-CUBE *locate( x0, y0, z0, xf, yf, zf )
-float x0, y0, z0, xf, yf, zf;
-{
- bool interior();
- CUBE *in, *out, *middle, *bitmap(), *bisect();
-
- in = bitmap( x0, y0, z0 );
- out = bitmap( xf, yf, zf );
-
- if ( ! (interior( in ) && !interior( out ) ) )
- {
- printf( "The given segment does not span the boundary. ");
- exit( 1 );
- }
- do {
- middle = bisect( in, out );
- if (interior( middle ) )
- {
- cfree( in );
- in = middle;
- }
- else {
- cfree( out );
- out = middle;
- }
- } while ( abs( in->x - out->x ) > 1 ||
- abs( in->y - out->y ) > 1 ||
- abs( in->z - out->z ) > 1 );
- return( in );
-}
-
-bool converges( p )
-CUBE *p;
-{
- float x, y, z;
- int index;
-
- index = 1;
- cartesian( p, &x, &y, &z );
- iterate->r = x;
- iterate->i = y;
-
- while( magnitude( iterate ) < BIG && index++ <= MAX_ITERATION )
- {
- mult( iterate, iterate, hold );
- subtract( iterate, hold, iterate );
- scale( iterate, z );
- }
-
- return( index > MAX_ITERATION );
-}
End.fractal.c
echo x - fractal.ms
sed s/.// >fractal.ms <<\End.fractal.ms
-.sp 2
-.ce 2
-Computer Science 175: Computer Graphics
-Term Project, Spring 1984
-.sp 2
-.ce
-Author: R. Martin Chavez
-.sp 5
-.SH
-General Introduction
-.PP
-Fractals are defined as geometrical objects for which the Hausdorff
-dimension, a measure of linear complexity, exceeds the topological
-dimension. Fractals may be given either a stochastic or a purely
-analytic representation. In any case, B. Mandelbrot has shown that
-fractals can mirror the complexities of natural phenomena in a
-rigorously mathematical fashion. See [Mandelbrot 1981] for details.
-.PP
-In the present case, fractals are derived from the iteration of
-a quadratic function over the ring of complex (or in the general
-case, quaternionic) numbers. More specifically, the transformation
-.DS
- T: z -> lambda*z(1-z)
-which may also be written (by a change of variables)
- T': x -> x*x - mu
-.DE
-is iterated for various starting values of z to obtain a sequence
-T(z), T(T(z)), T(T(T(z))), and so on. Depending on the values
-of z and the parameter lambda/mu, the sequence may either
-(a) diverge to infinity; (b) converge to a limit point; or
-(c) converge to a stable limit cycle of period w such that
-.DS
- lim T (z) = z
- n->inf nw mu
-.DE
-Mandelbrot has shown how to determine parameters lambda and
-mu so as to generate fractal surfaces of a particular kind.
-For the present, it suffices to note that the roots of the
-equation
-.DS
- T (0) = 0
- w
-.DE
-define values of mu for which all z iterates either diverge
-to infinity or else give rise to a limit cycle of period w.
-Neighboring values of mu (those that lie within the same
-"molecule" in Mandelbrot's mu-maps) also produce limit cycles
-of the same period, but the fractal separators S that partition
-the set of convergent and divergent points may look quite different.
-.SH
-Three-dimensional Fractal Displays
-.PP
-There are two ways in which to generalize Mandelbrot's analytical
-fractals to three dimensions. In the first case, x and y
-coordinates for graphic display are composed into a single
-complex number z; the z display coordinate is then used to
-represent an axis along which lambda varies continuously.
-The resulting figure is a 3d rendering of a parametric family
-of self-squared curves.
-The pixel files drape.col and drape.img in ~chavez/fractal
-define a view towards the origin from y = +infinity.
-In addition, I have studies parametric families given
-by variation of the phase of the parameter lambda
-while keeping its magnitude constant and equal to one.
-The files belt.{col,img} record the results.
-.SH
-Iteration over the Quaternions
-.PP
-In the second case, the quadratic transformation T is iterated
-over the quaternions.
-The files alien.col and alien.img define a view (in standard
-left-handed coordinates) of the fractal separator for
-points attracted to a cycle of length two.
-In brief, the following properties of
-the quaternions are needed to understand the present implementation:
-.IP [1]
-The quaternions constitute a four-dimensional division algebra.
-.IP [2]
-Quaternions are generalized complex numbers, with four components:
-real, i, j, and k.
-.IP [3]
-The multiplication of quaternions is non-commutative. When
-forming the product of two quaternions, the order of the factors must
-be preserved. The quantities i, j, and k multiply as follows:
-.DS
- i x i = -1 j x j = -1 k x k = -1
- i x j = k k x i = j j x k = i
- j x i = -k i x k = -j k x j = -i
-.DE
-.PP
-For the purposes of graphic display, the x, y, and z coordinates of
-each point in the model domain are mapped onto the three-dimensional
-subspace spanned by the quaternions 1, i, and j. In addition,
-I have discovered (mostly by luck) that the quaternionic iteration
-i*(z^2+mu) exhibits the same mu-dependency as the two-dimensional
-complex iteration. More specifically, a value of mu that produces
-an attractive cycle of length 4 in the complex iteration z^2-mu will
-also produce the same cycle
-when iterated over the quaternions in the form given above.
-.SH
-Image Computation and the Fractal Problem
-.PP
-Various source routines (qfractal.c, fractal.c, dfractal.c, wind.c)
-are defined for the generation of fractal surfaces. In each case,
-fractal surfaces are traced out in a 3d grid with integer coordinates
-onto which is mapped a subspace of Cartesian space. The programs
-differ only in the boolean function convergence( p ), which
-calculates the quadaratic iteration at grid point p. Dfractal.c
-has by far the most sophisticated convergence algorithm; cycles
-of arbitrary length are efficiently detected. The notion of
-convergence pertains to the fractal definition, not to
-the tracking algorithm; all that is required for successful
-determination of the boundary is a complete and unique mapping
-of all grid points into two disjoint sets that I have arbitrarily
-labeled "exterior" and "interior."
-.PP
-In general, four
-bits of information must be maintained for each grid point: one (TESTED)
-to tell whether the grid point has been tested for convergence, one (STATUS)
-to tell whether a TESTED point converges or diverges, one (INSPECTED) to
-tell whether a point has been inspected as a possible boundary point, and
-one (BOUNDARY) to indicate those grid points that lie along the connected
-boundary of the fractal domain. There is no limit on the number of grid
-points that may be used for tracking fractal surfaces. A multiply-buffered
-virtual addressing scheme allows for rapid access of grid points stored
-on a direct-access device; up to 26 buffers may be maintained in main
-memory at a given time, thereby minimizing page thrashing.
-.SH
-Boundary-Tracking Algorithm
-.PP
-The boundary of a fractal surface is determined in the following manner:
-.PP
-The user must first specify a line segment that crosses the fractal
-boundary. The tracer bisects the line segment until the boundary
-intersection is located. (When limit cycles are present, all
-interior points are those that are attracted to the same cycle
-as the interior point of the starting line segment.)
-For each newly discovered boundary
-points, all neighbors are listed according to a specified criterion
-of stringency. (Depending on the nature of the figure, more or less
-restrictive notions of adjacency are appropriate. For most of the surfaces
-in this project, cubes that share at least a corner are considered to
-be neighbors.) Now examine each neighbor in turn; those neighbors
-that lie in the interior of the fractal domain and yet adjoin at least
-one point NOT in the domain are then added to the list of newly discovered
-boundary points. The process continues until all neighbors of all the
-discovered boundary points have been examined.
-.PP
-Clearly, the algorithm lends itself easily to a recursive implementation.
-Nonetheless, computation of fractal surfaces at very high resolutions
-will often require extreme levels (~10K) of recursion. So as to minimize
-system overhead, I have implemented a LIFO stack on which newly discovered
-boundary points are kept until all neighboring cubes have been investigated.
-The search procedure terminates when the LIFO stack empties. Inasmuch
-as high-resolution images may require the generation of huge stacks,
-economy of space is of the utmost importance. Each grid-point is
-encoded into one long integer and stored into a stack that has been
-implemented as a simple linear array. In the current implementation,
-the stack must fit into virtual memory; that requirement effectively
-limits the number of stacked boundary points to 1,000,000. Grid
-resolutions that would cause the stack to explode would most likely
-take forever to compute anyway.
-.SH
-Performance Summary
-.PP
-Although I have optimized the search procedure to the best of my ability,
-it remains disappointingly slow. Tracking of a quaternionic surface on
-a 400x400x400 grid is likely to take up to 10 hours of CPU time on a
-VAX-11/780 with floating-point accelerator. The problem lies in the
-complexity of the task itself, and sheer processing speed is the absolute
-limiting factor. The convergence test has been carefully scrutinized
-so as to limit multiplication of complex numbers and quaternions to
-an absolute minimum, and the boundary-tracking method usually requires
-the application of a convergence test on only five times as many points as
-eventually end up on the boundary. (In the optimum case, only boundary
-points would be tested for convergence.) Inasmuch as the
-tracking algorithm tends to hover in
-a localized area while computing the boundary, storage of the grid on
-disk incurs very little overhead; performance summaries generated by
-the computation algorithms indicate that most references to points on
-the evaluation grid do not incur a page fault, e.g. the virtual page-management
-scheme successfully optimizes disk I/O.
-.PP
-The fractal generator makes up for its enormous consumption of CPU cycles
-by displaying robust behavior over a wide domain of tracking problems.
-The only absolute limits on the final image complexity are determined
-by the operating system's page- and disk-space requirements: refinement
-of the image grid requires no additional user intervention. Given the
-innate complexity of fractal domains, no other boundary-tracking scheme
-would seem to offer equally satisfactory results. Most important, the
-algorithm has been used -- with only trivial modifications -- to study
-both parametric and quaternionic fractal surfaces. Coupled with an
-algebraic root-finder like MACSYMA, the fractal.c programs provide a
-general-purpose package for fractal generation.
-.SH
-Raster-Scan Representation of Fractal Domains
-.PP
-The output of fractal.c is a file that contains the page-mapped model
-grid. (Note that only one of four informational bits, the BOUNDARY
-bit, is now necessary. Since the output of fractal.c is only temporary,
-there does not seem to be much of a need for storage compaction.)
-In the next stage of fractal image generation, a specialized display
-processor traverses the model grid and computes a pixel file and
-color table suitable for viewing on the JUPITER system. The display.c
-program takes as its only argument the name of the model grid; the
-user must then specify the dimension of the model grid, the desired
-view, and the output image file. The 3d fractal grid is taken to
-define a domain in left-handed Cartesian coordinates, where each
-covariant coordinate increases from zero to the grid resolution.
-There are six possible views, each of which is generated with an
-orthogonal projection. The integral values of the viewing parameter
-are defined as follows:
-.DS
-
- View: Description:
- ----- ------------
- 1 viewer at z = -infinity
- 2 viewer at x = +infinity
- 3 viewer at z = +infinity
- 4 viewer at x = -infinity
- 5 viewer at y = +infinity
- 6 viewer at y = -infinity
-.DE
-.PP
-In each case, the center of interest is the center of the fractal
-grid. Though the possible orthogonal projections are severely
-constrained, the resulting speed-up in display processing is worth
-the sacrifice.
-.PP
-The fractal display primitive is a surface cube; classification of
-a cube as a fractal boundary point implies that some component of
-the fractal domain lies within the cube's boundaries. Each cube
-maps onto one point in raster-scan space, depending on the orthogonal
-projection. A z-buffer with r^2 two-byte elements, where r is equal
-to the resolution of the model grid, records the results of the
-projection; only nearest points are kept. To simplify the processing
-and to avoid the need for multiple z-buffer projections, the light
-source is assumed to be located at the eyepoint.
-.PP
-By definition,
-fractals are infinitely discontinuous; as a consequence, surface
-normals do not exist. In order to compute a useful surface normal
-for Lambert shading, the neighbors of each pixel in the z-buffer
-are examined. Vectors connecting each pixel with neighboring
-pixel pairs are determined, and cross products are computed; the
-average z-component of the normalized cross products is taken to
-indicate the cosine of the angle between the surface normal and
-the line of projection. A simple Lambert shader incorporates
-the angle-of-incidence and depth information into an intensity
-calculation and generates a new entry into the color table if
-necessary. The shading is as smooth as is possible under the
-circumstances; insofar as fractal surfaces are never smooth,
-the shading method I have described provides an adequate
-means for computing the pixel intensity. In order to test the
-shader, I computed a simple spherical surface with my boundary
-tracer and then shaded it with display.c. The resulting Moire
-patterns are an artifact of the boundary calculation (which
-is guaranteed to generate no more than a biased sampling of
-boundary points, not a mathematically smooth surface) and in
-any case do not detract from the overall effect.
-.SH
-Acknowledgements
-.PP
-This term project is by far the most complex computing assignment
-I have ever undertaken. I couldn't have done it without the
-constant support of my friends, who brought me mountains of food
-and insisted that I sleep at least occasionally. (I'd especially
-like to thank Alyssa Eisenacher and Rod Fernandez for Herrell's
-incomparable chocolate-pudding ice cream and Babe's pizza.) My
-appreciation goes also to Steve Sistare, fellow sufferer, with
-whom I struggled through many long hours in the art-deco basement
-of Aiken Computation Lab. Finally, I'd like to thank the
-one-and-only (three-and-only?) bobrowns for their many efforts on
-my behalf throughout the course of the project. Without exclusive
-use of the hu-midgard VAX-11/750, I would never have gotten anywhere.
End.fractal.ms
echo x - graphics.h
sed s/.// >graphics.h <<\End.graphics.h
-/* general constant declarations */
-#define ERROR -1
-#define OK 0
-#define PI 3.1415926
-#define TRUE 1
-#define FALSE 0
-
-
-
-/* coordinate axis symbols */
-#define X 7
-#define Y 8
-#define Z 9
-
-/* some dimensions and limits */
-#define SCANLINES 484
-#define SCANWIDTH 640
-#define T_DIM 4
-
-
-
-/* data structures for basic graphic elements */
-
-typedef struct color
- {
- float red, green, blue;
- } COLOR;
-
-typedef struct point
- {
- float x, y, z, w;
- } VERTEX;
-
-typedef float TRANSFORM[T_DIM][T_DIM];
-
-
-struct polygon {
- int nr_vertices;
- VERTEX *vertex;
- COLOR *clr;
- };
-
-struct polylist {
- int nr_polygons;
- struct polygon *poly;
- };
-
-
-
-/* some crucial global variables */
-
-char *FILENAME; /* name of file current picture to be put in */
-
-struct /* characteristics of viewing environment */
- {
- float scr_dist;
- float scr_width;
- float scr_height;
- float scr_depth;
- TRANSFORM *view_tform;
- float ambient;
- COLOR background;
- } VIEW_COND;
-
-
-/* functions which return non-integer entities */
-
-TRANSFORM *compose(), *translate(), *rotate(), *scale();
-float *matmult();
-VERTEX *vectmult();
-struct polylist *polyclip ();
-
End.graphics.h
echo x - ifractal.c
sed s/.// >ifractal.c <<\End.ifractal.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The algorithm is by no means complete, and there
- are some imperfections in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- for CS175 (Computer Graphics) in the spring of 1984.
-
- For complete documentation, please refer to the companion file
- ~chavez/fractal/fractal.doc on hu-midgard on harvard.
-
- This program, "dfractal.c", computes the fractal domains that
- define convergent points for quadratic iterations over the
- quaternions. The program correctly recognizes the presence
- of a limit cycle when given a parameter mu that gives rise
- to such cycles and the corresponding cycle length. (Consult
- the primary documentation file for a brief summary of the
- methods for determining mu.) Interior points are those points
- that are attracted to the same cycle as the user-supplied
- interior point. This program represents a generalization
- (to cycles of length greater than two) of the quaternionic
- fractal-tracer qfractal.c.
-
- The program should be compiled with the csh command
-
- cc dfractal.c quaternion.o -lm -o dfractal
-
- and executed in the background at a "niced" priority.
- The program takes one command-line argument, the name
- of the fractal grid file. The file is written in a
- format suitable for processing with "display.c".
-
-*/
-
-#include "quaternion.h"
-#include <stdio.h>
-#include <ctype.h>
-
-#define MAX_ITERATION 100 /* don't loop forever in the convergence test */
-#define EPSILON 1.e-4 /* we're satisfied that we've reached the limit
- if we get this close to it */
-
-#define bool char /* standard boolean quantities */
-#define TRUE (1)
-#define FALSE (0)
-#define CORNER 0 /* mnemonics for neighborhood calculation */
-#define EDGE 1
-#define FACE 2
-
-#define TESTED 0 /* data bits in the data file */
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-#define NBUF 26 /* number of buffers for virtual scheme */
-
-/* have we reached the limit cycle yet? */
-#define at_limit(iterate,cycle) ((fabs( iterate->r-cycle->r ) < EPSILON) &&\
- (fabs( iterate->i-cycle->i ) < EPSILON) &&\
- (fabs( iterate->j-cycle->j ) < EPSILON) &&\
- (fabs( iterate->k-cycle->k ) < EPSILON) )
-#define round(x) ((int) ((x)+0.5))
-/* pack up a grid point for storage on the stack */
-#define pack(p) ((long) (res*(res*p->z+p->y)+p->x))
-#define mod(x,y) ((x)-((x)/(y)*(y)))
-#define EMPTY (CUBE *) 0
-
-/* the display primitive */
-typedef struct point {
- int x, y, z;
-} CUBE;
-
-/* the stack for storage of newly-discovered boundary points */
-long *stack;
-
-/* the grid virtual-memory file */
-FILE *grid, *fopen();
-/* grid resolution, its square, and the sort of neighbors we want */
-int res, res2, stringency;
-/* the size (in bytes) of the grid */
-long size;
-/* the mapping onto Cartesian space of the computation grid */
-float xlow, xmax, ylow, ymax, zlow, zmax, xinc, yinc, zinc;
-
-char *malloc();
-
-/* data areas for the grid buffering scheme */
-char *buffer_ring[NBUF];
-bool altered[NBUF];
-int block_ptr[NBUF];
-int current_buffer;
-
-/* quaternionic variables for the convergence test */
-quaternion *mu, *ifactor, *iterate, *cycle;
-/* the length of the mu-cycle and the limit point we're looking for */
-int cycle_length, limit_point;
-/* have we found the limit cycle yet? */
-bool initialized = FALSE;
-
-/* statistical information */
-int memrefs = 0, faults = 0, bpts = 0, ctests = 0,
- stack_ptr = 0, max_stack_size = 0;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- FILE *clear();
-
- CUBE *locate(), p;
- bool toggled();
-
- float x0, y0, z0, xf, yf, zf;
- int i;
-
- /* initialize the quaternionic variables */
- mu = qzero();
- iterate = qzero();
- ifactor = qzero();
- ifactor->i = 1.0;
- cycle = qzero();
-
- /* ask for the mu parameter */
- printf( "Mu?\n" );
- scanf( "%f %f", &mu->r, &mu->i );
- /* and the cycle length */
- printf( "Cycle length?\n" );
- scanf( "%d", &cycle_length );
-
- /* ask for resolution, stack size, a line segment that spans
- the boundary, and the stringency for adjacency computations */
- printf( "Resolution?\n" );
- scanf( "%d", &res );
- printf( "Maximum stack size?\n" );
- scanf( "%d", &max_stack_size );
- printf( "Enter the x, y, and z limits.\n" );
- scanf( "%f %f %f %f %f %f", &xlow, &xmax, &ylow, &ymax, &zlow, &zmax );
- printf( "Enter the coordinates of a segment that crosses the boundary.\n");
- scanf( "%f %f %f %f %f %f", &x0, &y0, &z0, &xf, &yf, &zf );
- printf( "Stringency? (0 -- corner, 1 -- edge, 2 -- face)\n" );
- scanf( "%d", &stringency );
-
- /* set up the mapping of Cartesian space onto our grid */
- xinc = (xmax-xlow)/res;
- yinc = (ymax-ylow)/res;
- zinc = (zmax-zlow)/res;
- res2 = res*res;
- size = res*res2/2 + 1;
-
- /* open the grid file and initialize the stack */
- grid = clear( argv[1], size );
- stack = (long *) calloc( max_stack_size, sizeof( long ) );
- if ( stack == NULL ) {
- printf( "fractal: couldn't create stack\n" );
- exit( 1 );
- }
-
- /* locate the intersection of the spanning line and the fractal domain;
- push that point onto the boundary stack and delve */
- assert( locate( x0, y0, z0, xf, yf, zf) );
-
- /* write out buffers that we still have lying around */
- for ( i = 0; i < NBUF; i++ )
- if (buffer_ring[i] != NULL && altered[i])
- {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
-
- /* print out some cross sections if we have enough room */
- printf( "Boundary points:\n\n" );
- if ( res < 80 )
- for ( p.z = 0; p.z < res; p.z++ )
- {
- printf( "\nz = %d\n", p.z );
- for ( p.y = res-1; p.y > 0; p.y-- )
- {
- for ( p.x = 0; p.x < res; p.x++ )
- {
- if (toggled( BOUNDARY, &p )) putchar('1');
- else putchar('0');
- }
- putchar('\n');
- }
- }
-
- /* generate a brief statistical report */
- printf( "\nStatistics:\n----------\n\n" );
- printf( "\tGrid size: %dx%dx%d = %d\n",
- res, res, res, res*res*res );
- printf( "\tMemory references:\t%d\n\tPage faults:\t\t%d\n",
- memrefs, faults );
- printf( "\tBoundary points:\t%d\n\tConvergence tests:\t%d\n",
- bpts, ctests );
- printf( "\tBitmap size (bytes):\t%d\n",
- size );
- fclose( grid );
-}
-
-/* open the grid data file and initialize to its maximum length */
-FILE *clear( fname, size )
-char *fname;
-long size;
-{
- FILE *buffer;
- char *empty;
- int i;
-
- empty = malloc( 1024 );
- buffer = fopen( fname, "w+" );
- if ( buffer == NULL ) {
- printf( "fractal: can't open %s\n", fname );
- exit( 1 );
- }
-
- for ( i = 1; i <= size/1024 + 1; i++ )
- fwrite( empty, 1024, 1, buffer );
-
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ )
- {
- buffer_ring[i] = NULL;
- altered[i] = FALSE;
- block_ptr[i] = -1;
- }
-
- cfree( empty );
- return( buffer );
-}
-
-/* map a grid point onto Cartesian space */
-cartesian( p, x, y, z )
-CUBE *p;
-float *x, *y, *z;
-{
- *x = xlow + p->x*xinc;
- *y = ylow + p->y*yinc;
- *z = zlow + p->z*zinc;
-}
-
-/* check the status of a grid point */
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
-
- /* compute the nibble, byte and block addresses;
- each grid point requires four bits for storage */
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
-
- /* figure out the item we want and compute a mask */
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- /* check to see if the block we need is already in core;
- if so, get the appropriate byte and & it with the mask */
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- /* no luck; we've got a page fault */
- ++faults;
-
- /* find another buffer */
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- /* write the old buffer out if necessary */
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- /* allocate some space if we haven't done so already */
- else buffer_ring[i] = malloc( BUFSIZ );
-
- /* read in the new block, then grab the byte we need */
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
-
-/* set informational bits for fractal grid points */
-set( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
- if ( item == BOUNDARY ) {
- ++bpts;
- if ( (bpts > 0) && (bpts-(bpts/1000)*1000 == 0) )
- printf( "\n%d boundary points processed.\n", bpts );
- }
-
- /* much of this came right out of toggle */
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- {
- buffer_ring[i][byte] |= mask; /* OR in the bit mask */
- altered[i] = TRUE;
- return;
- }
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- buffer_ring[i][byte] |= mask;
- block_ptr[i] = block;
- altered[i] = TRUE;
-
- return;
-}
-
-/* is this grid point on the inside or outside of the domain? */
-bool interior( p )
-CUBE *p;
-{
- bool hit, converges(), toggled(), out_of_range();
-
- /* if we're off the grid, assume we're out of the domain */
- if ( out_of_range( p ) ) return( FALSE );
-
- /* if we've already looked at this point, don't do it again */
- if ( toggled( TESTED, p ) ) return( toggled( STATUS, p ) );
- else {
- ++ctests;
- set( TESTED, p );
- /* we're testing now; go check for convergence */
- if ( (hit=converges( p )) ) set( STATUS, p );
- }
-
- return( hit );
-}
-
-/* assert that point p is on the boundary */
-assert( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
-
- int index;
- long nibble;
-
- /* keep looping until we've cleared the stack */
- while ( stack_ptr >= 0 )
- {
- index = 1;
- /* list all the neighbors; the interior ones are
- boundary candidates */
- while( n = neighbor( p, &index ) )
- if ( interior( n ) && !toggled( INSPECTED, n ) ) examine( n );
- else cfree( n );
-
- if ( stack_ptr > 0 ) {
- /* pop a point off the stack */
- nibble = stack[--stack_ptr];
- p->z = nibble/res2;
- p->y = nibble/res - res*p->z;
- p->x = nibble - (res*(res*p->z + p->y));
- }
- else --stack_ptr;
- }
-
- cfree( p );
-}
-
-/* examine a boundary candidate */
-examine( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
- int index;
-
- index = 1;
- /* if at least one of our neighbors is outside than we
- have to be on the boundary */
- while( n = neighbor( p, &index ) ) {
- if( !interior( n ) ) {
- /* push p onto the stack */
- stack[stack_ptr++] = pack( p );
- /* we've found that p is on the boundary, for sure */
- set( INSPECTED, p );
- set( BOUNDARY, p );
- cfree( n );
- break;
- }
- cfree( n );
- }
-
- set( INSPECTED, p );
-
- if ( n != EMPTY ) cfree( n );
- cfree( p );
-}
-
-/* enumerate the neighbors of p */
-CUBE *neighbor( p, index )
-CUBE *p;
-int *index;
-{
- bool out_of_range();
- CUBE *n;
-
- n = (CUBE *) malloc( sizeof( CUBE ) );
- n->x= p->x; n->y = p->y; n->z = p->z;
-
- if ( (stringency == FACE && *index > 6) ||
- (stringency == EDGE && *index > 18) ||
- (stringency == CORNER && *index > 26) )
- return( EMPTY );
-
- switch( *index ) {
- case 1: { --n->x; break; }
- case 2: { ++n->x; break; }
- case 3: { --n->y; break; }
- case 4: { ++n->y; break; }
- case 5: { --n->z; break; }
- case 6: { ++n->z; break; }
-
- case 7: { --n->z; --n->y; break; }
- case 8: { ++n->z; --n->y; break; }
- case 9: { --n->x; --n->y; break; }
- case 10: { ++n->x; --n->y; break; }
- case 11: { --n->z; --n->x; break; }
- case 12: { ++n->z; --n->x; break; }
- case 13: { --n->z; ++n->x; break; }
- case 14: { ++n->z; ++n->x; break; }
- case 15: { --n->z; ++n->y; break; }
- case 16: { ++n->z; ++n->y; break; }
- case 17: { --n->x; ++n->y; break; }
- case 18: { ++n->x; ++n->y; break; }
-
- case 19: { ++n->x; ++n->y; ++n->z; break; }
- case 20: { ++n->x; ++n->y; --n->z; break; }
- case 21: { ++n->x; --n->y; ++n->z; break; }
- case 22: { ++n->x; --n->y; --n->z; break; }
- case 23: { --n->x; ++n->y; ++n->z; break; }
- case 24: { --n->x; ++n->y; --n->z; break; }
- case 25: { --n->x; --n->y; ++n->z; break; }
- case 26: { --n->x; --n->y; --n->z; break; }
- }
-
- *index += 1;
- if ( out_of_range( n ) )
- {
- cfree( n );
- return( neighbor( p, index ) );
- }
-
- return( n );
-}
-
-/* check for grid points that are off the map */
-bool out_of_range( p )
-CUBE *p;
-{
- return ( (p->x < 0) || (p->x > res) ||
- (p->y < 0) || (p->y > res) ||
- (p->z < 0) || (p->z > res) );
-}
-
-/* turn Cartesian coordinates into a grid point */
-CUBE *bitmap( x, y, z )
-float x, y, z;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = round( (x-xlow)/xinc );
- new->y = round( (y-ylow)/yinc );
- new->z = round( (z-zlow)/zinc );
-
- return( new );
-}
-
-/* return the point halfway between the arguments */
-CUBE *bisect( left, right )
-CUBE *left, *right;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = floor( (left->x + right->x)/2. );
- new->y = floor( (left->y + right->y)/2. );
- new->z = floor( (left->z + right->z)/2. );
-
- return( new );
-}
-
-/* keep bisecting the starting line segment
- until we find ourselves a boundary point */
-CUBE *locate( x0, y0, z0, xf, yf, zf )
-float x0, y0, z0, xf, yf, zf;
-{
- bool interior();
- CUBE *in, *out, *middle, *bitmap(), *bisect();
-
- in = bitmap( x0, y0, z0 );
- out = bitmap( xf, yf, zf );
-
- /* we're looking for simply connected domains */
- if ( ! (interior( in ) && !interior( out ) ) )
- {
- printf( "The given segment does not span the boundary. ");
- exit( 1 );
- }
- do {
- middle = bisect( in, out );
- if (interior( middle ) )
- {
- cfree( in );
- in = middle;
- }
- else {
- cfree( out );
- out = middle;
- }
- } while ( abs( in->x - out->x ) > 1 ||
- abs( in->y - out->y ) > 1 ||
- abs( in->z - out->z ) > 1 );
- return( in );
-}
-
-/* perform the convergence test */
-bool converges( p )
-CUBE *p;
-{
- float x, y, z;
- bool convergence;
- int n;
-
- n = 1;
- convergence = FALSE;
- cartesian( p, &x, &y, &z );
-
- /* map x, y, and z onto the three-dimensional subspace of the
- quaternions spanned by i, j, and k */
- iterate->r = x;
- iterate->i = y;
- iterate->j = z;
- iterate->k = 0.0;
-
- /* keep looping until we explode, or we run out of permitted
- iterations, or we reach the right point on a stable limit cycle */
- while( qmagnitude( iterate ) < BIG &&
- n < MAX_ITERATION &&
- !( convergence = (!initialized? FALSE :
- ((mod(n, cycle_length) == limit_point) &&
- at_limit( iterate, cycle )))) )
- {
- /* apply the quadratic transformation */
- qmult( iterate, iterate, iterate );
- qadd( iterate, mu, iterate );
- qmult( iterate, ifactor, iterate );
- n++;
- }
-
- /* keep track of the limit point after the first time through */
- if (!initialized) {
- qmove( iterate, cycle );
- limit_point = mod(n, cycle_length);
- convergence = TRUE;
- initialized = TRUE;
- }
-
- return( convergence );
-}
End.ifractal.c
echo x - l400.sh
sed s/.// >l400.sh <<\End.l400.sh
-#! /bin/csh
-cd ~/fractal
-fractal l400.f <l400in >&l400.out
-display l400.f <dl
End.l400.sh
echo x - l400in
sed s/.// >l400in <<\End.l400in
-400
-500000
--2. 2. -2. 2. 1. 4.
-0. 0. 1. 100. 100. 1.
-0
End.l400in
echo x - lambda.c
sed s/.// >lambda.c <<\End.lambda.c
-#include <math.h>
-typedef struct complex_number {
- float r, i;
-} complex;
-
-complex *zero();
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- int n = 0;
- complex *seed;
-
- seed = zero();
- printf( "\nSeed? " );
- scanf( "%f %f", &seed->r, &seed->i );
- while (1) {
- printf( "\nIteration? " );
- scanf( "%d", &n );
- iterate( n, seed );
- printf( "\nResult: %.10f + %.10fi\n", seed->r, seed->i );
- }
-}
-
-complex *zero( )
-{
- complex *value;
-
- value = (complex *) malloc( sizeof( complex ) );
-
- value->r = 0.0;
- value->i = 0.0;
-
- return( value );
-}
-
-mult( a, b, c )
-complex *a, *b, *c;
-{
- float x, y;
-
- x = a->r*b->r - a->i*b->i;
- y = a->i*b->r + a->r*b->i;
-
- c->r = x;
- c->i = y;
-}
-
-scale( z, n )
-complex *z;
-int n;
-{
- z->r = n*z->r;
- z->i = n*z->i;
-}
-
-add( a, b, c )
-complex *a, *b, *c;
-{
- c->r = a->r + b->r;
- c->i = a->i + b->i;
-}
-
-subtract( a, b, c )
-complex *a, *b, *c;
-{
- c->r = a->r - b->r;
- c->i = a->i - b->i;
-}
-
-move( a, b )
-complex *a, *b;
-{
- b->r = a->r;
- b->i = a->i;
-}
-
-iterate( n, seed )
-int n;
-complex *seed;
-{
- int i;
- double pi, theta;
- complex *lambda, *one_z;
-
- one_z = zero();
- lambda = zero();
- pi = 4.*atan(1.);
- theta = -(2./5.)*pi;
- lambda->r = 1.;
-
- printf( "%f %f", lambda->r, lambda->i );
- for ( i=1; i<=n; i++ ) {
- one_z->r = 1.-seed->r;
- one_z->i = -seed->i;
- mult( seed, lambda, seed );
- mult( seed, one_z, seed );
- }
-}
End.lambda.c
echo x - locate.c
sed s/.// >locate.c <<\End.locate.c
-#include "complex.h"
-
-#define MAX_ITERATION 100
-#define bool char
-
-complex *lambda;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- int n = 0;
- complex *seed;
- float x1, y1, x2, y2, x, y, inc;
- bool converges();
-
- seed = zero();
- lambda = zero();
-
- printf( "\nLambda? ");
- scanf( "%f %f", &lambda->r, &lambda->i );
-
- printf( "\nx,y limits? " );
- scanf( "%f %f %f %f", &x1, &y1, &x2, &y2 );
- printf( "Increment? " );
- scanf( "%f", &inc );
- for ( x = x1; x <= x2; x += inc )
- for ( y = y1; y <= y2; y += inc )
- if ( converges( x, y ) ) printf( "%f + %fi\n", x, y );
-}
-
-bool converges( x, y )
-float x, y;
-{
- complex *iterate, *hold;
- int index;
-
- index = 1;
- iterate = zero();
- hold = zero();
- iterate->r = x;
- iterate->i = y;
-
- while( magnitude( iterate ) < BIG && index++ <= MAX_ITERATION )
- {
- mult( iterate, iterate, hold );
- subtract( iterate, hold, iterate );
- mult( iterate, lambda, iterate );
- }
-
- cfree( hold );
- cfree( iterate );
-
- return( index > MAX_ITERATION );
-}
End.locate.c
echo x - q2.sh
sed s/.// >q2.sh <<\End.q2.sh
-#! /bin/csh
-cd ~/fractal
-ifractal q2.f <q400in >&q400.out
-display q2.f <dq
End.q2.sh
echo x - q400in
sed s/.// >q400in <<\End.q400in
-1 0
-2
-50
-500000
--2. 2. -2. 2. -2. 2.
-0. 1. 0. 100. 100. 0.
-0
End.q400in
echo x - qdragon.c
sed s/.// >qdragon.c <<\End.qdragon.c
-#include "quaternion.h"
-
-quaternion *mu;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- int n = 0;
- quaternion *seed;
-
- seed = qzero();
- mu = qzero();
- mu->i = 1.;
-
- printf( "\nSeed? " );
- scanf( "%f %f %f", &seed->r, &seed->i, &seed->j );
- while (1) {
- printf( "\nIteration? " );
- scanf( "%d", &n );
- iterate( n, seed );
- printf( "\nResult: %.10f + %.10fi + %.10fj + %.10fk\n",
- seed->r, seed->i, seed->j, seed->k );
- }
-}
-
-iterate( n, seed )
-int n;
-quaternion *seed;
-{
- int i;
- double pi, theta;
-
- for ( i=1; i<=n; i++ ) {
- qmult( seed, seed, seed );
- seed->r += 1.;
- qmult( seed, mu, seed );
- }
-}
End.qdragon.c
echo x - qfractal.c
sed s/.// >qfractal.c <<\End.qfractal.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The information is by no means complete, and there
- are some imperfections in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- for CS175 (Computer Graphics) in the spring of 1984.
-
- The program "qfractal.f" traces one component (of which the origin
- is a member) of the domain of points attracted to a cycle of length
- two. This program is a specific case of the more general convergence
- algorithm in "dfractal.c", which correctly calculates any component
- of a limit cycle with arbitrary length. The operation of "qfractal.c"
- is exactly equivalent to "dfractal.c" with mu = 1.0 (purely real)
- and a cycle length of 2. The special-purpose "qfractal.c" exists
- only because I wrote it before I fully understood the intricacies
- of mu.
-
- Compile the program with the csh command
-
- cc qfractal.c quaternion.o -lm -o qfractal
-
- There is a single command-line argument, which is the
- name of the fractal grid file. The file is written in a
- format suitable for processing with "display.c".
-
-*/
-#include "quaternion.h"
-#include <stdio.h>
-#include <ctype.h>
-
-#define MAX_ITERATION 100
-#define EPSILON 1.e-4
-
-#define bool char
-#define TRUE (1)
-#define FALSE (0)
-#define CORNER 0
-#define EDGE 1
-#define FACE 2
-
-#define TESTED 0
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-#define NBUF 26
-
-#define round(x) ((int) ((x)+0.5))
-#define pack(p) ((long) (res*(res*p->z+p->y)+p->x))
-#define mod(x,y) ((x)-((x)/(y)*(y)))
-#define EMPTY (CUBE *) 0
-
-typedef struct point {
- int x, y, z;
-} CUBE;
-
-long *stack;
-
-FILE *grid, *fopen();
-int res, res2, stringency;
-long size;
-float xlow, xmax, ylow, ymax, zlow, zmax, xinc, yinc, zinc;
-
-char *malloc();
-
-char *buffer_ring[NBUF];
-bool altered[NBUF];
-int block_ptr[NBUF];
-int current_buffer;
-
-quaternion *mu, *iterate;
-
-int memrefs = 0, faults = 0, bpts = 0, ctests = 0,
- stack_ptr = 0, max_stack_size = 0;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- FILE *clear();
-
- CUBE *locate(), p;
- bool toggled();
-
- float x0, y0, z0, xf, yf, zf;
- int i;
-
- printf( "Resolution?\n" );
- scanf( "%d", &res );
- printf( "Maximum stack size?\n" );
- scanf( "%d", &max_stack_size );
- printf( "Enter the x, y, and z limits.\n" );
- scanf( "%f %f %f %f %f %f", &xlow, &xmax, &ylow, &ymax, &zlow, &zmax );
- printf( "Enter the coordinates of a segment that crosses the boundary.\n");
- scanf( "%f %f %f %f %f %f", &x0, &y0, &z0, &xf, &yf, &zf );
- printf( "Stringency? (0 -- corner, 1 -- edge, 2 -- face)\n" );
- scanf( "%d", &stringency );
-
- xinc = (xmax-xlow)/res;
- yinc = (ymax-ylow)/res;
- zinc = (zmax-zlow)/res;
- res2 = res*res;
- size = res*res2/2 + 1;
- grid = clear( argv[1], size );
- stack = (long *) calloc( max_stack_size, sizeof( long ) );
- if ( stack == NULL ) {
- printf( "fractal: couldn't create stack\n" );
- exit( 1 );
- }
-
- mu = qzero();
- mu->i = 1.;
- iterate = qzero();
-
- assert( locate( x0, y0, z0, xf, yf, zf) );
- printf( "Boundary points:\n\n" );
-
- for ( i = 0; i < NBUF; i++ )
- if (buffer_ring[i] != NULL && altered[i])
- {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
-
- if ( res < 80 )
- for ( p.z = 0; p.z < res; p.z++ )
- {
- printf( "\nz = %d\n", p.z );
- for ( p.y = res-1; p.y > 0; p.y-- )
- {
- for ( p.x = 0; p.x < res; p.x++ )
- {
- if (toggled( BOUNDARY, &p )) putchar('1');
- else putchar('0');
- }
- putchar('\n');
- }
- }
-
- printf( "\nStatistics:\n----------\n\n" );
- printf( "\tGrid size: %dx%dx%d = %d\n",
- res, res, res, res*res*res );
- printf( "\tMemory references:\t%d\n\tPage faults:\t\t%d\n",
- memrefs, faults );
- printf( "\tBoundary points:\t%d\n\tConvergence tests:\t%d\n",
- bpts, ctests );
- printf( "\tBitmap size (bytes):\t%d\n",
- size );
- fclose( grid );
-}
-
-FILE *clear( fname, size )
-char *fname;
-long size;
-{
- FILE *buffer;
- char *empty;
- int i;
-
- empty = malloc( 1024 );
- buffer = fopen( fname, "w+" );
- if ( buffer == NULL ) {
- printf( "fractal: can't open %s\n", fname );
- exit( 1 );
- }
-
- for ( i = 1; i <= size/1024 + 1; i++ )
- fwrite( empty, 1024, 1, buffer );
-
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ )
- {
- buffer_ring[i] = NULL;
- altered[i] = FALSE;
- block_ptr[i] = -1;
- }
-
- cfree( empty );
- return( buffer );
-}
-
-cartesian( p, x, y, z )
-CUBE *p;
-float *x, *y, *z;
-{
- *x = xlow + p->x*xinc;
- *y = ylow + p->y*yinc;
- *z = zlow + p->z*zinc;
-}
-
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
-
-set( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
- if ( item == BOUNDARY ) {
- ++bpts;
- if ( (bpts > 0) && (bpts-(bpts/1000)*1000 == 0) )
- printf( "\n%d boundary points processed.\n", bpts );
- }
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- {
- buffer_ring[i][byte] |= mask;
- altered[i] = TRUE;
- return;
- }
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- buffer_ring[i][byte] |= mask;
- block_ptr[i] = block;
- altered[i] = TRUE;
-
- return;
-}
-
-bool interior( p )
-CUBE *p;
-{
- bool hit, converges(), toggled(), out_of_range();
-
- if ( out_of_range( p ) ) return( FALSE );
- if ( toggled( TESTED, p ) ) return( toggled( STATUS, p ) );
- else {
- ++ctests;
- set( TESTED, p );
- if ( (hit=converges( p )) ) set( STATUS, p );
- }
-
- return( hit );
-}
-
-assert( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
-
- int index;
- long nibble;
-
- while ( stack_ptr >= 0 )
- {
- index = 1;
- while( n = neighbor( p, &index ) )
- if ( interior( n ) && !toggled( INSPECTED, n ) ) examine( n );
- else cfree( n );
-
- if ( stack_ptr > 0 ) {
- nibble = stack[--stack_ptr];
- p->z = nibble/res2;
- p->y = nibble/res - res*p->z;
- p->x = nibble - (res*(res*p->z + p->y));
- }
- else --stack_ptr;
- }
-
- cfree( p );
-}
-
-examine( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
- int index;
-
- index = 1;
- while( n = neighbor( p, &index ) ) {
- if( !interior( n ) ) {
- stack[stack_ptr++] = pack( p );
- set( INSPECTED, p );
- set( BOUNDARY, p );
- cfree( n );
- break;
- }
- cfree( n );
- }
-
- set( INSPECTED, p );
-
- if ( n != EMPTY ) cfree( n );
- cfree( p );
-}
-
-CUBE *neighbor( p, index )
-CUBE *p;
-int *index;
-{
- bool out_of_range();
- CUBE *n;
-
- n = (CUBE *) malloc( sizeof( CUBE ) );
- n->x= p->x; n->y = p->y; n->z = p->z;
-
- if ( (stringency == FACE && *index > 6) ||
- (stringency == EDGE && *index > 18) ||
- (stringency == CORNER && *index > 26) )
- return( EMPTY );
-
- switch( *index ) {
- case 1: { --n->x; break; }
- case 2: { ++n->x; break; }
- case 3: { --n->y; break; }
- case 4: { ++n->y; break; }
- case 5: { --n->z; break; }
- case 6: { ++n->z; break; }
-
- case 7: { --n->z; --n->y; break; }
- case 8: { ++n->z; --n->y; break; }
- case 9: { --n->x; --n->y; break; }
- case 10: { ++n->x; --n->y; break; }
- case 11: { --n->z; --n->x; break; }
- case 12: { ++n->z; --n->x; break; }
- case 13: { --n->z; ++n->x; break; }
- case 14: { ++n->z; ++n->x; break; }
- case 15: { --n->z; ++n->y; break; }
- case 16: { ++n->z; ++n->y; break; }
- case 17: { --n->x; ++n->y; break; }
- case 18: { ++n->x; ++n->y; break; }
-
- case 19: { ++n->x; ++n->y; ++n->z; break; }
- case 20: { ++n->x; ++n->y; --n->z; break; }
- case 21: { ++n->x; --n->y; ++n->z; break; }
- case 22: { ++n->x; --n->y; --n->z; break; }
- case 23: { --n->x; ++n->y; ++n->z; break; }
- case 24: { --n->x; ++n->y; --n->z; break; }
- case 25: { --n->x; --n->y; ++n->z; break; }
- case 26: { --n->x; --n->y; --n->z; break; }
- }
-
- *index += 1;
- if ( out_of_range( n ) )
- {
- cfree( n );
- return( neighbor( p, index ) );
- }
-
- return( n );
-}
-
-bool out_of_range( p )
-CUBE *p;
-{
- return ( (p->x < 0) || (p->x > res) ||
- (p->y < 0) || (p->y > res) ||
- (p->z < 0) || (p->z > res) );
-}
-
-bool converges( p )
-CUBE *p;
-{
- float x, y, z, mag;
- int index;
-
- index = 1;
- cartesian( p, &x, &y, &z );
- iterate->r = x;
- iterate->i = y;
- iterate->j = z;
- iterate->k = 0.0;
-
- while( (index++ < MAX_ITERATION) &&
- ((mag = qmagnitude( iterate )) < BIG) &&
- ( (fabs( iterate->i ) > EPSILON) ||
- (fabs( iterate->r ) > EPSILON) ||
- (fabs( iterate->j ) > EPSILON) ||
- (fabs( iterate->k ) > EPSILON) ) )
- {
- qmult( iterate, iterate, iterate );
- iterate->r += 1.;
- qmult( mu, iterate, iterate );
- }
-
- return( (mod( index, 2 ) == 0) && (mag < BIG) );
-}
-
-CUBE *bitmap( x, y, z )
-float x, y, z;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = round( (x-xlow)/xinc );
- new->y = round( (y-ylow)/yinc );
- new->z = round( (z-zlow)/zinc );
-
- return( new );
-}
-
-CUBE *bisect( left, right )
-CUBE *left, *right;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = floor( (left->x + right->x)/2. );
- new->y = floor( (left->y + right->y)/2. );
- new->z = floor( (left->z + right->z)/2. );
-
- return( new );
-}
-
-CUBE *locate( x0, y0, z0, xf, yf, zf )
-float x0, y0, z0, xf, yf, zf;
-{
- bool interior();
- CUBE *in, *out, *middle, *bitmap(), *bisect();
-
- in = bitmap( x0, y0, z0 );
- out = bitmap( xf, yf, zf );
-
- if ( ! (interior( in ) && !interior( out ) ) )
- {
- printf( "The given segment does not span the boundary. ");
- exit( 1 );
- }
- do {
- middle = bisect( in, out );
- if (interior( middle ) )
- {
- cfree( in );
- in = middle;
- }
- else {
- cfree( out );
- out = middle;
- }
- } while ( abs( in->x - out->x ) > 1 ||
- abs( in->y - out->y ) > 1 ||
- abs( in->z - out->z ) > 1 );
- return( in );
-}
End.qfractal.c
echo x - qin
sed s/.// >qin <<\End.qin
-50 50 50
--2. 2. -2. 2. -2. 2.
-0. 0. 0. 100. 100. 0.
-2
End.qin
echo x - quaternion.c
sed s/.// >quaternion.c <<\End.quaternion.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The algorithm is by no means complete, and there
- are many known bugs in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- in CS175 (Computer Graphics) in the spring of 1984.
-
- For complete documentation, please refer to the companion file
- ~chavez/fractal/fractal.doc on hu-midgard on harvard.
-
- The routines in "quaternion.c" provide a general scheme for
- manipulation of the algebraic quantities known as quaternions.
- To use the package, be sure to "#include quaternion.h"; in
- the C compilation, reference the object file quaternion.o.
-
-*/
-
-#include <math.h>
-
-#define BIG 1.e18
-
-typedef struct quaternionic_number {
- float r, i, j, k;
-} quaternion;
-
-float qmagnitude( z )
-quaternion *z;
-{
- if ( fabs( z->r ) < BIG && fabs( z->i ) < BIG &&
- fabs( z->j ) < BIG && fabs( z->k ) < BIG )
- return( sqrt( z->r*z->r + z->i*z->i +
- z->j*z->j + z->k*z->k ) );
- else return( BIG );
-}
-
-quaternion *qzero( )
-{
- quaternion *value;
-
- value = (quaternion *) malloc( sizeof( quaternion ) );
-
- value->r = 0.0;
- value->i = 0.0;
- value->j = 0.0;
- value->k = 0.0;
-
- return( value );
-}
-
-qmult( a, b, c )
-quaternion *a, *b, *c;
-{
- float r, i, j, k;
-
-
- r = a->r*b->r - a->i*b->i - a->j*b->j - a->k*b->k;
- i = a->r*b->i + a->i*b->r + a->j*b->k - a->k*b->j;
- j = a->r*b->j - a->i*b->k + a->j*b->r + a->k*b->i;
- k = a->r*b->k + a->i*b->j - a->j*b->i + a->k*b->r;
-
- c->r = r;
- c->i = i;
- c->j = j;
- c->k = k;
-}
-
-qscale( z, n )
-quaternion *z;
-float n;
-{
- z->r = n*z->r;
- z->i = n*z->i;
- z->j = n*z->j;
- z->k = n*z->k;
-}
-
-qadd( a, b, c )
-quaternion *a, *b, *c;
-{
- c->r = a->r + b->r;
- c->i = a->i + b->i;
- c->j = a->j + b->j;
- c->k = a->k + b->k;
-}
-
-qsubtract( a, b, c )
-quaternion *a, *b, *c;
-{
- c->r = a->r - b->r;
- c->i = a->i - b->i;
- c->j = a->j - b->j;
- c->k = a->k - b->k;
-}
-
-qmove( a, b )
-quaternion *a, *b;
-{
- b->r = a->r;
- b->i = a->i;
- b->j = a->j;
- b->k = a->k;
-}
End.quaternion.c
echo x - quaternion.h
sed s/.// >quaternion.h <<\End.quaternion.h
-#include <math.h>
-
-#define BIG 1.0e18
-
-typedef struct quaternionic_number {
- float r, i, j, k;
-} quaternion;
-
-float qmagnitude();
-quaternion *qzero();
End.quaternion.h
echo x - quaternion.sh
sed s/.// >quaternion.sh <<\End.quaternion.sh
-#! /bin/csh
-cd ~/fractal
-qfractal qfractal.f <q400in >&q400.out
-display qfractal.f <dq
End.quaternion.sh
echo x - reflect.c
sed s/.// >reflect.c <<\End.reflect.c
-#include <stdio.h>
-#include "/toe/graphics/lib/source/ngs/graphics.h"
-
-char pixel[SCANLINES*SCANWIDTH];
-
-main(argc, argv)
-char **argc;
-int argv;
-{
- FILE *p, *q;
- int i;
-
- p = fopen("qfractal.img", "r");
- q = fopen("qr.img", "w");
-
- fread( pixel, sizeof( char ), SCANLINES*SCANWIDTH*2, p );
- for ( i = SCANLINES*SCANWIDTH/2; i < SCANLINES*SCANWIDTH; i++ )
- pixel[i]=pixel[SCANLINES*SCANWIDTH-i];
-
- fwrite( pixel, sizeof( char ), SCANLINES*SCANWIDTH, q );
-}
End.reflect.c
echo x - wind.c
sed s/.// >wind.c <<\End.wind.c
-/*
- Computer Science 175: Computer Graphics
- Three-Dimensional Fractal Computation
-
- Programmer: R. Martin Chavez '84
- Implementation date: 5/8/84
-
-
- Target systems:
- harvard.ARPA and hu-midgard.ARPA
-
- The package of routines in ~chavez/fractal provide a fairly comprehensive
- set of algorithms for the generation and display of three-dimensional
- fractal figures. The information is by no means complete, and there
- are some imperfections in the raster display driver. Even so, these routines
- offer a reasonable first implementation, given time and resource
- limitations. The work represented herein was submitted as a term project
- in CS175 (Computer Graphics) in the spring of 1984.
-
- This file, "wind.c", contains a convergence checker that generates
- a parametric family of fractal curves where the phase of lambda
- takes on values from zlow to zmax (in radians) and the magnitude
- remains constant.
-
-*/
-
-#include "complex.h"
-#include <stdio.h>
-#include <ctype.h>
-
-#define MAX_ITERATION 100
-#define EPSILON 1.e-4
-
-#define bool char
-#define TRUE (1)
-#define FALSE (0)
-#define CORNER 0
-#define EDGE 1
-#define FACE 2
-
-#define TESTED 0
-#define STATUS 1
-#define INSPECTED 2
-#define BOUNDARY 3
-#define NBUF 26
-
-#define round(x) ((int) ((x)+0.5))
-#define pack(p) ((long) (res*(res*p->z+p->y)+p->x))
-#define mod(x,y) ((x)-((x)/(y)*(y)))
-#define EMPTY (CUBE *) 0
-
-typedef struct point {
- int x, y, z;
-} CUBE;
-
-long *stack;
-
-FILE *grid, *fopen();
-int res, res2, stringency;
-long size;
-float xlow, xmax, ylow, ymax, zlow, zmax, xinc, yinc, zinc;
-double *cs, *sn;
-
-char *malloc();
-
-char *buffer_ring[NBUF];
-bool altered[NBUF];
-int block_ptr[NBUF];
-int current_buffer;
-
-complex *hold, *iterate;
-
-int memrefs = 0, faults = 0, bpts = 0, ctests = 0,
- stack_ptr = 0, max_stack_size = 0;
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- FILE *clear();
-
- CUBE *locate(), p;
- bool toggled();
-
- float x0, y0, z0, xf, yf, zf;
- int i;
-
- printf( "Resolution?\n" );
- scanf( "%d", &res );
- printf( "Maximum stack size?\n" );
- scanf( "%d", &max_stack_size );
- printf( "Enter the x, y, and z limits.\n" );
- scanf( "%f %f %f %f %f %f", &xlow, &xmax, &ylow, &ymax, &zlow, &zmax );
- printf( "Enter the coordinates of a segment that crosses the boundary.\n");
- scanf( "%f %f %f %f %f %f", &x0, &y0, &z0, &xf, &yf, &zf );
- printf( "Stringency? (0 -- corner, 1 -- edge, 2 -- face)\n" );
- scanf( "%d", &stringency );
-
- xinc = (xmax-xlow)/res;
- yinc = (ymax-ylow)/res;
- zinc = (zmax-zlow)/res;
- res2 = res*res;
- size = res*res2/2 + 1;
- grid = clear( argv[1], size );
- stack = (long *) calloc( max_stack_size, sizeof( long ) );
- if ( stack == NULL ) {
- printf( "fractal: couldn't create stack\n" );
- exit( 1 );
- }
-
- hold = zero();
- iterate = zero();
- cs = (double *) calloc( res, sizeof( double ) );
- sn = (double *) calloc( res, sizeof( double ) );
- for ( i=0; i<res; i++ ) {
- cs[i] = cos(i*zinc + zlow);
- sn[i] = sin(i*zinc + zlow);
- }
-
- assert( locate( x0, y0, z0, xf, yf, zf) );
- printf( "Boundary points:\n\n" );
-
- for ( i = 0; i < NBUF; i++ )
- if (buffer_ring[i] != NULL && altered[i])
- {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
-
- if ( res < 80 )
- for ( p.z = 0; p.z < res; p.z++ )
- {
- printf( "\nz = %d\n", p.z );
- for ( p.y = res-1; p.y > 0; p.y-- )
- {
- for ( p.x = 0; p.x < res; p.x++ )
- {
- if (toggled( BOUNDARY, &p )) putchar('1');
- else putchar('0');
- }
- putchar('\n');
- }
- }
-
- printf( "\nStatistics:\n----------\n\n" );
- printf( "\tGrid size: %dx%dx%d = %d\n",
- res, res, res, res*res*res );
- printf( "\tMemory references:\t%d\n\tPage faults:\t\t%d\n",
- memrefs, faults );
- printf( "\tBoundary points:\t%d\n\tConvergence tests:\t%d\n",
- bpts, ctests );
- printf( "\tBitmap size (bytes):\t%d\n",
- size );
- fclose( grid );
-}
-
-FILE *clear( fname, size )
-char *fname;
-long size;
-{
- FILE *buffer;
- char *empty;
- int i;
-
- empty = malloc( 1024 );
- buffer = fopen( fname, "w+" );
- if ( buffer == NULL ) {
- printf( "fractal: can't open %s\n", fname );
- exit( 1 );
- }
-
- for ( i = 1; i <= size/1024 + 1; i++ )
- fwrite( empty, 1024, 1, buffer );
-
- current_buffer = 0;
- for ( i = 0; i < NBUF; i++ )
- {
- buffer_ring[i] = NULL;
- altered[i] = FALSE;
- block_ptr[i] = -1;
- }
-
- cfree( empty );
- return( buffer );
-}
-
-cartesian( p, x, y, z )
-CUBE *p;
-float *x, *y, *z;
-{
- *x = xlow + p->x*xinc;
- *y = ylow + p->y*yinc;
- *z = zlow + p->z*zinc;
-}
-
-bool toggled( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- return( buffer_ring[i][byte] & mask );
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- block_ptr[i] = block;
- altered[i] = FALSE;
-
- return( buffer_ring[i][byte] & mask );
-}
-
-set( item, p )
-int item;
-CUBE *p;
-{
- long block, nibble, byte;
- int i;
- int mask;
-
- ++memrefs;
- if ( item == BOUNDARY ) {
- ++bpts;
- if ( (bpts > 0) && (bpts-(bpts/1000)*1000 == 0) )
- printf( "\n%d boundary points processed.\n", bpts );
- }
-
- nibble = res*(res*p->z + p->y) + p->x;
- block = nibble/(2*BUFSIZ);
- byte = nibble/2 - (block*BUFSIZ);
- if ( nibble & 1 ) mask = (1 << (item+4));
- else mask = (1 << item);
-
- for ( i = 0; i < NBUF; i++ )
- if ( block_ptr[i] == block )
- {
- buffer_ring[i][byte] |= mask;
- altered[i] = TRUE;
- return;
- }
-
- ++faults;
- if (current_buffer < NBUF-1) ++current_buffer;
- else current_buffer = 0;
-
- i = current_buffer;
- if (buffer_ring[i] != NULL)
- {
- if ( altered[i] ) {
- fseek( grid, block_ptr[i]*BUFSIZ, 0 );
- fwrite( buffer_ring[i], BUFSIZ, 1, grid );
- }
- }
- else buffer_ring[i] = malloc( BUFSIZ );
-
- fseek( grid, block*BUFSIZ, 0 );
- fread( buffer_ring[i], BUFSIZ, 1, grid );
- buffer_ring[i][byte] |= mask;
- block_ptr[i] = block;
- altered[i] = TRUE;
-
- return;
-}
-
-bool interior( p )
-CUBE *p;
-{
- bool hit, converges(), toggled(), out_of_range();
-
- if ( out_of_range( p ) ) return( FALSE );
- if ( toggled( TESTED, p ) ) return( toggled( STATUS, p ) );
- else {
- ++ctests;
- set( TESTED, p );
- if ( (hit=converges( p )) ) set( STATUS, p );
- }
-
- return( hit );
-}
-
-assert( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
-
- int index;
- long nibble;
-
- while ( stack_ptr >= 0 )
- {
- index = 1;
- while( n = neighbor( p, &index ) )
- if ( interior( n ) && !toggled( INSPECTED, n ) ) examine( n );
- else cfree( n );
-
- if ( stack_ptr > 0 ) {
- nibble = stack[--stack_ptr];
- p->z = nibble/res2;
- p->y = nibble/res - res*p->z;
- p->x = nibble - (res*(res*p->z + p->y));
- }
- else --stack_ptr;
- }
-
- cfree( p );
-}
-
-examine( p )
-CUBE *p;
-{
- CUBE *n, *neighbor();
- bool interior();
- int index;
-
- index = 1;
- while( n = neighbor( p, &index ) ) {
- if( !interior( n ) ) {
- stack[stack_ptr++] = pack( p );
- set( INSPECTED, p );
- set( BOUNDARY, p );
- cfree( n );
- break;
- }
- cfree( n );
- }
-
- set( INSPECTED, p );
-
- if ( n != EMPTY ) cfree( n );
- cfree( p );
-}
-
-CUBE *neighbor( p, index )
-CUBE *p;
-int *index;
-{
- bool out_of_range();
- CUBE *n;
-
- n = (CUBE *) malloc( sizeof( CUBE ) );
- n->x= p->x; n->y = p->y; n->z = p->z;
-
- if ( (stringency == FACE && *index > 6) ||
- (stringency == EDGE && *index > 18) ||
- (stringency == CORNER && *index > 26) )
- return( EMPTY );
-
- switch( *index ) {
- case 1: { --n->x; break; }
- case 2: { ++n->x; break; }
- case 3: { --n->y; break; }
- case 4: { ++n->y; break; }
- case 5: { --n->z; break; }
- case 6: { ++n->z; break; }
-
- case 7: { --n->z; --n->y; break; }
- case 8: { ++n->z; --n->y; break; }
- case 9: { --n->x; --n->y; break; }
- case 10: { ++n->x; --n->y; break; }
- case 11: { --n->z; --n->x; break; }
- case 12: { ++n->z; --n->x; break; }
- case 13: { --n->z; ++n->x; break; }
- case 14: { ++n->z; ++n->x; break; }
- case 15: { --n->z; ++n->y; break; }
- case 16: { ++n->z; ++n->y; break; }
- case 17: { --n->x; ++n->y; break; }
- case 18: { ++n->x; ++n->y; break; }
-
- case 19: { ++n->x; ++n->y; ++n->z; break; }
- case 20: { ++n->x; ++n->y; --n->z; break; }
- case 21: { ++n->x; --n->y; ++n->z; break; }
- case 22: { ++n->x; --n->y; --n->z; break; }
- case 23: { --n->x; ++n->y; ++n->z; break; }
- case 24: { --n->x; ++n->y; --n->z; break; }
- case 25: { --n->x; --n->y; ++n->z; break; }
- case 26: { --n->x; --n->y; --n->z; break; }
- }
-
- *index += 1;
- if ( out_of_range( n ) )
- {
- cfree( n );
- return( neighbor( p, index ) );
- }
-
- return( n );
-}
-
-bool out_of_range( p )
-CUBE *p;
-{
- return ( (p->x < 0) || (p->x > res) ||
- (p->y < 0) || (p->y > res) ||
- (p->z < 0) || (p->z > res) );
-}
-
-CUBE *bitmap( x, y, z )
-float x, y, z;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = round( (x-xlow)/xinc );
- new->y = round( (y-ylow)/yinc );
- new->z = round( (z-zlow)/zinc );
-
- return( new );
-}
-
-CUBE *bisect( left, right )
-CUBE *left, *right;
-{
- CUBE *new;
-
- new = (CUBE *) malloc( sizeof( CUBE ) );
- new->x = floor( (left->x + right->x)/2. );
- new->y = floor( (left->y + right->y)/2. );
- new->z = floor( (left->z + right->z)/2. );
-
- return( new );
-}
-
-CUBE *locate( x0, y0, z0, xf, yf, zf )
-float x0, y0, z0, xf, yf, zf;
-{
- bool interior();
- CUBE *in, *out, *middle, *bitmap(), *bisect();
-
- in = bitmap( x0, y0, z0 );
- out = bitmap( xf, yf, zf );
-
- if ( ! (interior( in ) && !interior( out ) ) )
- {
- printf( "The given segment does not span the boundary. ");
- exit( 1 );
- }
- do {
- middle = bisect( in, out );
- if (interior( middle ) )
- {
- cfree( in );
- in = middle;
- }
- else {
- cfree( out );
- out = middle;
- }
- } while ( abs( in->x - out->x ) > 1 ||
- abs( in->y - out->y ) > 1 ||
- abs( in->z - out->z ) > 1 );
- return( in );
-}
-
-bool converges( p )
-CUBE *p;
-{
- complex *lambda, *iterate, *hold;
- float x, y, z;
- int index;
-
- index = 1;
- iterate = zero();
- hold = zero();
- lambda = zero();
- cartesian( p, &x, &y, &z );
-
- lambda->r = cs[p->z];
- lambda->i = sn[p->z];
- iterate->r = x;
- iterate->i = y;
-
- while( magnitude( iterate ) < BIG && index++ <= MAX_ITERATION )
- {
- mult( iterate, iterate, hold );
- subtract( iterate, hold, iterate );
- mult( iterate, lambda, iterate );
- }
-
- cfree( hold );
- cfree( iterate );
- cfree( lambda );
-
- return( index > MAX_ITERATION );
-}
End.wind.c