This article describes how to render a ridicuously large 500 Megapixel Buddhabrot with iteration depths in the range of millions. Lisp source included. Scroll down for the final image.
Why would you want to do this? To have it printed and on your wall, of course (see pictures below).
I will begin by outlining the general idea of the Buddhabrot (see the Wikipedia article for more), then talk about my optimizations and rendering choices.
The Buddhabrot
You need a maximum iteration depth n (sometimes called the Dwell Limit) because many samples will be within the Mandelbrot Set and you would run the iteration indefinitely for those. Whenever the maximum iteration depth is reached, the sample is considered to be inside the Mandelbrot set - although it possibly isn't and after another 1000 Iterations it would have escaped. There is no way to be sure.
We now consider only those samples "good" whose iteration escaped before reaching the maximum iteration depth, i.e. those who are definitely outside the Mandelbrot Set. Each good sample gives a set of points whose magnitude is equal to the number of iteration steps i that it took for us to find out that it was outside the Mandelbrot Set. Sometimes this is a very small number like 2 or 3, but samples nearer to the border of the Mandelbrot Set take much longer. A sample can take thousands or millions of iterations before escaping, thus leaving a much larger footprint in the resulting point cloud. It seems that for each Natural Number n there is a sample that takes at least n iteration steps, although for higher values of n it becomes harder and harder to find one.
Deep Iteration
It has been observed by many people that just rendering the Buddhabrot like this - using all the escaping samples' visited points - gives a nebula-like image but that a more interesting image can be made by considering only those samples as "interesting" whose iteration time i exceeds some minimum iteration time k with k < n. For example, you could use k = 1000 and n = 10000, i.e. keep only those samples that do not escape before their 1000th iteration but do escape before their 10000th iteration. If you look at one of those more deeply iterated good samples' visited path in the complex plane you will note some pattern: it forms some kind of orbit, sometimes some kind of circle or a set of circles or other shapes. I have no idea why the orbits are shaped this way and i don't know if anyone else does. Here is a part of the orbit of the sample -0.9114822433916717 + 0.2523940788714053i, for example (containing nearly five million points):
Incidentally, all interesting samples like this one are located very very close to the border of the Mandelbrot Set and make good targets to zoom into with your favourite fractal explorer. Try it.
Since the longer orbits look nice we will try to find only "interesting" samples that take a long time to finally escape. There have been several approches to find them. The most obvious one is of course just choosing a sample randomly and seeing how it goes, which means throwing away something like 99.999% of them. You could also use a regularly spaced grid of samples which basically means the same. The interesting samples are all located near the Mandelbrot Set's border, so there must be some spatial pattern to their distribution which might help us find them. If you ever played around with the Set a bit you probably know that the border is infinitely long and calling it "complex" is a bit of an understatement, but surely we can do better than just randomly sample the plane?
Optimization
- Rejecting all samples from within the two biggest parts of the Mandelbrot Set, the "Main Cardioid" and the circular "Period-2-Bulb", is a must-have. This simple check removes a large amount of all the points in the Set which are very costly because they will iterate until the dwell limit is reached. I just used the formula from the Wikipedia article. The areas removed by this test are shown in red on this image of the Mandelbrot Set (my output images are rotated clockwise by 90 degrees, don't get confused):
- People have also used the Metropolis-Hastings Algorithm, considering the mapping c -> i a propability density function. The method is especially useful for rendering a small detail of the bigger image, but does not work too well for the whole image with very big values of k. I kept finding many interesting samples with it, but their orbits turned out to be very similar because several good samples that are very close to each other produce orbits that are almost identical - giving an image of a single blurred orbit. I increased the mutation distance enough to produce sufficiently distinctive orbits, but then the speedup was gone. In the end i didn't use this method. Probably it is more useful for rendering the "normal" Buddhabrot or a detail of it without the minimum iteration constraint.
- To partly remove some of the costly larger areas (like the non-circular bulbs labeled "3" and up in the image above) from within the Set i divided the plane into rectangular cells and chose random samples only from those cells that contained some part of the border of the Mandelbrot Set. Depending on the cell grid resolution and the iteration depth this gave a speedup of 3 to 8.
Since we still choose our samples randomly, every rendering will look different - depending on which orbits your random search finds. There is no "the" Buddhabrot (see Melinda's comment below - rendering for an infinite amount of time yields "the" Buddhabrot). This also means that the image is not symmetic on the small scale which is a good thing if you ask me. You could of course get a symmetric image of the same quality in half the time by mirroring it, but i think it's worth waiting for the asymmetric one. On the right is a part of the final rendering that exhibits superficial symmetry, but is composed of asymmetric orbits. Click it to view this detail in the final resolution.
If i had had my renderer run for much longer (maybe a week or a month), more and more orbits would have shown up, occluding each other and finally blurring the whole thing into the standard nebula-like buddhabrot. The key is to render just long enough until a sufficient amount of orbits have appeared and then stop.
Lisp Implementation
The renderer is written in Common Lisp which is not your classic number crunching language. If the goal was to do it as fast as possible, you'd probably end up doing it in C with inline assembly or on the GPU. On the other hand, i had never really optimized any lisp code to the limit and wanted to do just that: see how far you can get. One long night of low-level optimization at entropia with the help of some friends later there was is only a factor of 2 left compared to the C implementation that one entropian made. Which is quite okay, considering the cumulative speedup of about 2000 compared to my first implementation. Trading of sbcl's builtin complex numbers for traditional pairs of floats, declaring all types, working in-place instead of on the stack, considering cache efficiency, prompting usage of inline arithmetic instead of jumps to costly sbcl functions, etc. gained us a factor of about 3 to 5 while the big rest was due to the algorithmic improvements.
The renderer uses double-floats throughout. I think i can safely say that the orbits' shapes are not due to cumulative floating-point errors (which i had suspected at some time) because tiny variations in the starting value usually do not produce significantly different orbits. Orbits are not chaotic systems.
Since tracing the few good samples is quite cheap compared to finding them, the renderer first multithreadedly collects some good samples by randomly sampling from the grid cells (see optimizations above) and only afterwards processes them to keep track of their jumping around on the complex plane. Their path is recorded in a big 16bit greyscale image which is repeatedly written to disk.
Postprocessing, Color and Noise
Converting this greyscale HDR image to a LDR image with color has to be done by hand which is a nontrivial task for most image manipulation programs because 20000x25000 is not exactly your standard image resolution. Gimp and Cinepaint both crash or freeze (and Gimp can't handle 16-bit images), but Photoshop worked for me. I also cheated a bit by gently surface-blurring the image, removing some noise.
The noise is the only reason you need the really large iteration depths in the range of millions that i used. If you increase both the minimum iteration time and the maximum iteration limit, your orbits will consist of more and more points, allowing you to choose a higher resolution. If, on the other hand, you get less than, say, five hits per pixel in the darker regions of the orbit you will start to notice that some pixels may get only two hits and their neighbour pixel may get ten hits which will be visible as noise in the output image. If you iterate more deeply, this effect evens out and the signal-to-noise ratio increases. To remove noise, you can always downscale the output image or, as i did, selectively blur the image a bit to improve the look of the noisy areas.
For the output image that i consider the most successful one i used a minimum iteration time of 1000000 and a maximum iteration depth of 5000000. it took about 16 hours on my university's 8-core Xeon machine which i abused for this task overnight.
Results
Here is a scaled version of about 5000x6000 pixels for your firefox's non-crashing convenience.
I am quite satisfied with the digitally printed one as well, although some of the darker reds lose their saturation.
Do not open the large image in your browser. It is a 93MiB JPG file of 20000x25000 pixels which consumes 1.5 GiB of RAM when opened. Don't try this at home. I warned you.
You can also have a look at the Lisp source.
Johann, I love what you've done. My only small issue is with your statement "There is no 'the' Buddhabrot." I understand what you mean and I agree that it makes sense given the aesthetic choices that you made but there are other choices in which this concept does make sense. If you let your renderer run long enough, all those lovely and distinct trajectories will blend together such that further computation will not change the image at all after histogram equalization. Given this requirement, images can be stable, converged, and otherwise "fully cooked". Doing this for your image may take 100 times as many processing cycles and be therefore out of reach with your current resources, but that doesn't mean that such canonical b-brot images do not exist.
Hi Melinda
You are of course right: i'm stopping very early in the approximation process towards something that probably is "the" Buddhabrot. Forgive the simplification - i wanted to make it very clear that the structures in my images are not set in mathematical stone, as e.g. the Mandelbrot Set's are.
Anyway, is it known whether the deeply iterated Buddhabrot converges towards the same image as the normal rendering without the minimum iteration constraint?
Hello Johann,
There is very little in this field that is mathematically known (I.E. proved) but two points are pretty clear. 1) Images generated as I described eventually converge, and 2) Choice of maximum iteration threshold can have a profound effect on the image converged upon.
You are asking about the effects of an additional *minimum* iteration threshold, and my experience suggests that adding such a constraint does not greatly affect deeply iterated b-brot images because the contributions from fast escapers tend to get drowned out by the slow escapers. There's a bit of a tug-of-war between the fast escapers which are more common but individually make less of a contribution, and the slow escapers which are rare but make large individual contributions. With low maximum thresholds the fast escapers dominate, whereas with high maximums the slow escapers win even though they are relatively rare. implemented an option to specify minimum thresholds but I quickly stopped using them for this reason.
This also shows why early experimenters missed discovering the b-brot because without filtering out the non-escapers, the resulting images are dominated by them resulting in anti-b-brot images which do contain the b-brot but completely drowned out. It's like trying to take a photo of the sun's corona which is essentially impossible without blocking out the blinding light from the disk of the sun.
Gimp crashing? I have several 15000x15000 images Gimp managed a year or so ago, and a poster with a Julia set in my office which (if I remember correctly, as I am not at home now) was 30kx30k and worked okay (but slooow).
Good work you have done! If you find c's such that the n-th iteration is eventually mapped to zero, you get points that are inside the M set but accumulate eventually at the boundary (check a poster I presented, you can find it at www.maia.ub.es/~ruben/ but the result was folks knowledge). Maybe you can use this to speed-up?
Ruben
I now utilise a Xeon 8-core for geologic tasks, geographic information systems projects, massive multi-tasking(doc+PDF+GoogleEarth+streaming+data-access+photoshop+photofiltre+firewall+antivirus+presentation+spreadsheet+more...). The system is amazing.
That said, I remember trying to compute Mandelbrot sets back in the 1990s, and early 2000s with a Pentium 200Mhz, a Xeon single-core 2.4 Ghz, and then somehow forgot all about it, until having recently assembled a Xeon 8-core, 10GB-RAM system with a 1333Mhz FSSB. Now I am trying to understand what the most recent multi-threading fractal software for non-programmers is...I've taken basic (Java, C++) programming at the entry level, but cannot afford the time required to write new code just to make fractals. Is there a compiled program to work with the complexity described on this page? I can easily afford 16 hours on an 8-core 2.4 Ghz Xeon machine with 10GB RAM, about twice a week....but I don't know how to generate this kind of imagery. I would gladly run fractal models with my machine, make them available to whomever wanted the results, and be satisfied with a neat, 50000x50000 JPEG at the end of it all....anyone interested in working together?
I Love it. And The Large file is nö Problem for Photoshop. Which Lears me to my question. is it possible to get Otter Backgrounds. or The file As 32 Bit png without Background. I'd like to print it for my Home but dark ist not so fine.
Best regards
Harry
This is a beautiful rendering. I love the coloring and the really high iteration counts give it a very unique look compared to the typical cloudy Buddhabrot. Bravo!
-William Milberry
How much costed this printed banner?
Hello. And Bye.