remesh 0.0.5

// SPDX-License-Identifier: MIT OR Apache-2.0
// Copyright (c) 2025 lacklustr@protonmail.com https://github.com/eadf

use crate::common::VertexIndex;
use crate::common::remesh_error::RemeshError;
use crate::common::sealed::{IndexType, ScalarType};
use std::fmt::Debug;
use std::marker::PhantomData;
use vector_traits::num_traits::{AsPrimitive, Float};
use vector_traits::prelude::GenericVector3;

#[derive(Debug, Clone)]
pub enum CollapseStrategy {
    /// Don't collapse edges
    Disabled,
    /// Use dihedral angle as quality measurement
    DihedralAngle,
    /// Use displacement as quality measurement
    Displacement,
    /// Use qem as quality measurement
    Qem,
}

#[derive(Debug, Clone)]
pub enum SplitStrategy {
    /// Don't collapse edges
    Disabled,
    /// Use dihedral angle as quality measurement
    DihedralAngle,
    /// Use displacement as quality measurement
    Displacement,
}

#[derive(Clone)]
pub enum FlipStrategy<S: ScalarType> {
    /// Don't flip edges
    Disabled,
    /// Only flip to improve valence distribution
    Valence,
    /// Prioritize valence improvement, then weight by triangle aspect ratio quality
    WeightedQuality { quality_threshold: S },
}

impl<S: ScalarType + Default> FlipStrategy<S>
where
    f64: AsPrimitive<S>,
{
    /// Create a new `FlipStrategy::Valence`
    pub fn default_valence() -> Self {
        FlipStrategy::Valence
    }

    /// Create a new `FlipStrategy::Weighted` with default values
    pub fn default_quality() -> Self {
        FlipStrategy::WeightedQuality {
            quality_threshold: DEFAULT_EDGE_FLIP_QUALITY_THRESHOLD.as_(),
        }
    }

    /// Create a new `FlipStrategy::WeightedQuality` with custom threshold
    pub fn quality(quality_threshold: S) -> Self {
        FlipStrategy::WeightedQuality { quality_threshold }
    }
}

#[derive(Debug, Clone)]
#[allow(dead_code)]
/// Remeshing parameters
pub(crate) struct RemeshParams<S: ScalarType> {
    /// Target edge length squared
    pub(super) target_edge_length: S,
    pub(super) min_area_threshold_sq: S,
    /// Number of iterations
    pub(super) iterations: u32,
    /// Split edges longer than this factor * `target_edge_length`
    split_multiplier: S,
    pub(super) split_strategy: SplitStrategy,
    /// Collapse edges shorter than this factor * `target_edge_length`
    collapse_multiplier: S,
    pub(super) collapse_strategy: CollapseStrategy,
    pub(super) flip_strategy: FlipStrategy<S>,
    pub(super) corner_table_defragmentation_ratio: Option<S>,
    pub(super) smooth_weight: Option<S>,
    pub(super) smooth_normal_threshold: S,
    /// Cosine threshold for coplanarity: dot product of normals above this value
    /// indicates triangles are nearly coplanar (parallel, same direction). Range: `[0.5..1]`
    pub(super) coplanar_threshold: S,
    /// coplanar_threshold^2
    pub(super) coplanar_threshold_sq: S,
    /// Cosine threshold for inversion: dot product of normals below this value
    /// indicates triangles are inverted (parallel, opposite directions). Range: `[-1..-0.5]`
    pub(super) crease_limit_threshold: S,
    /// crease_limit_threshold^2
    pub(super) crease_limit_threshold_sq: S,
    /// The acceptable normal dot product (negative) between two triangles already in the mesh.
    /// Below this value the triangles are considered inverted.
    pub(super) inversion_validation_threshold: S,
    /// `inversion_validation_threshold^2`
    pub(super) inversion_validation_threshold_sq: S,
    pub(super) max_valence: i16,
    pub(super) collapse_threshold_sq: S,
    pub(super) collapse_qem_threshold: S,
    pub(super) collapse_qem_threshold_sq: f64,
    /// the largest squared distance a qem solution is allowed to land at (from either endpoint)
    pub(super) collapse_qem_max_dist_sq: f64,
    /// The smallest squared distance a qem solution is allowed to land at (from either endpoint)
    /// If smaller than this the endpoint will be used instead
    pub(super) collapse_qem_min_dist_sq: f64,
    pub(super) split_threshold_sq: S,
    pub(super) max_projection_distance_sq: S,
    pub(super) fix_non_manifold: bool,
    /// The vertex_limit of the `print_input_status()` method
    pub(super) print_stats: Option<usize>,
}

