alice-sdf 5.1.0

//! SDF to mesh conversion using Marching Cubes (Deep Fried Edition)
//!
//! # Deep Fried Optimizations
//! - **Z-Slab Parallelization**: Process each Z-layer independently with Rayon.
//! - **Merge Sub-meshes**: Avoid mutex contention by merging results after parallel phase.
//! - **Forced Inlining**: `#[inline(always)]` on hot-path helpers.
//!
//! Author: Moroya Sakamoto

use crate::compiled::{eval_compiled, eval_compiled_batch_simd_parallel, CompiledSdf};
use crate::eval::{eval, eval_material, normal};
use crate::mesh::Vertex;
use crate::SdfNode;
use glam::{Vec2, Vec3, Vec4};
use rayon::prelude::*;

/// Configuration for marching cubes
#[derive(Debug, Clone, Copy, serde::Serialize, serde::Deserialize)]
pub struct MarchingCubesConfig {
    /// Grid resolution along each axis
    pub resolution: usize,
    /// Iso-level (usually 0 for SDF surface)
    pub iso_level: f32,
    /// Whether to compute vertex normals
    pub compute_normals: bool,
    /// Whether to generate triplanar UV coordinates
    pub compute_uvs: bool,
    /// UV tiling scale (world units per UV tile)
    pub uv_scale: f32,
    /// Whether to compute tangent vectors (requires normals and UVs)
    pub compute_tangents: bool,
    /// Whether to evaluate material IDs per vertex
    pub compute_materials: bool,
}

impl Default for MarchingCubesConfig {
    fn default() -> Self {
        Self {
            resolution: 64,
            iso_level: 0.0,
            compute_normals: true,
            compute_uvs: false,
            uv_scale: 1.0,
            compute_tangents: false,
            compute_materials: false,
        }
    }
}

/// Configuration preset with all AAA features enabled
impl MarchingCubesConfig {
    /// AAA preset: normals + triplanar UV + tangents + materials
    pub const fn aaa(resolution: usize) -> Self {
        Self {
            resolution,
            iso_level: 0.0,
            compute_normals: true,
            compute_uvs: true,
            uv_scale: 1.0,
            compute_tangents: true,
            compute_materials: true,
        }
    }
}

/// Compute triplanar UV coordinates from position and normal (Deep Fried)
///
/// Projects the vertex position onto the dominant axis plane determined
/// by the surface normal, producing seamless UVs for arbitrary SDF surfaces.
#[inline(always)]
fn triplanar_uv(position: Vec3, normal: Vec3, inv_scale: f32) -> Vec2 {
    let abs_n = normal.abs();

    if abs_n.x >= abs_n.y && abs_n.x >= abs_n.z {
        // X-dominant: project onto YZ plane
        Vec2::new(position.y * inv_scale, position.z * inv_scale)
    } else if abs_n.y >= abs_n.z {
        // Y-dominant: project onto XZ plane
        Vec2::new(position.x * inv_scale, position.z * inv_scale)
    } else {
        // Z-dominant: project onto XY plane
        Vec2::new(position.x * inv_scale, position.y * inv_scale)
    }
}

/// Compute tangent vector from normal (Deep Fried)
///
/// Generates a tangent perpendicular to the normal, aligned with the
/// triplanar UV projection direction. W component stores handedness.
#[inline(always)]
fn compute_tangent(normal: Vec3) -> Vec4 {
    let abs_n = normal.abs();

    let tangent = if abs_n.x >= abs_n.y && abs_n.x >= abs_n.z {
        // X-dominant: tangent along Y
        let t = Vec3::new(0.0, 1.0, 0.0);
        (t - normal * normal.dot(t)).normalize()
    } else if abs_n.y >= abs_n.z {
        // Y-dominant: tangent along X
        let t = Vec3::new(1.0, 0.0, 0.0);
        (t - normal * normal.dot(t)).normalize()
    } else {
        // Z-dominant: tangent along X
        let t = Vec3::new(1.0, 0.0, 0.0);
        (t - normal * normal.dot(t)).normalize()
    };

    Vec4::new(tangent.x, tangent.y, tangent.z, 1.0)
}

/// Compute MikkTSpace-compatible tangent vectors from UV gradients
///
/// For each triangle, computes the tangent and bitangent from position/UV deltas,
/// accumulates at shared vertices, then orthogonalizes against the normal and
/// computes handedness. Produces results compatible with MikkTSpace-based normal
/// map baking tools (xNormal, Substance, UE5, Unity, Godot).
fn compute_mikktspace_tangents(mesh: &mut Mesh) {
    let vert_count = mesh.vertices.len();
    if vert_count == 0 {
        return;
    }

    let mut tangents = vec![Vec3::ZERO; vert_count];
    let mut bitangents = vec![Vec3::ZERO; vert_count];

    // Phase 1: Accumulate tangent/bitangent per triangle
    for tri in mesh.indices.chunks_exact(3) {
        let i0 = tri[0] as usize;
        let i1 = tri[1] as usize;
        let i2 = tri[2] as usize;

        if i0 >= vert_count || i1 >= vert_count || i2 >= vert_count {
            continue;
        }

        let p0 = mesh.vertices[i0].position;
        let p1 = mesh.vertices[i1].position;
        let p2 = mesh.vertices[i2].position;

        let uv0 = mesh.vertices[i0].uv;
        let uv1 = mesh.vertices[i1].uv;
        let uv2 = mesh.vertices[i2].uv;

        let dp1 = p1 - p0;
        let dp2 = p2 - p0;
        let duv1 = uv1 - uv0;
        let duv2 = uv2 - uv0;

        let det = duv1.x * duv2.y + (-(duv1.y * duv2.x));
        if det.abs() < 1e-8 {
            continue; // Degenerate UV mapping
        }

        let inv_det = 1.0 / det;
        let t = Vec3::new(
            (dp1.x * duv2.y + (-(dp2.x * duv1.y))) * inv_det,
            (dp1.y * duv2.y + (-(dp2.y * duv1.y))) * inv_det,
            (dp1.z * duv2.y + (-(dp2.z * duv1.y))) * inv_det,
        );
        let b = Vec3::new(
            (dp2.x * duv1.x + (-(dp1.x * duv2.x))) * inv_det,
            (dp2.y * duv1.x + (-(dp1.y * duv2.x))) * inv_det,
            (dp2.z * duv1.x + (-(dp1.z * duv2.x))) * inv_det,
        );

        // Area-weighted accumulation (implicit via dp magnitude)
        tangents[i0] += t;
        tangents[i1] += t;
        tangents[i2] += t;
        bitangents[i0] += b;
        bitangents[i1] += b;
        bitangents[i2] += b;
    }

    // Phase 2: Orthogonalize and compute handedness per vertex
    for i in 0..vert_count {
        let n = mesh.vertices[i].normal;
        let t = tangents[i];
        let b = bitangents[i];

        if t.length_squared() < 1e-10 {
            // Fallback for vertices with no valid UV gradient
            mesh.vertices[i].tangent = compute_tangent(n);
            continue;
        }

        // Gram-Schmidt orthogonalize: T' = normalize(T - N * dot(N, T))
        let ortho_t = (t - n * n.dot(t)).normalize();

        // Handedness: sign of dot(cross(N, T'), B)
        let w = if n.cross(ortho_t).dot(b) < 0.0 {
            -1.0
        } else {
            1.0
        };

        mesh.vertices[i].tangent = Vec4::new(ortho_t.x, ortho_t.y, ortho_t.z, w);
    }
}

