alice-sdf 5.0.0

//! Mesh to SDF conversion (Deep Fried Edition v2)
//!
//! Converts polygon meshes to SDF representations with multiple strategies:
//! - **BVH-Accelerated Exact**: O(log n) queries using BVH
//! - **Capsule Approximation**: Fast but approximate (legacy mode)
//! - **Hybrid**: BVH for distance, capsule tree for SdfNode
//!
//! # Deep Fried Optimizations
//! - **Edge Deduplication**: Uses HashSet to eliminate duplicate capsules.
//! - **BVH Acceleration**: O(log n) distance queries for exact SDF.
//! - **Parallel Batch Evaluation**: Rayon-powered batch queries.
//! - **Forced Inlining**: `#[inline(always)]` on hot-path helpers.
//!
//! Author: Moroya Sakamoto

use crate::mesh::bvh::MeshBvh;
use crate::mesh::mesh_sign::{ExteriorField, MeshSignMode};
use crate::mesh::MeshInputError;
use crate::types::SdfNode;
use glam::Vec3;
use rayon::prelude::*;
use std::collections::HashSet;
use std::sync::Arc;

/// Cells along the longest bounding box axis used by
/// [`MeshSignMode::ExteriorFloodFill`] unless the caller overrides it.
pub const DEFAULT_SIGN_FLOOD_FILL_RESOLUTION: u32 = 64;

/// Conversion strategy for mesh_to_sdf
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum MeshToSdfStrategy {
    /// Use capsule approximation (fast, creates SdfNode tree)
    #[default]
    Capsule,
    /// Use BVH for exact distance (accurate, returns MeshSdf wrapper)
    BvhExact,
    /// Hybrid: create capsule tree but use BVH for evaluation
    Hybrid,
}

/// Configuration for mesh to SDF conversion
#[derive(Debug, Clone)]
pub struct MeshToSdfConfig {
    /// Conversion strategy
    pub strategy: MeshToSdfStrategy,
    /// Whether to use bounding volumes for acceleration
    pub use_bvh: bool,
    /// Maximum triangles per leaf node (if using BVH)
    pub max_triangles_per_leaf: usize,
    /// Capsule radius factor (multiplied by average edge length)
    pub capsule_radius_factor: f32,
    /// Whether to compute vertex normals for sign determination
    pub compute_normals: bool,
    /// How the sign of the resulting field is decided.
    ///
    /// Defaults to [`MeshSignMode::NearestFaceNormal`], which needs a closed,
    /// consistently wound mesh. Use [`MeshSignMode::ExteriorFloodFill`] for
    /// meshes you did not author — open surfaces, mixed winding and separate
    /// components all come out right.
    pub sign_mode: MeshSignMode,
    /// Cells along the longest bounding box axis for
    /// [`MeshSignMode::ExteriorFloodFill`]. Ignored by the other modes.
    ///
    /// The thinnest feature that keeps its interior spans a few cells, so raise
    /// this for thin-walled meshes.
    pub sign_flood_fill_resolution: u32,
}

impl Default for MeshToSdfConfig {
    fn default() -> Self {
        Self {
            strategy: MeshToSdfStrategy::Capsule,
            use_bvh: true,
            max_triangles_per_leaf: 4,
            capsule_radius_factor: 0.05,
            compute_normals: true,
            sign_mode: MeshSignMode::NearestFaceNormal,
            sign_flood_fill_resolution: DEFAULT_SIGN_FLOOD_FILL_RESOLUTION,
        }
    }
}

impl MeshToSdfConfig {
    /// Create config for fast capsule approximation
    pub fn fast() -> Self {
        Self {
            strategy: MeshToSdfStrategy::Capsule,
            use_bvh: false,
            ..Default::default()
        }
    }

    /// Create config for accurate BVH-based SDF.
    ///
    /// Uses [`MeshSignMode::ExteriorFloodFill`], so the sign is a function of the
    /// geometry rather than of the triangle winding: an open surface reports no
    /// interior, reversing the winding changes nothing, and each closed component
    /// gets its own interior. Pair with [`Self::sign_flood_fill_resolution`] when
    /// the mesh has features thinner than a few grid cells.
    ///
    /// For the pre-4.0 behaviour, set `sign_mode` back to
    /// [`MeshSignMode::NearestFaceNormal`].
    pub fn accurate() -> Self {
        Self {
            strategy: MeshToSdfStrategy::BvhExact,
            use_bvh: true,
            max_triangles_per_leaf: 4,
            sign_mode: MeshSignMode::ExteriorFloodFill,
            ..Default::default()
        }
    }

