curvo 0.1.88

pub mod constrained_triangulation;
pub mod trimmed_surface_constraints;
pub mod trimmed_surface_ext;

use std::collections::HashMap;

use super::adaptive_tessellation_node::AdaptiveTessellationNode;
use super::adaptive_tessellation_option::AdaptiveTessellationOptions;
use super::advancing_front::AdvancingFrontOptions;
use super::surface_tessellation::SurfaceTessellation;
use super::{ConstrainedTessellation, DividableDirection, Tessellation};
use itertools::Itertools;
use nalgebra::allocator::Allocator;
use nalgebra::{ComplexField, DefaultAllocator, DimName, OVector, Point2, Point3, Vector3};
use nalgebra::{Vector2, U4};
use spade::{HasPosition, SpadeNum, Triangulation};

use crate::curve::NurbsCurve2D;
use crate::misc::FloatingPoint;
use crate::misc::PolygonBoundary;
use crate::prelude::{Contains, SurfaceTessellation3D};
use crate::region::CompoundCurve2D;
use crate::surface::TrimmedSurface;
use crate::tessellation::trimmed_surface::trimmed_surface_ext::TrimmedSurfaceExt;
pub use constrained_triangulation::TrimmedSurfaceConstrainedTriangulation;
pub use trimmed_surface_constraints::*;

#[derive(Debug, Clone, Copy)]
pub struct Vertex<T: FloatingPoint> {
    point: Point3<T>,
    normal: Vector3<T>,
    uv: Vector2<T>,
}

impl<T: FloatingPoint> Vertex<T> {
    pub fn new(point: Point3<T>, normal: Vector3<T>, uv: Vector2<T>) -> Self {
        Self { point, normal, uv }
    }

    pub fn point(&self) -> Point3<T> {
        self.point
    }

    pub fn normal(&self) -> Vector3<T> {
        self.normal
    }

    pub fn uv(&self) -> Vector2<T> {
        self.uv
    }
}

impl<T: FloatingPoint + SpadeNum> HasPosition for Vertex<T> {
    type Scalar = T;

    fn position(&self) -> spade::Point2<Self::Scalar> {
        spade::Point2::from([self.uv.x, self.uv.y])
    }
}

impl<T: FloatingPoint + SpadeNum, F> Tessellation<Option<AdaptiveTessellationOptions<T, U4, F>>>
    for TrimmedSurface<T>
where
    F: Fn(&AdaptiveTessellationNode<T, U4>) -> Option<DividableDirection> + Copy,
{
    type Output = anyhow::Result<SurfaceTessellation3D<T>>;

    /// Tessellate a trimmed surface using an adaptive algorithm
    fn tessellate(&self, options: Option<AdaptiveTessellationOptions<T, U4, F>>) -> Self::Output {
        trimmed_surface_adaptive_tessellate(self, None, options)
    }
}

impl<T: FloatingPoint + SpadeNum, F>
    ConstrainedTessellation<Option<AdaptiveTessellationOptions<T, U4, F>>> for TrimmedSurface<T>
where
    F: Fn(&AdaptiveTessellationNode<T, U4>) -> Option<DividableDirection> + Copy,
{
    type Constraint = TrimmedSurfaceConstraints<T>;
    type Output = anyhow::Result<SurfaceTessellation3D<T>>;

    /// Tessellate a trimmed surface using an adaptive algorithm with constraints
    fn constrained_tessellate(
        &self,
        constraints: Self::Constraint,
        adaptive_options: Option<AdaptiveTessellationOptions<T, U4, F>>,
    ) -> Self::Output {
        trimmed_surface_adaptive_tessellate(self, Some(constraints), adaptive_options)
    }
}

impl<T: FloatingPoint + SpadeNum> Tessellation<AdvancingFrontOptions<T>> for TrimmedSurface<T> {
    type Output = anyhow::Result<SurfaceTessellation3D<T>>;

    fn tessellate(&self, options: AdvancingFrontOptions<T>) -> Self::Output {
        super::advancing_front::AdvancingFrontMesher::new(self, options).mesh()
    }
}

impl<T: FloatingPoint + SpadeNum> ConstrainedTessellation<AdvancingFrontOptions<T>>
    for TrimmedSurface<T>
{
    type Constraint = TrimmedSurfaceConstraints<T>;
    type Output = anyhow::Result<SurfaceTessellation3D<T>>;

    fn constrained_tessellate(
        &self,
        constraints: Self::Constraint,
        options: AdvancingFrontOptions<T>,
    ) -> Self::Output {
        super::advancing_front::AdvancingFrontMesher::new(self, options)
            .mesh_with_constraints(Some(constraints))
    }
}

/// Tessellate a trimmed surface using an adaptive algorithm with or without constraints
pub fn trimmed_surface_adaptive_tessellate<
    T: FloatingPoint + SpadeNum,
    S: TrimmedSurfaceExt<T, F>,
    F,
