High resolution Buddhabrot generator

Generates high resolution pictures of the Buddhabrot using standard modern C++ with very few dependencies.

16384 x 16384 version: large 95 MiB file

Requirements

png++

Building

On UNIX-like systems,

g++ -Ofast -march=native -lpng -lpthread -o buddhabrot buddhabrot.cpp

clang++ -O3 -march=native -lpng -lpthread -o buddhabrot buddhabrot.cpp

Optional: CubeHelix colouring

g++ -Ofast -march=native -lpng -o cubehelix cubehelix.cpp

clang++ -O3 -march=native -lpng -lpthread -o buddhabrot buddhabrot.cpp

Running

./buddhabrot image_size iterations num_threads max_samples_per_pixel

image_size is how big your image is. The output will always be a square iamge representing the complex plane from -2 to 2 on each axis.
iterations is the max iterations.
num_threads is the number of threads to use. Use the physical cores, not logical threads. For example my AMD Ryzen 9 3900X performs better with num_threads = 12 although it is hyperthreaded and has 24 logical threads. Please also note that memory usage scales with number of threads, so if you are running out of RAM, you may wish to use fewer threads.
max_samples_per_pixel is the maximum number of random samples per pixel. If this value is too low, the output may be grainy.

The program will automatically output a 16-bit grayscale PNG image that is brightness-normalized and gamma-corrected. If it appears washed out, you can adjust the constrast in your preferred image editing program. The 16-bit depth is much more than conventional 8-bit images so you have lots of leeway to adjust the image.

The program cubehelix will apply the CubeHelix colour palette to output a cool-looking colourful image, with an option to adjust the contrast as needed.

Theory

The Buddhabrot is the probability distribution over trajectories that escape the Mandelbrot fractal.

The pseudocode is the following:

function buddhabrot(size, max_iterations):
    image = zeros(size, size)
    for some random pixels c:
        z = complex(0, 0)
        trajectory = []
        for iterations in range(max_iterations):
            z = z * z + c
            trajectory.append(z)
            if |z| > escape_radius:
                for point in trajectory:
                    image[point.x, point.y]++
                break
    return image

However, naively choosing the random pixels will be very slow.

We notice that the vast majority of random pixels c don't do much because they either escape very quickly, or they never escape.

Hence, we use a crude version of importance sampling.

For a given region, we can estimate how likely it to be interesting by starting out with very few random samples. If these early samples all turn out to be boring (e.g. if all of them fail to escape), then we can terminate the sampling immediately.

However, if they do escape, we can look at how fast they escape. If they escape really quickly, then it is kind of boring too, so we can terminate after a few samples. But if some of them take a long time to escape, then these are interesting trajectories, so we will increase the number of samples accordingly.

Another note is that the most interesting samples tend to happen on the edges of the Mandelbrot set, so if we find that, for a given region, some points escape and some don't, then this is a good indicator that it is on the edge and therefore supremely interesting indeed.

Finally, after sampling, the pixel values along the trajectory should be divided by the number of samples.

Although my code is not well optimized, on my computer it is able to generate a 1000 iteration 16384 x 16384 image with up to 128 samples per pixel within 4 minutes.

g++ -Ofast -march=native -lpng -lpthread -o buddhabrot buddhabrot.cpp; and time ./buddhabrot 16384 1000 12 128

________________________________________________________
Executed in  225.85 secs   fish           external
   usr time   34.89 mins    0.00 micros   34.89 mins
   sys time    0.09 mins  418.00 micros    0.09 mins