    /// Create config for hybrid approach
    pub fn hybrid() -> Self {
        Self {
            strategy: MeshToSdfStrategy::Hybrid,
            use_bvh: true,
            max_triangles_per_leaf: 4,
            ..Default::default()
        }
    }

    /// Exact BVH distance with a sign that does not depend on triangle winding.
    ///
    /// Identical to [`Self::accurate`] — delegates to it so the two cannot
    /// drift — and spelled out for call sites where the reason for the choice is
    /// the mesh's provenance rather than the accuracy: scans, generative models,
    /// assemblies of independently wound parts. See
    /// [`MeshSignMode::ExteriorFloodFill`] for the resolution caveat.
    pub fn topology_robust() -> Self {
        Self::accurate()
    }
}

/// BVH-accelerated mesh SDF
///
/// Provides exact signed distance queries using BVH acceleration.
/// This is more accurate than the capsule approximation but cannot be
/// converted to shader code directly.
pub struct MeshSdf {
    bvh: Arc<MeshBvh>,
    bounds: crate::types::Aabb,
    sign_mode: MeshSignMode,
    /// `Some` exactly when `sign_mode` is [`MeshSignMode::ExteriorFloodFill`];
    /// construction fails rather than silently falling back to the winding rule.
    exterior: Option<Arc<ExteriorField>>,
}

impl MeshSdf {
    /// Create MeshSdf from vertices and indices.
    ///
    /// Returns `None` on the errors [`MeshSdf::try_new`] reports; use that form
    /// when you need to know which one.
    pub fn new(vertices: &[Vec3], indices: &[u32], config: &MeshToSdfConfig) -> Option<Self> {
        Self::try_new(vertices, indices, config).ok()
    }

    /// Create MeshSdf from vertices and indices.
    ///
    /// # Errors
    ///
    /// Empty input, a mesh with no bounding box, or — for
    /// [`MeshSignMode::ExteriorFloodFill`] — a degenerate extent or a
    /// `sign_flood_fill_resolution` over the cell budget.
    pub fn try_new(
        vertices: &[Vec3],
        indices: &[u32],
        config: &MeshToSdfConfig,
    ) -> Result<Self, MeshInputError> {
        if vertices.is_empty() || indices.is_empty() {
            return Err(MeshInputError {
                reason: "mesh_to_sdf needs at least one triangle",
            });
        }

        let bvh = MeshBvh::build(vertices, indices, config.max_triangles_per_leaf);
        let bounds = bvh.bounds().ok_or(MeshInputError {
            reason: "mesh has no bounding box",
        })?;

        let exterior = match config.sign_mode {
            MeshSignMode::ExteriorFloodFill => Some(Arc::new(ExteriorField::build(
                &bvh,
                config.sign_flood_fill_resolution,
            )?)),
            _ => None,
        };

        Ok(Self {
            bvh: Arc::new(bvh),
            bounds,
            sign_mode: config.sign_mode,
            exterior,
        })
    }

    /// Which rule decides the sign of [`Self::eval`].
    #[inline]
    pub const fn sign_mode(&self) -> MeshSignMode {
        self.sign_mode
    }

    /// Evaluate signed distance at a point
    #[inline]
    pub fn eval(&self, point: Vec3) -> f32 {
        self.exterior.as_ref().map_or_else(
            || self.bvh.signed_distance(point),
            |field| field.signed_distance(&self.bvh, point),
        )
    }

    /// Evaluate unsigned distance at a point
    #[inline]
    pub fn eval_unsigned(&self, point: Vec3) -> f32 {
        self.bvh.unsigned_distance(point)
    }

    /// Batch evaluate signed distances (parallel)
    pub fn eval_batch(&self, points: &[Vec3]) -> Vec<f32> {
        self.exterior.as_ref().map_or_else(
            || self.bvh.signed_distance_batch(points),
            |field| {
                points
                    .par_iter()
                    .map(|&p| field.signed_distance(&self.bvh, p))
                    .collect()
            },
        )
    }

    /// Batch evaluate unsigned distances (parallel)
    pub fn eval_unsigned_batch(&self, points: &[Vec3]) -> Vec<f32> {
        self.bvh.unsigned_distance_batch(points)
    }

