remesh 0.0.5

Isotropic remeshing library
Documentation
// 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::macros::{integrity_assert, integrity_println};
use crate::corner_table::CornerIndex;
use crate::isotropic_remesh::IsotropicRemeshAlgo;
#[allow(unused_imports)]
use rayon::iter::{IndexedParallelIterator, IntoParallelRefMutIterator, ParallelIterator};
use std::fmt::Debug;
use vector_traits::num_traits::AsPrimitive;
use vector_traits::prelude::{GenericVector3, SimdUpgradable};

impl<S, V, const ENABLE_UNSAFE: bool> IsotropicRemeshAlgo<S, V, ENABLE_UNSAFE>
where
    S: crate::common::sealed::ScalarType,
    f64: AsPrimitive<S>,
    V: Debug + Copy + From<[S; 3]> + Into<[S; 3]> + Sync + 'static,
{
    pub(crate) fn smooth_project_vertices(
        &mut self,
        damping_weight: S,
        projection_distance_limit_sq: S,
        projection_coplanar_threshold: S,
    ) -> bool {
        self.vertices_buffer
            .resize(self.vertices.len(), S::Vec3::ZERO);

        //#[cfg(feature = "integrity_check")]
        //self.dump_status("pre smooth_vertices()");

        // 'borrow' vertices_buffer
        let mut vertices_buffer = std::mem::take(&mut self.vertices_buffer);

        vertices_buffer
            .par_iter_mut()
            //.iter_mut()
            .enumerate()
            .map(|(i, v)| (VertexIndex(i as u32), v))
            .for_each(|(v_index, vertex)| {
                let corner = self.corner_table.corner(v_index);
                if !corner.is_valid() {
                    // this vertex was unused.
                    integrity_assert!(!self.vertex_pool.is_used(v_index));
                    //integrity_assert!(self.corner_table.is_deleted(corner.triangle()));
                    return;
                }
                integrity_assert!(self.vertex_pool.is_used(v_index));
                integrity_assert!(!self.corner_table.is_corner_deleted(corner));

                let original_pos = self.vertices[v_index.0 as usize];
                let original_pos_as_simd = original_pos.to_simd();

                // Apply Laplacian smoothing with damping
                if let Some((neighbor_avg, vertex_normal)) = self.check_planarity_and_neighbor_avg(
                    corner,
                    self.params.coplanar_threshold,
                    self.params.min_area_threshold_sq,
                ) {
                    let pos = original_pos_as_simd + (neighbor_avg - original_pos_as_simd) * damping_weight;
                    if let Some((projection_pos, _projection_normal)) = self.bvh.closest_point(S::Vec3::from_simd(pos), vertex_normal, projection_coplanar_threshold) {
                        integrity_assert!(_projection_normal.magnitude_sq() <= S::ONE);
                        integrity_assert!(vertex_normal.magnitude_sq() <= S::ONE);

                        if projection_pos.distance_sq(original_pos) > projection_distance_limit_sq {
                            integrity_println!(
                                "Found bad projection: Old:{:?} smoothed:{:?} projection:{:?} distance:{} limit:{}",
                                original_pos_as_simd.into(),
                                pos.into(),
                                projection_pos.into(),
                                original_pos.distance(projection_pos),
                                projection_distance_limit_sq.sqrt(),
                            );
                            *vertex = original_pos;
                        } else {
                            #[cfg(any(feature = "integrity_check", debug_assertions))]
                            let dot = _projection_normal.to_simd().dot(vertex_normal);
                            integrity_assert!(
                                 dot > projection_coplanar_threshold,
                                "Found bad projection (normal): Old:{:?} smoothed:{:?} projection:{:?} distance:{} limit:{} dot:{dot} limit:{projection_coplanar_threshold}",
                                original_pos_as_simd.into(),
                                pos.into(),
                                projection_pos.into(),
                                original_pos.distance(projection_pos),
                                projection_distance_limit_sq.sqrt(),
                            );
                            *vertex = projection_pos;
                            integrity_println!(
                                "Old:{:?} smoothed:{:?} projection:{:?} distance:{}",
                                original_pos_as_simd.into(),
                                pos.into(),
                                projection_pos.into(),
                                original_pos.distance(projection_pos)
                            );
                        }
                    } else {
                        integrity_println!(
                            "Old:{:?} smoothed:{:?} found no projection",
                            original_pos_as_simd.into(),
                            pos.into()
                        );
                        *vertex = S::Vec3::from_simd(original_pos_as_simd);
                    }
                } else {
                    integrity_println!("Old:{:?} bad coplanarity", original_pos_as_simd.into(),);
                    *vertex = original_pos;
                }
            });
        // put back vertices_buffer
        self.vertices_buffer = vertices_buffer;
        std::mem::swap(&mut self.vertices, &mut self.vertices_buffer);
        // smoothing always does something, so we can't use that as a break condition
        false
    }

    fn check_planarity_and_neighbor_avg(
        &self,
        start_corner: CornerIndex,
        planarity_threshold: S,
        min_area_threshold_sq: S,
    ) -> Option<(S::Vec3Simd, S::Vec3Simd)> {
        #[cfg(feature = "integrity_check")]
        let vertex_fan = self.corner_table.ccw_vertex_fan(start_corner);
        integrity_println!(
            "start_corner:{} surrounding vertices: {}",
            self.dbg_corner(start_corner),
            self.dbg_fan_prev(&vertex_fan)
        );

        let c0 = start_corner;
        let (c1, c2) = self.corner_table.next_prev(c0);

        let p0 = self.vertex_of_corner(c0).to_simd();
        let p1 = self.vertex_of_corner(c1).to_simd();
        let p2 = self.vertex_of_corner(c2).to_simd();

        integrity_println!(
            "checking vertex p0:{} {:?}",
            self.dbg_corner(c0),
            self.vertex_of_corner(c0).into()
        );
        integrity_println!(
            "checking vertex p1:{} {:?}",
            self.dbg_corner(c1),
            self.vertex_of_corner(c1).into()
        );
        integrity_println!(
            "checking vertex p2:{} {:?}",
            self.dbg_corner(c2),
            self.vertex_of_corner(c2).into()
        );

        let mut neighbor_sum = p1 + p2;
        let mut neighbor_count = 1;

        // Initial edge vectors
        let mut edge1 = p1 - p0;
        let mut edge2 = p2 - p0;

        let reference_cross = edge1.cross(edge2);
        let reference_area_sq = reference_cross.magnitude_sq();

        // Check if starting triangle is degenerate
        if reference_area_sq < min_area_threshold_sq {
            integrity_println!(
                "bail: degenerate reference triangle, area² = {:?}",
                reference_area_sq
            );
            return None;
        }

        let reference_normal = reference_cross.normalize();

        let mut corner = self.corner_table.swing_ccw(start_corner);

        while corner != start_corner {
            // Reuse: edge2 becomes edge1 for next triangle
            edge1 = edge2;

            // Fetch only the new vertex
            let corner_prev = self.corner_table.prev(corner);
            let p_prev = self.vertex_of_corner(corner_prev).to_simd();

            integrity_println!(
                "checking vertex prev {} {:?}",
                self.dbg_corner(corner_prev),
                p_prev.into()
            );

            neighbor_sum += p_prev;
            neighbor_count += 1;

            edge2 = p_prev - p0;

            // Check normal
            let new_cross = edge1.cross(edge2);
            let new_area_sq = new_cross.magnitude_sq();

            // Check for degenerate triangle
            if new_area_sq < min_area_threshold_sq {
                integrity_println!(
                    "bail: degenerate triangle at corner {}, area² = {:?}",
                    self.dbg_corner(corner),
                    new_area_sq
                );
                return None;
            }

            let new_normal = new_cross.normalize();

            if reference_normal.dot(new_normal) < planarity_threshold {
                integrity_println!(
                    "bail: not coplanar: {:?} >= {}",
                    reference_normal.dot(new_normal),
                    planarity_threshold
                );
                return None; // Early exit!
            }

            corner = self.corner_table.swing_ccw(corner);
        }

        #[cfg(feature = "integrity_check")]
        assert_eq!(neighbor_count, vertex_fan.len());
        integrity_println!(
            "smoothing: average = {:?}",
            (neighbor_sum / (neighbor_count as f64).as_()).into()
        );

        Some((
            neighbor_sum / ((neighbor_count as f64).as_()),
            reference_normal,
        ))
    }
}