/// Simple mesh structure
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct Mesh {
    /// Mesh vertices
    pub vertices: Vec<Vertex>,
    /// Triangle indices
    pub indices: Vec<u32>,
}

impl Mesh {
    /// Create an empty mesh
    pub const fn new() -> Self {
        Self {
            vertices: Vec::new(),
            indices: Vec::new(),
        }
    }

    /// Get the number of vertices
    pub fn vertex_count(&self) -> usize {
        self.vertices.len()
    }

    /// Get the number of triangles
    pub fn triangle_count(&self) -> usize {
        self.indices.len() / 3
    }
}

impl Default for Mesh {
    fn default() -> Self {
        Self::new()
    }
}

/// Convert an SDF to a mesh using marching cubes
///
/// Produces a closed mesh whose vertices are shared across cell boundaries:
/// see [`marching_cubes`]. When UVs are enabled, tangent vectors are computed
/// using the MikkTSpace-compatible UV gradient method for correct normal map
/// rendering in all major engines.
///
/// # Arguments
/// * `node` - The SDF tree
/// * `min` - Minimum corner of bounding box
/// * `max` - Maximum corner of bounding box
/// * `config` - Marching cubes configuration
///
/// # Returns
/// Generated mesh
pub fn sdf_to_mesh(node: &SdfNode, min: Vec3, max: Vec3, config: &MarchingCubesConfig) -> Mesh {
    let mut mesh = marching_cubes(node, min, max, config);

    if config.compute_tangents && config.compute_normals && config.compute_uvs {
        compute_mikktspace_tangents(&mut mesh);
    }

    mesh
}

/// Marching cubes on a `config.resolution`³ grid over `[min, max]`.
///
/// Every vertex belongs to one lattice edge (its lower lattice point and an
/// axis) and is computed once, from that edge's lower endpoint; the cells
/// around the edge refer to it by index. The mesh is therefore closed and
/// 2-manifold without any welding by position, and adjacent cells cannot
/// disagree on a vertex. A lattice point whose value equals `iso_level` is
/// outside (only `value < iso_level` is inside), on the edges and in the
/// cell classification alike.
///
/// Such a lattice point is reached by the crossing on every edge from it to
/// an inside neighbour, and each of those edges keeps its own vertex: the
/// mesh can hold **several vertices with the same position** (one per edge,
/// separate indices) and **zero-area triangles** between them. The index
/// buffer stays a closed 2-manifold; merging those vertices by position
/// would collapse their fan and can pinch the surface. A consumer that
/// identifies vertices by position (STL, a positional weld) sees the zero-area
/// triangles as degenerate faces. For a sphere of integer radius on an
/// integer lattice, `tests/test_mc_shared_vertex_oracle.rs` pins both counts
/// in closed form.
///
/// The vertex normal is the SDF gradient at the
/// vertex position (`compute_normals`), or the area-weighted mean of the
/// adjacent face normals otherwise.
///
/// Until 5.0 every triangle corner was a separate vertex and callers welded
/// by position, which left open or pinched edges wherever the copies
/// differed in the last bits or a lattice point sat on the surface.
pub fn marching_cubes(node: &SdfNode, min: Vec3, max: Vec3, config: &MarchingCubesConfig) -> Mesh {
    let resolution = config.resolution;
    let cell_size = (max - min) / resolution as f32;

    // Evaluate SDF on grid. The grid is the hot loop of mesh generation
    // ((res + 1)³ evaluations, 2.1 M at res 128): compile once and run the
    // 8-lane SIMD batch evaluator (tree ≡ compiled ≡ SIMD by the parity
    // corpus); trees the compiler rejects (Terrain, Triangle, Bezier) fall
    // back to the tree walker. Review SDF-R2-4: the tree path ran at
    // 6.7 M evals/s where the SIMD batch does ~4× that.
    let grid_size = resolution + 1;
    let total_points = grid_size * grid_size * grid_size;
    let grid_point = |i: usize| {
        let x = i % grid_size;
        let y = (i / grid_size) % grid_size;
        let z = i / (grid_size * grid_size);
        lattice_point(min, cell_size, [x, y, z])
    };
    let values: Vec<f32> = match CompiledSdf::try_compile(node) {
        Ok(compiled) => {
            let points: Vec<Vec3> = (0..total_points).into_par_iter().map(grid_point).collect();
            eval_compiled_batch_simd_parallel(&compiled, &points)
        }
        Err(_) => (0..total_points)
            .into_par_iter()
            .map(|i| eval(node, grid_point(i)))
            .collect(),
    };

    let (positions, indices) =
        extract_shared_vertices(&values, resolution, min, cell_size, config.iso_level);
    let material = |p: Vec3| eval_material(node, p);
    build_mesh(
        positions,
        indices,
        config,
        |p| normal(node, p, 0.001),
        config.compute_materials.then_some(&material),
    )
}

/// Lattice point `g` of the grid with corner `min` and spacing `cell`.
///
/// The one formula for a lattice position on the CPU path; the GPU path
/// uploads the per-axis coordinates computed with it, so both evaluate the
/// field at bit-identical points.
#[inline(always)]
pub fn lattice_point(min: Vec3, cell: Vec3, g: [usize; 3]) -> Vec3 {
    min + Vec3::new(
        g[0] as f32 * cell.x,
        g[1] as f32 * cell.y,
        g[2] as f32 * cell.z,
    )
}