    /// Compute gradient (approximate normal) at a point
    pub fn gradient(&self, point: Vec3, epsilon: f32) -> Vec3 {
        let dx = self.eval(point + Vec3::X * epsilon) - self.eval(point - Vec3::X * epsilon);
        let dy = self.eval(point + Vec3::Y * epsilon) - self.eval(point - Vec3::Y * epsilon);
        let dz = self.eval(point + Vec3::Z * epsilon) - self.eval(point - Vec3::Z * epsilon);
        Vec3::new(dx, dy, dz).normalize_or_zero()
    }

    /// Get mesh bounds
    pub const fn bounds(&self) -> (Vec3, Vec3) {
        (self.bounds.min, self.bounds.max)
    }

    /// Get triangle count
    pub fn triangle_count(&self) -> usize {
        self.bvh.triangle_count()
    }

    /// Get reference to underlying BVH
    pub fn bvh(&self) -> &MeshBvh {
        &self.bvh
    }

    /// Convert to approximate SdfNode tree using capsule method
    ///
    /// This allows the MeshSdf to be used where SdfNode is required,
    /// but with reduced accuracy.
    pub fn to_sdf_node(&self, radius_factor: f32) -> SdfNode {
        // Extract vertices and edges from triangles
        let mut vertices: Vec<Vec3> = Vec::new();
        let mut vertex_map: std::collections::HashMap<[u32; 3], usize> =
            std::collections::HashMap::new();
        let mut unique_edges: HashSet<(usize, usize)> = HashSet::new();

        for tri in &self.bvh.triangles {
            let mut get_or_add = |v: Vec3| -> usize {
                let key = [v.x.to_bits(), v.y.to_bits(), v.z.to_bits()];
                if let Some(&idx) = vertex_map.get(&key) {
                    idx
                } else {
                    let idx = vertices.len();
                    vertices.push(v);
                    vertex_map.insert(key, idx);
                    idx
                }
            };

            let i0 = get_or_add(tri.v0);
            let i1 = get_or_add(tri.v1);
            let i2 = get_or_add(tri.v2);

            unique_edges.insert(edge_key_usize(i0, i1));
            unique_edges.insert(edge_key_usize(i1, i2));
            unique_edges.insert(edge_key_usize(i2, i0));
        }

        // Compute average edge length
        let avg_edge_len: f32 = if unique_edges.is_empty() {
            0.1
        } else {
            let total: f32 = unique_edges
                .iter()
                .map(|&(i0, i1)| (vertices[i1] - vertices[i0]).length())
                .sum();
            total / unique_edges.len() as f32
        };

        let radius = avg_edge_len * radius_factor;

        // Build capsule nodes
        let nodes: Vec<SdfNode> = unique_edges
            .into_iter()
            .map(|(i0, i1)| SdfNode::capsule(vertices[i0], vertices[i1], radius))
            .collect();

        reduce_to_union(nodes)
    }
}

/// Convert a polygon mesh to an SDF representation
///
/// Returns an SdfNode tree that approximates the mesh using capsules.
/// For more accurate SDF evaluation, use `mesh_to_sdf_exact` instead.
///
/// # Arguments
/// * `vertices` - Mesh vertex positions
/// * `indices` - Triangle indices (3 per triangle)
/// * `config` - Conversion configuration
///
/// # Returns
/// SDF tree approximating the mesh
pub fn mesh_to_sdf(vertices: &[Vec3], indices: &[u32], config: &MeshToSdfConfig) -> SdfNode {
    if vertices.is_empty() || indices.is_empty() {
        return SdfNode::sphere(0.001); // Degenerate case
    }

    // Compute average edge length for capsule radius
    let mut total_edge_len = 0.0;
    let mut edge_count = 0;

    for chunk in indices.chunks(3) {
        if chunk.len() == 3 {
            let v0 = vertices[chunk[0] as usize];
            let v1 = vertices[chunk[1] as usize];
            let v2 = vertices[chunk[2] as usize];

            total_edge_len += (v1 - v0).length();
            total_edge_len += (v2 - v1).length();
            total_edge_len += (v0 - v2).length();
            edge_count += 3;
        }
    }

    let avg_edge_len = if edge_count > 0 {
        total_edge_len / edge_count as f32
    } else {
        0.1
    };

    // Use configured fraction of edge length as capsule radius
    let radius = avg_edge_len * config.capsule_radius_factor;

    // Edge deduplication using HashSet
    // Key: (min_idx, max_idx) ensures (a, b) == (b, a)
    let mut unique_edges: HashSet<(u32, u32)> = HashSet::new();

    for chunk in indices.chunks(3) {
        if chunk.len() == 3 {
            let i0 = chunk[0];
            let i1 = chunk[1];
            let i2 = chunk[2];

            // Insert edges with canonical ordering (min, max)
            unique_edges.insert(edge_key(i0, i1));
            unique_edges.insert(edge_key(i1, i2));
            unique_edges.insert(edge_key(i2, i0));
        }
    }

    // Build SDF from unique edges only
    let nodes: Vec<SdfNode> = unique_edges
        .into_iter()
        .map(|(i0, i1)| {
            let v0 = vertices[i0 as usize];
            let v1 = vertices[i1 as usize];
            SdfNode::capsule(v0, v1, radius)
        })
        .collect();

    // Build balanced binary tree of unions
    reduce_to_union(nodes)
}

