Fractal Drawing Styles 3.0

 About this file
        This file discusses five different ways to color fractals. See
"The Science of Fractal Images" for the mathematical background on
Binary Decomposition, Continuous Potential, and the Distance Estimator.
All five have been implemented in a program called Mandella for the
Macintosh.
	
 The fractal formula
        Almost all fractals are generated using complex numbers. Every
complex number can be represented by a point on the plane. To generate a
fractal we iterate something like this z := z^2 + c
        Where z is a point in the plane and c is a complex constant.
When this is iterated one of two things can happen: z can blow up to
infinity or z can become periodic. If z never goes to infinity it is
inside the fractal.
        To simplify things we use a bailout value to see if z is headed
for infinity and a number called maxIter. If the number of iterations
exceeds maxIter we assume the point is periodic and therefore in the
fractal.
        After iterating a point we can color it based on the number of
iterations or on the final z. Most fractal programs use the number of
iterations to assign colors to points. This is called the "Level Set"
method.

 Variables
        z       starts as a point in the complex plane and is then iterated
	orbits	an array that saves the value of z for each iteration
        count   the number of iterations
	hiCount	the maximum number of iterations

1) Level Sets
		This method uses the count value to index into the palette:
			index := (NoColors*count)/hiCount + FirstColor
        However this will produce poor results. There are two problems
        with this approach: 1) If the image has been zoomed the lowest
        count in the image could be quite large. This means that a
        possibly large portion of the palette is never used. 2) The high
        count values cluster close to the set. This causes a large
        amount of the palette to be dedicated to a very small area next
        to the fractal. This problem will become worse if the hiCount is
        increased.
                To alleviate these problems we can use three cutoffs:
        high, middle, and low. The high cutoff will be used for the
        hiCount. If the count is between the low and middle cutoffs we
        use the count value to index into the palette. if the count is
        greater then the middle cutoff we use the last color in the
        palette. If the count is less then the low cutoff use white.
        Putting it all together we have:
                        IF count < low THEN
				index := White;
			ELSIF count >= high THEN
				index := Black;
			ELSIF count > middle THEN
				index := LastColor;
			ELSE
				range := LONG(middle - low);
				delta := LONG(count - low);
				index := (delta*NoColors) DIV range + FirstColor;
			END;
                This is much better but we have only partially solved
        the problem with the counts clustering close to the fractal. To
        bring out the most detail we need to somehow use the bulk of the
        palette for the lower count values. One possibility is to use a
        formula like this:
                        newCount := hiCount*(count/hiCount)^(1/root)
        This formula will take a count from 0 to hiCount and return a
        new count value in the same range if root > 0. If root = 1 then
        newCount = count. If root > 1 then the low counts account for
        the majority of newCounts. This is what we want since the low
        counts will then be assigned to most of the palette. If root < 1
        most of the newCounts come from the higher count values. This
        can sometimes be interesting. To integrate this with cutoffs we
        can write:
                                range := LONG(middle - low);
				delta := LONG(count - low);
				new := range*(delta/range)^(1/root);
				index := (new*NoColors) DIV range + FirstColor;

                Another idea is to allow the user to somehow draw a
        graph to control the count to palette mapping. An easy and
        powerful method is to use Bezier curves. The December 1986 Byte
        has a good article on Bezier curves.
	
        The rest of the drawing styles return real numbers (x). To index
into the palette the real number must be somehow scaled into a number
from 0 to 1. Then the palette index is simply: index :=
(NoColors-1)*scale + FirstColor. To do the scaling we need to find the
largest x. Then scale := x/maxX. In an interactive program maxX can
sometimes be hard to find.

        Many programs using the Level Set method store each count value
in an array. This array can then be used to redraw the image (perhaps
after a change in the root scaling). To integrate this with the other
drawing styles we can multiply scale by some large number Unity. Then
0..Unity will be inside the palette. And Unity+1 will be inside the set.
	
