2010-02-26

And now for something completely different! The Mandelbrot Set!
This is the basis for all those paisley computer fractal pictures we see
everywhere, especially when they were popular in the 1990s. For your viewing
pleasure, **a Mandlebrot explorer app is included**, with (Flash) source,
**right down there in this blog post! Woo!**

In this article, I'm going to talk about how the images are generated, and use that as an example of something I find I do quite commonly: remapping number ranges. High excitement, there.

I do apologize in advance that this post is a bit *mathy*. Now,
I was never particularly good with math, so I figure that if I'm
grasping it, it must be pretty straightforward (the actual manipulations
never get beyond elementary
algebra.) But it's not as easy reading as, say, *Captain
Blood*.** If you don't care for math, I encourage you to just
jump, guilt-free, down to the pseudocode and the app!** Enjoy!

First, let's talk about the Set itself. What is it? It is this: the
set of complex
numbers that, when run through a iterative function over and over,
never grow above a certain value. If you subject some complex number to
this treatment, and keeps getting bigger, that number is *not* in
the Mandelbrot Set. If the number remains forever small, then it
*is* in the Set.

(An analogy is how if you keep squaring the result of squaring a number over and over, you'd expect that number to keep getting larger forever, e.g. 2, 4, 16, 256, etc. And that's true, unless the number is less than one, in which case it never gets large and goes to zero e.g. 0.2, 0.04, 0.0016, 0.00000256, etc. The number either grows to infinity or doesn't, depending on what number it starts with. It's the same with the Mandlebrot Set, except we're doing it with a slightly more complicated equation, and we're using complex numbers.)

I'm glad you asked! A complex number is a number that is made up of both a real part and an imaginary part. These parts are added together and that makes up a single complex number. The imaginary part is multiplied by

i, the imaginary unit number (which is defined to be the square root of -1, whatever that means.) So you have a complex number likex:

x= 12 +i3In that case, the real part of

xis 12, and the imaginary part isi3.Why do you need such a thing? Well, there are plenty of applications for complex numbers (similar to the way, as Wikipedia points out, there are plenty of uses for negative numbers—even though you cannot posses a negative number of, say, aardvarks.) In this case, we'll use them to create pretty pictures.

One more thing: notice that if you put the real number on a horizontal axis, and the imaginary number on a vertical axis, you can plot complex numbers on what is called the

complex plane. This will become apparently important in just a moment.

So what is this magical series that we apply to complex numbers to
see if they're in the Mandlebrot Set or not? Well, for a given complex
number, *c*, we watch the values of the complex number *z* in
the following equation:

z_{0} = 0 (complex zero)

z_{n+1} = z_{n}2 + c

What we do is calculate *z _{n}* up to a
fairly large

You can skip this sidebar, if you want—it shows how to go from the above equation to the math the algorithm uses. If you hate that kind of stuff, then you can just jump ahead to the algorithm.

How do we do that math with the complex numbers? It's pretty straightforward, like elementary algebra. Let's say you have these two complex numbers

zand_{n}c, stolen from the above Mandelbrot equation (the values ofx,y,a, andbare arbitrary for this example):

z_{n}=x+iy

c=a+ibAnd we want to compute the next number in the sequence

z_{n+1}:

z_{n+1}=z_{n}^{2}+c(Mandelbrot equation)Then we just expand it like with regular algebra and go. Let's start with the z

_{n}^{2}part:

z_{n}^{2}= (x+iy)^{2}[from definition ofz_{n}, above]=

x^{2}+i2xy+i^{2}y^{2}[expand]=

x^{2}+i2xy-y^{2}[sincei^{2}= -1 by definition]= (

x^{2}-y^{2}) +i2xy[consolidate real and imaginary parts]That's a complex number, see? It has real part (

x^{2}-y^{2}) and imaginary parti2xy. We can finish up here by adding onc, which we do by just adding in the real and imaginary components:= (

x^{2}-y^{2}) +i2xy+c= (

x^{2}-y^{2}) +i2xy+a+ib[from the definition ofc, above]= (

x^{2}-y^{2}+a) +i(2xy+b) [consolidate real and imaginary parts]The real and imaginary parts of this final equation are clearly visible in the inner

