Kontakt
Kontakt

The Mandelbrot Set: A WebAssembly Example

by Christoph Dähne on 01.01.2022

During some quiet time at the end of the year I finally played around with WebAssembly. Writing an online Mandelbrot Set ("Apfelmännchen" in German) has been on my mind for quiet some time now. Those two things are the perfect match.

Since there are so many and so good resources about WebAssembly and the Mandelbrot Set around already (including the video), I skip the basics here. You can try it out and use the code as an example.

Disclaimer: This is by no means a full Mandelbrot set visualizer. There is room for image improvement (ie the use of colors) and the algorithm is not optimized at all. The code is meant to be a more complex "Hello World" WebAssembly program. After all: we must deal with complex numbers.

In this example code we run f(z) = z^2 + c a defined number of times for each c. If the end result is larger than the defined threshold, we consider f(z) = z^2 + c to be unbounded.

Note that the computation runs entirely on your machine and there is no protection against insanely high parameter-values. Memory consumption is constant, since you cannot change the image resolution. Execution time increases with the number of calculation iterations. The default values run a few seconds.

The image generator








(loading…)



Reading Suggestions

Those sources helped me the most during the implementation:

The Source Code

runMe.sh

#!/usr/bin/env bash
 
pushd $(dirname $0)
 
echo "compiling WebAssembly"
npx wat2wasm mandelbrot.wat -o mandelbrot.wasm
 
echo "starting server, open http://localhost:8080/mandelbrot.html"
php -S localhost:8080

mandelbrot.html

<!DOCTYPE html>
<html>
<head>
    <title>Hello Mandelbrot</title>
</head>
<body>
    <h1>Hello Mandelbrot</h1>
    <form id="my-form" action="#">
        <label for="x-axis-start">X-Axis Start:</label>
        <input type="text" id="x-axis-start" value="-2"><br>
        <label for="x-axis-end">X-Axis End:</label>
        <input type="text" id="x-axis-end" value="0.5"><br>
        <label for="y-axis-start">Y-Axis Start:</label>
        <input type="text" id="y-axis-start" value="-1.25"><br>
        <label for="y-axis-end">Y-Axis End:</label>
        <input type="text" id="y-axis-end" value="1.25"><br>
        <br>
        <label for="calculation-threshold">Calculation Threshold:</label>
        <input type="text" id="calculation-threshold" value="0.005"><br>
        <label for="calculation-iterations">Calculation Iterations:</label>
        <input type="text" id="calculation-iterations" value="75"><br>
 
        <p id="my-hint">(loading…)</p>
        <input id="my-submit" type="submit" value="Render Image">
    </form>
    <br><br>
    <canvas id="my-canvas" width="660" height="660" style="border: 1px solid black;"/>
    <script src="mandelbrot.js"></script>
</body>
</html>

mandelbrot.js

console.log('loading…');
 
console.log('canvas setup…');
const canvas = document.getElementById('my-canvas');
const width = Number(canvas.attributes.width.value);
const height = Number(canvas.attributes.height.value);
const canvasCtx = canvas.getContext('2d');
canvasCtx.clearRect(0, 0, width, height);
 
/**
 * shared memory setup
 * 
 * memory layout: color(0,0), color(0,1), …, color(x,y), …, color(x_max, y_max)
 * color(x,y) = rgba (4 bytes)
 */
console.log('shared memory setup…');
const pixelCount = width * height;
const bytesPerPixel = 4; // rgba
const pixelBytesCount = pixelCount * bytesPerPixel;
const pageSize = 65536; // 64kB
const pageCount = Math.ceil(pixelBytesCount / pageSize);
if (pageCount > 32) {
    throw "insufficient memory, working on 32 * 64kB"
}
const memory = new WebAssembly.Memory({initial: 32 /* must be in sync with WAT file */});
// shared memory as image data
const imageData = new ImageData(
    new Uint8ClampedArray(memory.buffer, 0, pixelBytesCount),
    width,
    height);
 
/**
 * form setup
 */
const xAxisStartField = document.getElementById('x-axis-start');
const xAxisEndField = document.getElementById('x-axis-end');
const yAxisStartField = document.getElementById('y-axis-start');
const yAxisEndField = document.getElementById('y-axis-end');
const thresholdField = document.getElementById('calculation-threshold');
const iterationsField = document.getElementById('calculation-iterations');
const hintField = document.getElementById('my-hint');
const myForm = document.getElementById('my-form');
 
function updateHint() {
    const overallIterations = width
      * height
      * iterationsField.value;
    hintField.innerHTML = `Required executions of f(z) = z^2 + c:<br>${overallIterations}`;
    return false;
}
updateHint();
myForm.onkeyup = updateHint
 
/**
 * load WebAssembly
 */
