plozone 0.1.0

//! Terrain mesh export (feature `terrain`).
//!
//! The octree point cloud encodes occupancy: occupied voxels are solid, empty
//! ones are air. This module provides the dependency-free pieces — a [`TriMesh`]
//! type, OBJ/PLY serializers, and a per-voxel density export
//! ([`marching_cubes_export_density`]) — plus an in-process Marching Cubes
//! surface extractor ([`marching_cubes`]) backed by the `mcubes` crate.

use std::fmt::Write as _;

use crate::octree::OctreeNode;

/// A triangle mesh: positions plus triangle indices into `vertices`.
#[derive(Debug, Clone, Default, PartialEq)]
pub struct TriMesh {
    pub vertices: Vec<[f32; 3]>,
    pub indices: Vec<[u32; 3]>,
}

/// Serialize a mesh to Wavefront OBJ text (1-based indices).
pub fn to_obj(mesh: &TriMesh) -> String {
    let mut out = String::new();
    for v in &mesh.vertices {
        let _ = writeln!(out, "v {} {} {}", v[0], v[1], v[2]);
    }
    for tri in &mesh.indices {
        let _ = writeln!(out, "f {} {} {}", tri[0] + 1, tri[1] + 1, tri[2] + 1);
    }
    out
}

/// Serialize a mesh to binary little-endian PLY.
///
/// Layout: ASCII header, then each vertex as 3 little-endian `f32` (12 bytes),
/// then each face as a `uchar` count of 3 followed by 3 little-endian `u32`.
pub fn to_ply(mesh: &TriMesh) -> Vec<u8> {
    let header = format!(
        "ply\nformat binary_little_endian 1.0\n\
         element vertex {}\nproperty float x\nproperty float y\nproperty float z\n\
         element face {}\nproperty list uchar uint vertex_indices\nend_header\n",
        mesh.vertices.len(),
        mesh.indices.len()
    );
    let mut out = header.into_bytes();
    for v in &mesh.vertices {
        for &f in v {
            out.extend_from_slice(&f.to_le_bytes());
        }
    }
    for tri in &mesh.indices {
        out.push(3u8);
        for &i in tri {
            out.extend_from_slice(&i.to_le_bytes());
        }
    }
    out
}

/// Export `(center, density)` for every occupied voxel materialised at `depth`.
///
/// `density` is `min(point_count / 8, 1.0)` — a crude occupancy estimate to feed
/// an isosurface extractor. Note that only voxels the octree actually
/// subdivided to are visited (empty regions never materialise), so this returns
/// occupied cells, not the full grid.
pub fn marching_cubes_export_density(octree: &OctreeNode, depth: u8) -> Vec<([f32; 3], f32)> {
    octree
        .nodes_at_depth(depth)
        .into_iter()
        .filter(|n| !n.points.is_empty())
        .map(|n| {
            let c = [n.center[0] as f32, n.center[1] as f32, n.center[2] as f32];
            let d = (n.points.len() as f32 / 8.0).min(1.0);
            (c, d)
        })
        .collect()
}

