Willamette Web Weavers
Yet Another Fractal -- Source Code
[*THERE IS A FRACTAL GRAPHIC HERE*]
The Algorithm
The algorithm for this fractal is similar to the classic set popularized by
Benoit Mandelbrot (discovered by others, earlier). But there are some
crucial differences:
The regular Mandelbrot set algorithm is:
Z = Z * Z + c
where 'c' is a complex "constant" which is used a little more like a
variable.
That algorithm is repeatedly run on each point in the image. Before starting
the iterations, 'Z' is zeroed, and 'c' is initialized with some scaled
version of the x and y pixel coordinates.
The iterations continue until either the absolute value of the complex
answer 'Z' goes above some arbitrary limit, or the routine stops because of
some set limit to the number of iterations. (You can't loop forever and get
a picture done.) The point is then assigned a color based on how many times
it was exponentiated before it went above the limit.
Some points never go above the limit. They comprise the real Mandelbrot set,
and have that classic cartoid shape, with all of the littler cartoids
hanging off of it. They are traditionally colored black or white, while the
fringe points that are almost in the set get all of the pretty colors.
* The algorithm for this fractal is:
Z = (Z + c) ** 2
Note that we add the constant 'c' before exponentiating. Another big
difference is that 'c' really is a constant: it never changes through the
whole process of generating an image. Also, we initialize 'Z' with the pixel
coordinate information. The other major difference is that we graph the
absolute value of the point when it went out of bounds, not paying any
attention to how many iterations were done before it went out.
This amounts to a crude way of graphing the derivative of the point (the
slope, or rate of increase) as the point is exponentiated. Some points will
slowly increase in value, and will just crawl over the line as they go out
of bounds, and will have a low resulting value. Some other points will be
skyrocketing towards infinity as they are repeatedly squared, and they will
jump way beyond the limit on their last iteration, and will have a very high
value. They will get colored black (in these images), and will be considered
to be outside of the fractal set.
And again, some points never go out. They are colored white in these images.
As an added twist, you can either assign colors based on the absolute value
of the point, or on the absolute value minus the upper limit, to get a
little more contrast.
This fractal is much larger than the Mandelbrot fractal. Mandelbrot's
favorite is approximately +-2 vertically (imaginary) and +-3 horizontally
(real). This thing appears to be infinitely large. The more you zoom out,
the more outer concentric circles you find, until you can't see them any
more because they are so thin, when viewed from that far away, that they
fall into the cracks between the pixels on the screen, and break up into
vague dot patterns. It will take some time to explore infinity, and see if
these things really do go on forever in both directions (microscopically and
macroscopically).
David Joyce, Professor of Math. & Comp. Sci. at Clark University, Worcester,
MA,
djoyce@black.clarku.edu
http://aleph0.clarku.edu/~djoyce/home.html
commented:
Terry, I'm not sure, but it looks like a Julia set to me. When finding the
Mandelbrot set, z starts out at 0, c is complex, and you iterate. When
finding a Julia set, c is fixed and z is complex, and you iterate. Actually,
when Julia first studied these things in the 20's, he considered any
polynomial and iterated it, not just quadratic polynomials. I'm not sure,
but I seem to remember that he restricted his attention to real numbers
instead of complex numbers.
-- Dave
By the way, Dave has the best fractal web pages I've ever seen. When you are
done with this page, check out The Fractal Explorer on the previous page.
The source code for generating these things follows.
[Image]
How to make an image:
These images are all GIF files. I used "RAW2GIF.EXE" to make GIF files out
of the raw binary files that came out of the following program. Raw2gif is
easily fetched from oak.oakland.edu using anonymous FTP. It is in
/SimTel/msdos/gif, and is hidden in a screen-grabbing package called
"GRABSC11.ZIP".
Raw2gif wants a pallette file for input, as well as raw image data. A
pallette file is simply a list of numbers for color number (1 to 256), and
Red, Green, and Blue values (0 to 255) for the colors. A pallette file is
just plain ASCII text, and looks like this:
1 255 255 255
2 137 0 0
3 136 0 0
4 134 0 0
5 132 0 0
...
252 0 0 127
253 0 0 130
254 0 0 134
255 0 0 137
256 0 0 0
For best results, adhere to the rigid format of 3 columns per number,
padding with leading spaces where needed, and putting two spaces between the
numbers' columns. By the way, movecols (part of the Columns Toolkit; see The
Willamette Web Weavers home page) is handy for swapping columns of numbers,
so you can easily exchange colors in your pallette as you try out different
color schemes.
You can make up you own pallette, or use a default one. There is no such
thing as the "right" pallette to use. The assignment of colors to numbers is
strictly arbitrary, and is a matter of one's own personal artistic tastes.
If you save the ".RAW" data file, you can try it out with several different
pallettes, to find something that is pleasing to the eye. If you use the VGA
graphics screen-grabber TSR, and snapshot some screen you like, you can
steal the pallette from that, too, because the screen-grabber TSR writes out
both a binary image data file and the (ascii) pallette file for it. It just
sucks the color table right out of the VGA board.
The pallette I used in most all of these images is just a modified rainbow,
that goes "white, pink, red, orange, yellow, green, blue, indigo, violet,
ultra-violet (black)" as you go from zero to 255.
You make a GIF image file with a command like:
raw2gif -s 640 480 -p pall006.pal hfract.raw >hfracte3.gif
where 640x480 is the size of the image, pall006.pal is the name of the
pallette file, hfract.raw is the raw binary image data, and hfracte3.gif is
the name of the desired output file.
Once you have converted the data to a GIF image, you can use any GIF viewer
to view the image. I like "Lview" a lot. It's also easily fetched from out
there on the Net.
Web Page About LView Pro
[Image]
The Source Code
If you tell your Web browser to save this page, "Save As" then you will have
the following source code in a text file, and can chop out the separate
files. They are small.
Or, you can download a ZIP file of the separate source files.
Note that most Web browsers will make some of the following code look a
little strange, particulary the "#include " statements. The
names of the files will be missing from the screen. Don't worry, they are
still in the source code. It's just that your Web browswer doesn't
understand what it thinks are commands to it like , so it doesn't
display them on the screen.
I decided not to change the code to make it display better: you would have
to undo all of the "fixes" to get good compilable code.
/* file = hfracts.c
*
* writes out a raw binary file of a weird fractal.
* Does not display on screen.
*
* Copyright (C) Terrance E. Hodgins 1996
* All rights reserved.
*
* Permission is hereby granted to the general
* public for people to copy this source code,
* modify it, play with it, and give copies to
* others for non-commercial applications, provided
* that these copyright notices are retained in all
* such copies. Thank you. Have fun.
*/
#include
#include
#include
#include
#include
#include "cmpxmath.h"
/* sets up screen coordinates for a nominal 640x480x256 VGA.
* change if you have something like 1024x768x256 VGA,
* or any other kind of screen. This code should work okay
* on any machine, like a 1280x1024 Sun.
*/
int maxx = 640;
int maxy = 480;
FILE * Fp2 ;
FILE * Fp3 ;
int fract_it( int, int );
main( int argc, char ** argv)
{
int i, ic ;
unsigned int color, x, y, xend, yend, mx2, my2 ;
/* this file is for the image data */
Fp2 = fopen( "hfract.raw", "wb" );
/* this file is for logging the point values, so that you
* can debug things if you wander into some unknown space and
* get lost. Actual writing to it is commented out down in
* "fract_it()".
*/
Fp3 = fopen( "hfract.txt", "w" );
xend = maxx ;
yend = maxy ;
mx2 = maxx / 2 ;
my2 = maxy / 2 ;
for( y = 0 ; y < yend ; y++ )
{
for( x = 0 ; x < xend ; x++ )
{
/* shift coordinates:
* the VGA screen coordinates
* 0 <= x < 640 and 0 <= y < 480
* become
* -320 <= x < 320 and -240 <= y < 240
* for mathematical convenience, to center
* the image. Note that the y coordinate
* is also inverted, so that the image will
* not be (mathematically speaking) upside
* down.
*/
color = fract_it( x - mx2, my2 - y );
/* write out data to file */
fputc( (unsigned char) color, Fp2 );
/* allow bailing out */
if( kbhit() )
/* check that it's really a 'quit' command,
* and not somebody accidentally bumping
* the keyboard...
*/
if(( ic = getch()) == 'q' || ic == 'Q' ||
ic == 27 )
goto done;
}
printf( "Done line %i.\n", y );
}
done:
fflush( Fp2 );
fclose( Fp2 );
fflush( Fp3 );
fclose( Fp3 );
exit(0);
}
/*
* the fract_it() routine does all the work of iterating on a
* point, and returns a scaled color value related to the rate
* of increase (derivative) of the absolute value of the
* complex value of the point as it is repeatedly squared, and
* goes above some arbitrary limit.
* Returns 0 if the point never goes above the limit,
* and 255 for points way above the limit.
*/
/* fract_it()
*
* color = fract_it( x, y );
*/
/* this scales screen coordinates to the complex fractal plane.
* If you change the image size from 640x480 to something else,
* you will want to fudge this too.
*/
#define FX_SCALE (double) 250.0
/* that 250 is for zooming way in. For the big picture, try 3.0 */
/* For way far out pictures, try 0.3 or 0.03, but there, you will be
* hitting the limits of screen resolution, and the outer circles
* will fall into the cracks between the pixels on the screen.
*/
/* For extreme zooming in, increase it much higher, but first,
* make sure that you are centered on something to see. There is
* a lot of black, empty space in there.
*/
/* these offsets center the image.
* Change these to pan left or right, up or down.
* Note that the center of the fractal is not at {0,0},
* it is closer to the following numbers.
*/
#define X_OFFSET (double) -3.750
#define Y_OFFSET (double) -1.3636
/* The limit of number of iterations on a point:
* Only if you have a Cray or a Sun:
* #define FRAC_LIMIT 25000
* better for a PC:
*/
#define FRAC_LIMIT 500
/* or 1000 if you feel more patient...
* But note that you don't seem to gain anything by greatly
* increasing this limit. Sometimes you will get more detail
* in the image by lowering it.
*/
/* the upper limit of absolute value of exponentiated point: */
#define FRAC_BOUNDS (double) 2000.0
/* could be higher, much higher */
/* scales a point's absolute value to a color (0 to 255).
* if you change FRAC_BOUNDS, you will need to proportionally
* fudge this, too.
*/
//#define FRAC_VAL_SCALE (double) 40.0
#define FRAC_VAL_SCALE (double) 20.0
/* You will shift what points get what colors by changing this. */
/* the complex value to add with each iteration */
CMPLX incre = { 2.75, 2.75 };
/* also try { 1.414, 1.414 } or { 1.2, 1.414 }
* If the real and imaginary numbers are not equal, it adds a
* warp to the circles and ovals, and increases the thickness
* of the detailed fringe on the inside edge of each circle and
* oval.
* I haven't had the time to try a lot of combinations of
* these. I would guess that you can produce some pretty wild
* results if you pick the right numbers here.
*/
int fract_it( int x, int y )
{
unsigned int color, i;
double fx, fy, xc, absval;
union
{
CMPLX cpx; /* my complex numbers */
struct complex complx; /* Microsoft's structure */
} cmpx;
CMPLX cx1, cx2 ;
/* convert screen x,y coordinates to complex plane coords. */
fx = (((double) x ) / FX_SCALE ) + X_OFFSET;
fy = (((double) y ) / FX_SCALE ) + Y_OFFSET;
cx1.real = fx; /* initialize complex value for point */
cx1.imag = fy;
for( i = 0 ; i < FRAC_LIMIT ; i++ )
{
/* note that we add first, before exponentiation.
* This is the reverse of the algorithm for the
* more familiar fractal set popularized by Mandelbrot.
*/
cpx_add( &cx1, &incre, &cx2 ) ; /* answer is in cx2 */
cpx_mul( &cx2, &cx2, &cmpx ); /* answer is in cmpx */
absval = cabs( cmpx.complx );
if( absval > FRAC_BOUNDS )
{
/* translate absval to a color */
/* subtracting FRAC_BOUNDS enhances the contrast
* a bit. Is optional.
*/
xc = absval - FRAC_BOUNDS ;
if(( color = (int)((( xc ) / FRAC_VAL_SCALE )))
> 255 )
color = 255;
goto ending ;
}
/* else, continue */
cx1.real = cmpx.cpx.real;
cx1.imag = cmpx.cpx.imag;
}
color = 0; /* we never went above the limit. */
ending:
/* log data values to a file for debugging. Commented out now,
* because it slows the program down. And it produces a huge
* output file, like 4 MB.
*/
// fprintf( Fp3, "%3i %3i %e %e %e\n", x, y, fx, fy, absval );
return( color );
}
#include "cmpxmath.c"
/* end of file */
--------------------------------------------------------------------------
/* FILE = cmpxmath.c
*
* Description: contains routines to do complex addition,
* subtraction, multiplication, division, etc., of complex
* variables.
*
* These routines are for double-precision floating point
* numbers. It is assumed that the complex variables are made
* up of two double-precision floating point variables: one for
* the real part, and one for the imaginary part.
*
* Copyright (C) Terrance E. Hodgins 1996
* All rights reserved.
*/
#include
#include
#ifndef CPX_XDEF
#include "cmpxmath.h"
#endif
/* function = cpx_mul()
*
* Description: does a multiplication of two complex
* values.
*
* Returns: a dummy status int.
* indirectly returns the answer in cpx_out.
*/
int cpx_mul( CMPLX * cpx_in, CMPLX * cpx_2in, CMPLX * cpx_out )
{
double avar ; /* real half of 1st complex variable */
double bivar ; /* imag half of 1st complex variable */
double cvar ; /* real half of 2nd complex variable */
double divar ; /* imag half of 2nd complex variable */
double ac_prod, ad_prod, bc_prod, bd_prod ;
avar = cpx_in->real ;
bivar = cpx_in->imag ;
cvar = cpx_2in->real ;
divar = cpx_2in->imag ;
/* do binomial multiplication:
* (a + bi) * (c + di) = ac + adi + bci + bd(i*i)
*
* = ac - bd + (ad + bc)i
*/
ac_prod = avar * cvar ;
ad_prod = avar * divar ;
bc_prod = bivar * cvar ;
bd_prod = bivar * divar ;
cpx_out->real = ac_prod - bd_prod ;
cpx_out->imag = ad_prod + bc_prod ;
return( 1 );
}
/* **************** END OF FUNCTION ***************** */
/* function = cpx_mcc()
*
* Description: multiplies one complex variable by the
* complex conjugate of another.
*
* Returns: a status int.
* indirectly returns the answer in cpx_out.
*
*/
int cpx_mccjg( CMPLX * cpx_in, CMPLX * cpx_2in, CMPLX * cpx_out )
{
int status ;
double divar ; /* imag half of 2nd complex variable */
/* make complex conjugate out of 2nd var, then multiply */
divar = cpx_2in->imag ;
cpx_2in->imag = divar * (double) -1.0 ;
status = cpx_mul( cpx_in, cpx_2in, cpx_out ) ;
return( status );
}
/* **************** END OF FUNCTION ***************** */
/* function = cpx_exp()
*
* Description: does complex exponentiation, base e
*
* Returns: a dummy status int.
* indirectly returns e**( cpx_in ) as cpx_out.
*
*/
int cpx_exp( CMPLX * cpx_in, CMPLX * cpx_out )
{
double expdv_a;
/* use the algorithm:
* e**(a + bi) = (e**a) * ( e**bi )
* = (e**a) * ( cos(b) + sin(b)i )
* = (e**a) * ( cos(b)) + (e**a)*(sin(b))i
*/
expdv_a = exp( cpx_in->real );
cpx_out->real = ( expdv_a * cos( cpx_in->imag ));
cpx_out->imag = ( expdv_a * sin( cpx_in->imag ));
return( 1 );
}
/* **************** END OF FUNCTION ***************** */
/* function = cpx_add()
*
* Description: adds two complex variables
*
* Returns: a dummy status int.
* indirectly returns sum as cpx_out.
*
*/
int cpx_add( CMPLX * cpx_in, CMPLX * cpx_2in, CMPLX * cpx_out )
{
cpx_out->real = cpx_in->real + cpx_2in->real ;
cpx_out->imag = cpx_in->imag + cpx_2in->imag ;
return( 1 );
}
/* **************** END OF FUNCTION ***************** */
/* function = cpx_sub()
*
* Description: subtracts two complex variables
*
* Returns: a dummy status int.
* indirectly returns difference as cpx_out.
*
*/
int cpx_sub( CMPLX * cpx_in, CMPLX * cpx_2in, CMPLX * cpx_out )
{
cpx_out->real = cpx_in->real - cpx_2in->real ;
cpx_out->imag = cpx_in->imag - cpx_2in->imag ;
return( 1 );
}
/* **************** END OF FUNCTION ***************** */
/* **************** END OF FILE ***************** */
--------------------------------------------------------------------------
/* FILE = cmpxmath.h
*
* Description: contains definitions for double-precision
* complex math routines.
*
* Written by Terrance E. Hodgins 1988-1996
* Copyright (C) Terrance E. Hodgins 1996
* All rights reserved.
*/
#ifndef CPX_XDEF
#define CPX_XDEF 1
/* define the double-precision complex variable */
typedef struct dcmplex
{
double real ;
double imag ;
} CMPLX ;
/* define argument types for routines, for strong type-checking */
/* from file cmpxmath.c : */
int cdecl cpx_mul( CMPLX *, CMPLX *, CMPLX * ) ;
int cdecl cpx_mccjg( CMPLX *, CMPLX *, CMPLX * ) ;
int cdecl cpx_exp( CMPLX *, CMPLX * ) ;
int cdecl cpx_sub( CMPLX *, CMPLX *, CMPLX * ) ;
int cdecl cpx_add( CMPLX *, CMPLX *, CMPLX * ) ;
#endif
#ifndef PI
#define PI (double)3.141592653589793238462643383279502
#endif
/* **************** END OF FILE ***************** */
--------------------------------------------------------------------------
# makefile for hfracts, using MSC 5 or 6. Assumes 80286 or above, with
# math coprocessor. Change as needed.
#
LIB1 = C:\msc\lib
#
hfracts.obj: hfracts.c cmpxmath.c cmpxmath.h
cl -AS -Gs -G2 -On -W3 -Zp4 -FPi87 -c hfracts.c >hfracts.erl
hfracts.exe: hfracts.obj cmpxmath.c cmpxmath.h
link hfracts,hfracts,hfracts,$(LIB1)\slibc7;
--------------------------------------------------------------------------
[Image]
Return to the Fractal page
[Image]
This Web page was written and made by Terrance Hodgins, Willamette Web
Weavers of Portland.
Copyright © 1997 Willamette Web Weavers.