AOH :: DRAW30.TXT Fractal Drawing Styles 3.0
```
Fractal Drawing Styles 3.0

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.