/// Build a triangle mesh from octree point density using Marching Cubes.
///
/// Samples the octree on a regular `(resolution+1)³` grid, computes a density
/// value at each sample point (fraction of nearby octree points), and extracts
/// the isosurface at `isolevel` using the `mcubes` crate.
///
/// `resolution` controls grid resolution per axis (e.g. 16 → 17³ = 4913 samples).
/// Higher = finer mesh, slower extraction. The octree's world bounds map to the
/// full grid extent.
///
/// Returns an empty mesh if the octree has no points or the isosurface doesn't
/// intersect the density field.
pub fn marching_cubes(octree: &OctreeNode, resolution: usize, isolevel: f32) -> TriMesh {
    if resolution == 0 || octree.points.is_empty() && octree.children.is_none() {
        return TriMesh::default();
    }

    let n = resolution + 1;
    let total = n * n * n;
    if total > 5_000_000 {
        return TriMesh::default();
    }

    let half = octree.half_len;
    let world_size = half * 2.0;
    let step = world_size / resolution as f64;
    let origin = [octree.center[0] - half, octree.center[1] - half, octree.center[2] - half];

    let voxel = step;

    let mut values = vec![0.0f32; total];
    for ix in 0..n {
        for iy in 0..n {
            for iz in 0..n {
                let x = origin[0] + ix as f64 * step;
                let y = origin[1] + iy as f64 * step;
                let z = origin[2] + iz as f64 * step;

                let pts = octree.range_query(
                    [x - voxel, y - voxel, z - voxel],
                    [x + voxel, y + voxel, z + voxel],
                );
                let density = (pts.len() as f32 / 8.0).min(1.0);
                values[ix * n * n + iy * n + iz] = density;
            }
        }
    }

    let dims = (n, n, n);
    let size_f = (
        world_size as f32,
        world_size as f32,
        world_size as f32,
    );
    let interval = (step as f32, step as f32, step as f32);
    let offset = lin_alg::f32::Vec3::new(
        origin[0] as f32,
        origin[1] as f32,
        origin[2] as f32,
    );

    let mc = match mcubes::MarchingCubes::new(
        dims,
        size_f,
        interval,
        offset,
        values,
        isolevel,
    ) {
        Ok(mc) => mc,
        Err(_) => return TriMesh::default(),
    };

    let mesh = mc.generate(mcubes::MeshSide::Both);

    let vertices: Vec<[f32; 3]> = mesh
        .vertices
        .iter()
        .map(|v| [v.posit.x, v.posit.y, v.posit.z])
        .collect();

    let indices: Vec<[u32; 3]> = mesh
        .indices
        .chunks_exact(3)
        .map(|tri| [tri[0] as u32, tri[1] as u32, tri[2] as u32])
        .collect();

    TriMesh { vertices, indices }
}

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

    fn single_triangle() -> TriMesh {
        TriMesh {
            vertices: vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
            indices: vec![[0, 1, 2]],
        }
    }

    #[test]
    fn obj_has_vertices_and_one_based_faces() {
        let obj = to_obj(&single_triangle());
        assert!(obj.contains("v 0 0 0"));
        assert!(obj.contains("v 1 0 0"));
        // Indices are emitted 1-based.
        assert!(obj.contains("f 1 2 3"), "got:\n{obj}");
        assert_eq!(obj.lines().filter(|l| l.starts_with("v ")).count(), 3);
        assert_eq!(obj.lines().filter(|l| l.starts_with("f ")).count(), 1);
    }

    #[test]
    fn ply_header_and_byte_layout() {
        let mesh = single_triangle();
        let bytes = to_ply(&mesh);
        let header_end = b"end_header\n";
        // Find the header boundary.
        let pos = bytes
            .windows(header_end.len())
            .position(|w| w == header_end)
            .map(|p| p + header_end.len())
            .expect("header present");
        let header = std::str::from_utf8(&bytes[..pos]).unwrap();
        assert!(header.contains("element vertex 3"));
        assert!(header.contains("element face 1"));
        assert!(header.contains("binary_little_endian"));

        // Body: 3 vertices × 12 bytes + 1 face × (1 + 12) bytes.
        let body = bytes.len() - pos;
        assert_eq!(body, 3 * 12 + (1 + 12));

        // First vertex decodes back to (0,0,0).
        let x = f32::from_le_bytes(bytes[pos..pos + 4].try_into().unwrap());
        assert_eq!(x, 0.0);
        // The face block starts with the count byte 3.
        assert_eq!(bytes[pos + 3 * 12], 3u8);
    }

    #[test]
    fn density_export_reports_occupied_voxels() {
        // High threshold keeps points at the root (depth 0) without subdividing.
        let mut octree = OctreeNode::new([0.0; 3], 50.0);
        octree.insert([1.0, 1.0, 1.0], 100);
        octree.insert([2.0, 2.0, 2.0], 100);
        let cells = marching_cubes_export_density(&octree, 0);
        assert_eq!(cells.len(), 1, "root holds both points");
        assert_eq!(cells[0].0, [0.0, 0.0, 0.0], "root center");
        assert!((cells[0].1 - 2.0 / 8.0).abs() < 1e-6, "density = 2/8");
    }

    #[test]
    fn density_export_after_subdivision() {
        let mut octree = OctreeNode::new([0.0; 3], 50.0);
        // Force subdivision: spread points across octants past a low threshold.
        for p in [
            [10.0, 10.0, 10.0],
            [-10.0, 10.0, 10.0],
            [10.0, -10.0, 10.0],
            [-10.0, -10.0, -10.0],
        ] {
            octree.insert(p, 1);
        }
        let cells = marching_cubes_export_density(&octree, 1);
        assert!(!cells.is_empty(), "occupied depth-1 voxels exist after subdivision");
        // Every reported center lies within the root cube.
        for (c, _) in &cells {
            for axis in c {
                assert!(axis.abs() <= 50.0, "center within world: {c:?}");
            }
        }
    }

    #[test]
    fn marching_cubes_produces_mesh_from_points() {
        let mut octree = OctreeNode::new([0.0; 3], 50.0);
        // Insert a cluster of points to create a density blob.
        for x in -5..=5 {
            for y in -5..=5 {
                for z in -5..=5 {
                    octree.insert([x as f64, y as f64, z as f64], 100);
                }
            }
        }
        let mesh = marching_cubes(&octree, 16, 0.05);
        assert!(!mesh.vertices.is_empty(), "MC should produce vertices from point cluster");
        assert!(!mesh.indices.is_empty(), "MC should produce triangles from point cluster");
    }

    #[test]
    fn marching_cubes_empty_octree_returns_empty() {
        let octree = OctreeNode::new([0.0; 3], 50.0);
        let mesh = marching_cubes(&octree, 8, 0.5);
        assert!(mesh.vertices.is_empty());
        assert!(mesh.indices.is_empty());
    }

    #[test]
    fn marching_cubes_zero_resolution_returns_empty() {
        let mut octree = OctreeNode::new([0.0; 3], 50.0);
        octree.insert([1.0, 2.0, 3.0], 8);
        let mesh = marching_cubes(&octree, 0, 0.5);
        assert!(mesh.vertices.is_empty());
    }
}