Fractals are among the most visually stunning phenomena in mathematics, and rendering them on the GPU with GLSL is a perfect way to explore both the beauty of complex dynamics and the power of parallel computation. In this post, we will build a complete, animated fractal shader from scratch — a Julia set that smoothly morphs over time with vibrant, psychedelic coloring — and then break down every piece of the math and technique behind it.
At the heart of fractal rendering lies complex number iteration. The Mandelbrot set and Julia sets both derive from the same recurrence relation: z = z² + c. The difference is simple: for the Mandelbrot set, c varies across the plane and z₀ = 0; for a Julia set, c is fixed and z₀ varies across the plane.
Below is a fully self-contained WebGL1 fragment shader. It renders an animated Julia set with smooth iteration counting, orbit trap coloring, and a vibrant multi-hue palette.
precision mediump float;
uniform vec2 iResolution;
uniform float iTime;
const int MAX_ITER = 256;
const float BAILOUT = 256.0;
vec3 hsv2rgb(float h, float s, float v) {
vec3 c = vec3(h, s, v);
vec3 rgb = clamp(
abs(mod(c.x * 6.0 + vec3(0.0, 4.0, 2.0), 6.0) - 3.0) - 1.0,
0.0, 1.0
);
return c.z * mix(vec3(1.0), rgb, c.y);
}
vec3 fractalPalette(float t) {
float r = 0.5 + 0.5 * cos(3.0 + t * 4.2 + 0.0);
float g = 0.5 + 0.5 * cos(3.0 + t * 4.2 + 2.1);
float b = 0.5 + 0.5 * cos(3.0 + t * 4.2 + 4.2);
return vec3(r, g, b);
}
void main() {
vec2 uv = (gl_FragCoord.xy - 0.5 * iResolution.xy) / iResolution.y;
float zoom = 1.5 + 0.3 * sin(iTime * 0.15);
uv *= zoom;
vec2 z = uv;
vec2 c = vec2(
0.355 + 0.15 * cos(iTime * 0.23) + 0.05 * sin(iTime * 0.73),
-0.02 + 0.15 * sin(iTime * 0.23) + 0.05 * cos(iTime * 0.51)
);
float iter = 0.0;
float min_dist = 1e10;
for (int i = 0; i < 256; i++) {
z = vec2(
z.x * z.x - z.y * z.y,
2.0 * z.x * z.y
) + c;
float d = dot(z, z);
min_dist = min(min_dist, d);
if (d > BAILOUT) break;
iter += 1.0;
}
float smooth_iter = iter;
if (iter < 256.0) {
float log_zn = log(dot(z, z)) * 0.5;
float nu = log(log_zn / log(2.0)) / log(2.0);
smooth_iter = iter + 1.0 - nu;
}
vec3 color;
if (iter >= 256.0) {
float trap = sqrt(min_dist);
color = hsv2rgb(0.65 + 0.1 * trap, 0.8, trap * 0.4);
} else {
float t = smooth_iter / 256.0;
vec3 primary = fractalPalette(t * 5.0 + iTime * 0.08);
float trap_influence = 1.0 - exp(-min_dist * 2.0);
vec3 secondary = hsv2rgb(
fract(t * 3.0 + iTime * 0.05),
0.7 + 0.3 * trap_influence,
0.9
);
color = mix(primary, secondary, 0.3 + 0.2 * sin(t * 20.0));
float glow = 1.0 - t;
color *= 0.7 + 0.6 * pow(glow, 0.4);
}
vec2 vign_uv = gl_FragCoord.xy / iResolution.xy;
float vignette = 1.0 - 0.3 * length(vign_uv - 0.5);
color *= vignette;
color = pow(color, vec3(0.9));
gl_FragColor = vec4(color, 1.0);
}
Every pixel on screen corresponds to a point on the complex plane. We iterate f(z) = z² + c and ask: does the sequence remain bounded, or does it escape to infinity? The number of iterations before escape determines the color.
// Complex squaring in GLSL using vec2
vec2 z_squared = vec2(
z.x * z.x - z.y * z.y, // real: a² - b²
2.0 * z.x * z.y // imaginary: 2ab
);
The Mandelbrot set and Julia sets are intimately connected. Each point c in the Mandelbrot set corresponds to a unique Julia set. If c lies inside the Mandelbrot set, the Julia set is connected (one piece). If c lies outside, the Julia set is a Cantor dust (totally disconnected). The most visually interesting Julia sets come from c values near the boundary of the Mandelbrot set.
// Animate c through parameter space near the Mandelbrot boundary
vec2 c = vec2(
0.355 + 0.15 * cos(iTime * 0.23) + 0.05 * sin(iTime * 0.73),
-0.02 + 0.15 * sin(iTime * 0.23) + 0.05 * cos(iTime * 0.51)
);
Naive fractal renderers assign color based on the raw integer iteration count, producing ugly discrete color bands. The normalized (smooth) iteration count uses a double logarithm to get a continuous value that smoothly interpolates between integer iteration counts.
float log_zn = log(dot(z, z)) * 0.5; // = log(|z|) float nu = log(log_zn / log(2.0)) / log(2.0); float smooth_iter = iter + 1.0 - nu;
Our shader combines cosine palettes (phase-shifted cosine waves for each RGB channel), orbit trap coloring (tracking minimum distance the orbit comes to the origin), and HSV color space for the fractal interior. Adding time offsets creates slow, mesmerizing hue rotation.
Try implementing a Burning Ship fractal by taking absolute values before squaring. Experiment with higher-power polynomials like z³ + c or z⁴ + c for different symmetries. Add deep zoom by exponentially decreasing the view scale over time. Use distance estimation to render crisp fractal boundaries without oversampling, or use the fractal as a heightmap for a 3D raymarched terrain.