whileloop in the pseudocode, below.

So what is the magic bounding range that determines if a point is in the Set or not? It is this: a circle of radius 2 on the complex plane...

That is, take a complex number *c* and use it to calculate a
bunch of iterations of *z* out to *z*_{1000}. If *z*_{1000} falls
on the complex plane within a circle of radius 2 (from the origin), then
*c* is in the Mandelbrot Set! (Well, it *probably* is. Maybe
it would have broken out on iteration 1001, but we have to give up
sometime, and at some point we just assume that it's never getting
bigger.)

Let's take an example where you calculate *z*_{1000} for a starting value of *c* and it comes
up with a complex number with real part *x* and imaginary part
*y*:

*z*_{1000} = *x* +
i*y*

And remember that the equation for a circle is:

*x*^{2} + *y*^{2} = *r*^{2}

where *r* is the circle radius. That means that for our
starting value of *c* to be in the Mandlebrot Set, the point
*z*_{1000} on the complex plane must fall
within a circle of radius 2, so this must be true:

*x*^{2} + *y*^{2} ≤ 2^{2}

If that becomes false at any point during our calculation, we can
bail out and be certain that *c* was not in the Mandlebrot Set.
This condition is clearly visible in the pseudocode's `while`
loop, below.

Which brings me to an important point. Some values of *c* lead
us out past the radius 2 circle *faster than others*. We can count
how many iterations it took for *z* to break the radius 2 barrier,
and *use that count to choose a color for the pixel* that
represents our number *c*. This coloring scheme is called the *escape
time algorithm*; there are other
algorithms that produce different-looking results.

If the number never breaks radius 2, it is in the Set, and historically these pixels are colored black. This means that the Set is uniform black, and the periphery of the set is what makes all those squirrely lines and colors!

Finally, it's **pseudocode** time for the basic algorithm! Let's
do it!

Pseudocodefor x,y in (all pixels on screen): # First, we need to convert from screen coordinates to the # coordinates on the complex plane they represent. We want to store # this result in our number "c". We'll map x to the real part of c, # and y to the imaginary part of c: cReal =computeRealPartFromX(x) # [these implemented below] cImg =computeImaginaryPartFromY(y) # And we'll start z at 0: zReal = 0 zImg = 0 # Remember our count: count = 0 MAXCOUNT = 1000 # how high we're willing to go # We want to count up until we break out of the radius 2 barrier, or # until we hit our count limit, which I've chosen to be 1000 in this # example: while (zReal*zReal + zImg*zImg <= 2*2) AND (count < MAXCOUNT): # compute the next iteration's values: nextZReal = zReal*zReal - zImg*zImg + cReal nextZImg = 2 * zReal * zImg + cImg zReal = nextZReal zImg = nextZImg count = count + 1 # Out of the loop, let's check the results. If we hit MAXCOUNT, # then it's IN the Set, and we'll color it black: if (count == MAXCOUNT): # in the set color = black else: # NOT in the set color =mapCountToColor(count) # [discussed below] # Finally, draw the pixel: setPixel(x, y, color)

And that's all there is to it. That's the basis for all those pretty fractal pictures!

