Tony Finch - Smooth colouring is the key to the Mandelbrot set

dotatfanf wrote
on 10th November 2010 at 15:54
Previous Entry Share Next Entry

Smooth colouring is the key to the Mandelbrot set

At Cambridge Geek Night 6 on Monday I gave a brief talk about what has been occupying some of my spare time in recent weeks. Since Benoît Mandelbrot died last month I have been playing around with making pretty pictures of Mandelbrot sets, culminating (so far) in a movie of a deep zoom targeted near z = -0.15 + 1.04i. The software and materials for my talk are at http://dotat.at/prog/mandelbrot/.

One of the nicest renderings I did when I was a teenager was a 300dpi laser-printed binary decomposition of the level sets of the exterior of the Mandelbrot set. It took a whole afternoon to generate using BBC BASIC on my Archimedes, which was the best tool for numeric code I had available. I wanted to be able to produce smoothly coloured renderings like the ones in the fractal art books by Heinz-Otto Peitgen and his collaborators, but a few things prevented me. Firstly the lack of a computer that could display true colour images :-) and secondly the difficulty of the mathematics. If I remember correctly, I just about grasped enough of the distance estimation method to understand that it requires significantly more work per iteration, which also made it unattractive on the slow computers I had to hand.

Revisiting the Mandelbrot set now, I have the luxury of computers that are thousands of times faster and thousands of times more colourful and millions of times more accurate with their arithmetic. The Internet also provides more accessible resources on the mathematics of the Mandelbrot set. A page by Linas Vepstas was particularly helpful with the smooth colouring problem. Even better, it turns out that his simple recipe arises from the mathematics that gives us our deepest understanding of the Mandelbrot set.

First some basic notation to make the formulae less cluttered.

To generate a Julia set, c is a fixed parameter and z varies across the complex plane. For the Mandelbrot set, c varies across the plane and z starts its iteration at 0. In both cases we iterate until z escapes to infinity, in which case it is outside the set, or remains bounded, in which case it is inside. Because we are most interested in the behaviour of z as it is iterated, we usually leave the parameter c implicit.

The simplest Julia set, when c = 0, is just the unit circle. In this case we can write down a direct formula for zn:

There is something called a Böttcher map which transforms the exterior of Julia set c = 0 into the exterior of other connected Julia sets and back again. It's like a Riemann map, but it operates on the complement of the set.

The definition of the map is as follows. Observe that Φ0(z) is the identity function.

Adrien Douady and John Hubbard showed that the Böttcher map Φc(0) is a conformal mapping from the exterior of the Mandelbrot set M to the exterior of the unit disk, and therefore the Mandelbrot set is connected.

The Böttcher map allows us to treat the exterior of the set as an electrostatic field. Imagine an infinitely long crinkly cylinder whose cross-section is the set. If you electrically charge the cylinder it creates a field in the surrounding space. The Böttcher map tells you how the simple circular equipotential curves and radial field lines of the Julia set c = 0 correspond to those surrounding the Mandelbrot set and the other connected Julia sets.

This is a powerful tool for understanding the Mandelbrot set. The binary decomposition makes a lot of the field structure visible, since the boundaries between the black and white areas are formed from segments of field lines and equipotential curves. If you express the angle of a field line as a binary fraction of a whole turn, you get the sequence of black and white segments that the field line crosses in the binary decomposition. Field lines with rational angles correspond to pinch points on the Mandelbrot set.

For our purposes we are less interested in field lines and more interested in the Douady-Hubbard potential function. We use the Böttcher map to warp the space around the set into its simplest form. This allows us to use the standard formula for the electrostatic potential around an infinite wire, which is the log of the distance from the wire. The base of the logarithm is related to how highly charged the wire is; we can choose it for our convenience.

Since we don't have time to count to infinity, in practice we define a bail-out radius R and stop iterating points that get further than that from the origin. The iteration count n is used for the flat blocky level-set colouring scheme.

The way to get a smooth colouring is to ask what fractional iteration value ν would make the point land exactly on the bail-out radius.

We can use the potential function to formalize this idea. If R is large, varying it will not make much difference to the potential we calculate for a point. So to derive ν we assume that we get the same result from the two expressions for the potential in the first line below. (I think this is what causes the small errors and discontinuities discussed by Linas.)

The final term in that equation is a constant, since R is a fixed parameter of our implementation. We might as well throw it away since it does not add any useful information.

And thus we get the equation to use for smoothly colouring of the exterior of the Mandelbrot set.


(Leave a comment)
From:ewx
Date:2010-11-10 16:24 (UTC)
(Link)
Good stuff. I foresee an attempt at 32.96 fixed point log in my future l-)
(Reply) (Thread)
From:crazyscot
Date:2010-11-10 17:08 (UTC)
(Link)
... how deep does that let you go, and what will you find when you get there?