/// The inside test of every marching cubes path: strictly below the
/// iso-level. A value equal to the iso-level is outside, so a lattice point
/// on the surface belongs to exactly one side and the edge and cell
/// classifications cannot disagree.
#[inline(always)]
fn is_inside(value: f32, iso_level: f32) -> bool {
    value < iso_level
}

/// For each of the 12 cube edges (table numbering): the offset of its lower
/// corner and its axis (0 = x, 1 = y, 2 = z), derived from
/// [`CORNER_OFFSETS`] and [`EDGE_CONNECTIONS`].
const EDGE_LOWER_AXIS: [([usize; 3], usize); 12] = edge_lower_axis();

const fn edge_lower_axis() -> [([usize; 3], usize); 12] {
    let mut out = [([0usize; 3], 0usize); 12];
    let mut e = 0;
    while e < 12 {
        let a = CORNER_OFFSETS[EDGE_CONNECTIONS[e][0]];
        let b = CORNER_OFFSETS[EDGE_CONNECTIONS[e][1]];
        let mut lo = [0usize; 3];
        let mut axis = 0usize;
        let mut k = 0;
        while k < 3 {
            lo[k] = if a[k] < b[k] { a[k] } else { b[k] };
            if a[k] != b[k] {
                axis = k;
            }
            k += 1;
        }
        out[e] = (lo, axis);
        e += 1;
    }
    out
}

/// Shared-vertex marching cubes on a grid of `(res + 1)³` values (x fastest,
/// then y, then z).
///
/// Pass 1 visits every lattice edge once, z-layer by z-layer in parallel,
/// and emits a vertex for each sign change, interpolated from the edge's
/// lower endpoint. The vertex order (lattice point, then axis) makes the
/// edge ids ascending, so pass 2 — the cells, z-slab by z-slab — finds a
/// cell edge's vertex by binary search. Both passes concatenate their slabs
/// in z order, so the output is the same for any thread count.
///
/// Returns `(positions, triangle indices)`.
fn extract_shared_vertices(
    values: &[f32],
    res: usize,
    min: Vec3,
    cell: Vec3,
    iso_level: f32,
) -> (Vec<Vec3>, Vec<u32>) {
    let g = res + 1;
    let stride = [1, g, g * g];
    debug_assert_eq!(values.len(), g * g * g);

    let layers: Vec<(Vec<u64>, Vec<Vec3>)> = (0..g)
        .into_par_iter()
        .map(|z| {
            let mut ids = Vec::new();
            let mut pos = Vec::new();
            for y in 0..g {
                for x in 0..g {
                    let c = [x, y, z];
                    let i = x + y * g + z * g * g;
                    let va = values[i];
                    let inside_a = is_inside(va, iso_level);
                    for axis in 0..3 {
                        if c[axis] == res {
                            continue;
                        }
                        let vb = values[i + stride[axis]];
                        if inside_a == is_inside(vb, iso_level) {
                            continue;
                        }
                        let mut cb = c;
                        cb[axis] += 1;
                        let (p, _, _) = interpolate_edge(
                            lattice_point(min, cell, c),
                            lattice_point(min, cell, cb),
                            va,
                            vb,
                            iso_level,
                        );
                        ids.push(i as u64 * 3 + axis as u64);
                        pos.push(p);
                    }
                }
            }
            (ids, pos)
        })
        .collect();
    let total: usize = layers.iter().map(|l| l.0.len()).sum();
    let mut ids = Vec::with_capacity(total);
    let mut positions = Vec::with_capacity(total);
    for (i, p) in layers {
        ids.extend(i);
        positions.extend(p);
    }

    let slabs: Vec<Vec<u32>> = (0..res)
        .into_par_iter()
        .map(|z| {
            let mut out = Vec::new();
            for y in 0..res {
                for x in 0..res {
                    let mut cube_index = 0usize;
                    for (k, off) in CORNER_OFFSETS.iter().enumerate() {
                        let i = (x + off[0]) + (y + off[1]) * g + (z + off[2]) * g * g;
                        if is_inside(values[i], iso_level) {
                            cube_index |= 1 << k;
                        }
                    }
                    if EDGE_TABLE[cube_index] == 0 {
                        continue;
                    }
                    let tri = &TRI_TABLE[cube_index];
                    let mut t = 0;
                    while tri[t] != -1 {
                        let (lo, axis) = EDGE_LOWER_AXIS[tri[t] as usize];
                        let i = (x + lo[0]) + (y + lo[1]) * g + (z + lo[2]) * g * g;
                        let id = i as u64 * 3 + axis as u64;
                        let v = ids.binary_search(&id).unwrap_or_else(|_| {
                            unreachable!("cell edge {id} has no vertex: edge and cell classifications disagree")
                        });
                        out.push(v as u32);
                        t += 1;
                    }
                }
            }
            out
        })
        .collect();
    let indices = slabs.concat();
    (positions, indices)
}

/// Attach the per-vertex attributes `config` asks for to a shared-vertex
/// marching cubes surface. Each attribute is a function of the vertex
/// position alone (or, without `compute_normals`, of the adjacent faces), so
/// a vertex has one value whichever cell refers to it.
fn build_mesh(
    positions: Vec<Vec3>,
    indices: Vec<u32>,
    config: &MarchingCubesConfig,
    gradient_normal: impl Fn(Vec3) -> Vec3 + Sync,
    material: Option<&(dyn Fn(Vec3) -> u32 + Sync)>,
) -> Mesh {
    let normals: Vec<Vec3> = if config.compute_normals {
        positions.par_iter().map(|&p| gradient_normal(p)).collect()
    } else {
        let mut acc = vec![Vec3::ZERO; positions.len()];
        for t in indices.chunks_exact(3) {
            let [a, b, c] = [t[0], t[1], t[2]].map(|i| i as usize);
            let n = (positions[b] - positions[a]).cross(positions[c] - positions[a]);
            acc[a] += n;
            acc[b] += n;
            acc[c] += n;
        }
        acc.into_iter()
            .map(|n| n.try_normalize().unwrap_or(Vec3::Y))
            .collect()
    };
    let inv_uv = 1.0 / config.uv_scale;
    let vertices: Vec<Vertex> = positions
        .par_iter()
        .zip(normals.par_iter())
        .map(|(&p, &n)| {
            let mut v = Vertex::new(p, n);
            if config.compute_uvs {
                v.uv = triplanar_uv(p, n, inv_uv);
            }
            // With UVs, tangents come from the MikkTSpace pass in sdf_to_mesh.
            if config.compute_tangents && config.compute_normals && !config.compute_uvs {
                v.tangent = compute_tangent(n);
            }
            if let Some(m) = material {
                v.material_id = m(p);
            }
            v
        })
        .collect();
    Mesh { vertices, indices }
}