Of course, I glossed over a couple functions there, and this is what is leading to Part II of this post. (I know, you're thinking, "Gaahhh! That was only PART ONE?" Well, this part is much shorter and straightforward.

Nevertheless, let's have an intermission! Check this thing out! Click on "Restart" to get it going, then click on the bitmap to zoom in, or the "Undo Zoom" button to zoom back out. (Click here for a 1000-pixel-wide version of this app.)

Things to note:

- I use Adam7 interlacing to present a rough image as quickly as possible. This ties in to my earlier post on interlacing. You can interrupt the render at any time and keep zooming in.
- I use cooperative multitasking during the render, also the topic of a previous blog entry. The app looks how long a "frame" is, and makes sure it doesn't use more time than that during the render phase. If it runs out of time, it cooperatively yields control, then picks up again on the next frame where it left off. (This is extra-necessary, since Flash is single-threaded.)
- If you zoom in really really far, you'll eventually come up
against the resolution of the floating point numbers that are used
in the calculations. This proves that humans are better than
computers! [
*You just keep believing that, Beej —Computer*] - It doesn't take long to zoom in really far, and there's a lot of space in there. By the time you start seeing floating-point artifacts, you've blown up the original image so much that it would stretch from your desk to somewhere out past Jupiter.
- This could be done in HTML5. I wanted it to run on most
platforms, though, and Flash more easily embeds in the blog. That
being said, a port to using the canvas element would be pretty
straightforward, and I'd love to hear
if someone does it [
*Update*: Kevin did it!] You can see my previous blog entry on pixel manipulations with the canvas. And my Flash source for this app is in the public domain... **Here is a link to the source code**. It's straight ActionScript compiled with mxmlc, but you could use FlexBuilder, I imagine. Also see my other blog entry on getting started with Flash.

Welcome back! Now the reason I want to talk about mapping one number range to another number range is because I do so much of it in the Mandlebrot explorer app.

Let me first give you a straightfoward example:

If you have a variable that ranges from 0.0 to 1.0, and you want to remap it to a number that ranges from 0.0 to 212.0, you can just multiply it by 212, right? You've just remapped [0.0,1.0] to [0.0,212.0]. Easy.

One more:

If you have a variable that ranges from 0.0 to 3490.0 and you want to remap it to a number that ranges from 0.0 to 1.0, you can just divide it by 3490.0, right? You've just remapped [0.0,3490.0] to [0.0,1.0]. Easy, again.

What this means is that you can map from any range zero-to-*n*
to the range zero-to-one. And you can map the range zero-to-one to the
range zero-to-*m*. Therefore you can map any range zero-to-*n*
to zero-to-*m*. Any number to any number. (If you then add a
number, say *x*, to the result, then you've mapped to
*x*-to-(*m*+*x*).)

It's a simple matter of scaling and ratios, when you look at it, right? I find that thinking about the mapping passing through the range zero-to-one helps me see how it's working. Of course, you can roll all this into a single equation when you code it up.

This is used all the time when you have a random number generator that returns a value in the range 0-to-1. You see code like this, where

random()is a function that returns a random number between 0.0 and 1.0:// remap a random number from 0..1 to 0..10: f = random() * 10.0; // remap random number from 0..1 to 15..22: g = (random() * 7.0) + 15.0;

In the Mandelbrot generator, we want to map a pixel coordinate on the screen to some real and imaginary coordinates on the complex plane. Which complex coordinates in particular depend on which part of the Set we're "zoomed in" on.

To make this happen, we'll do some remapping, like so:

Pseudocode# The complexPlaneWidth is how much of the real axis we can see # in our bitmap viewport; making this number smaller would mean # more "zoom". And changing the complexPlaneLeftEdgeCoord would # mean panning left and right: complexPlaneWidth = 3 complexPlaneLeftEdgeCoord = -1 # And here is our bitmap width in pixels: bitmapWidth = 400 # The following function converts an x pixel coordinate on the # bitmap plane into the corresponding real coordinate on the # complex plane: computeRealPartFromX(x): # x was in the range 0-to-bitmapWidth; # we want to remap it to the range 0-to-1: f = x / bitmapWidth # floating-point calculation # Now f is in the range 0-to-1; remap it to the range # 0-to-complexPlaneWidth: realPart = f * complexPlaneWidth # Finally, we need to shift it over from the range # 0-to-complexPlaneWidth # # to # # complexPlaneLeftEdgeCoord-to- # (complexPlaneLeftEdgeCoord + complexPlaneWidth) realPart = realPart + complexPlaneLeftEdgeCoord return realPart

(Note that when we divide `x` by `bitmapWidth`, we need
this to be a floating point calculation, not an integer calculation,
since the integer calculation would be zero.)

It's a very similar situation with mapping the *y* coordinate to
the imaginary part of the complex number.

And in fact, it's a very similar situation when we map from an escape value to a color. The escape value goes from 0 to 1000, but we need the red, green, and blue values to go from 0 to 255. (The actual mapping I used for the app wasn't a straight linear mapping, though—for that, I just messed with the numbers until I got something visually trippy. There are, of course, an infinite number of ways to "theme" the output.)