/// Convert a polygon mesh to an exact BVH-accelerated SDF
///
/// Returns a MeshSdf that provides accurate signed distance queries.
/// This cannot be converted to shader code directly, but is useful for:
/// - Accurate distance queries
/// - Mesh processing algorithms
/// - Ray marching on CPU
///
/// # Arguments
/// * `vertices` - Mesh vertex positions
/// * `indices` - Triangle indices (3 per triangle)
/// * `config` - Conversion configuration
///
/// # Returns
/// `Option<MeshSdf>` - None if mesh is empty
pub fn mesh_to_sdf_exact(
    vertices: &[Vec3],
    indices: &[u32],
    config: &MeshToSdfConfig,
) -> Option<MeshSdf> {
    MeshSdf::new(vertices, indices, config)
}

/// Create canonical edge key (min, max) for HashSet deduplication
#[inline(always)]
const fn edge_key(a: u32, b: u32) -> (u32, u32) {
    if a < b {
        (a, b)
    } else {
        (b, a)
    }
}

/// Create canonical edge key for usize indices
#[inline(always)]
const fn edge_key_usize(a: usize, b: usize) -> (usize, usize) {
    if a < b {
        (a, b)
    } else {
        (b, a)
    }
}

/// Reduce a list of nodes to a single union using binary tree structure (Deep Fried)
#[inline(always)]
fn reduce_to_union(mut nodes: Vec<SdfNode>) -> SdfNode {
    if nodes.is_empty() {
        return SdfNode::sphere(0.001);
    }

    while nodes.len() > 1 {
        let mut next = Vec::with_capacity(nodes.len().div_ceil(2));

        for chunk in nodes.chunks(2) {
            if chunk.len() == 2 {
                next.push(chunk[0].clone().union(chunk[1].clone()));
            } else {
                next.push(chunk[0].clone());
            }
        }

        nodes = next;
    }

    nodes.pop().unwrap()
}

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

    #[test]
    fn test_mesh_to_sdf_single_triangle() {
        let vertices = vec![
            Vec3::new(0.0, 0.0, 0.0),
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(0.5, 1.0, 0.0),
        ];
        let indices = vec![0, 1, 2];

        let sdf = mesh_to_sdf(&vertices, &indices, &MeshToSdfConfig::default());

        // Should be a valid SDF tree
        assert!(sdf.node_count() > 0);
    }

    #[test]
    fn test_mesh_to_sdf_cube() {
        // Simple cube mesh
        let vertices = vec![
            Vec3::new(-1.0, -1.0, -1.0),
            Vec3::new(1.0, -1.0, -1.0),
            Vec3::new(1.0, 1.0, -1.0),
            Vec3::new(-1.0, 1.0, -1.0),
            Vec3::new(-1.0, -1.0, 1.0),
            Vec3::new(1.0, -1.0, 1.0),
            Vec3::new(1.0, 1.0, 1.0),
            Vec3::new(-1.0, 1.0, 1.0),
        ];
        // Front face, back face (simplified - just 2 triangles for test)
        let indices = vec![
            0, 1, 2, 0, 2, 3, // Front
            4, 6, 5, 4, 7, 6, // Back
        ];

        let sdf = mesh_to_sdf(&vertices, &indices, &MeshToSdfConfig::default());
        assert!(sdf.node_count() > 0);
    }

    #[test]
    fn test_mesh_to_sdf_empty() {
        let vertices: Vec<Vec3> = vec![];
        let indices: Vec<u32> = vec![];

        let sdf = mesh_to_sdf(&vertices, &indices, &MeshToSdfConfig::default());
        // Should return degenerate sphere
        assert_eq!(sdf.node_count(), 1);
    }

    #[test]
    fn test_reduce_to_union() {
        let nodes = vec![
            SdfNode::sphere(1.0),
            SdfNode::sphere(1.0),
            SdfNode::sphere(1.0),
            SdfNode::sphere(1.0),
        ];

        let union = reduce_to_union(nodes);
        // 4 spheres + 2 unions at level 1 + 1 union at root = 7
        assert_eq!(union.node_count(), 7);
    }

    #[test]
    fn test_mesh_sdf_exact() {
        let vertices = vec![
            Vec3::new(-1.0, -1.0, 0.0),
            Vec3::new(1.0, -1.0, 0.0),
            Vec3::new(1.0, 1.0, 0.0),
            Vec3::new(-1.0, 1.0, 0.0),
        ];
        let indices = vec![0, 1, 2, 0, 2, 3];

        let mesh_sdf = mesh_to_sdf_exact(&vertices, &indices, &MeshToSdfConfig::accurate())
            .expect("Should create MeshSdf");

        // Point on surface
        let d_surface = mesh_sdf.eval(Vec3::ZERO);
        assert!(d_surface.abs() < 0.01, "Expected ~0, got {}", d_surface);

        // Point above. The quad is an open sheet, so it encloses nothing and the
        // signed distance is +1 on *both* sides — this used to be asserted on
        // `.abs()` because the winding rule made the sign arbitrary.
        let d_above = mesh_sdf.eval(Vec3::new(0.0, 0.0, 1.0));
        assert!(
            (d_above - 1.0).abs() < 0.01,
            "Expected ~+1 above an open sheet, got {}",
            d_above
        );
        let d_below = mesh_sdf.eval(Vec3::new(0.0, 0.0, -1.0));
        assert!(
            (d_below - 1.0).abs() < 0.01,
            "Expected ~+1 below an open sheet too, got {}",
            d_below
        );

        // Batch evaluation
        let points = vec![
            Vec3::ZERO,
            Vec3::new(0.0, 0.0, 1.0),
            Vec3::new(0.0, 0.0, -1.0),
        ];
        let distances = mesh_sdf.eval_batch(&points);
        assert_eq!(distances.len(), 3);
    }

    #[test]
    fn test_mesh_sdf_to_node() {
        let vertices = vec![
            Vec3::new(0.0, 0.0, 0.0),
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(0.5, 1.0, 0.0),
        ];
        let indices = vec![0, 1, 2];

        let mesh_sdf = mesh_to_sdf_exact(&vertices, &indices, &MeshToSdfConfig::accurate())
            .expect("Should create MeshSdf");

        let node = mesh_sdf.to_sdf_node(0.05);
        assert!(node.node_count() > 0);
    }

    #[test]
    fn test_config_presets() {
        let fast = MeshToSdfConfig::fast();
        assert_eq!(fast.strategy, MeshToSdfStrategy::Capsule);
        assert!(!fast.use_bvh);

        let accurate = MeshToSdfConfig::accurate();
        assert_eq!(accurate.strategy, MeshToSdfStrategy::BvhExact);
        assert!(accurate.use_bvh);
        assert_eq!(accurate.sign_mode, MeshSignMode::ExteriorFloodFill);

        // documented as an alias, so pin that rather than repeating the fields
        assert_eq!(
            MeshToSdfConfig::topology_robust().sign_mode,
            accurate.sign_mode
        );
        assert_eq!(
            MeshToSdfConfig::topology_robust().strategy,
            accurate.strategy
        );

        // the default stays on the winding rule so existing callers are unaffected
        assert_eq!(
            MeshToSdfConfig::default().sign_mode,
            MeshSignMode::NearestFaceNormal
        );

        let hybrid = MeshToSdfConfig::hybrid();
        assert_eq!(hybrid.strategy, MeshToSdfStrategy::Hybrid);
        assert!(hybrid.use_bvh);
    }
}