/// Merge multiple sub-meshes into a single mesh (Deep Fried)
#[inline(always)]
fn merge_meshes(sub_meshes: Vec<Mesh>) -> Mesh {
    // Pre-calculate total capacity
    let total_vertices: usize = sub_meshes.iter().map(|m| m.vertices.len()).sum();
    let total_indices: usize = sub_meshes.iter().map(|m| m.indices.len()).sum();

    let mut merged = Mesh {
        vertices: Vec::with_capacity(total_vertices),
        indices: Vec::with_capacity(total_indices),
    };

    for sub in sub_meshes {
        let base_idx = merged.vertices.len() as u32;
        merged.vertices.extend(sub.vertices);
        merged
            .indices
            .extend(sub.indices.iter().map(|i| i + base_idx));
    }

    merged
}

/// Linear interpolation along an edge (Deep Fried)
///
/// Branchless NaN-safe: copysign ensures |denom| >= epsilon without branching.
#[inline(always)]
fn interpolate_vertex(p0: Vec3, p1: Vec3, v0: f32, v1: f32, iso_level: f32) -> Vec3 {
    interpolate_edge(p0, p1, v0, v1, iso_level).0
}

/// [`interpolate_vertex`] that also returns the edge parameter `t` and
/// whether the endpoints were swapped into canonical order, so a caller can
/// blend per-endpoint attributes (normals) with the same `t` the position
/// used: `(position, t, swapped)`, position `= q0 + (q1 - q0) * t` where
/// `(q0, q1)` is `(p0, p1)` or, when `swapped`, `(p1, p0)`.
#[inline(always)]
pub fn interpolate_edge(p0: Vec3, p1: Vec3, v0: f32, v1: f32, iso_level: f32) -> (Vec3, f32, bool) {
    // A grid edge is shared by up to four cells, each of which reaches it
    // through a different local edge (and endpoint order). Evaluating from a
    // canonical endpoint — the lexicographically smaller corner — makes the
    // result bit-identical in every cell. Until 1.10.3 the endpoint order
    // was the local one, the last bits differed between neighbours and
    // `sdf_to_mesh` left 16–64 open edges at res ≥ 64 (external review).
    let swapped = (p1.x, p1.y, p1.z) < (p0.x, p0.y, p0.z);
    let (p0, p1, v0, v1) = if swapped {
        (p1, p0, v1, v0)
    } else {
        (p0, p1, v0, v1)
    };
    let denom = v1 - v0;
    let safe_denom = f32::copysign(denom.abs().max(1e-10), denom);
    let t = ((iso_level - v0) / safe_denom).clamp(0.0, 1.0);
    (p0 + (p1 - p0) * t, t, swapped)
}

// ---------------------------------------------------------------------------
// Compiled VM variants — use eval_compiled instead of interpreted eval
// ---------------------------------------------------------------------------

/// Compute surface normal using compiled evaluator (6 eval calls per vertex)
#[inline(always)]
fn normal_compiled(sdf: &CompiledSdf, point: Vec3, epsilon: f32) -> Vec3 {
    let ex = Vec3::new(epsilon, 0.0, 0.0);
    let ey = Vec3::new(0.0, epsilon, 0.0);
    let ez = Vec3::new(0.0, 0.0, epsilon);

    let grad = Vec3::new(
        eval_compiled(sdf, point + ex) - eval_compiled(sdf, point - ex),
        eval_compiled(sdf, point + ey) - eval_compiled(sdf, point - ey),
        eval_compiled(sdf, point + ez) - eval_compiled(sdf, point - ez),
    );

    let len_sq = grad.length_squared();
    if len_sq < 1e-20 {
        return Vec3::Y;
    }
    grad / len_sq.sqrt()
}

/// Convert compiled SDF to mesh (compiled VM path)
///
/// Uses `eval_compiled` for grid evaluation and normal computation,
/// providing 2-5x speedup over the interpreted `sdf_to_mesh`.
pub fn sdf_to_mesh_compiled(
    sdf: &CompiledSdf,
    min: Vec3,
    max: Vec3,
    config: &MarchingCubesConfig,
) -> Mesh {
    let mut mesh = marching_cubes_compiled(sdf, min, max, config);

    if config.compute_tangents && config.compute_normals && config.compute_uvs {
        compute_mikktspace_tangents(&mut mesh);
    }

    mesh
}

/// Marching cubes using the compiled VM evaluator.
///
/// The grid is evaluated with the SIMD batch path; the surface is extracted
/// as in [`marching_cubes`] (shared vertices, one per lattice edge). The
/// vertex normal is the central-difference gradient of the compiled field
/// at the vertex position. Until 5.0 it was the grid gradient at the local
/// edge's first endpoint, which differs between the cells sharing an edge
/// (measured: 10,240 of 15,366 sphere vertices at res 64 got two normals, up
/// to 103° apart at box corners), and the welded mesh stayed open.
/// Material IDs need the tree ([`marching_cubes`]); this path leaves them 0.
pub fn marching_cubes_compiled(
    sdf: &CompiledSdf,
    min: Vec3,
    max: Vec3,
    config: &MarchingCubesConfig,
) -> Mesh {
    let resolution = config.resolution;
    let cell_size = (max - min) / resolution as f32;

    let grid_size = resolution + 1;
    let total_points = grid_size * grid_size * grid_size;

    // Build all grid points, then SIMD batch evaluate
    let points: Vec<Vec3> = (0..total_points)
        .map(|i| {
            let x = i % grid_size;
            let y = (i / grid_size) % grid_size;
            let z = i / (grid_size * grid_size);
            lattice_point(min, cell_size, [x, y, z])
        })
        .collect();

    let values = eval_compiled_batch_simd_parallel(sdf, &points);

    let (positions, indices) =
        extract_shared_vertices(&values, resolution, min, cell_size, config.iso_level);
    build_mesh(
        positions,
        indices,
        config,
        |p| normal_compiled(sdf, p, 0.001),
        None,
    )
}

