Rendering a Julia set zoom animation in TypeScript

Fractals are beautiful mathematical constructs: infinitely repeating, self-similar, awe-inspiring artefacts of beauty. Lately I have been yearning to render one myself, for no particular reason other than to enjoy the act of creating such a thing. My target goal for this post will thus be to illustrate how I went from an empty project to a full (albeit low-quality) animation:

The Julia set

For this project, I chose to render the Julia set. Or, more specifically, a Julia set, given there are infinitely many. If you snoop on Wikipedia, you will find a very rigorous and likely very obscure definition. Given that my inuitive understanding of complex algebra is relegated to a deep memory of my univeristy days, I went straight away for the more practical construction. Here is what I got from it, in a more digestible form.

Let’s say we have a function f(z) that can be simply expressed as a polynomial in the complex domain over the input variable z. We can similarly define these functions:

• f₂(z) = f(f(z))
• f₃(z) = f(f(f(z)))
• …

And so on. In general, f_n(z) is called the n-th iterate of f and is simply the result of repeatedly applying f to its own result n times. The iterations of f build a sequence. The Julia set of f is the set of all those complex numbers that make such sequence converge. In more mundane terms, no matter how many times you keep re-applying f onto itself, its result will always be bound by a finite constant.

There is a famous complex polynomial that is traditionally used to illustrate the Julia set - and its more distant cousing, the Mandelbrot set. Such polynomial is written as:

f_c(z) = z² + c

where c is an arbitrarily chosen constant. Since the behaviour of f_c depends on the value of the constant, this actually defines a broad family of quadratic polynomials, each of which has its own particular Julia set. In the following paragraphs, I will refer to the Julia set of f_c simply as “the Julia set”, implying that we are working on the specific quadric polynomial mentioned above, and that we have a fixed value of c.

Colouring points

Starting from the empirical construction given above, it follows that any point in the complex domain either belongs to the Julia set or not. So, how can its visualization contain all those sparkly colour gradients? The answer lies in the rendering implementation. Traditionally, computing whether a point belongs to the set is done through a loop. The number of iterations completed before exiting yields an estimation of how deep into the sequence we got before discovering it does not actually converge.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

class Complex {
    public constructor (public re: number, public im: number) { }
}

/**
 * Represents the JuliaSet for the given value of c passed to its constructor.
 */
class JuliaSet {
    private static ITERATION_LIMIT: number = 128;
    private static VALUE_THRESHOLD: number = 4;

    public constructor (private c: Complex) { }

    /**
     * Compute the convergence speed of the given point z with respect to the Julia set quadratic form.
     * @returns A real number between 0 and 1. Low numbers near 0 indicate that z causes the series
     *     to diverge really fast. High numbers near 1 indicate that z causes the series to either
     *     diverge really slowly or not at all (although it is not possible to know which clause
     *     is true in a finite number of steps).
     */
    public computeConvergence (z: Complex) : number {
        let i = 0;
        let c = this.c;

        let zrsqr = z.re * z.re;
        let zisqr = z.im * z.im;

        // This is an optimized algorithm to compute z^2 + c that operates directly on the
        // real an imaginary parts of z and c using the least amount of sums and multiplications.
        // It's useful to avoid the usage of Math.pow as it is a lot slower than this technique

        while (zrsqr + zisqr <= JuliaSet.VALUE_THRESHOLD) {
            z.im = z.re * z.im;
            z.im += z.im; // Multiply by two
            z.im += c.im;
            z.re = zrsqr - zisqr + c.re;
            zrsqr = z.re * z.re;
            zisqr = z.im * z.im;

            if (++i > JuliaSet.ITERATION_LIMIT)
                break;
        }

        return Math.min(i / JuliaSet.ITERATION_LIMIT, 1.0);
    }
}

export { Complex, JuliaSet };

In literature, if f_c is always less than 2 for a “good number” of iteration, we can assume that the current point belongs to the Julia set. Of course this is just an approximation that programs can use for faster rendering. Notice that in the code above I used 4 as my limit because the condition is applied to the squared magnitude: another trick to get rid of a very expensive square root operation that isn’t really necessary.

In the end, computeConvergence returns a number between 0 and and 1. In broad terms we can say that this value is a kind of “fuzzy boolean” representing the truth value of the assertion “z belongs to the Julia set”. By defining a palette of various colours ordered in a continuous gradient, it is thus very easy to assign a different colour to each point. You can see how I implemented my colour palette here.

Rendering a complete image

In order to finally render a snapshot of the Julia set onto an image, we just need a few more things:

• the value of the constant c, of course;
• the complex viewport, i.e. which area of the complex plane we are currently displaying;
• and the resulting image resolution for the output.

In my approach, I decomposed the second requirement into two distinct parts that are simpler to assign for the user:

• the complex coordinate of the “focal point”, i.e. the complex number that is going to be rendered exactly in the center of the image. This makes zooming much easier because we just keep the center constant and we shrink the viewport around it symmetrically;
• the level of zoom currently applied, which is simply a number between 1 and… well, infinity, potentially - but for the sake of arguments, let’s say between 1 and 10⁷. The limitation is mainly due to the fact that native floating point arithmetic isn’t suitable for such fine-grained computation. The precision barrier soon makes further zooming impossible without resorting to custom infinite-precision libraries.

In the following excerpt from my implementation, I use a node module called Jimp to easily create and draw images in memory.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

import * as Jimp from 'jimp';
import * as path from 'path';

import { Complex, JuliaSet } from './julia';
import { Palette } from './palette';

class Size {
    public constructor (public width: number, public height: number) {
        if ((width <= 0) || (height <= 0)) {
            throw new Error('Expected both with and height to be positive non-zero numbers');
        }
    }
}

/**
 * This is used to render a julia set onto an image with the given palette.
 */
class JuliaSetVisualizer {
    public constructor (private juliaSet: JuliaSet, private palette: Palette) { }

    private async createEmptyImage (width: number, height: number) : Promise<any> {
        return new Promise((resolve, reject) => {
            new Jimp(width, height, (err: any, image: any) => err ? reject(err) : resolve(image));
        });
    }

    /**
     * Generate a visualization of the julia set with the palette given in the constructor.
     * @param center Complex point that is rendered at the center of the image. This indirectly defines where 
     *     visualizer is looking within the fractal. Using different centers you can render different sections of it.
     * @param imageSize Resolution of the final image. Virtually unbounded: just remember it has to fit in your working memory.
     * @param zoom Level of zoom to use for rendering. Sensible values range from 1 to 10^7. Beyond the latter,
     *     default arithmetic implementation is not precise enough. If you want to go depper, you'll need to replace
     *     javascript numbers with an alternative implementation of arbitrary-precision floating points (which is generally very slow).
     * @param outputFileName Path or file name in which to save the output. Must contain a valid image extension (png or jpg).
     */
    public async generateImage (center: Complex, imageSize: Size, zoom: number, outputFileName: string) : Promise<any> {
        if (zoom < 1) {
            throw new Error('Zoom level must be at least 1');
        }

        if (!outputFileName.endsWith('.jpg') && !outputFileName.endsWith('.png')) {
            throw new Error('Invalid output file name extension. Allowed extensions are: .png, .jpg');
        }

        const image = await this.createEmptyImage(imageSize.width, imageSize.height);

        const DELTA = 1 / zoom;
        const ASPECT_RATIO = imageSize.width / imageSize.height;
        const R_MIN = center.re - DELTA * ASPECT_RATIO;
        const R_MAX = center.re + DELTA * ASPECT_RATIO;
        const I_MIN = center.im - DELTA;
        const I_MAX = center.im + DELTA;

        for (let x = 0; x < imageSize.width; ++x) {
            for (let y = 0; y < imageSize.height; ++y) {
                const r = R_MIN + (R_MAX - R_MIN) * (x / imageSize.width);
                const i = I_MIN + (I_MAX - I_MIN) * (y / imageSize.height);
                const convergence = this.juliaSet.computeConvergence(new Complex(r, i));
                const color = this.palette.interpolate(Math.sqrt(convergence));
                image.setPixelColor(Jimp.rgbaToInt(color.r, color.g, color.b, 255, undefined), x, y);
            }
        }

        return image.write(outputFileName);
    }
}

Notice that I use the square root of the converge value to interpolate a color: this is just an utility trick to obtain a more aestethically-pleasing result.

Finally, by keeping center and imageSize fixed, and by incrementing the zoom factor exponentially over time, it is possible to obtain a sequence of images that can be converted into a video. In my case, I simply used ffmpeg:

ffmpeg -f image2  -framerate 30 -i sequence/out%04d.jpg -vcodec mpeg4 -b 7M output.avi

Limitations and further improvements

The most obvious limitation is of course the JavaScript floating point airthmetic. It is possible to go around it by using other libraries for infinite precision arithemtic, but they have a huge performance impact on the program. A more refined approach would be to use the superfractal math paper, published some years ago now, to dramatically speed up computation of deep zooms; and to use native libraries for 128, 256 and 512 bit floating point operations. On top of that, note that this is a trivially parallelizable problem.

Additionally, my color scheme is very naive and gets stuck into a monochrome gradient after a short while. There are many solutions to the palette limitation, like using the smooth iteration count for better gradients. There is another article I read some time ago about how to create an infinite-recurring palette at any level of zoom using logarithmic level curves, but I can no longer find it 😢

Anyway, I hope you can play around with my code and enjoy some fractal exploration!