From: mrob@world.std.com (Robert P Munafo)
Subject: Re: Fractal Witchcraft speed?
Summary: I was successful in duplicating its speed; here's how:
Keywords: Mandelbrot algorithm witchcraft iteration
Organization: The World Public Access UNIX, Brookline, MA
Date: Thu, 24 Dec 1992 07:23:23 GMT


This is a response to Tim Wegner's description of his efforts to
duplicate the Fractal Witchcraft "Synchronous Orbit Iteration"
algorithm (which I will call SOI for brevity).
 
It's long, so here's the summary: well-implemented SOI uses
curve-fitting and curve-based interpolation.  mu-molecules (a.k.a.
islands or midgets) cause the distortion; skip to the bottom if you
want to know why.
 
For those of you that don't get _Amygdala_ (and why don't you?? :-)
here are some details from Steven Stoft's article that I feel I
should mention. Stoft (who is one of _Witchcraft_'s creators)
describes "the four corners of a square that more than covers the
entire [desired] image" and instructs the reader to iterate the
corners until "they begin to deviate slightly from a perfect square".
Stoft suggests that an error amounting to 1/10 of the distance
between adjecent pixels is a reasonable limit, beyond which one
should break the square up into "16 smaller squares".  The
coordinates of Steven's sample image are center_real =
-0.1049775673880424, center_imag = 0.9270413637809322, size =
0.000000000012 = 1.2e-11, max_iterations = 10000.
 
I read this back in mid-September and immediately set out
implementing the algorithm.  As per Stoft's description I measured
the edge lengths of the iterated corners and went until the ratio of
longest-edge to shortest-edge got bigger than (1 + (1/16 *
1/pixel_width)) where pixel_width is the width in pixels of my image.
(I decided to cut off at 1/16 pixel rather than 1/10.)  This turns
out to be just a first approximation and I shall call it the "Rhombus
Method" because you are essentially tracking a rectangular array of
iterated points and watching to see when it stops being a rhombus.
 
As you found, Tim, this method is woefully inadaquate.  The rhombus
is a very poor model of the iterated point set.  Iterating the sample
image, I got 270 iterations before having to split the first time, a
total of 343 iterations before the second split, and 613 before the
third.  (This is just a 160x160 pixel image; expect much fewer
iterations for 640x480).  By comparison, the average number of
iterations for the pixels in the final image is 5453.8.
 
The next method I tried is the "Quadrilateral Method".  I removed
assumptions about edge lengths but required that interior points be
closely approximated by linear interpolation.  To get an arbitrary
point in the interior you must first interpolate along two edges to
generate two intermediate points, then interpolate between these two
points to get your desired point.  To use the Quadrilateral Method
you must now iterate about 8 points: the four corners ("key points"),
plus about four "test points": on each iteration you try to guess the
values of the test points by interpolation and see how close your
guess is to the actual iterates of the test points.
 
The Quadrilateral Method is quite an improvement, it turns out.  This
time I got 541 iterations before the first split, then 1082, then
2160, 4320, 5130, 5385.  I split into 4 pieces each time, so the
last step had 1024 pieces.
 
The next method involves the use of a curve-fitting function
(quadratic, cubic, spline, take your pick) to model the four edges of
the iterated point set.  The type of curve you pick must lend itself
well to interpolating an arbitrary point inside the deformed set. 
You must now iterate at least eight "key points" because each edge
needs at least three key points to serve as the parameters for curve
fitting, and you should iterate an equal number of "test points".
 
I did not implement a Curve Method because at this point in my work I
discovered an extremely beautiful feature (nested imbedded Julia
Sets) and have been creating QuickTime movies of them ever since. (An
example of one of these, which I'd like to call the "Munafo Midget",
is center_real = -1.769110375463767385, center_imag =
+0.009020388228023440, size = 2.0e-17, max_iterations = 10000.  Zoom
out to size = 1.0e-12 to see one imbedded Julia Set, and zoom out to
size = 5.2e-10 to see the other. You really need to use the Distance
Estimator method to do justice to these images.)
 
-> The math behind SOI (very informal: don't complain, sci.math folks!)
 
SOI is a lot of work for the computer.  The formulas used to implement
SOI are essentially equivalent to the deriviative (or differential?)
of the Mandelbrot iteration function.  For curve-fitting methods you
are effectively iterating the second deriviative as well.  The set of
points you are tracking bounces around the complex plane; on each
iteration it is rotated and the edge furthest from the origin is
stretched a little; at the same time a slight nonlinear curvature is
introduced.  After P_min iterations (P_min being the period of the
smallest-period mu-atom contained in the original rectangle) your
iterated set will straddle the origin; when this happens you know
that the next iteration will introduce a singularity.  (For example,
the iterations-before-splitting in the Quadrilateral Method test
above, 541, is the period of a mu-molecule located near, but not in,
the image.
 
Because of this effect of mu-atoms on the iterate set distortion, SOI
does not perform as well on images that contain low-period mu-atoms.
This is the main reason why SOI works so poorly on low-magnification
images.

Robert P. Munafo, mrob@world.std.com   (lost my signature file!)
Newsgroups: sci.fractals