// Corner offsets for the 8 corners of a cube, in the numbering of the
// Bourke / Lorensen tables below: corners 0–3 walk the y = 0 face
// (x, then z), corners 4–7 the y = 1 face. Until 1.10.3 the second and
// third axes were swapped (0–3 on the z = 0 face), a mirror image of the
// table's cube, so every triangle came out wound inward (negative signed
// volume, STL facets contradicting their normals).
pub const CORNER_OFFSETS: [[usize; 3]; 8] = [
    [0, 0, 0],
    [1, 0, 0],
    [1, 0, 1],
    [0, 0, 1],
    [0, 1, 0],
    [1, 1, 0],
    [1, 1, 1],
    [0, 1, 1],
];

// Edge connections (which corners each edge connects)
pub const EDGE_CONNECTIONS: [[usize; 2]; 12] = [
    [0, 1],
    [1, 2],
    [2, 3],
    [3, 0],
    [4, 5],
    [5, 6],
    [6, 7],
    [7, 4],
    [0, 4],
    [1, 5],
    [2, 6],
    [3, 7],
];

// Edge table - which edges are intersected for each cube configuration
pub const EDGE_TABLE: [u16; 256] = [
    0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03,
    0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f,
    0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0xaa, 0x7a6,
    0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569,
    0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69,
    0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6,
    0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, 0xe5c,
    0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf,
    0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9, 0xac3,
    0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a,
    0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
    0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65,
    0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3, 0xfaa,
    0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30,
    0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33,
    0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f,
    0x596, 0x29a, 0x393, 0x99, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0,
];