const DEFAULT_EDGE_LENGTH: f64 = 1.0;
const DEFAULT_SPLIT_MULTIPLIER: f64 = 1.4;
const DEFAULT_COLLAPSE_MULTIPLIER: f64 = 0.7;
const DEFAULT_SMOOTH_WEIGHT: f64 = 0.1;
const DEFAULT_EDGE_FLIP_QUALITY_THRESHOLD: f64 = 1.1;
const DEFAULT_SMOOTH_NORMAL_THRESHOLD: f64 = 0.90;
const DEFAULT_CORNER_DEFRAGMENTATION_LIMIT: f64 = 0.4;
const DEFAULT_COPLANAR_THRESHOLD: f64 = 0.9961946980917455; // cos(5°)
const MIN_COPLANAR_THRESHOLD: f64 = 0.5; // cos(60°)
const MAX_COPLANAR_THRESHOLD: f64 = 1.0; // cos(0°)
const DEFAULT_COLLAPSE_INVERSION_THRESHOLD: f64 = -0.984807753012208; // cos(-170°)
const MIN_COLLAPSE_INVERSION_THRESHOLD: f64 = -0.1736481776669303; // cos(100)
const MAX_COLLAPSE_INVERSION_THRESHOLD: f64 = -1.0; // cos(180)
const DEFAULT_VALIDATION_INVERSION_THRESHOLD: f64 = -0.999_390_827_019_095_8; // cos(-178°)
const DEFAULT_COLLAPSE_QEM_THRESHOLD: f64 = 0.05;
const DEFAULT_ITERATIONS: u32 = 10;
const MIN_ITERATIONS: u32 = 1;
const MAX_ITERATIONS: u32 = 1000;
const MIN_TARGET_EDGE_LENGTH: f64 = 0.00001;

use crate::corner_table::DEFAULT_MAX_VALENCE;
use crate::isotropic_remesh::IsotropicRemeshAlgo;
use crate::prelude::IsotropicRemesh;

impl<S> Default for RemeshParams<S>
where
    S: ScalarType,
    f64: AsPrimitive<S>,
{
    fn default() -> Self {
        let target_edge_length = DEFAULT_EDGE_LENGTH.as_();
        let target_edge_length_sq = Float::powi(target_edge_length, 2);
        let coplanar_threshold: S = DEFAULT_COPLANAR_THRESHOLD.as_();
        let inversion_collapse_threshold: S = DEFAULT_COLLAPSE_INVERSION_THRESHOLD.as_();
        let inversion_validation_threshold: S = DEFAULT_VALIDATION_INVERSION_THRESHOLD.as_();
        let min_area_threshold_sq = target_edge_length_sq * Float::powi(0.001.as_(), 2);
        let max_projection_distance_sq = target_edge_length_sq * Float::powi(0.5.as_(), 2);
        let collapse_threshold_sq =
            Float::powi(target_edge_length * DEFAULT_COLLAPSE_MULTIPLIER.as_(), 2);
        Self {
            target_edge_length,
            min_area_threshold_sq,
            iterations: DEFAULT_ITERATIONS,
            split_multiplier: DEFAULT_SPLIT_MULTIPLIER.as_(),
            collapse_multiplier: DEFAULT_COLLAPSE_MULTIPLIER.as_(),
            smooth_weight: None,
            smooth_normal_threshold: DEFAULT_SMOOTH_NORMAL_THRESHOLD.as_(),
            flip_strategy: FlipStrategy::Disabled,
            collapse_strategy: CollapseStrategy::DihedralAngle,
            split_strategy: SplitStrategy::DihedralAngle,
            coplanar_threshold,
            coplanar_threshold_sq: Float::powi(coplanar_threshold, 2),
            crease_limit_threshold: inversion_collapse_threshold,
            crease_limit_threshold_sq: Float::powi(inversion_collapse_threshold, 2),
            inversion_validation_threshold,
            inversion_validation_threshold_sq: Float::powi(inversion_validation_threshold, 2),
            max_valence: DEFAULT_MAX_VALENCE,
            corner_table_defragmentation_ratio: Some(DEFAULT_CORNER_DEFRAGMENTATION_LIMIT.as_()),
            split_threshold_sq: Float::powi(target_edge_length * DEFAULT_SPLIT_MULTIPLIER.as_(), 2),
            collapse_threshold_sq,
            collapse_qem_threshold: DEFAULT_COLLAPSE_QEM_THRESHOLD.as_(),
            collapse_qem_threshold_sq: Float::powi(
                DEFAULT_EDGE_LENGTH * DEFAULT_COLLAPSE_MULTIPLIER * DEFAULT_COLLAPSE_QEM_THRESHOLD,
                2,
            ),
            collapse_qem_max_dist_sq: (collapse_threshold_sq * S::TWO).into(),
            collapse_qem_min_dist_sq: collapse_threshold_sq.into() * 0.0001, // 0.1% squared
            max_projection_distance_sq,
            fix_non_manifold: false,
            print_stats: None,
        }
    }
}

