Home

Julia Sets

This is an interactive visualization of a set of points in the complex plane, where a complex number zz may be represented in (x,y)(x, y) coordinates by the decomposition z=x+iyz = x + iy Specifically, the Julia set(s) that we visualize on this page are the points that stay close to the origin when repeatedly mapped according to

zz2+c;[xy][x2y2+cx2xy+cy]\begin{equation} \label{eqn-1}\tag{1} z \rightarrow z^{2} + c \quad ; \quad \begin{bmatrix} x \\ y \end{bmatrix} \mapsto \begin{bmatrix} x^{2} - y^{2} + c_{x} \\ 2xy + c_{y} \end{bmatrix} \end{equation}

where cc is a point determined by the coordinates of your mouse. Each value of cc thus defines a unique set, which we visualize by mapping how long an initial point (a pixel in this image) stays close to the origin using color. Go ahead, move your mouse over the image!

View fullscreen

You can click on the image as well! Doing so animates a single iteration of the map Eq. (1). Note that, in the complex plane, multiplication looks like a combination of rotation and a scaling, while addition looks like a translation.

Notice that, because the points in the set must approximately stay in this set (lest they soon fly off to infinity), the Julia set describes fixed points and cycles of the mapping (Eq. (1)).

Also note that, if you could record all of the mouse positions for which the origin of the image is in the Julia set, you would map the Mandelbrot set! In fact, the Julia set(s) and the Mandelbrot set are different two-dimensional slices of the same four-dimensional object, where the Julia sets fix values of cc and the Mandelbrot set fixes x=0x = 0 . We can visualize the Mandelbrot set in the background of the set above. Pay attention to the origin, noting that the origin is in the Julia set if and only if your mouse is in the Mandelbrot set.

Show Mandelbrot set above

Fragment Shader:

Using twgl.js behind the scenes, here is the code to generate the Julia sets above:

precision mediump float;

// default uniforms (built in)
uniform vec2 resolution;
uniform float min_dim;
uniform bool mouseover;
uniform bool mousedown;
uniform vec2 mouse;
uniform vec2 mouse_click;
uniform float click_elapsed;

vec2 c = mouse;
float p = max(0., 1000. - click_elapsed) / 1000.;

const int iter = 130;
const float fiter = 130.;

/*
 * map interval [0, 1] to RGB color
 * use YCoCg color space for simplicity
 * https://en.wikipedia.org/wiki/YCoCg
 */
vec3 colormap(float v) {
    float Y = pow(1. - pow(v - 1., 2.), 2.); // logistic on v [0, 1] -> [0, 1]
    float cr = (1. - cos(6.283 * v)) / 2.;   // radius in "hue space"
    float Cg = cr * sin(4. * Y + 3.);        // in [-1, 1]
    float Co = cr * cos(4. * Y + 3.);        // in [-1, 1]
    float R = Y + Co - Cg;
    float G = Y + Cg;
    float B = Y - Co - Cg;

    return vec3(R, G, B);
}

/* multiply two complex numbers */
vec2 c_mul(vec2 x, vec2 y) {
    vec2 z;

    z.x = x.x * y.x - x.y * y.y;
    z.y = x.x * y.y + x.y * y.x;

    return z;
}

/* raise complex number to power e */
const int MAX_EXPONENT = 10;
vec2 c_pow(vec2 x, int e) {

    vec2 y = vec2(1.0, 0.0);

    for(int i=0; i < MAX_EXPONENT; i++) {

        if (i >= e) { break; }

        y = c_mul(y, x);

    }

    return y;

}

/* get magnitude-squared of 2-vector */
float mag2(vec2 z) {
    return (z.x * z.x + z.y * z.y);
}

float iterate(vec2 z) {

    /* Compute map on z and c up to max iteration */

    vec2 w;
    float v;
    float r2;

    // iterate map
    for(int i=0; i < iter; i++) {

        w = c_pow(z, 2) + c;

        v = float(i);
        r2 = mag2(w);

        // bailout radius
        if(r2 > 5.0) {
            break;
        }
        z.x = w.x;
        z.y = w.y;
    }

    // "real" iteration number using potential function
    if (v == (fiter - 1.)) {
        v = v / fiter;
    } else {
        v = (v - log2( log(r2) / log(5.))) / fiter;
    }

    return v;

}

float single_point(vec2 z) {

    /*
     * Compute value for single point, allowing
     * for parameterization of mapping with cut
     *
     * Call multiple times for anti-aliasing
     */

    float r2;
    float a;
    float theta;

    float theta2;
    vec2 zcut;
    float v = 0.;

    float rc2 = (z.x - c.x) * (z.x - c.x) + (z.y - c.y) * (z.y - c.y);

    // display c
    if (rc2 < 0.0001) {
        return 0.0;
    }

    // apply parameterized, partial inverse mapping
    // so we can visualize the map
    if (p > 0.) {

        z.x = z.x + (p - 1.) * c.x;
        z.y = z.y + (p - 1.) * c.y;

        a = 1.0 - 0.5 * (1. - p);
        r2 = z.x * z.x + z.y * z.y;

        theta = a * atan(z.y, z.x);
        z.y = pow(r2, a / 2.) * sin(theta);
        z.x = pow(r2, a / 2.) * cos(theta);

        theta2 = theta + sign(z.y) * 3.14159 * (1. - p);
        zcut.y = pow(r2, a / 2.) * sin(theta2);
        zcut.x = pow(r2, a / 2.) * cos(theta2);

        if (sign(zcut.y) != sign(z.y)) {
            v = iterate(zcut);
        }

    }

    return max(v, iterate(z));

}

void main() {

    float v = 0.;

    // 4-Rook Antialiasing
    v += single_point(
        vec2(
            ((2. * (gl_FragCoord.x + 0.125)) - resolution.x) / min_dim,
            ((2. * (gl_FragCoord.y + 0.375)) - resolution.y) / min_dim
        )
    );
    v += single_point(
        vec2(
            ((2. * (gl_FragCoord.x + 0.375)) - resolution.x) / min_dim,
            ((2. * (gl_FragCoord.y - 0.125)) - resolution.y) / min_dim
        )
    );
    v += single_point(
        vec2(
            ((2. * (gl_FragCoord.x - 0.125)) - resolution.x) / min_dim,
            ((2. * (gl_FragCoord.y - 0.375)) - resolution.y) / min_dim
        )
    );
    v += single_point(
        vec2(
            ((2. * (gl_FragCoord.x - 0.375)) - resolution.x) / min_dim,
            ((2. * (gl_FragCoord.y + 0.125)) - resolution.y) / min_dim
        )
    );

    vec4 value = vec4(colormap(v / 4.0), v / 4.0);

    gl_FragColor = value;
}