// Triangle table - which edges form triangles for each cube configuration
pub const TRI_TABLE: [[i8; 16]; 256] = [
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
    [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
    [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
    [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
    [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
    [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
    [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
    [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
    [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
    [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
    [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
    [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
    [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
    [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
    [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
    [5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1],
    [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1],
    [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
    [0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
    [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1],
    [10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1],
    [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1],
    [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1],
    [5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1],
    [9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1],
    [0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1],
    [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1],
    [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1],
    [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1],
    [2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1],
    [7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1],
    [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1],
    [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1],
    [11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1],
    [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1],
    [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1],
    [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1],
    [11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
    [1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1],
    [9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1],
    [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1],
    [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
    [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1],
    [6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1],
    [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1],
    [6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1],
    [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1],
    [1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
    [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1],
    [6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1],
    [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1],
    [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
    [3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
    [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1],
    [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1],
    [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1],
    [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1],
    [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1],
    [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
    [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1],
    [10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1],
    [10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1],
    [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1],
    [1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1],
    [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1],
    [0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1],
    [10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1],
    [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1],
    [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1],
    [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1],
    [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
    [3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1],
    [6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1],
    [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1],
    [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1],
    [10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1],
    [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
    [7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1],
    [7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1],
    [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1],
    [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1],
    [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1],
    [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
    [0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1],
    [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
    [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
    [2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
    [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1],
    [7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1],
    [2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1],
    [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1],
    [10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1],
    [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1],
    [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1],
    [7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1],
    [6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1],
    [8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1],
    [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1],
    [6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1],
    [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1],
    [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1],
    [8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1],
    [0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1],
    [1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1],
    [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1],
    [10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1],
    [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
    [10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
    [5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
    [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1],
    [9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
    [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1],
    [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1],
    [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
    [7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1],
    [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1],
    [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1],
    [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
    [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1],
    [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
    [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1],
    [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1],
    [6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1],
    [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1],
    [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1],
    [6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1],
    [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1],
    [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1],
    [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1],
    [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1],
    [9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1],
    [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
    [1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
    [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1],
    [0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1],
    [5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1],
    [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1],
    [11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1],
    [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1],
    [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
    [2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1],
    [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1],
    [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1],
    [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
    [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1],
    [9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1],
    [9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1],
    [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1],
    [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1],
    [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1],
    [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1],
    [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1],
    [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
    [9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1],
    [5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1],
    [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
    [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1],
    [8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1],
    [0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1],
    [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1],
    [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1],
    [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1],
    [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1],
    [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1],
    [11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1],
    [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1],
    [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1],
    [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
    [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1],
    [1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1],
    [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1],
    [4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1],
    [3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1],
    [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1],
    [0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1],
    [9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1],
    [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
];

/// Configuration for adaptive marching cubes
#[derive(Debug, Clone)]
pub struct AdaptiveConfig {
    /// Maximum octree depth (resolution = 2^max_depth per axis)
    pub max_depth: u32,
    /// Minimum octree depth (prevents overly coarse cells)
    pub min_depth: u32,
    /// Distance threshold: cells further than this from surface stay coarse
    pub surface_threshold: f32,
    /// Iso-level
    pub iso_level: f32,
    /// Compute vertex normals
    pub compute_normals: bool,
    /// Compute triplanar UVs
    pub compute_uvs: bool,
    /// UV scale
    pub uv_scale: f32,
    /// Compute tangents
    pub compute_tangents: bool,
    /// Compute materials
    pub compute_materials: bool,
}

impl Default for AdaptiveConfig {
    fn default() -> Self {
        Self {
            max_depth: 6, // 64 effective resolution
            min_depth: 2, // 4 minimum
            surface_threshold: 0.0,
            iso_level: 0.0,
            compute_normals: true,
            compute_uvs: false,
            uv_scale: 1.0,
            compute_tangents: false,
            compute_materials: false,
        }
    }
}

impl AdaptiveConfig {
    /// AAA preset: high resolution near surface, efficient far away
    pub const fn aaa(max_depth: u32) -> Self {
        Self {
            max_depth,
            min_depth: 2,
            surface_threshold: 0.0,
            iso_level: 0.0,
            compute_normals: true,
            compute_uvs: true,
            uv_scale: 1.0,
            compute_tangents: true,
            compute_materials: true,
        }
    }
}

/// Octree node for adaptive subdivision
struct OctreeCell {
    min: Vec3,
    max: Vec3,
    depth: u32,
}

/// Convert SDF to mesh using adaptive octree marching cubes
///
/// Subdivides space using an octree: cells near the surface are refined
/// to max_depth resolution while cells far from the surface stay coarse.
/// Typically produces 60-80% fewer triangles than uniform grid for the
/// same effective resolution near the surface.
pub fn adaptive_marching_cubes(
    node: &SdfNode,
    bounds_min: Vec3,
    bounds_max: Vec3,
    config: &AdaptiveConfig,
) -> Mesh {
    // Build leaf cells via octree subdivision
    let root = OctreeCell {
        min: bounds_min,
        max: bounds_max,
        depth: 0,
    };

    let mut leaf_cells = Vec::new();
    subdivide_octree(node, &root, config, &mut leaf_cells);

    // Process leaf cells in parallel
    let mc_config = MarchingCubesConfig {
        resolution: 1, // Each leaf is a single cell
        iso_level: config.iso_level,
        compute_normals: config.compute_normals,
        compute_uvs: config.compute_uvs,
        uv_scale: config.uv_scale,
        compute_tangents: false, // Computed in post-processing
        compute_materials: config.compute_materials,
    };

    let sub_meshes: Vec<Mesh> = leaf_cells
        .par_iter()
        .map(|cell| process_adaptive_cell(node, cell, &mc_config))
        .collect();

    let mut mesh = merge_meshes(sub_meshes);

    // Deduplicate for watertight output
    super::optimize::deduplicate_vertices(&mut mesh);

    // Post-process tangents
    if config.compute_tangents && config.compute_normals && config.compute_uvs {
        compute_mikktspace_tangents(&mut mesh);
    }

    mesh
}

/// Recursively subdivide octree based on SDF distance
#[allow(clippy::suspicious_operation_groupings)]
fn subdivide_octree(
    node: &SdfNode,
    cell: &OctreeCell,
    config: &AdaptiveConfig,
    leaves: &mut Vec<OctreeCell>,
) {
    let center = (cell.min + cell.max) * 0.5;
    let half_diag = (cell.max - cell.min).length() * 0.5;

    // Sample SDF at cell center
    let dist = eval(node, center).abs();

    // Determine threshold: auto-compute from cell diagonal if user set 0
    let threshold = if config.surface_threshold > 0.0 {
        config.surface_threshold
    } else {
        half_diag * 1.2
    };

    // Subdivide if:
    // 1. Not at max depth yet
    // 2. Cell is near the surface (distance < threshold)
    // 3. Above minimum depth
    let should_subdivide =
        cell.depth < config.max_depth && (dist < threshold || cell.depth < config.min_depth);

    if should_subdivide {
        // Split into 8 octants
        for octant in 0..8 {
            let ox = (octant & 1) as f32;
            let oy = ((octant >> 1) & 1) as f32;
            let oz = ((octant >> 2) & 1) as f32;

            let child_min = Vec3::new(
                cell.min.x + ox * (center.x - cell.min.x),
                cell.min.y + oy * (center.y - cell.min.y),
                cell.min.z + oz * (center.z - cell.min.z),
            );
            let child_max = Vec3::new(
                center.x + ox * (cell.max.x - center.x),
                center.y + oy * (cell.max.y - center.y),
                center.z + oz * (cell.max.z - center.z),
            );

            let child = OctreeCell {
                min: child_min,
                max: child_max,
                depth: cell.depth + 1,
            };

            subdivide_octree(node, &child, config, leaves);
        }
    } else if dist < half_diag * 2.0 {
        // Only emit leaf if it could contain surface
        leaves.push(OctreeCell {
            min: cell.min,
            max: cell.max,
            depth: cell.depth,
        });
    }
}

/// Process a single adaptive cell using marching cubes
fn process_adaptive_cell(node: &SdfNode, cell: &OctreeCell, config: &MarchingCubesConfig) -> Mesh {
    let mut mesh = Mesh::new();
    let cell_size = cell.max - cell.min;

    // Evaluate 8 corners
    let mut corner_values = [0.0f32; 8];
    let mut corner_positions = [Vec3::ZERO; 8];

    #[allow(clippy::needless_range_loop)]
    for i in 0..8 {
        let dx = CORNER_OFFSETS[i][0] as f32;
        let dy = CORNER_OFFSETS[i][1] as f32;
        let dz = CORNER_OFFSETS[i][2] as f32;

        let pos = cell.min + Vec3::new(dx * cell_size.x, dy * cell_size.y, dz * cell_size.z);
        corner_positions[i] = pos;
        corner_values[i] = eval(node, pos);
    }

    // Compute cube index
    let mut cube_index = 0;
    #[allow(clippy::needless_range_loop)]
    for i in 0..8 {
        if corner_values[i] < config.iso_level {
            cube_index |= 1 << i;
        }
    }

    if EDGE_TABLE[cube_index] == 0 {
        return mesh;
    }

    // Compute edge vertices
    let mut edge_vertices = [Vec3::ZERO; 12];
    for i in 0..12 {
        if EDGE_TABLE[cube_index] & (1 << i) != 0 {
            let e0 = EDGE_CONNECTIONS[i][0];
            let e1 = EDGE_CONNECTIONS[i][1];
            edge_vertices[i] = interpolate_vertex(
                corner_positions[e0],
                corner_positions[e1],
                corner_values[e0],
                corner_values[e1],
                config.iso_level,
            );
        }
    }

    // Generate triangles
    let mut i = 0;
    while TRI_TABLE[cube_index][i] != -1 {
        let v0 = edge_vertices[TRI_TABLE[cube_index][i] as usize];
        let v1 = edge_vertices[TRI_TABLE[cube_index][i + 1] as usize];
        let v2 = edge_vertices[TRI_TABLE[cube_index][i + 2] as usize];

        let base_idx = mesh.vertices.len() as u32;

        let (n0, n1, n2) = if config.compute_normals {
            (
                normal(node, v0, 0.001),
                normal(node, v1, 0.001),
                normal(node, v2, 0.001),
            )
        } else {
            let face_normal = (v1 - v0).cross(v2 - v0).normalize();
            (face_normal, face_normal, face_normal)
        };

        let mut vert0 = Vertex::new(v0, n0);
        let mut vert1 = Vertex::new(v1, n1);
        let mut vert2 = Vertex::new(v2, n2);

        if config.compute_uvs {
            // Division Exorcism: precompute reciprocal once per triangle
            let inv_uv = 1.0 / config.uv_scale;
            vert0.uv = triplanar_uv(v0, n0, inv_uv);
            vert1.uv = triplanar_uv(v1, n1, inv_uv);
            vert2.uv = triplanar_uv(v2, n2, inv_uv);
        }

        if config.compute_materials {
            vert0.material_id = eval_material(node, v0);
            vert1.material_id = eval_material(node, v1);
            vert2.material_id = eval_material(node, v2);
        }

        mesh.vertices.push(vert0);
        mesh.vertices.push(vert1);
        mesh.vertices.push(vert2);

        mesh.indices.push(base_idx);
        mesh.indices.push(base_idx + 1);
        mesh.indices.push(base_idx + 2);

        i += 3;
    }

    mesh
}

// ---------------------------------------------------------------------------
// Adaptive Marching Cubes — Compiled VM variant
// ---------------------------------------------------------------------------

/// Adaptive marching cubes using compiled VM evaluator
///
/// Same octree subdivision logic as `adaptive_marching_cubes`, but uses
/// `eval_compiled` for distance queries and `normal_compiled` for normals.
/// Typically 2-5x faster than the interpreted variant.
pub fn adaptive_marching_cubes_compiled(
    sdf: &CompiledSdf,
    bounds_min: Vec3,
    bounds_max: Vec3,
    config: &AdaptiveConfig,
) -> Mesh {
    let root = OctreeCell {
        min: bounds_min,
        max: bounds_max,
        depth: 0,
    };

    let mut leaf_cells = Vec::new();
    subdivide_octree_compiled(sdf, &root, config, &mut leaf_cells);

    let mc_config = MarchingCubesConfig {
        resolution: 1,
        iso_level: config.iso_level,
        compute_normals: config.compute_normals,
        compute_uvs: config.compute_uvs,
        uv_scale: config.uv_scale,
        compute_tangents: false,
        compute_materials: false, // Compiled path doesn't support materials
    };

    let sub_meshes: Vec<Mesh> = leaf_cells
        .par_iter()
        .map(|cell| process_adaptive_cell_compiled(sdf, cell, &mc_config))
        .collect();

    let mut mesh = merge_meshes(sub_meshes);

    super::optimize::deduplicate_vertices(&mut mesh);

    if config.compute_tangents && config.compute_normals && config.compute_uvs {
        compute_mikktspace_tangents(&mut mesh);
    }

    mesh
}

/// Recursively subdivide octree using compiled evaluator
#[allow(clippy::suspicious_operation_groupings)]
fn subdivide_octree_compiled(
    sdf: &CompiledSdf,
    cell: &OctreeCell,
    config: &AdaptiveConfig,
    leaves: &mut Vec<OctreeCell>,
) {
    let center = (cell.min + cell.max) * 0.5;
    let half_diag = (cell.max - cell.min).length() * 0.5;

    let dist = eval_compiled(sdf, center).abs();

    let threshold = if config.surface_threshold > 0.0 {
        config.surface_threshold
    } else {
        half_diag * 1.2
    };

    let should_subdivide =
        cell.depth < config.max_depth && (dist < threshold || cell.depth < config.min_depth);

    if should_subdivide {
        for octant in 0..8 {
            let ox = (octant & 1) as f32;
            let oy = ((octant >> 1) & 1) as f32;
            let oz = ((octant >> 2) & 1) as f32;

            let child_min = Vec3::new(
                cell.min.x + ox * (center.x - cell.min.x),
                cell.min.y + oy * (center.y - cell.min.y),
                cell.min.z + oz * (center.z - cell.min.z),
            );
            let child_max = Vec3::new(
                center.x + ox * (cell.max.x - center.x),
                center.y + oy * (cell.max.y - center.y),
                center.z + oz * (cell.max.z - center.z),
            );

            let child = OctreeCell {
                min: child_min,
                max: child_max,
                depth: cell.depth + 1,
            };

            subdivide_octree_compiled(sdf, &child, config, leaves);
        }
    } else if dist < half_diag * 2.0 {
        leaves.push(OctreeCell {
            min: cell.min,
            max: cell.max,
            depth: cell.depth,
        });
    }
}

/// Process a single adaptive cell using compiled evaluator
fn process_adaptive_cell_compiled(
    sdf: &CompiledSdf,
    cell: &OctreeCell,
    config: &MarchingCubesConfig,
) -> Mesh {
    let mut mesh = Mesh::new();
    let cell_size = cell.max - cell.min;

    let mut corner_values = [0.0f32; 8];
    let mut corner_positions = [Vec3::ZERO; 8];

    #[allow(clippy::needless_range_loop)]
    for i in 0..8 {
        let dx = CORNER_OFFSETS[i][0] as f32;
        let dy = CORNER_OFFSETS[i][1] as f32;
        let dz = CORNER_OFFSETS[i][2] as f32;

        let pos = cell.min + Vec3::new(dx * cell_size.x, dy * cell_size.y, dz * cell_size.z);
        corner_positions[i] = pos;
        corner_values[i] = eval_compiled(sdf, pos);
    }

    let mut cube_index = 0;
    #[allow(clippy::needless_range_loop)]
    for i in 0..8 {
        if corner_values[i] < config.iso_level {
            cube_index |= 1 << i;
        }
    }

    if EDGE_TABLE[cube_index] == 0 {
        return mesh;
    }

    let mut edge_vertices = [Vec3::ZERO; 12];
    for i in 0..12 {
        if EDGE_TABLE[cube_index] & (1 << i) != 0 {
            let e0 = EDGE_CONNECTIONS[i][0];
            let e1 = EDGE_CONNECTIONS[i][1];
            edge_vertices[i] = interpolate_vertex(
                corner_positions[e0],
                corner_positions[e1],
                corner_values[e0],
                corner_values[e1],
                config.iso_level,
            );
        }
    }

    let mut i = 0;
    while TRI_TABLE[cube_index][i] != -1 {
        let v0 = edge_vertices[TRI_TABLE[cube_index][i] as usize];
        let v1 = edge_vertices[TRI_TABLE[cube_index][i + 1] as usize];
        let v2 = edge_vertices[TRI_TABLE[cube_index][i + 2] as usize];

        let base_idx = mesh.vertices.len() as u32;

        let (n0, n1, n2) = if config.compute_normals {
            (
                normal_compiled(sdf, v0, 0.001),
                normal_compiled(sdf, v1, 0.001),
                normal_compiled(sdf, v2, 0.001),
            )
        } else {
            let face_normal = (v1 - v0).cross(v2 - v0).normalize();
            (face_normal, face_normal, face_normal)
        };

        let mut vert0 = Vertex::new(v0, n0);
        let mut vert1 = Vertex::new(v1, n1);
        let mut vert2 = Vertex::new(v2, n2);

        if config.compute_uvs {
            // Division Exorcism: precompute reciprocal once per triangle
            let inv_uv = 1.0 / config.uv_scale;
            vert0.uv = triplanar_uv(v0, n0, inv_uv);
            vert1.uv = triplanar_uv(v1, n1, inv_uv);
            vert2.uv = triplanar_uv(v2, n2, inv_uv);
        }

        mesh.vertices.push(vert0);
        mesh.vertices.push(vert1);
        mesh.vertices.push(vert2);

        mesh.indices.push(base_idx);
        mesh.indices.push(base_idx + 1);
        mesh.indices.push(base_idx + 2);

        i += 3;
    }

    mesh
}

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

    #[test]
    fn test_sdf_to_mesh_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let config = MarchingCubesConfig {
            resolution: 16,
            iso_level: 0.0,
            compute_normals: true,
            ..Default::default()
        };

        let mesh = sdf_to_mesh(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        assert!(mesh.vertex_count() > 0);
        assert!(mesh.triangle_count() > 0);
    }

    #[test]
    fn test_sdf_to_mesh_box() {
        let box3d = SdfNode::box3d(2.0, 2.0, 2.0);
        let config = MarchingCubesConfig {
            resolution: 16,
            iso_level: 0.0,
            compute_normals: false,
            ..Default::default()
        };

        let mesh = sdf_to_mesh(&box3d, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        assert!(mesh.vertex_count() > 0);
        assert!(mesh.triangle_count() > 0);
    }

    #[test]
    fn test_mesh_vertex_normals() {
        let sphere = SdfNode::sphere(1.0);
        let config = MarchingCubesConfig {
            resolution: 8,
            iso_level: 0.0,
            compute_normals: true,
            ..Default::default()
        };

        let mesh = sdf_to_mesh(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        // Check that normals are normalized
        for vertex in &mesh.vertices {
            let len = vertex.normal.length();
            assert!((len - 1.0).abs() < 0.1);
        }
    }

    #[test]
    fn test_watertight_mesh() {
        // Verify that sdf_to_mesh produces a watertight mesh (shared vertices)
        let sphere = SdfNode::sphere(1.0);
        let config = MarchingCubesConfig {
            resolution: 8,
            iso_level: 0.0,
            compute_normals: true,
            ..Default::default()
        };

        let mesh = sdf_to_mesh(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        // After dedup, vertex count should be less than 3 * triangle_count
        // (shared vertices means fewer vertices than the per-triangle case)
        let tri_count = mesh.triangle_count();
        assert!(tri_count > 0);
        assert!(
            mesh.vertex_count() < tri_count * 3,
            "Expected shared vertices: {} verts < {} (3 * {} tris)",
            mesh.vertex_count(),
            tri_count * 3,
            tri_count
        );

        // Verify manifold: every edge should be shared by exactly 2 triangles
        use std::collections::HashMap;
        let mut edge_count: HashMap<(u32, u32), u32> = HashMap::new();
        for tri in mesh.indices.chunks_exact(3) {
            let edges = [
                (tri[0].min(tri[1]), tri[0].max(tri[1])),
                (tri[1].min(tri[2]), tri[1].max(tri[2])),
                (tri[0].min(tri[2]), tri[0].max(tri[2])),
            ];
            for &edge in &edges {
                *edge_count.entry(edge).or_insert(0) += 1;
            }
        }

        // Allow some boundary edges (count == 1) but no edge should have > 2 triangles
        let bad_edges = edge_count.values().filter(|&&c| c > 2).count();
        assert_eq!(
            bad_edges, 0,
            "Found {} non-manifold edges (shared by >2 triangles)",
            bad_edges
        );
    }

    #[test]
    fn test_mikktspace_tangents() {
        let sphere = SdfNode::sphere(1.0);
        let config = MarchingCubesConfig::aaa(16);

        let mesh = sdf_to_mesh(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        for v in &mesh.vertices {
            let t = Vec3::new(v.tangent.x, v.tangent.y, v.tangent.z);
            let n = v.normal;

            // Tangent should be normalized
            let t_len = t.length();
            assert!(
                (t_len - 1.0).abs() < 0.1,
                "Tangent not normalized: length = {}",
                t_len
            );

            // Tangent should be perpendicular to normal
            let dot = t.dot(n).abs();
            assert!(
                dot < 0.15,
                "Tangent not perpendicular to normal: dot = {}",
                dot
            );

            // Handedness should be -1 or +1
            let w = v.tangent.w;
            assert!(
                (w.abs() - 1.0).abs() < 0.01,
                "Invalid handedness: w = {}",
                w
            );
        }
    }

    #[test]
    fn test_adaptive_marching_cubes_basic() {
        let sphere = SdfNode::sphere(1.0);
        let config = AdaptiveConfig::default();

        let mesh = adaptive_marching_cubes(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        assert!(
            mesh.vertex_count() > 0,
            "Adaptive MC should produce vertices"
        );
        assert!(
            mesh.triangle_count() > 0,
            "Adaptive MC should produce triangles"
        );
    }

    #[test]
    fn test_adaptive_fewer_triangles() {
        // Small sphere in a large bounding box — most volume is empty,
        // so adaptive skips far-away cells while uniform evaluates everything.
        let sphere = SdfNode::sphere(1.0);

        // Uniform at resolution 64 over a large domain
        let uniform_config = MarchingCubesConfig {
            resolution: 64,
            iso_level: 0.0,
            compute_normals: true,
            ..Default::default()
        };
        let uniform_mesh = sdf_to_mesh(
            &sphere,
            Vec3::splat(-5.0),
            Vec3::splat(5.0),
            &uniform_config,
        );

        // Adaptive: same effective max resolution but tight surface threshold
        let adaptive_config = AdaptiveConfig {
            max_depth: 6,
            min_depth: 2,
            surface_threshold: 0.5,
            compute_normals: true,
            ..Default::default()
        };
        let adaptive_mesh = adaptive_marching_cubes(
            &sphere,
            Vec3::splat(-5.0),
            Vec3::splat(5.0),
            &adaptive_config,
        );

        // Adaptive should produce fewer or equal triangles
        assert!(
            adaptive_mesh.triangle_count() <= uniform_mesh.triangle_count(),
            "Adaptive ({}) should have <= triangles than uniform ({})",
            adaptive_mesh.triangle_count(),
            uniform_mesh.triangle_count()
        );
    }

    #[test]
    fn test_adaptive_aaa_preset() {
        let shape = SdfNode::sphere(1.0).smooth_union(SdfNode::box3d(0.5, 0.5, 0.5), 0.2);
        let config = AdaptiveConfig::aaa(5);

        let mesh = adaptive_marching_cubes(&shape, Vec3::splat(-2.0), Vec3::splat(2.0), &config);

        assert!(mesh.vertex_count() > 0);
        // Verify normals are normalized
        for v in &mesh.vertices {
            let len = v.normal.length();
            assert!((len - 1.0).abs() < 0.2, "Normal not normalized: {}", len);
        }
    }
}