impl<S> RemeshParams<S>
where
    S: ScalarType,
    f64: AsPrimitive<S>,
{
    pub(crate) fn validate(&mut self) -> Result<(), RemeshError> {
        if self.split_multiplier < 0.1.as_() {
            Err(RemeshError(
                "The split_multiplier was too low: {split_multiplier}".into(),
            ))?;
        }
        self.split_threshold_sq = Float::powi(self.target_edge_length * self.split_multiplier, 2);

        if self.collapse_multiplier < 0.1.as_() {
            Err(RemeshError(
                "The collapse_multiplier was too low: {collapse_multiplier}".into(),
            ))?;
        }
        self.collapse_threshold_sq =
            Float::powi(self.target_edge_length * self.collapse_multiplier, 2);

        let collapse_threshold = self.collapse_threshold_sq.sqrt();
        let split_threshold = self.split_threshold_sq.sqrt();
        if collapse_threshold >= split_threshold {
            Err(RemeshError(format!(
                "The split and collapse length parameters overlap collapse:{collapse_threshold} split:{split_threshold}",
            )))?;
        }
        self.collapse_qem_threshold_sq =
            Float::powi(self.collapse_threshold_sq * self.collapse_qem_threshold, 2).as_();

        Ok(())
    }

    pub(crate) fn print_essentials(&self) {
        println!("### RemeshParams ###");
        println!("iterations:{}", self.iterations);
        println!("target_edge_length:{:?}", self.target_edge_length);
        println!("split_multiplier:{:?}", self.split_multiplier);
        println!("collapse_multiplier:{:?}", self.collapse_multiplier);
        match self.collapse_strategy {
            CollapseStrategy::Qem => println!(
                "collapse_edges mode:{:?} {}% max_dist_sq:{}",
                self.collapse_strategy,
                self.collapse_qem_threshold * 100.0.into(),
                self.collapse_qem_max_dist_sq
            ),
            _ => println!("collapse_edges mode:{:?}", self.collapse_strategy),
        }
        println!("flip_edges mode:{:?}", self.flip_strategy);
        println!("smooth_weight:{:?}", self.smooth_weight);
        println!(
            "coplanar_threshold(cos):{:.4?}-> accept triangle modifications with normal diff below or equal {:.4?}°",
            self.coplanar_threshold,
            Float::to_degrees(self.coplanar_threshold.acos())
        );
        println!(
            "inversion_collapse_threshold(cos):{:.4?}-> accept adjacent triangles with normals that stay above {:.4?}° diff",
            self.crease_limit_threshold,
            -Float::to_degrees(self.crease_limit_threshold.acos())
        );
        println!(
            "inversion_validation_threshold(cos):{:.4?}-> accept adjacent triangles, already in the mesh, with normals that stay above {:.4?}° diff",
            self.inversion_validation_threshold,
            -Float::to_degrees(self.inversion_validation_threshold.acos())
        );
        println!("Fix non-manifold mesh: {}", self.fix_non_manifold);
        println!();
    }

    /// Set target edge length
    pub fn with_target_edge_length(&mut self, target_edge_length: S) -> Result<(), RemeshError> {
        if target_edge_length < MIN_TARGET_EDGE_LENGTH.as_() {
            Err(RemeshError(format!(
                "The target_edge_length was too small: {target_edge_length} limit:{}",
                MIN_TARGET_EDGE_LENGTH
            )))?;
        }
        self.target_edge_length = target_edge_length;
        Ok(())
    }
}