>(
    s: &S,
    constraints: Option<TrimmedSurfaceConstraints<T>>,
    options: Option<AdaptiveTessellationOptions<T, U4, F>>,
) -> anyhow::Result<SurfaceTessellation3D<T>>
where
    F: Fn(&AdaptiveTessellationNode<T, U4>) -> Option<DividableDirection> + Copy,
{
    let TrimmedSurfaceConstrainedTriangulation {
        cdt,
        exterior,
        interiors,
    } = TrimmedSurfaceConstrainedTriangulation::try_new(s, constraints, options)?;

    let vmap: HashMap<_, _> = cdt
        .vertices()
        .enumerate()
        .map(|(i, v)| (v.fix(), i))
        .collect();

    let uv_exterior_boundary =
        exterior.map(|c| PolygonBoundary::new(c.iter().map(|v| v.uv.into()).collect()));

    let uv_interior_boundaries = interiors
        .into_iter()
        .map(|c| PolygonBoundary::new(c.iter().map(|v| v.uv.into()).collect()))
        .collect_vec();

    let inv_3 = T::from_f64(1. / 3.).unwrap();

    let faces = cdt
        .inner_faces()
        .filter_map(|f| {
            let vs = f.vertices();
            let tri = vs.iter().map(|v| v.as_ref()).map(|p| p.uv).collect_vec();

            let (a, b) = (tri[1] - tri[0], tri[2] - tri[1]);
            let area = a.x * b.y - a.y * b.x;
            if ComplexField::abs(area) < T::default_epsilon() {
                None
            } else {
                let center: Point2<T> = ((tri[0] + tri[1] + tri[2]) * inv_3).into();
                let exterior_contains = uv_exterior_boundary
                    .as_ref()
                    .map(|exterior| exterior.contains(&center, ()).unwrap_or(false))
                    .unwrap_or(true);
                let interior_contains = !uv_interior_boundaries.is_empty()
                    && uv_interior_boundaries
                        .iter()
                        .any(|interior| interior.contains(&center, ()).unwrap_or(false));
                if exterior_contains && !interior_contains {
                    let a = vmap[&vs[0].fix()];
                    let b = vmap[&vs[1].fix()];
                    let c = vmap[&vs[2].fix()];
                    Some([a, b, c])
                } else {
                    None
                }
            }
        })
        .collect_vec();

    // filter out isolated vertices
    let mut remap: HashMap<usize, usize> = HashMap::new();
    let mut vertices = vec![];
    let vs = cdt.vertices().collect_vec();
    let remapped_faces = faces
        .iter()
        .filter_map(|face| {
            face.iter()
                .map(|v| {
                    *remap.entry(*v).or_insert_with(|| {
                        let i = vertices.len();
                        vertices.push(*vs[*v].as_ref());
                        i
                    })
                })
                .collect_array::<3>()
        })
        .collect_vec();

    let mut points = vec![];
    let mut normals = vec![];
    let mut uvs = vec![];
    vertices.iter().for_each(|v| {
        points.push(v.point);
        normals.push(v.normal);
        uvs.push(v.uv);
    });

    Ok(SurfaceTessellation {
        faces: remapped_faces,
        points,
        normals,
        uvs,
    })
}

/// Tessellate the compound curve using an adaptive algorithm recursively
fn tessellate_uv_compound_curve_adaptive<T: FloatingPoint, S: TrimmedSurfaceExt<T, F>, F>(
    curve: &CompoundCurve2D<T>,
    surface: &S,
    angle_tolerance: T,
) -> Vec<Vertex<T>> {
    curve
        .spans()
        .iter()
        .enumerate()
        .flat_map(|(i, span)| {
            let mut vertices = tessellate_uv_curve_adaptive(span, surface, angle_tolerance);
            if i > 0 {
                vertices.remove(0); // Skip the first vertex for spans after the first
            }
            vertices
        })
        .collect_vec()
}

