Fix/visibility (#30)

* test: add pan-handle test

* refactor: add triangulation trait

* feat: add triangulation to dev cli

* refactor: replace visibility algorithm

* test: update visibility tests

* test: enable rdp tests

* refactor: add type conversions from spade to geo

* refactor: separate extension traits from types

* fix: failing split across index 0

* fix: underflow for null polygon

* chore: clean-up

* chore: add license

* chore: formatting

Author: NiallOswald

src/algorithms/simplify_charshape.rs

@@ -18,10 +18,11 @@
 
 // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan
 
-use geo::{Coord, CoordsIter, GeoFloat, LineString, Polygon};
+use crate::extensions::triangulate::Triangulate;
+use geo::{Coord, GeoFloat, LineString, Polygon};
 use hashbrown::HashSet;
 use spade::handles::{DirectedEdgeHandle, VertexHandle};
-use spade::{CdtEdge, ConstrainedDelaunayTriangulation, Point2, SpadeNum, Triangulation};
+use spade::{CdtEdge, Point2, SpadeNum, Triangulation};
 use std::cmp::Ordering;
 use std::collections::BinaryHeap;
 use std::hash::Hash;
@@ -98,28 +99,7 @@ where
 
     let eps_2 = eps * eps;
 
-    // Construct Delaunay triangulation
-    let num_vertices = orig.exterior().0.len() - 1;
-
-    let vertices = orig
-        .exterior_coords_iter()
-        .take(num_vertices) // duplicate points are removed
-        .map(|c| Point2::new(c.x, c.y))
-        .collect::<Vec<_>>();
-
-    let edges = (0..num_vertices)
-        .map(|i| {
-            if i == 0 {
-                [vertices.len() - 1, i]
-            } else {
-                [i - 1, i]
-            }
-        })
-        .collect::<Vec<_>>();
-
-    let tri =
-        ConstrainedDelaunayTriangulation::<Point2<T>>::bulk_load_cdt(vertices, edges).unwrap();
-
+    let tri = orig.triangulate();
     let boundary_edges = tri.convex_hull().map(|edge| edge.rev()).collect::<Vec<_>>();
     let mut boundary_nodes: HashSet<_> =
         HashSet::from_iter(boundary_edges.iter().map(|&edge| BoundaryNode(edge.from())));
@@ -169,7 +149,6 @@ fn recompute_boundary<'a, T>(
 ) where
     T: GeoFloat + SpadeNum,
 {
-    //
     let choices = [edge.prev(), edge.next()];
     for new_edge in choices {
         let e = CharScore {

src/algorithms/simplify_rdp.rs

@@ -18,96 +18,139 @@
 
 // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan
 
-use crate::algorithms::visibility::visiblity_polygon;
-use crate::extensions::segments::{FromSegments, HullSegments};
-use geo::{Coord, Distance, Euclidean, GeoFloat, Line, LineString, Polygon};
-use rayon::prelude::*;
+use crate::algorithms::visibility::visibility_intersection;
+use crate::extensions::conversions::IntoCoord;
+use crate::extensions::segments::FromSegments;
+use crate::extensions::triangulate::Triangulate;
+use geo::{Distance, Euclidean, GeoFloat, Line, LineString, Polygon};
+use spade::handles::{FixedVertexHandle, VertexHandle};
+use spade::{CdtEdge, ConstrainedDelaunayTriangulation, Point2, SpadeNum, Triangulation};
+
+struct CircularIterator<'a, T: SpadeNum> {
+    current: VertexHandle<'a, Point2<T>, (), CdtEdge<()>>,
+    until: VertexHandle<'a, Point2<T>, (), CdtEdge<()>>,
+    cdt: &'a ConstrainedDelaunayTriangulation<Point2<T>>,
+}
+
+impl<'a, T: SpadeNum> CircularIterator<'a, T> {
+    fn new(
+        from: VertexHandle<'a, Point2<T>, (), CdtEdge<()>>,
+        until: VertexHandle<'a, Point2<T>, (), CdtEdge<()>>,
+        cdt: &'a ConstrainedDelaunayTriangulation<Point2<T>>,
+    ) -> Self {
+        CircularIterator {
+            current: from,
+            until,
+            cdt,
+        }
+    }
+}
+
+impl<'a, T: SpadeNum> Iterator for CircularIterator<'a, T> {
+    type Item = VertexHandle<'a, Point2<T>, (), CdtEdge<()>>;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        if self.current == self.until {
+            return None;
+        }
+
+        let vertex = self.current;
+
+        self.current = {
+            let handle =
+                FixedVertexHandle::from_index((self.current.index() + 1) % self.cdt.num_vertices());
+            self.cdt.get_vertex(handle).unwrap()
+        };
 
-fn rdp_preserve<T>(ls: &[Coord<T>], eps: T) -> Vec<Coord<T>>
+        Some(vertex)
+    }
+}
+
+fn rdp_preserve<T>(
+    from: VertexHandle<'_, Point2<T>, (), CdtEdge<()>>,
+    to: VertexHandle<'_, Point2<T>, (), CdtEdge<()>>,
+    cdt: &ConstrainedDelaunayTriangulation<Point2<T>>,
+    eps: T,
+) -> Vec<Point2<T>>
 where
