• Wikipedia on the Buddhabrot
• Wikipedia on the Metropolis - Hastings algorithm
• Buddhabrot discoverer Melinda Green's site.
• Another fast zooming technique.

- The Mandelbrot Set -

The Mandelbrot set is very probably the most recognizable fractal. It is visualized by picking a complex number and iterating it. Each successive iteration is a point in an "orbit". Some orbits escape to infinity. That is, the total area covered by it will keep getting larger. Other orbits never escape and just hang around forever. These two sets of orbits form the interior and exterior of the set.

Pseudocode might be:

int Mandelbrot(complex c)
{
	complex z = c
	int iterations = 0
	while( z.real2 + z.imaginary2 < 4.f AND iterations < max_iters)
	{
		z = z2 + c
		iterations++
	}

	return iterations
}

- The Buddhabrot -

About two years ago, I found a page describing something called the Buddhabrot, which is just the Mandelbrot set viewed in a novel way. The Buddhabrot computes orbits as normal, but instead of colouring each pixel based upon how long it took for the corresponding complex number to escape to infinity, it saves the orbit and increments every pixel that the orbit passes through.

The technique produces some interesting images but has a fundamental shortcoming that limits any real "exploration". As you start zooming in, the probability that any given orbit passes through our area of interest falls rapidly. Unless you're willing to wait a very, very long time you're limited to viewing the fractal as a whole.

It doesn't have to be that way though. You just need to be careful about which complex numbers you iterate. Similar complex numbers tend to produce similar orbits so what's needed is some means to identify important complex numbers and then explore the complex numbers around them.

- The Metropolis-Hastings Algorithm -

The Metropolis-Hastings algorithm was developed in 1953 and can effectively sample traditionally difficult probability distributions. All that's required is the ability to sample the distribution at x. The algorithm works by proposing small changes to x, called mutations, and either accepting or rejecting them based upon a transition probability.

Pretend we have some function F that has a very difficult to sample probability distribution. We start with some value of x and propose small changes to it. Two tentative samples at F(x) and F(x') are drawn, and the "better" sample is chosen.

Pseudocode might look something like this:

It's a bit more complicated than that, but you get the idea. The transition probability is, oddly enough, the probability that one sample will get chosen over another.

- Making It Work -

But fractals don't have probability distributions, right? Not per se, however the area we're trying to render does. Pretend that we're interested in part of the fractal from (-0.5, 0.5) to (-0.4, 0.4). Every orbit that passes through or starts inside the square has a probability of one and everything else has a probability of zero.

What about the transition probability? It's not possible to compute just how likely one orbit is to be chosen over another. We can instead construct a transition probability around how "interesting" each orbit is. Orbits that take longer to escape are usually more interesting than those that escape quickly. Without a good transition probability, the images generated would have little or no contrast because orbits that affect many pixels would be given the same computation time as orbits that affect only a single pixel.

- A Mutation Strategy -

A good mutation strategy is critical to producing clean images within a reasonable timeframe. I've found that two very simple mutation types work quite well. The first is to offset the current complex number by a small amount. The second is to choose an entirely new complex number. Small mutations can be perturbed based upon an exponential distribution. The minimum and maximum size of the distribution should be proportional to the current level of magnification.

It's important that the values of r1 and r2 not be too large or too small. If they're too large, the ratio of accepted mutations will go down, if they're too small, they won't efficiently sample the entire space we're interested in. The ideal values are different for every area of the fractal but I've found that one tenth and one ten thousandth of the magnification works well. If one were feeling saucy, the values of r1 and r2 could themselves be chosen using the Metropolis-Hastings algorithm.

The second mutation type needs no explanation.

I've found that giving a probability of 4/5 to mutation one and 1/5 to mutation two works fairly well.

And that's pretty much it.

- Pretty Pictures -

For comparison, here is an image centered at (-0.0443594, -0.986749) with a magnification of 88.2 rendered in two minutes with the naive method.

Here is the same area in two minutes rendered with the Metropolis-Hastings method.

- Source Code -

The source for a very simply program can be found here. All of the images on this page were created with it. It uses SDL to create the window and should compile on any platform that supports it with minimal changes.

- Addendum -

The algorithm as described on this page will not translate into three dimensions - the so called "Buddhagram" - without some small modifications.

Extending the described mutation types into three dimensions is trivial, but you'll quickly find that the Metropolis-Hastings algorithm performs worse than a naive implementation when attempting to view the fractal from certain regions.

The contribution function is too blame because it only looks at how many times an orbit deposits pixels on the screen. Without knowing what's in the foreground and what's in the background, the algorithm will constantly assign too high importance to distant details.

The solution is to add the z value of each point in the orbit into the contribution function:

There is also a 3d specific mutation type that will increase the performance. The mutation simply picks a new point in front of the camera with an exponential distribution in z. This will help the algorithm quickly locate structures that are very close to the camera. Unlike the other two mutations, this one shouldn't be used that often. Maybe once every hundred mutations or so. Otherwise, distant structures will be poorly sampled.

A view of a minibrot on the main stem with the main body of the fractal in the background. Rendered with a fisheye lens.