I'm currently seeing precision artefacts at the same place they crop up in fanf's video - 2^52 zoom or so - which tallies with scraping the bottom of double, and I'm wondering whether there's a way to reliably detect that that is happening and disallow further zoom. Possibly (L)DBL_EPSILON is going to be my friend...
(Reply) (Parent) (Thread)
From:ewx
Date:2010-11-10 17:19 (UTC)
(Link)

All the interesting stuff has magnitudes between about 0.25 and 2, so with double your pixels can't go much below 2-53 or so and the 'floating' in floating point is of little use - fixed point is what you need. So 96 fractional bits should get you around 40 2x zooms beyond the double version.

32 bits for the integer part is a bit excessive but it makes the C implementation easier to put the units column at the bottom of a word. (That might not be true for assembler though, I've not thought this through. Another 16 bits wouldn't hurt.)

(Reply) (Parent) (Thread)
From:fanf
Date:2010-11-10 17:22 (UTC)
(Link)
Just calculate the size of a pixel and put a lower bound of about 2^-51 on it. Although you can in principle zoom in further if you are closer to the origin, that buys you only three or so binary orders of magnitude so it's probably not worth it.

32.96 would allow me to make a 6 minute movie zooming at a factor of 2 every 4 seconds (i.e. 1-2^-100 per frame).
(Reply) (Parent) (Thread)
From:fanf
Date:2010-11-10 23:21 (UTC)
(Link)
The code is slightly simpler if you implement log to the base two - see my addendum below.
(Reply) (Parent) (Thread)
From:ewx
Date:2010-11-14 18:14 (UTC)
(Link)
I’ve added smooth coloring to mine, though currently it cheats a bit and does the logs in double. Still, a fixed point log2 sounds like a fun exercise sometime.
(Reply) (Parent) (Thread)
From:nonameyet
Date:2010-11-10 22:32 (UTC)
(Link)
Was the target point chosen "at random", after a some manual searching or based on some clever maths ?

Did you just pick the most appropriate compression algorithm available in MPEG-4 or did you analyze the special features of this video ? I'm fairly sure that it must have coherencies that make a custom compression algorithm incredibly good (though not so much use unless everyone has the decompression algorithm).
(Reply) (Thread)
From:fanf
Date:2010-11-10 22:49 (UTC)
(Link)
Manual search, I'm afraid :-) The precise co-ordinates are buried in mandelbrot.c.

I initially tried making the movie with ffmpeg's own h.264 compressor and the result was disastrous. I installed x264 (which required a new binutils so that the assembler could support SSE3 instructions!) and recreated the movie and the result looked fine and was also smaller.
(Reply) (Parent) (Thread)
From:fanf
Date:2010-11-10 23:17 (UTC)
(Link)

Mandelbrot set implementations coded without a complex number abstraction usually deal with the square of the modulus rn2 since that saves the cost of a square root in each iteration. When you use r 2 in the formula for the smoothed iteration count, the log allows you to replace the square root with a division:

n' = n - log2(log(r 2) / 2)

This division can be eliminated for the same reason as the log(R) term:

n' = n - log2(log(r 2)) + log2(2)
n'' = n - log2(log(r 2))

If you only have natural logarithms,

n'' = n - log(log(r 2)) / log(2)

If you only have logs to the base 2, then you can play the constant offset game again to get

n''' = n - log2(log2(r 2))


Edited at 2010-11-10 23:20 (UTC)
(Reply) (Thread)
From:dwmalone
Date:2010-11-12 22:57 (UTC)
(Link)
I produced some animations years ago on the Amiga and then on an IBM RS6000 that I got access to in college. I still have them, but they are in Amiga ANIM format, so I'm not sure if even xanim will play them any more.

http://www.maths.tcd.ie/~dwmalone/maths/zoomanim.mov
http://www.maths.tcd.ie/~dwmalone/maths/cricanim.mov
http://www.maths.tcd.ie/~dwmalone/maths/FractAnim.mov
http://www.maths.tcd.ie/~dwmalone/maths/Baranim.mov

The first is a straight forward zoom. The second morphs a Circle (i.e. Julia set with constant 0) into dust by changing the constant value. The third morphs the regular Mandelbrot set into the z -> z^3 Mandelbrot set by varying the power between 2 and 3. Because you have some choices about which branch of the power to take for non-int powers, you get some interesting edges. The final one morphs the z -> conj(z^2) Mandelbrot set into the cubic one in a similar way. The Mandelbrot set with the conjugate Mandelbrot set is quite interesting 'cos it seems to have both smooth and fractal bits.
(Reply) (Thread)
From:fanf
Date:2010-11-14 21:18 (UTC)
(Link)
Nice. I managed to play them using a Java program called MultiShow http://www.randelshofer.ch/multishow/
(Reply) (Parent) (Thread)
From:ext_329527
Date:2010-11-21 17:07 (UTC)

Nicely colored visualization of the Mandebrot Set

(Link)
There's a nicely colored visualization of the Mandebrot Set here: http://mathcadworksheets.blogspot.com/2010/11/generate-colorful-images-of-mandelbrot.html
(Reply) (Thread)

(Leave a comment)

Powered by LiveJournal.com