And it's somewhat similar when we map a pixel number to an Adam7 interlaced coordinate.

And when we map a complex number *c* on the screen to a
Mandelbrot result *z*_{1000}, too.

Sure, not all these are linear mappings, but look at how the data flows as we start with pixel 0 and run through all the pixels in the image:

- Take the number of the next pixel to process, call it
*n*. - Use the Adam7 mapping to convert
*n*to an*x*,*y*screen coordinate. - Convert the
*x*,*y*screen coordinate to a complex plane coordinate's real and imaginary parts, represented by a complex number*c*. - Run the Mandelbrot iteration process on
*c*, producing subsequent values of the complex number*z*. Keep track of the iteration number in*t*. - If the point
*z*gets outside a circle of radius 2 on the complex plane, map the iteration number*t*to a color*q*. Else set*q*to black. - Paint the pixel at
*x*,*y*the color*q*. - GOTO step 1

It's a whole lot of munging and mapping from one representation to another. A computer's work is never done!

Nice article about Mandelbrot algorithms. I have a few comments.

I'm not really sure about how your mapCountToColor function really works. I've seen some other renderings, all with very different color patterns, even though the black part is the same. Perhaps you can tell us more about coloring algorithms? After all, the color is a very important part of the fractal rendering.

From your pictures, it looks like colors are somewhat resembling gradients, but not true gradients as you can still see contour lines between different shades of a color. Is this what's supposed to happen?

@luckytoilet That's a good idea to add more on that topic, since I glossed over it pretty much, huh. :)

A simple coloring algorithm might take the escape number (the iteration count) and remap it to the red channel (that is, remap the 0..1000 count to a 0..255 color value, assuming 1000 is your max escape number):

red = 255 * count / 1000

green = 0

blue = 0

That would make a red image, generally dark.

The Beej Coloring Technique used above is:

red = (count * 8) MOD 255

green = (count * 9) MOD 255

blue = (255 - count * 4) MOD 255

I didn't really plan any of that. I wanted it to look zany, and so I just threw down some numbers and then tweaked them. 8 and 9 are relatively prime, so my gut tells me this would be good for making interference patterns (but I have no idea if that's really true). And I made blue walk backward just 'cause. Making the multipliers larger packs the contour lines on the final image closer together.

Yes, as for those lines, I implemented the discrete time escape algorithm, as opposed to the continuous one, so it is "supposed" to be that way. But it certainly doesn't have to be. See the comparison here:

http://en.wikipedia.org/wiki/Mandelbrot_set#Continuous_.28smooth.29_coloring

The continuous one isn't actually much more complex at all, but it's one more thing to explain in the (already lengthy) post, and I wanted the demo to have the "classic" look.

I wrote a Mandelbrot and Julia set viewer a few days ago using the HTML5 canvas tag. You can see it here:

http://www.scale18.com/canvas2.html

at last i have found an understandable explanation of the iteration of the mandelbrot set

many thanks !

I am thinking about complex numbers , that their existence is somewhat arbitrary definition .

Could i similarly regard them simply as numbers that are 'positions' relative to a point in a 2D space ? ie (a+ ib) is simply (x,y) in the 'conventional' x,y , plane .

Am I missing something ?

Does the fact that i=square root(-1) matter ?

@michael hilton In the case of the Mandelbrot set, for viewing purposes, we do simply regard them as that—it's a mapping from the "complex plane" to the XY Cartesian plane so we can see it. It's just a representation we've chosen, just like you can choose to plot, say, ages of fish on either a scatterplot or a histogram.

But the nature of the Set is tied up in the mathematical behavior of complex numbers, so it's a disservice to completely ignore the

i. :-)