impl<S, V, const ENABLE_UNSAFE: bool> IsotropicRemesh<S, V, ENABLE_UNSAFE>
where
    S: ScalarType,
    f64: AsPrimitive<S>,
    V: Debug + Copy + From<[S; 3]> + Into<[S; 3]> + Sync + 'static,
{
    /// Create a new isotropic remesher from a triangle mesh.
    ///
    /// Takes a vertex and index iterators describing a triangulated mesh.
    /// The indices are interpreted as consecutive triplets, where each group of
    /// three indices defines one triangle face.
    ///
    /// # Parameters
    ///
    /// - `vertices`: Iterator yielding vertex position references in any format convertible to/from `[S; 3]`
    /// - `indices`: Iterator yielding triangle index references, interpreted as groups of 3 forming triangles
    ///
    /// # Supported Index Types
    ///
    /// The index type can be any of:
    /// - `usize`, `u64` : Runtime range-checked conversion
    /// - `u32`, `u16`: Zero-cost conversion (recommended)
    /// - `i32`: Runtime range-checked conversion (convenient for small examples)
    ///
    /// **Note**: The output indices from `run()` are currently returned as `Vec<u32>`,
    /// regardless of input type.
    ///
    /// # Vertex Type Flexibility
    ///
    /// The vertex type `V` can be any type that implements `Copy + From<[S; 3]> + Into<[S; 3]>`:
    /// - Array notation: `[f32; 3]`, `[f64; 3]`
    /// - Math libraries: `glam::Vec3`, `nalgebra::Vector3<f64>`
    /// - Custom types: Implement the `From`/`Into` traits for `[S; 3]`
    ///
    ///
    /// # Errors
    ///
    /// Returns `RemeshError` if:
    /// - The index iterator length is not divisible by 3
    /// - Any index is out of bounds for the vertex buffer
    /// - Any vertex contains non-finite floats
    /// - Index conversion fails (e.g., negative `i32` values)
    ///
    /// # Examples
    ///
    /// ```rust,ignore
    /// // With f32 arrays and u32 indices
    /// let vertices = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
    /// let indices: Vec<u32> = vec![0, 1, 2];
    /// let remesher = IsotropicRemesh::<f32, _>::new(&vertices, &indices)?;
    ///
    /// // With glam and usize indices
    /// let vertices = vec![
    ///     glam::Vec3::ZERO,
    ///     glam::Vec3::X,
    ///     glam::Vec3::Y,
    /// ];
    /// let indices = vec![0_usize, 1, 2];
    /// let remesher = IsotropicRemesh::<f32, _>::new(&vertices, &indices)?;
    ///
    /// // With i32 indices (convenient for examples)
    /// let indices = vec![0, 1, 2];
    /// let remesher = IsotropicRemesh::<f32, _>::new(&vertices, &indices)?;
    /// ```
    pub fn new<'a, VI, II, I>(vertices: VI, indices: II) -> Result<Self, RemeshError>
    where
        VI: IntoIterator<Item = &'a V>,
        II: IntoIterator<Item = &'a I>,
        V: 'a,
        I: IndexType + 'a,
        I::Error: Debug,
    {
        let vertices: Vec<_> = vertices
            .into_iter()
            .map(|v| {
                let vertex = Self::to_glam(*v);
                if vertex.is_finite() {
                    Ok(vertex)
                } else {
                    Err(RemeshError(format!(
                        "Non-finite vertex detected: {:?}",
                        vertex
                    )))
                }
            })
            .collect::<Result<Vec<_>, RemeshError>>()?;

        let mut max_index = 0;
        let indices: Vec<_> = indices
            .into_iter()
            .map(|&i| {
                TryInto::<u32>::try_into(i)
                    .map_err(|e| RemeshError(format!("{:?}", e)))
                    .and_then(|i| {
                        max_index = std::cmp::max(max_index, i);
                        if i < vertices
                            .len()
                            .try_into()
                            .map_err(|e| RemeshError(format!("{:?}", e)))?
                        {
                            Ok(VertexIndex(i))
                        } else {
                            Err(RemeshError(format!("Index out of range {:?}", i)))
                        }
                    })
            })
            .collect::<Result<Vec<_>, _>>()?;

        // Check triangle count after collection
        if !indices.len().is_multiple_of(3) {
            return Err(RemeshError(
                "Indices must represent triangles. (was not multiple of 3)".to_string(),
            ));
        }

        Ok(Self {
            vertices,
            indices,
            max_index,
            params: RemeshParams::default(),
            _pd: PhantomData,
        })
    }

    /// Set target edge length
    pub fn with_target_edge_length(mut self, target_edge_length: S) -> Result<Self, RemeshError> {
        self.params.with_target_edge_length(target_edge_length)?;
        Ok(self)
    }

    /// Sets the policy of the flip_edges phase
    ///
    /// **Default**: Disabled (if not called, equivalent to passing `FlipStrategy::Disabled`)
    pub fn with_flip_edges(mut self, flip_strategy: FlipStrategy<S>) -> Result<Self, RemeshError> {
        if let FlipStrategy::WeightedQuality { quality_threshold } = flip_strategy {
            if quality_threshold > 1.5.as_() || quality_threshold < 1.0.as_() {
                return Err(RemeshError(
                    "FlipStrategy::WeightedQuality should be within [1.0..1.5] range".into(),
                ));
            }
        }

        self.params.flip_strategy = flip_strategy;
        Ok(self)
    }

    pub fn with_default_flip_edges(self) -> Result<Self, RemeshError> {
        self.with_flip_edges(FlipStrategy::default_quality())
    }

    pub fn without_flip_edges(self) -> Result<Self, RemeshError> {
        self.with_flip_edges(FlipStrategy::Disabled)
    }

    /// Set smooth weight.
    ///
    /// Controls the smoothing operation.
    ///
    pub fn with_smooth_weight(mut self, weight: S) -> Result<Self, RemeshError> {
        self.params.smooth_weight = Some(weight);
        Ok(self)
    }

    pub fn with_default_smooth_weight(self) -> Result<Self, RemeshError> {
        self.with_smooth_weight(DEFAULT_SMOOTH_WEIGHT.as_())
    }

    pub fn without_smooth_weight(mut self) -> Result<Self, RemeshError> {
        self.params.smooth_weight = None;
        Ok(self)
    }

    /// Configure edge splitting behavior.
    ///
    /// Controls the maximum edge length threshold before edges are subdivided.
    /// Edges longer than `target_edge_length * split_multiplier` may be split.
    ///
    pub fn with_split_multiplier(mut self, split_multiplier: S) -> Result<Self, RemeshError> {
        self.params.split_multiplier = split_multiplier;
        Ok(self)
    }

    pub fn with_default_split_multiplier(self) -> Result<Self, RemeshError> {
        self.with_split_multiplier(DEFAULT_SPLIT_MULTIPLIER.as_())
    }

    pub fn with_split_edges(mut self, strategy: SplitStrategy) -> Result<Self, RemeshError> {
        self.params.split_strategy = strategy;
        Ok(self)
    }

    pub fn without_split_edges(mut self) -> Result<Self, RemeshError> {
        self.params.split_strategy = SplitStrategy::Disabled;
        Ok(self)
    }

    /// Configure edge collapsing behavior.
    ///
    /// Controls the minimum edge length threshold before edges are merged.
    /// Edges shorter than `target_edge_length * collapse_multiplier` may be collapsed.
    ///
    pub fn with_collapse_multiplier(mut self, collapse_threshold: S) -> Result<Self, RemeshError> {
        self.params.collapse_multiplier = collapse_threshold;
        Ok(self)
    }

    pub fn with_default_collapse_multiplier(self) -> Result<Self, RemeshError> {
        self.with_collapse_multiplier(DEFAULT_COLLAPSE_MULTIPLIER.as_())
    }

    /// Selects the collapse variant to use, will use the default parameters of each variant.
    pub fn with_collapse_edges(mut self, strategy: CollapseStrategy) -> Result<Self, RemeshError> {
        self.params.collapse_strategy = strategy;
        Ok(self)
    }

    pub fn with_collapse_qem_threshold(mut self, qem_threshold: S) -> Result<Self, RemeshError> {
        self.params.collapse_strategy = CollapseStrategy::Qem;
        self.params.collapse_qem_threshold = qem_threshold;
        Ok(self)
    }

    pub fn with_default_collapse_edges(self) -> Result<Self, RemeshError> {
        self.with_collapse_edges(CollapseStrategy::DihedralAngle)
    }

    pub fn without_collapse_edges(self) -> Result<Self, RemeshError> {
        self.with_collapse_edges(CollapseStrategy::Disabled)
    }

    /// Configure the coplanarity threshold for feature preservation using cosine directly.
    ///
    /// This is a lower-level alternative to `with_coplanar_angle_threshold` that uses
    /// the cosine of the dihedral angle directly. Higher values preserve sharper features.
    ///
    /// Most users should prefer `with_coplanar_angle_threshold` for better ergonomics.
    ///
    /// # Arguments
    /// * `coplanar_threshold_cos` - Cosine of the maximum dihedral angle (range: -1.0 to 1.0)
    ///   - `cos(15°) ≈ 0.966` - Sharp feature preservation
    ///   - `cos(30°) ≈ 0.866` - Moderate preservation
    ///   - `cos(45°) ≈ 0.707` - Soft preservation
    pub fn with_coplanar_threshold(
        mut self,
        coplanar_threshold_cos: S,
    ) -> Result<Self, RemeshError> {
        if coplanar_threshold_cos < MIN_COPLANAR_THRESHOLD.as_() {
            let min_c = MIN_COPLANAR_THRESHOLD.as_();
            return Err(RemeshError(format!(
                "Coplanar threshold too low: {coplanar_threshold_cos}({}°). Minimum is {min_c}({}°)",
                Float::to_degrees(coplanar_threshold_cos.acos()),
                Float::to_degrees(min_c.acos()),
            )));
        }

        if coplanar_threshold_cos > MAX_COPLANAR_THRESHOLD.as_() {
            let max_c = MAX_COPLANAR_THRESHOLD.as_();
            return Err(RemeshError(format!(
                "Coplanar threshold too high: {coplanar_threshold_cos}({}°). Maximum is {max_c}({}°)",
                Float::to_degrees(coplanar_threshold_cos.acos()),
                Float::to_degrees(max_c.acos()),
            )));
        }

        self.params.coplanar_threshold = coplanar_threshold_cos;
        self.params.coplanar_threshold_sq = Float::powi(coplanar_threshold_cos, 2);
        Ok(self)
    }

    /// Configure the coplanarity threshold for feature preservation using an angle in degrees.
    ///
    /// Controls the maximum dihedral angle between adjacent triangles that are considered
    /// coplanar. Smaller angles preserve sharper features, while larger angles allow more
    /// aggressive remeshing across surface discontinuities.
    ///
    /// # Arguments
    /// * `angle_degrees` - Maximum angle in degrees (0° = perfectly coplanar, 180° = opposite facing)
    ///
    /// # Recommended values
    /// * `5-15°` - Preserve sharp features and hard edges
    /// * `15-30°` - Balance between feature preservation and smooth remeshing
    /// * `30-45°` - Allow more aggressive smoothing, softer feature preservation
    ///
    pub fn with_coplanar_angle_threshold(self, angle_degrees: S) -> Result<Self, RemeshError> {
        if angle_degrees < S::zero() || angle_degrees > 90.0.into() {
            return Err(RemeshError(format!(
                "Coplanar angle threshold must be between 0° and 90° {angle_degrees}°"
            )));
        }
        self.with_coplanar_threshold(Float::to_radians(angle_degrees).cos())
    }

    pub fn with_default_coplanar_threshold(self) -> Result<Self, RemeshError> {
        self.with_coplanar_threshold(DEFAULT_COPLANAR_THRESHOLD.as_())
    }

    pub fn with_crease_threshold(mut self, crease_threshold_cos: S) -> Result<Self, RemeshError> {
        if crease_threshold_cos < MAX_COLLAPSE_INVERSION_THRESHOLD.as_() {
            let min_c = MAX_COLLAPSE_INVERSION_THRESHOLD.as_();
            return Err(RemeshError(format!(
                "Crease threshold too low: {crease_threshold_cos}({}°). Minimum is {min_c}({}°)",
                Float::to_degrees(crease_threshold_cos.acos()),
                Float::to_degrees(min_c.acos()),
            )));
        }

        if crease_threshold_cos > MIN_COLLAPSE_INVERSION_THRESHOLD.as_() {
            let max_c = MIN_COLLAPSE_INVERSION_THRESHOLD.as_();
            return Err(RemeshError(format!(
                "Crease threshold too high: {crease_threshold_cos}({}°). Maximum is {max_c}({}°)",
                Float::to_degrees(crease_threshold_cos.acos()),
                Float::to_degrees(max_c.acos()),
            )));
        }

        self.params.crease_limit_threshold = crease_threshold_cos;
        self.params.crease_limit_threshold_sq = Float::powi(crease_threshold_cos, 2);
        Ok(self)
    }

    /// Configure the sharp crease threshold for feature preservation using an angle in degrees.
    ///
    pub fn with_crease_angle_threshold(self, angle_degrees: S) -> Result<Self, RemeshError> {
        let angle_degrees = Float::abs(angle_degrees);
        if angle_degrees < 100.0.into() || angle_degrees > 180.0.into() {
            return Err(RemeshError(format!(
                "Crease angle threshold must be between 100° and 180° : {angle_degrees}"
            )));
        }
        self.with_crease_threshold(Float::to_radians(angle_degrees).cos())
    }

    pub fn with_default_crease_threshold(self) -> Result<Self, RemeshError> {
        self.with_crease_threshold(DEFAULT_COLLAPSE_INVERSION_THRESHOLD.as_())
    }

    /// Enable the 'fix non-manifold mesh' functionality
    pub fn with_fix_non_manifold(mut self) -> Result<Self, RemeshError> {
        self.params.fix_non_manifold = true;
        Ok(self)
    }

    pub fn with_print_stats(mut self, print_stats_limit: usize) -> Result<Self, RemeshError> {
        self.params.print_stats = Some(print_stats_limit);
        Ok(self)
    }
}