console.log('WASM setup…');
(async () => {
    const wasm = await WebAssembly.instantiateStreaming(
        fetch('mandelbrot.wasm'),
        {env: {memory}});
    const { render } = wasm.instance.exports;
 
    myForm.onsubmit = () => {
        // d - domain, x - pixels
        const x_min = 0;
        const x_max = width;
        const dx_min = xAxisStartField.value;
        const dx_max = xAxisEndField.value;
        const dx_scale = (dx_max - dx_min) / (x_max - x_min);
        const dx_offset = dx_min - dx_scale * x_min;
        // same for y-axis (but keep x/y ration)
        const y_min = height; // inverted x-axis !!
        const y_max = 0;
        const dy_min = yAxisStartField.value;
        const dy_max = yAxisEndField.value;
        const dy_scale = (dy_max - dy_min) / (y_max - y_min);
        const dy_offset = dy_min - dy_scale * y_min;
 
        const before = Date.now();
        render(
            0, // (param $x_start i32)
            width, // (param $x_count i32)
            dx_scale, // (param $dx_scale f32)
            dx_offset, // (param $dx_offset f32)
            0, // (param $y_start i32)
            height, // (param $y_count i32)
            dy_scale, // (param $dy_scale f32)
            dy_offset, // (param $dy_offset f32)
            thresholdField.value, // (local $d_threshold f32)
            iterationsField.value // (param $iterations i32)
        );
        const afterRender = Date.now();
        canvasCtx.putImageData(imageData, 0 ,0);
        const afterDraw = Date.now();
        console.log('done, timing in ms', {
            rendering: afterRender - before,
            drawing: afterDraw - afterRender,
            all: afterDraw - before
        });
        // stay on page
        return false;
    };
    console.log('ready to start rendering');
})();

mandelbrot.wat

