2017-10-13

The Mandelbrot Set

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.

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 \(\sqrt{-1}\), whatever that means.)  So you have a complex number like \(z\):

A complex number \(z\) plotted on the complex plane.

$$z=3+4i$$

In that case, the real part of \(z\) is \(3\), and the imaginary part is \(4i\).

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$$ $$z_{n+1} = z_{n}^2+c$$

where \(z_0\) is a complex zero: \(0+0i\).

What we do is calculate \(z_n\) 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 \(z_n\) and \(c\), stolen from the above Mandelbrot equation (the values of \(x\), \(y\), \(a\), and \(b\) are arbitrary for this example):

$$\begin{align} z_n&=x+yi \\ c&=a+bi \end{align}$$

And we want to compute the next number in the sequence \(z_{n+1}\) from the Mandelbrot equation:

$$z_{n+1} = z_n^2 + c$$

Then we just expand it like with regular algebra and go. Let's start with the \(z_n^2\) part:

$$\begin{align} z_n^2 &= (x + yi)^2 && \small{\text{from definition of }z_n\text{, above}} \\ &= x^2 + 2xyi + i^2y^2 && \small{\text{expand}} \\ &= x^2 + 2xyi - y^2 && \small{\text{since }i^2=-1\text{ by definition}} \\ &=(x^2 - y^2) - 2xyi && \small{\text{consolidate real and imaginary parts}} \end{align}$$

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

$$\begin{align} &=(x^2-y^2) + 2xyi + c && \\ &=(x^2-y^2) + 2xyi + a + bi && \small{\text{from the definition of }c\text{, above}} \\ &=(x^2-y^+a) + (2xy + b)i && \small{\text{consolidate real and imaginary parts}} \end{align}$$

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 \(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 + yi$$

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 \le 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!

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)      # 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.)

Things to note about the app at the top of the page:

• 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 renders until 30 milliseconds have elapsed and then cooperatively yields control to the JavaScript interpreter. Then it picks up again on the next frame where it left off. (This is extra-necessary, since JavaScript 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.

• Here is a link to the source code on GitHub.

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:

One more:

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,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:

# 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:

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!

References

• Mandelbrot Explorer JS source on GitHub

• A C-based Mandelbrot set renderer

• Wikipedia

Historic Comments

Comments