impl<S, V, const ENABLE_UNSAFE: bool> IsotropicRemeshAlgo<S, V, ENABLE_UNSAFE>
where
    S: ScalarType,
    f64: AsPrimitive<S>,
    V: Debug + Copy + From<[S; 3]> + Into<[S; 3]> + Sync + 'static,
{
    /// Set maximum number of iterations. If the previous iteration didn't produce any changes
    /// the iteration is terminated.
    pub(super) fn check_iterations<I>(&mut self, raw_iterations: I) -> Result<u32, RemeshError>
    where
        I: TryInto<u32>,
        <I as TryInto<u32>>::Error: Debug,
    {
        let iterations = raw_iterations
            .try_into()
            .map_err(|e| RemeshError(format!("Invalid iterations value: {:?}", e)))?;

        if iterations < MIN_ITERATIONS {
            return Err(RemeshError(format!(
                "Iterations too low: {}. Minimum is {}",
                iterations, MIN_ITERATIONS
            )));
        }

        if iterations > MAX_ITERATIONS {
            return Err(RemeshError(format!(
                "Iterations too high: {}. Maximum is {}",
                iterations, MAX_ITERATIONS
            )));
        }

        Ok(iterations)
    }

    #[cfg(test)]
    /// Set target edge length
    pub(crate) fn with_target_edge_length(
        &mut self,
        target_edge_length: S,
    ) -> Result<(), RemeshError> {
        self.params.with_target_edge_length(target_edge_length)?;
        Ok(())
    }
}