2) ColorSet
                Here we color the fractal and leave everything outside
        the fractal black. To get a color we use the final z value. This
        style works well with Julia type fractals. The code looks like
        this:

                                IF count >= hiCount THEN
					mag := z.r*z.r + z.i*z.i;
					scale := mag - TRUNC(mag);	(* ignore integer part *)
					count := Unity*scale;
				ELSE
					count := -1;					(* use background color *)				
				END;
					
	
3) Binary Decomposition
                In this style we treat z as a point in the plane and use
        the angle z makes with the x axis to index into the palette:

                                IF count >= hiCount THEN
					count := Unity+1;
				ELSE
					theta := ArcTan(z.i/z.r) + PI/2;
					scale := theta/PI;
					count := Unity*scale;
				END;
				
4) Continuous Potential
                Here we use the magnitude of z and the count value to
        get a smoother image. The potential is calculated like this: pot
        := Ln(Mag(z))/(2^count). To scale the potential into the range
        [0..1] we need to find the largest potential: maxPot :=
        Ln(bailout)/(2^loCount). Then scale := pot/maxPot. As count
        increases scale will go rapidly to zero. To smooth things out we
        will use a variable called slope. Using a large bailout value
        will also help matters (500 works well for the classic
        Mandelbrot fractal). The code follows:

                                IF count >= hiCount THEN
					count := Unity+1;
				ELSE
					pot := Ln(Mag(z))/(2^count);
					maxPot := Ln(bailout)/(2^low);
					scale := pot/maxPot;
					scale := scale^(1/slope);
					scale := 1 - scale;				(* shift high counts to the right of palette *)
					count := Unity*scale;
				END;
				
5) Distance Estimator
                This method is often used to create monochrome images.
        It uses orbits to determine how close a pixel is to the set. If
        the pixel is less then delta it is colored black. Delta is
        threshold times dx. The user should input threshold. For the
        entire Mandelbrot set 0.6 works well. As the user zooms in
        threshold should be decreased because the pixels should be
        closer to the set. Dx is the distance between pixels. See pages
        196-198 of the the Science of Fractal Images for more details.

                                overflow := 1.0E32;
				delta := threshold*dx;
				IF count >= hiCount THEN
					count := Unity+1;
				ELSE
					zPrime.r := 0;
					zPrime.i := 0;
					k := 1;
					WHILE (k <= count) AND (ABS(zPrime.r) < overflow) AND (ABS(zPrime.i) < overflow) DO
						zPrime := Mult(orbits[k], zPrime);
						zPrime.r := 2*zPrime.r + 1;
						zPrime.i := 2*zPrime.i;
						INC(k)
					END;
					IF k <= count THEN
						count := Unity+1;					(* very close to set *)
					ELSE
						mag1 := Sqrt(Mag(z));
						mag2 := Sqrt(Mag(zPrime));
						dist := Ln(mag1)*(mag1/mag2);
						IF dist < delta THEN
							count := Unity+1;				(* close to set *)
						ELSE
							count := 1;						(* far from set *)
						END;
					END;
				END;
		
                Instead of just using white when the distance is greater
        then delta we can use distance to index into the palette. Since
        we do not know the maximum distance we'll cheat and use a real
        number version of the modulus function:

                                dist := dist/delta;                                     (* scale dist *)
				scale := ABS(Remainder(dist);		(* scale is in [0, 1] *)
				scale := ABS(scale - 1);			(* invert scale *)
				count := Unity*scale;
				
		The Remainder function looks like this:
				PROCEDURE Remainder (x, y: REAL): REAL;
					VAR quotient: REAL;
				BEGIN
					quotient := Rint(x/y);		(* round to integer *)
					RETURN x - (y*quotient);
				END Remainder;
	

*** Change History ***

Version 3.0
	1) Instead of using the Remainder function ColorSet uses the fractional
		part of Mag(z).
        2) Continuous Potential uses a dedicated variable (slope) to do
        root scaling.
        3) The Distance Estimator uses a slightly different method for indexing
		into the palette.
	4) The max constant has been renamed Unity.
	
Version 2.0
	1) First version posted.
	
	
	
Address comments and suggestions to:
	Jesse Jones 			
	
	Phone:	(206) 328-3966
	BIX:		JesJones
	CServe:	73627,152