Beej's Bit Bucket

 ⚡ Tech and Programming Fun

The Mandelbrot Set

2010-02-26
One of the "copies" of the Mandlebrot Set floating somewhere around its outside edge.

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 like x:

x = 12 + i3

A complex number z plotted on the complex plane.

In that case, the real part of x is 12, and the imaginary part is i3.

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:

z0 = 0  (complex zero)

zn+1 = zn2 + c

What we do is calculate zn up to a fairly large n (like 1000, for example) using a loop. If z stays bounded below a certain magic range, then congratulations, c is in the Mandlebrot Set. Otherwise c is not. I know "below a magic range" really isn't giving you much to go on, but I will elaborate momentarily.

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 zn and c, stolen from the above Mandelbrot equation (the values of x, y, a, and b are arbitrary for this example):

zn = x + iy

c = a + ib

And we want to compute the next number in the sequence zn+1:

zn+1 = zn2 + c (Mandelbrot equation)

Then we just expand it like with regular algebra and go. Let's start with the zn2 part:

zn2 = (x + iy)2 [from definition of zn, above]

= x2 + i2xy + i2y2 [expand]

= x2 + i2xy - y2 [since i2 = -1 by definition]

= (x2 - y2) + i2xy [consolidate real and imaginary parts]

That's a complex number,  see? It has real part (x2-y2) and imaginary part i2xy. We can finish up here by adding on c, which we do by just adding in the real and imaginary components:

= (x2 - y2) + i2xy + c

= (x2 - y2) + i2xy + a + ib [from the definition of c, above]

= (x2 - y2 + a) + i(2xy + b) [consolidate real and imaginary parts]

The real and imaginary parts of this final equation are clearly visible in the inner while loop 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 z1000. If z1000 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 z1000 for a starting value of c and it comes up with a complex number with real part x and imaginary part y:

z1000 = x + iy

And remember that the equation for a circle is:

x2 + y2 = r2

Points are colored differently depending on how quickly the corresponding complex number breaks past the radius-2 barrier. The black pixels in the center are in the Mandlebrot Set.

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

x2 + y2 ≤ 22

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!

Pseudocode
for 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.)

Either scripts and active content are not permitted to run or Adobe Flash Player version 10.0.0 or greater is not installed.

Get Adobe Flash Player

Things to note:

  1. 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.
  2. 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.)
  3. 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]
  4. 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.
  5. 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...
  6. 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 z1000, 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:

  1. Take the number of the next pixel to process, call it n.
  2. Use the Adam7 mapping to convert n to an x,y screen coordinate.
  3. Convert the x,y screen coordinate to a complex plane coordinate's real and imaginary parts, represented by a complex number c.
  4. Run the Mandelbrot iteration process on c, producing subsequent values of the complex number z. Keep track of the iteration number in t.
  5. 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.
  6. Paint the pixel at x,y the color q.
  7. GOTO step 1

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

Share me!

Historic Comments

 luckytoilet 2010-02-27 16:11:10

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?

 beej 2010-02-27 19:24:51

@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.

 Kevin 2010-03-11 22:21:47

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

 michael hilton 2010-05-10 09:43:09

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 ?

 beej 2010-05-14 17:04:40

@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. :-)

Comments

blog comments powered by Disqus
Blog  ⚡  Email beej@beej.us  ⚡  Home page