(module
  (;
    The main and only exported functions is "render".
    It write the color values of pixels directly in the shared linear memory.
    This memory can be painted on a canvas afterwards (via JS).
 
    Note that we operate on complex numbers at some places: z = x + yi
    To avoid working on x and y separately everywhere, which would be messy
    especially when it comes to return values of functions like f(z) = z^2 + c
    we encode x and y into one variable.
    Since x and y are f32 (32bit floats) we concat the bits to an i64
    (64bit integer) and reverse the operations whenever we need to perform
    calculations.
    I am not sure about the performance penalty yet.
  ;)
 
  ;; load shared memory, number of pages must be in sync with JS file
  (import "env" "memory" (memory 32))
  ;; main rendering function
  ;; updates the pixels in the shared memory location
  ;; renders the "Apfelmännchen" (mandelbrot set)
  ;; by executing f(z) = z^2 + c (z, c are complex numbers)
  (func (export "render")
        (param $x_start i32) ;; minimal x-coordinate
        (param $x_count i32) ;; number of columns in the pixel matrix
        ;; we transform the x-coordinate to the domain value by
        ;; cx = dx_scale * x + dx_offset
        ;; since it is a linear transformation from the pixel axis to the domain axis
        (param $dx_scale f32) ;; 
        (param $dx_offset f32)
        ;; the y parameters are similar to the x ones
        (param $y_start i32)
        (param $y_count i32)
        (param $dy_scale f32)
        (param $dy_offset f32)
        ;; complex vectors (x,y) which length is greater than this threshold
        ;; are considered "infinite"
        (param $d_threshold f32)
        ;; number of repeated executions of f(z) = z^2 + c for each c
        (param $iterations i32)
 
        ;; increasing x in outer loop
        (local $x i32)
        ;; increasing y in inner loop
        (local $y i32)
        ;; real component x of c = x + yi
        (local $cx f32)
        ;; imaginary component y of c = x + yi
        (local $cy f32)
        ;; length of the current result (in innermost loop)
        (local $length f32)
        ;; resuling color of the curren pixel (in innermost loop)
        (local $color i32)
    ;; for (x = 0; x < x_count; x++)
    (local.set $x (i32.const 0))
    (block $for_x_break (loop $for_x_loop
      ;; break unless x < x_count
      (br_if $for_x_break
        (i32.ge_s (local.get $x) (local.get $x_count)))
 
      ;; transform pixel coordinate x to real component cx of c
      ;; cx = dx_scale * x + dx_offset
      (local.set $cx
        (f32.add
          (f32.mul
            (local.get $dx_scale)
            (f32.convert_i32_s (local.get $x)))
          (local.get $dx_offset)))
 
      ;; for (y = 0; y < y_count; y++)
      (local.set $y (i32.const 0))
      (block $for_y_break (loop $for_y_loop
        ;; break unless y < y_count
        (br_if $for_y_break
          (i32.ge_s (local.get $y) (local.get $y_count)))
 
        ;; transform pixel coordinate y to real component cy of c
        ;; cy = dy_scale * y + dy_offset
        (local.set $cy
          (f32.add
            (f32.mul
              (local.get $dy_scale)
              (f32.convert_i32_s (local.get $y)))
            (local.get $dy_offset)))
 
        ;; execute f(z) = z^2 + c multiple times
        ;; and meaure length of resulting vector
        (local.set $length
          (call $c_length
            (call $mandelbrot_fn
              (call $c_combine (local.get $cx) (local.get $cy))
              (local.get $iterations))))
 
        ;; determine color
        (if (f32.gt (local.get $length) (local.get $d_threshold))
          (then
            (local.set $color (i32.const 0xffffffff)))
          (else
            (local.set $color (i32.const 0xff000000)))
        )
        ;; set color
        (call $set_pixel_color
          (local.get $x)
          (local.get $x_count)
          (local.get $y)
          (local.get $color))
 
        ;; y++
        (local.set $y
          (i32.add (local.get $y) (i32.const 1)))
        ;; repeat
        br $for_y_loop
      ))
      ;; x++
      (local.set $x
        (i32.add (local.get $x) (i32.const 1)))
      ;; repeat
      br $for_x_loop
    ))
  )
  ;; helper function to set the color of one pixel
  (func $set_pixel_color
        ;; x coordinate to et
        (param $x i32)
        ;; number or colums in pixel matrix
        (param $width i32)
        ;; y coordinate
        (param $y i32)
        ;; RGBA color code in reverse order (due to endianess)
        (param $abgr i32)
    ;; index = 4 * (y * width + x)
    ;; stack: -
    local.get $y
    ;; stack: y
    local.get $width
    ;; stack: y, width
    i32.mul
    ;; stack: y * width
    local.get $x
    ;; stack: x, y * width
    i32.add
    ;; stack: x + y * width
    i32.const 4
    ;; stack: 4, x + y * width
    i32.mul
    ;; stack: 4 * (x + y * width)
    ;; memory[index] = abgr
    local.get $abgr
    ;; stack: index, 4 * (x + y * width)
    i32.store
    ;; stack: -
  )
  ;; executes mandelbrot_f many times: f(f(f(…f(z))))
  (func $mandelbrot_fn
        (param $c i64)
        (param $n i32)
        (result i64)
        (local $z i64)
    (local.set $z
      (call $c_combine (f32.const 0) (f32.const 0)))
    ;; for (n > 0; n--)
    (block $for_n_break (loop $for_n_loop
      ;; break unless n > 0
      (br_if $for_n_break
        (i32.le_s (local.get $n) (i32.const 0)))
      ;; call f(z)
      (local.set $z 
        (call $mandelbrot_f (local.get $z) (local.get $c)))
      ;; n--
      (local.set $n
        (i32.sub (local.get $n) (i32.const 1)))
      ;; repeat
      br $for_n_loop
    ))
    local.get $z
    return
  )
  ;; f(z) = z^2 + c (z and c are complex numbers)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $mandelbrot_f
        (param $z i64)
        (param $c i64)
        (result i64)
    ;; f(z, c) = z^2 + c
    (call $c_add
      (call $c_square (local.get $z))
      (local.get $c))
    return
  )
  ;; f(z) = z^2 (z is a complex number)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_square
        (param $z i64)
        (result i64)
        (local $x f32)
        (local $y f32)
    ;; z^2 = (x + yi)^2 = x^2 - y^2 + 2xyi
    (local.set $x (call $c_real_of (local.get $z)))
    (local.set $y (call $c_imaginary_of (local.get $z)))
    (call $c_combine
      (f32.sub
        (f32.mul
          (local.get $x)
          (local.get $x))
        (f32.mul
          (local.get $y)
          (local.get $y)))
      (f32.mul
        (f32.const 2)
        (f32.mul
          (local.get $x)
          (local.get $y))))
    return
  )
  ;; f(a, b) = a + b (a and b are complex numbers)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_add
        (param $a i64)
        (param $b i64)
        (result i64)
    ;; a + b = (x + yi) + (u + vi) = (x + u) + (y + vi)
    (call $c_combine
      (f32.add
        (call $c_real_of (local.get $a))
        (call $c_real_of (local.get $b)))
      (f32.add
        (call $c_imaginary_of (local.get $a))
        (call $c_imaginary_of (local.get $b))))
    return
  )
  ;; f(z) = |z| (z is a complex number)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_length
        (param $z i64)
        (result f32)
        (local $x f32)
        (local $y f32)
    ;; |z| = sqrt(x^2 + y^2)
    (local.set $x (call $c_real_of (local.get $z)))
    (local.set $y (call $c_imaginary_of (local.get $z)))
    (f32.sqrt
      (f32.add
        (f32.mul (local.get $x) (local.get $x))
        (f32.mul (local.get $y) (local.get $y))))
    return
  )
  ;; f(z) = f(x + yi) = x (z is a complex number)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_real_of
        (param $z i64)
        (result f32)
    (i64.shr_u (local.get $z) (i64.const 32))
    i32.wrap_i64
    f32.reinterpret_i32
    return
  )
  ;; f(z) = f(x + yi) = y (z is a complex number)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_imaginary_of
        (param $z i64)
        (result f32)
    local.get $z
    i32.wrap_i64
    f32.reinterpret_i32
    return
  )
  ;; f(x, y) = x + yi = z (z is a complex number)
  ;; z[64bit] = real[32bit] || imaginary[32bit]
  (func $c_combine
        (param $x f32)
        (param $y f32)
        (result i64)
    (i64.extend_i32_u (i32.reinterpret_f32 (local.get $x)))
    (i64.shl (i64.const 32))
    (i64.extend_i32_u (i32.reinterpret_f32 (local.get $y)))
    i64.or
    return
  )
)

Thanks for reading and feel free to ask questions or give feedback :-).