/// Tessellate the curve using an adaptive algorithm recursively
fn tessellate_uv_curve_adaptive<T: FloatingPoint, S: TrimmedSurfaceExt<T, F>, F>(
    curve: &NurbsCurve2D<T>,
    surface: &S,
    angle_tolerance: T,
) -> Vec<Vertex<T>> {
    let degree = curve.degree();
    match degree {
        1 => {
            // if the curve is a linear curve, should start tessellation from the knots
            let knots = curve.knots();
            let n = knots.len();
            let (min, max) = curve.knots_domain();
            // let min = min + T::default_epsilon();
            // let max = max - T::default_epsilon();

            let pts = (1..n - 2)
                .filter_map(|i| {
                    let start = knots[i].clamp(min, max);
                    let end = knots[i + 1].clamp(min, max);
                    if start == end {
                        return None;
                    }
                    let evaluated =
                        iterate_uv_curve_tessellation(curve, surface, start, end, angle_tolerance);
                    #[allow(clippy::iter_skip_zero)]
                    if i == 1 {
                        Some(evaluated.into_iter().skip(0))
                    } else {
                        Some(evaluated.into_iter().skip(1))
                    }
                })
                .flatten();

            pts.map(|uv| {
                let p = surface.point_at(uv.x, uv.y);
                let n = surface.normal_at(uv.x, uv.y);
                Vertex::new(p, n, uv.coords)
            })
            .collect_vec()
        }
        _ => {
            let (start, end) = curve.knots_domain();
            let pts = iterate_uv_curve_tessellation(curve, surface, start, end, angle_tolerance);
            pts.into_iter()
                .map(|uv| {
                    let p = surface.point_at(uv.x, uv.y);
                    let n = surface.normal_at(uv.x, uv.y);
                    Vertex::new(p, n, uv.coords)
                })
                .collect_vec()
        }
    }
}

/// Adaptively tessellate a 2D curve in the surface's UV domain.
///
/// `angle_tolerance` is in **radians** and bounds the angle traversed in each
/// half-segment by both:
///   * the curve's UV-tangent direction (catches curved boundaries on flat
///     surfaces — the surface-normal criterion below is constant there)
///   * the surface normal evaluated at the curve sample (catches highly
///     curved base surfaces)
fn iterate_uv_curve_tessellation<T: FloatingPoint, S: TrimmedSurfaceExt<T, F>, F>(
    curve: &NurbsCurve2D<T>,
    surface: &S,
    start: T,
    end: T,
    angle_tolerance: T,
) -> Vec<Point2<T>> {
    let (u_domain, v_domain) = surface.knots_domain();
    let eps = T::default_epsilon();
    let min = Point2::new(u_domain.0 + eps, v_domain.0 + eps);
    let max = Point2::new(u_domain.1 - eps, v_domain.1 - eps);

    // Clamp end to avoid NaN at exact knot domain boundary
    let curve_domain = curve.knots_domain();
    let safe_start = start.max(curve_domain.0 + eps);
    let safe_end = end.min(curve_domain.1 - eps);

    let (p1, n1) = curve.point_tangent_at(safe_start);
    let delta = safe_end - safe_start;
    if delta < T::from_f64(1e-8).unwrap() {
        return vec![Point2::new(
            p1.x.clamp(min.x, max.x),
            p1.y.clamp(min.y, max.y),
        )];
    }

    let exact_mid = safe_start + (safe_end - safe_start) * T::from_f64(0.5).unwrap();
    let (_, n2) = curve.point_tangent_at(exact_mid);
    let (p3, n3) = curve.point_tangent_at(safe_end);

    let flag = {
        if curve.degree() == 1 {
            // Linear curves are straight in UV; their tangent is constant so
            // there is no angular contribution to detect.
            false
        } else {
            // Closed curves have nn1 ≈ nn3, so the half-segment checks via the
            // midpoint are required to detect them at the top recursion level.
            let a1 = angle_between(&n1, &n2, eps);
            let a2 = angle_between(&n2, &n3, eps);
            a1 > angle_tolerance || a2 > angle_tolerance
        }
    } || {
        // Surface normal across the segment. Closed-curve detection is already
        // handled by the tangent branch above, so a single start-to-end check
        // is sufficient and gives a clean "total normal change ≤ tolerance"
        // semantic for the leaf segment.
        let sn1 = surface.normal_at(p1.x, p1.y);
        let sn3 = surface.normal_at(p3.x, p3.y);
        angle_between(&sn1, &sn3, eps) > angle_tolerance
    };
    if flag {
        let mut left_pts =
            iterate_uv_curve_tessellation(curve, surface, start, exact_mid, angle_tolerance);
        let right_pts =
            iterate_uv_curve_tessellation(curve, surface, exact_mid, end, angle_tolerance);
        left_pts.pop();
        [left_pts, right_pts].concat()
    } else {
        vec![
            Point2::new(p1.x.clamp(min.x, max.x), p1.y.clamp(min.y, max.y)),
            Point2::new(p3.x.clamp(min.x, max.x), p3.y.clamp(min.y, max.y)),
        ]
    }
}

/// Angle (radians) between two vectors of arbitrary dimension. Returns 0 if
/// either vector is below `eps` so callers can treat degenerate samples as
/// "no angular deviation" instead of triggering subdivision.
fn angle_between<T, D>(a: &OVector<T, D>, b: &OVector<T, D>, eps: T) -> T
where
    T: FloatingPoint,
    D: DimName,
    DefaultAllocator: Allocator<D>,
{
    match (a.clone().try_normalize(eps), b.clone().try_normalize(eps)) {
        (Some(na), Some(nb)) => na.dot(&nb).clamp(-T::one(), T::one()).acos(),
        _ => T::zero(),
    }
}