-    T: GeoFloat + Send + Sync,
+    T: SpadeNum + GeoFloat + Send + Sync,
 {
-    let (first, last) = match ls {
-        [] => return vec![],
-        &[only] => return vec![only],
-        &[first, last] => return vec![first, last],
-        &[first, .., last] => (first, last),
-    };
-
-    let visible = visiblity_polygon(ls);
-    let chord = Line::new(first, last);
-
-    let split_index = visible
-        .into_iter()
-        .skip(1)
-        .fold(
-            (0usize, T::zero()),
-            |(farthest_index, farthest_distance), (index, coord)| {
-                let distance = Euclidean.distance(coord, &chord);
-                if distance < farthest_distance {
-                    (farthest_index, farthest_distance)
-                } else {
-                    (index, distance)
-                }
-            },
-        )
-        .0;
-
-    if split_index == 0 || split_index == ls.len() - 1 {
-        println!("Failed to reduce. Skipping.");
-        return vec![first, last];
+    if cdt.exists_constraint(from.fix(), to.fix()) {
+        return vec![from.position(), to.position()];
     }
 
-    let farthest_distance = ls.iter().map(|&v| Euclidean.distance(v, &chord)).fold(
-        T::zero(),
-        |farthest_distance, distance| {
+    let chord = {
+        let [from, to] = [from, to].map(|v| v.position().into_coord());
+        Line::new(from, to)
+    };
+
+    let farthest_distance =
+        CircularIterator::new(from, to, cdt).fold(T::zero(), |farthest_distance, v| {
+            let distance = Euclidean.distance(v.position().into_coord(), &chord);
             if distance > farthest_distance {
                 distance
             } else {
                 farthest_distance
             }
+        });
+
+    if farthest_distance <= eps {
+        return vec![from.position(), to.position()];
+    }
+
+    let (split_vertex, _) = visibility_intersection(from, to, cdt).into_iter().fold(
+        (from, -T::one()), // Placeholder, should always be overwritten
+        |(farthest_vertex, farthest_distance), v| {
+            let distance = Euclidean.distance(v.position().into_coord(), &chord);
+            if distance > farthest_distance {
+                (v, distance)
+            } else {
+                (farthest_vertex, farthest_distance)
+            }
         },
     );
 
-    if farthest_distance > eps {
-        let (mut left, right) = rayon::join(
-            || rdp_preserve(&ls[..=split_index], eps),
-            || rdp_preserve(&ls[split_index..], eps),
-        );
+    // This should never occur
+    if split_vertex == from || split_vertex == to {
+        panic!("Attempted to split at endpoint");
+    }
 
-        left.pop();
-        left.extend_from_slice(&right);
+    let (mut left, right) = rayon::join(
+        || rdp_preserve(from, split_vertex, cdt, eps),
+        || rdp_preserve(split_vertex, to, cdt, eps),
+    );
 
-        return left;
-    }
+    left.pop();
+    left.extend_from_slice(&right);
 
-    vec![first, last]
+    left
 }
 
 pub trait SimplifyRDP<T, Epsilon = T> {
     fn simplify_rdp(&self, eps: Epsilon) -> Self;
 }
 
-impl<T> SimplifyRDP<T> for LineString<T>
-where
-    T: GeoFloat + Send + Sync,
-{
-    fn simplify_rdp(&self, eps: T) -> Self {
-        LineString::new(rdp_preserve(&self.0, eps))
-    }
-}
-
 impl<T> SimplifyRDP<T> for Polygon<T>
 where
-    T: GeoFloat + Send + Sync,
+    T: SpadeNum + GeoFloat + Send + Sync,
 {
     fn simplify_rdp(&self, eps: T) -> Self {
-        let reduced_segments = self
-            .hull_segments()
-            .into_par_iter() // parallelize with rayon
-            .map(|ls| ls.simplify_rdp(eps))
-            .collect::<Vec<_>>();
-
-        Polygon::from_segments(reduced_segments)
+        if self.exterior().0.len() < 3 {
+            return self.clone();
+        }
+
+        let cdt = self.triangulate();
+
+        let segments = cdt
+            .convex_hull()
+            .map(|edge| {
+                rdp_preserve(edge.from(), edge.to(), &cdt, eps)
+                    .into_iter()
+                    .map(|point| point.into_coord())
+                    .collect::<LineString<_>>()
+            })
+            .collect();
+
+        Polygon::from_segments(segments)
     }
 }

src/algorithms/simplify_vw.rs

@@ -18,7 +18,7 @@
 
 // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan
 
-use crate::extensions::ord_triangles::{OrdTriangle, OrdTriangles};
+use crate::types::ord_triangle::{OrdTriangle, OrdTriangles};
 
 use geo::algorithm::{Area, Intersects};
 use geo::geometry::{Coord, Line, LineString, Point, Polygon};