SDF level-set — marching squares

//! 2-D SDF level-set — marching squares.
//!
//! Walk a regular `w × h` grid of scalar values, emit one line segment per
//! cell that the iso-value `iso` crosses, using the 16-case marching-squares
//! lookup. Pure compute. V-class L0-closed.
//!
//! Output is a packed `Vec<f32>` of segment endpoints in pixel coordinates
//! relative to the grid: `[x0, y0, x1, y1, x0, y0, x1, y1, ...]`. The JS
//! side strokes them on a canvas.

use serde::{Deserialize, Serialize};

/// Per-frame report from the marching-squares pass.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct SdfKernelReport {
    pub grid_w: u32,
    pub grid_h: u32,
    pub segments: u32,
    pub bit_ops: u64,
}

/// Marching squares on a `w × h` scalar grid, iso-level `iso`. Appends
/// segment endpoints to `out` and returns the kernel report.
///
/// `cells` is row-major: `cells[y * w + x]`. `out` is cleared at entry.
pub fn marching_squares(cells: &[f32], w: u32, h: u32, iso: f32, out: &mut Vec<f32>) -> SdfKernelReport {
    out.clear();
    let w_us = w as usize;
    let h_us = h as usize;
    debug_assert_eq!(cells.len(), w_us * h_us);

    let mut segments = 0u32;

    for y in 0..(h_us - 1) {
        for x in 0..(w_us - 1) {
            let tl = cells[y * w_us + x];
            let tr = cells[y * w_us + (x + 1)];
            let bl = cells[(y + 1) * w_us + x];
            let br = cells[(y + 1) * w_us + (x + 1)];

            let mut idx = 0u8;
            if tl > iso { idx |= 1; }
            if tr > iso { idx |= 2; }
            if br > iso { idx |= 4; }
            if bl > iso { idx |= 8; }

            let lerp = |a: f32, b: f32| -> f32 {
                if (b - a).abs() < 1e-12 {
                    0.5
                } else {
                    ((iso - a) / (b - a)).clamp(0.0, 1.0)
                }
            };

            let fx = x as f32;
            let fy = y as f32;
            // Edge midpoints (lerped) at:
            //   top:    (fx + t_top, fy)
            //   right:  (fx + 1, fy + t_right)
            //   bottom: (fx + t_bot, fy + 1)
            //   left:   (fx, fy + t_left)
            let t_top = lerp(tl, tr);
            let t_right = lerp(tr, br);
            let t_bot = lerp(bl, br);
            let t_left = lerp(tl, bl);

            let p_top = (fx + t_top, fy);
            let p_right = (fx + 1.0, fy + t_right);
            let p_bot = (fx + t_bot, fy + 1.0);
            let p_left = (fx, fy + t_left);

            // Push helper.
            let mut push = |a: (f32, f32), b: (f32, f32)| {
                out.extend_from_slice(&[a.0, a.1, b.0, b.1]);
                segments += 1;
            };

            // 16-case lookup. Ambiguous cases (5, 10) resolved with
            // saddle-point disambiguation via average of corners.
            match idx {
                0 | 15 => {}
                1 => push(p_left, p_top),
                2 => push(p_top, p_right),
                3 => push(p_left, p_right),
                4 => push(p_right, p_bot),
                5 => {
                    // Ambiguous: tl and br above iso; tr and bl below.
                    let avg = 0.25 * (tl + tr + bl + br);
                    if avg > iso {
                        push(p_left, p_top);
                        push(p_right, p_bot);
                    } else {
                        push(p_left, p_bot);
                        push(p_top, p_right);
                    }
                }
                6 => push(p_top, p_bot),
                7 => push(p_left, p_bot),
                8 => push(p_left, p_bot),
                9 => push(p_top, p_bot),
                10 => {
                    let avg = 0.25 * (tl + tr + bl + br);
                    if avg > iso {
                        push(p_top, p_right);
                        push(p_left, p_bot);
                    } else {
                        push(p_left, p_top);
                        push(p_right, p_bot);
                    }
                }
                11 => push(p_right, p_bot),
                12 => push(p_left, p_right),
                13 => push(p_top, p_right),
                14 => push(p_left, p_top),
                _ => unreachable!(),
            }
        }
    }

    // Bit-ops estimate: per cell we do ~4 comparisons + up to 4 lerps
    // (~6 ops each) + a couple of writes. Round to ~40 f32 ops per cell,
    // 32 bit-ops per op.
    let cells_visited = ((w_us - 1) * (h_us - 1)) as u64;
    let bit_ops = cells_visited * 40 * 32;

    SdfKernelReport {
        grid_w: w,
        grid_h: h,
        segments,
        bit_ops,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Fill a grid with the SDF of a circle at center, radius r.
    /// Positive outside, negative inside (SDF convention).
    fn circle_sdf(w: u32, h: u32, cx: f32, cy: f32, r: f32) -> Vec<f32> {
        let mut v = Vec::with_capacity((w * h) as usize);
        for y in 0..h {
            for x in 0..w {
                let dx = x as f32 - cx;
                let dy = y as f32 - cy;
                v.push((dx * dx + dy * dy).sqrt() - r);
            }
        }
        v
    }

    #[test]
    fn empty_grid_under_iso_emits_no_segments() {
        let cells = vec![-1.0f32; 16 * 16];
        let mut out = Vec::new();
        let r = marching_squares(&cells, 16, 16, 0.0, &mut out);
        assert_eq!(r.segments, 0);
        assert!(out.is_empty());
    }

    #[test]
    fn circle_emits_a_closed_contour() {
        let cells = circle_sdf(64, 64, 32.0, 32.0, 16.0);
        let mut out = Vec::new();
        let r = marching_squares(&cells, 64, 64, 0.0, &mut out);
        // Circle of radius 16 → circumference ≈ 100. Expect ~80-120 segments.
        assert!(r.segments > 60 && r.segments < 200, "got {} segments", r.segments);
        assert_eq!(out.len() as u32, r.segments * 4);
    }

    #[test]
    fn higher_iso_means_smaller_circle() {
        let cells = circle_sdf(64, 64, 32.0, 32.0, 20.0);
        let mut out = Vec::new();
        let r0 = marching_squares(&cells, 64, 64, 0.0, &mut out);
        let mut out2 = Vec::new();
        let r1 = marching_squares(&cells, 64, 64, 5.0, &mut out2);
        // iso = 5 picks the contour at radius 25 — bigger, more segments.
        assert!(r1.segments > r0.segments, "{} vs {}", r1.segments, r0.segments);
    }

    #[test]
    fn bit_ops_scale_with_grid_area() {
        let cells_a = vec![0.0f32; 16 * 16];
        let cells_b = vec![0.0f32; 32 * 32];
        let mut out = Vec::new();
        let ra = marching_squares(&cells_a, 16, 16, 0.0, &mut out);
        out.clear();
        let rb = marching_squares(&cells_b, 32, 32, 0.0, &mut out);
        let ratio = rb.bit_ops as f64 / ra.bit_ops as f64;
        assert!(ratio > 3.5 && ratio < 5.0, "ratio was {ratio}");
    }
}

This is the exact Rust file compiled into the WASM module the page just loaded. Every line that runs on your device is here. The receipt is a function of this code, not a bespoke benchmark.