diff --git a/src/algorithms/simplify_charshape.rs b/src/algorithms/simplify_charshape.rs index 43ea68c..7a21dee 100644 --- a/src/algorithms/simplify_charshape.rs +++ b/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::>(); - - let edges = (0..num_vertices) - .map(|i| { - if i == 0 { - [vertices.len() - 1, i] - } else { - [i - 1, i] - } - }) - .collect::>(); - - let tri = - ConstrainedDelaunayTriangulation::>::bulk_load_cdt(vertices, edges).unwrap(); - + let tri = orig.triangulate(); let boundary_edges = tri.convex_hull().map(|edge| edge.rev()).collect::>(); 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 { diff --git a/src/algorithms/simplify_rdp.rs b/src/algorithms/simplify_rdp.rs index cec48ca..6b722f7 100644 --- a/src/algorithms/simplify_rdp.rs +++ b/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, (), CdtEdge<()>>, + until: VertexHandle<'a, Point2, (), CdtEdge<()>>, + cdt: &'a ConstrainedDelaunayTriangulation>, +} + +impl<'a, T: SpadeNum> CircularIterator<'a, T> { + fn new( + from: VertexHandle<'a, Point2, (), CdtEdge<()>>, + until: VertexHandle<'a, Point2, (), CdtEdge<()>>, + cdt: &'a ConstrainedDelaunayTriangulation>, + ) -> Self { + CircularIterator { + current: from, + until, + cdt, + } + } +} + +impl<'a, T: SpadeNum> Iterator for CircularIterator<'a, T> { + type Item = VertexHandle<'a, Point2, (), CdtEdge<()>>; + + fn next(&mut self) -> Option { + 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(ls: &[Coord], eps: T) -> Vec> + Some(vertex) + } +} + +fn rdp_preserve( + from: VertexHandle<'_, Point2, (), CdtEdge<()>>, + to: VertexHandle<'_, Point2, (), CdtEdge<()>>, + cdt: &ConstrainedDelaunayTriangulation>, + eps: T, +) -> Vec> 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 { fn simplify_rdp(&self, eps: Epsilon) -> Self; } -impl SimplifyRDP for LineString -where - T: GeoFloat + Send + Sync, -{ - fn simplify_rdp(&self, eps: T) -> Self { - LineString::new(rdp_preserve(&self.0, eps)) - } -} - impl SimplifyRDP for Polygon 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::>(); - - 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::>() + }) + .collect(); + + Polygon::from_segments(segments) } } diff --git a/src/algorithms/simplify_vw.rs b/src/algorithms/simplify_vw.rs index 2399e4e..38c8612 100644 --- a/src/algorithms/simplify_vw.rs +++ b/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}; diff --git a/src/algorithms/visibility.rs b/src/algorithms/visibility.rs index 6492aa3..2c4edce 100644 --- a/src/algorithms/visibility.rs +++ b/src/algorithms/visibility.rs @@ -18,170 +18,434 @@ // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan -use geo::{coord, Coord, GeoFloat, GeoNum, Kernel, Orientation}; -use std::ops::Sub; - -fn left_visibility_polygon( - mut points_iter: impl Iterator)>, -) -> Vec<(usize, Coord)> { - let first: (usize, Coord) = points_iter.next().unwrap(); - let second = points_iter.next().unwrap(); - let mut stack = Vec::from([first, second]); - - while let Some(v) = points_iter.next() { - if let Some(v) = if matches!( - T::Ker::orient2d(first.1, stack.last().unwrap().1, v.1), - Orientation::CounterClockwise +use crate::extensions::conversions::IntoCoord; +use geo::{GeoNum, Kernel, Orientation}; +use spade::handles::{DirectedEdgeHandle, FixedVertexHandle, VertexHandle}; +use spade::{CdtEdge, ConstrainedDelaunayTriangulation, Point2, SpadeNum, Triangulation}; +use std::collections::BinaryHeap; + +/// Take a given window (left, right) on an edge e and recurse on the visible edges +fn visit_edge<'a, T>( + source: Point2, + edge: DirectedEdgeHandle<'a, Point2, (), CdtEdge<()>, ()>, + left: VertexHandle, (), CdtEdge<()>, ()>, + right: VertexHandle, (), CdtEdge<()>, ()>, + visible: &mut BinaryHeap, (), CdtEdge<()>>>, +) where + T: GeoNum + SpadeNum, +{ + if edge.is_constraint_edge() || edge.is_outer_edge() { + return; + } + + // If coprime point is visible, push to stack + let coprime = edge.opposite_vertex().unwrap(); + // TODO: Replace with LineSideInfo + if matches!( + T::Ker::orient2d( + source.into_coord(), + left.position().into_coord(), + coprime.position().into_coord() + ), + Orientation::Clockwise | Orientation::Collinear + ) && matches!( + T::Ker::orient2d( + source.into_coord(), + right.position().into_coord(), + coprime.position().into_coord() + ), + Orientation::CounterClockwise | Orientation::Collinear + ) { + visible.push(coprime); + } + + // Iterate over new edges + for edge in [edge.prev(), edge.next()] { + // Update horizons + let new_left = if matches!( + T::Ker::orient2d( + source.into_coord(), + left.position().into_coord(), + edge.to().position().into_coord() + ), + Orientation::CounterClockwise, ) { - // Line extends visible region - Some(v) - } else if matches!( + // Vision is restricted by the new edge + left + } else { + edge.to() + }; + + let new_right = if matches!( T::Ker::orient2d( - stack.get(stack.len().sub(2)).unwrap().1, - stack.last().unwrap().1, - v.1 + source.into_coord(), + right.position().into_coord(), + edge.from().position().into_coord() ), - Orientation::CounterClockwise + Orientation::Clockwise, ) { - // Line lies inside the visible region (blocking) - reduce(v, first, &mut stack) // Drop all shadowed points + // Vision is restricted by the new edge + right } else { - // Line lies outside the visible region (shadowed) - skip(first, *stack.last().unwrap(), &mut points_iter) // Iterate until visible - } { - stack.push(v); + edge.from() + }; + + if matches!( + T::Ker::orient2d( + source.into_coord(), + new_left.position().into_coord(), + new_right.position().into_coord() + ), + Orientation::CounterClockwise + ) { + // Left and right have changed side: this edge is not visible + continue; } + + visit_edge(source, edge.rev(), new_left, new_right, visible); } - stack } -fn reduce( - v: (usize, Coord), - first: (usize, Coord), - stack: &mut Vec<(usize, Coord)>, -) -> Option<(usize, Coord)> { - let x = stack.pop().unwrap(); - while matches!( - T::Ker::orient2d(first.1, v.1, stack.last().unwrap().1), - Orientation::CounterClockwise - ) && matches!( - T::Ker::orient2d(x.1, v.1, stack.last().unwrap().1), - Orientation::CounterClockwise - ) { - stack.pop(); +fn visibility_vertex( + mut edge: DirectedEdgeHandle, (), CdtEdge<()>, ()>, + direction: Orientation, +) -> BinaryHeap, (), CdtEdge<()>, ()>> { + let mut visible = BinaryHeap::from([edge.to()]); + let source = edge.from().position(); + + if edge.is_outer_edge() && edge.is_constraint_edge() { + // There are no faces to iterate over + return visible; } - if matches!( - T::Ker::orient2d(first.1, stack.last().unwrap().1, v.1), - Orientation::CounterClockwise - ) { - Some(v) - } else { - None + if matches!(direction, Orientation::Clockwise) { + edge = edge.rev(); } -} -fn skip( - first: (usize, Coord), - last: (usize, Coord), - points_iter: &mut impl Iterator)>, -) -> Option<(usize, Coord)> { - points_iter.find(|&v| { - matches!( - T::Ker::orient2d(first.1, last.1, v.1), - Orientation::CounterClockwise - ) - }) + while let Some(coprime) = edge.opposite_vertex() { + visible.push(coprime); + + let window = match direction { + Orientation::CounterClockwise => edge.next().rev(), + Orientation::Clockwise => edge.prev().rev(), + Orientation::Collinear => panic!(), + }; + visit_edge(source, window, window.from(), window.to(), &mut visible); + + // Advance to next edge + edge = match direction { + Orientation::CounterClockwise => edge.prev().rev(), + Orientation::Clockwise => edge.next().rev(), + Orientation::Collinear => panic!(), + }; + + if edge.is_constraint_edge() { + break; + } + } + visible } -fn merge_walk(x: Vec<(S, T)>, y: Vec<(S, T)>) -> Vec<(S, T)> -where - S: PartialOrd + PartialEq, -{ - let x_iter = x.into_iter(); - let mut y_iter = y.into_iter(); +pub fn visibility_intersection<'a, T: GeoNum + SpadeNum>( + from: VertexHandle, (), CdtEdge<()>>, + to: VertexHandle, (), CdtEdge<()>>, + cdt: &'a ConstrainedDelaunayTriangulation>, +) -> Vec, (), CdtEdge<()>>> { + let from_edge = { + let next = FixedVertexHandle::from_index((from.index() + 1) % cdt.num_vertices()); + cdt.get_edge_from_neighbors(from.fix(), next).unwrap() + }; + let to_edge = { + let index = if to.index() == 0 { + cdt.num_vertices() - 1 + } else { + to.index() - 1 + }; + let prev = FixedVertexHandle::from_index(index); + cdt.get_edge_from_neighbors(to.fix(), prev).unwrap() + }; + + let from_vis = visibility_vertex(from_edge, Orientation::CounterClockwise); + let to_vis = visibility_vertex(to_edge, Orientation::Clockwise); + + binary_heap_intersection(from_vis, to_vis) + .into_sorted_vec() + .into_iter() + .filter(|v| { + if from.index() <= to.index() { + v.index() >= from.index() && v.index() <= to.index() + } else { + v.index() <= to.index() || v.index() >= from.index() + } + }) + .collect() +} - let mut intersection = Vec::new(); - let Some(mut other) = y_iter.next() else { +fn binary_heap_intersection(mut x: BinaryHeap, mut y: BinaryHeap) -> BinaryHeap { + let mut intersection = BinaryHeap::new(); + let Some(mut other) = y.pop() else { return intersection; }; - for item in x_iter { - while item.0 > other.0 { - other = match y_iter.next() { + while let Some(item) = x.pop() { + while item < other { + other = match y.pop() { Some(other) => other, None => return intersection, }; } - if item.0 == other.0 { + if item == other { intersection.push(item); } } intersection } -fn reverse_coords( - ls: impl DoubleEndedIterator)>, -) -> impl Iterator)> { - ls.into_iter().rev().map(|v| { - let (x, y) = v.1.x_y(); - (v.0, coord! {x: -x, y: y}) - }) -} +#[cfg(test)] +mod test { + use crate::algorithms::visibility::{visibility_intersection, visit_edge}; + use crate::extensions::triangulate::Triangulate; + use geo::polygon; + use spade::handles::FixedVertexHandle; + use spade::{Point2, Triangulation}; + use std::collections::BinaryHeap; + + #[test] + fn collinear_test() { + let poly = polygon![ + (x: 0.0, y: 1.0), + (x: 1.0, y: 1.0), + (x: 1.0, y: 0.0), + (x: 0.0, y: 0.0), + ]; -pub fn visiblity_polygon(ls: &[Coord]) -> Vec<(usize, Coord)> { - let iter = ls.iter().copied().enumerate(); - let left = left_visibility_polygon(iter); + let cdt = poly.triangulate(); - let rev_iter = reverse_coords(ls.iter().copied().enumerate()); - let rev_right = left_visibility_polygon(rev_iter); - let right = reverse_coords(rev_right.into_iter()).collect::>(); - merge_walk(left, right) -} + let from = { + let handle = FixedVertexHandle::from_index(0); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(1); + cdt.get_vertex(handle).unwrap() + }; -#[cfg(test)] -mod tests { - use super::visiblity_polygon; - use geo::coord; + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = Vec::::new(); + + assert_eq!(vis, correct); + } #[test] - fn visibility_test() { - let ls = vec![ - coord! { x: 0.0, y: 0.0 }, - coord! { x: 0.0, y: -1.0 }, - coord! { x: -2.0, y: -1.0 }, - coord! { x: -2.0, y: -3.0 }, - coord! { x: 2.0, y: -3.0 }, - coord! { x: 2.0, y: -2.0 }, - coord! { x: -1.0, y: -2.0 }, - coord! { x: 2.0, y: -1.0 }, - coord! { x: 2.0, y: 0.0 }, + fn constrained_test() { + let poly = polygon![ + (x: 0.0, y: 1.0), + (x: 0.5, y: 0.5), + (x: 1.0, y: 1.0), + (x: 1.0, y: 0.0), + (x: 0.0, y: 0.0), ]; - let correct = vec![ - (0, coord! { x: 0.0, y: 0.0 }), - (1, coord! { x: 0.0, y: -1.0 }), - (7, coord! { x: 2.0, y: -1.0 }), - (8, coord! { x: 2.0, y: 0.0 }), + + let cdt = poly.triangulate(); + + let from = { + let handle = FixedVertexHandle::from_index(0); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(2); + cdt.get_vertex(handle).unwrap() + }; + + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![1]; + + assert_eq!(vis, correct); + } + + #[test] + fn central_peak_test() { + let poly = polygon![ + (x: 0.0, y: 1.0), + (x: 0.5, y: 0.5), + (x: 1.0, y: 0.9), + (x: 1.5, y: 0.5), + (x: 2.0, y: 1.0), + (x: 2.0, y: 0.0), + (x: 0.0, y: 0.0), ]; - assert_eq!(visiblity_polygon(&ls), correct); + + let cdt = poly.triangulate(); + + let from = { + let handle = FixedVertexHandle::from_index(0); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(4); + cdt.get_vertex(handle).unwrap() + }; + + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![2]; + + assert_eq!(vis, correct); } #[test] - fn collinear_test() { - let ls = vec![ - coord! { x: 0.0, y: 0.0 }, - coord! { x: 1.0, y: -2.0 }, - coord! { x: 2.0, y: 0.0 }, - coord! { x: 3.0, y: 0.0 }, - coord! { x: 4.0, y: -2.0 }, - coord! { x: 5.0, y: 0.0 }, + fn sawtooth_test() { + let poly = polygon![ + (x: 0.0, y: 0.0), + (x: 1.0, y: -1.0), + (x: 0.0, y: -2.0), + (x: 1.0, y: -3.0), + (x: 2.0, y: -3.0), + (x: 1.0, y: -2.0), + (x: 2.0, y: -1.0), + (x: 2.0, y: 0.0), + (x: 3.0, y: 0.0), + (x: 3.0, y: -4.0), + (x: -1.0, y: -4.0), + (x: -1.0, y: 0.0), ]; - let correct = vec![ - (0, coord! { x: 0.0, y: 0.0 }), - (2, coord! { x: 2.0, y: 0.0 }), - (3, coord! { x: 3.0, y: 0.0 }), - (5, coord! { x: 5.0, y: 0.0 }), + + let cdt = poly.triangulate(); + + let source = Point2::new(1.0, 1.0); + let edge = { + let from = FixedVertexHandle::from_index(0); + let to = FixedVertexHandle::from_index(7); + cdt.get_edge_from_neighbors(from, to).unwrap().rev() + }; + + let mut vis = BinaryHeap::new(); + visit_edge(source, edge, edge.from(), edge.to(), &mut vis); + + let vis = vis + .into_sorted_vec() + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![1, 3, 5, 6]; + + assert_eq!(vis, correct); + } + + #[test] + fn snail_test() { + let poly = polygon![ + (x: 0.0, y: 0.0), + (x: 0.0, y: -3.0), + (x: -2.0, y: -3.0), + (x: -2.0, y: -2.0), + (x: -1.0, y: -2.0), + (x: -3.0, y: -1.0), + (x: -3.0, y: -4.0), + (x: 3.0, y: -4.0), + (x: 3.0, y: 0.0), + (x: 4.0, y: 0.0), + (x: 4.0, y: -5.0), + (x: -4.0, y: -5.0), + (x: -4.0, y: 0.0), + ]; + + let cdt = poly.triangulate(); + + let from = { + let handle = FixedVertexHandle::from_index(0); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(8); + cdt.get_vertex(handle).unwrap() + }; + + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![1, 7]; + + assert_eq!(vis, correct); + } + + #[test] + fn pan_handle_test() { + let poly = polygon![ + (x: -1.0, y: 0.0), + (x: 0.0, y: 0.0), + (x: 0.0, y: -1.0), + (x: 1.0, y: -1.0), + (x: 1.0, y: 0.0), + (x: 2.0, y: 0.0), + (x: 3.0, y: 1.0), + (x: 3.0, y: -2.0), + (x: -2.0, y: -2.0), + (x: -2.0, y: 1.0), + ]; + + let cdt = poly.triangulate(); + + let from = { + let handle = FixedVertexHandle::from_index(0); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(5); + cdt.get_vertex(handle).unwrap() + }; + + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![1, 4]; + + assert_eq!(vis, correct); + } + + #[test] + fn index_shift_test() { + let poly = polygon![ + (x: 3.0, y: 2.0), + (x: 4.0, y: 1.0), + (x: 5.0, y: 1.0), + (x: 5.0, y: 3.0), + (x: 6.0, y: 3.0), + (x: 6.0, y: 0.0), + (x: 0.0, y: 0.0), + (x: 0.0, y: 3.0), + (x: 1.0, y: 3.0), + (x: 1.0, y: 1.0), + (x: 2.0, y: 1.0), ]; - assert_eq!(visiblity_polygon(&ls), correct); + + let cdt = poly.triangulate(); + + let from = { + let handle = FixedVertexHandle::from_index(8); + cdt.get_vertex(handle).unwrap() + }; + let to = { + let handle = FixedVertexHandle::from_index(3); + cdt.get_vertex(handle).unwrap() + }; + + let vis = visibility_intersection(from, to, &cdt) + .into_iter() + .map(|v| v.index()) + .collect::>(); + let correct = vec![0, 2, 9]; + + assert_eq!(vis, correct); } } diff --git a/src/extensions/conversions.rs b/src/extensions/conversions.rs new file mode 100644 index 0000000..d638d46 --- /dev/null +++ b/src/extensions/conversions.rs @@ -0,0 +1,38 @@ +// Copyright 2025- European Centre for Medium-Range Weather Forecasts (ECMWF) + +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at + +// http://www.apache.org/licenses/LICENSE-2.0 + +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// In applying this licence, ECMWF does not waive the privileges and immunities +// granted to it by virtue of its status as an intergovernmental organisation nor +// does it submit to any jurisdiction. + +// Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan + +use geo::{Coord, CoordNum}; +use spade::{Point2, SpadeNum}; + +pub trait IntoCoord { + fn into_coord(self) -> Coord; +} + +impl IntoCoord for Point2 +where + T: CoordNum + SpadeNum, +{ + fn into_coord(self) -> Coord { + Coord { + x: self.x, + y: self.y, + } + } +} diff --git a/src/extensions/mod.rs b/src/extensions/mod.rs index 6063af5..e144f44 100644 --- a/src/extensions/mod.rs +++ b/src/extensions/mod.rs @@ -18,6 +18,7 @@ // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan -pub mod ord_triangles; +pub mod conversions; pub mod segments; +pub mod triangulate; pub mod validation; diff --git a/src/extensions/triangulate.rs b/src/extensions/triangulate.rs new file mode 100644 index 0000000..eb54a7e --- /dev/null +++ b/src/extensions/triangulate.rs @@ -0,0 +1,53 @@ +// Copyright 2025- European Centre for Medium-Range Weather Forecasts (ECMWF) + +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at + +// http://www.apache.org/licenses/LICENSE-2.0 + +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// In applying this licence, ECMWF does not waive the privileges and immunities +// granted to it by virtue of its status as an intergovernmental organisation nor +// does it submit to any jurisdiction. + +// Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan + +use geo::{CoordsIter, GeoNum, Polygon}; +use spade::{ConstrainedDelaunayTriangulation, Point2, SpadeNum}; + +pub trait Triangulate { + fn triangulate(&self) -> ConstrainedDelaunayTriangulation>; +} + +impl Triangulate for Polygon +where + T: SpadeNum + GeoNum, +{ + fn triangulate(&self) -> ConstrainedDelaunayTriangulation> { + let num_vertices = self.exterior().0.len() - 1; + + let vertices = self + .exterior_coords_iter() + .take(num_vertices) // duplicate points are removed + .map(|c| Point2::new(c.x, c.y)) + .collect(); + + let edges = (0..num_vertices) + .map(|i| { + if i == 0 { + [num_vertices - 1, 0] + } else { + [i - 1, i] + } + }) + .collect(); + + ConstrainedDelaunayTriangulation::>::bulk_load_cdt(vertices, edges).unwrap() + } +} diff --git a/src/lib.rs b/src/lib.rs index d3c9a70..a50725b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -18,6 +18,7 @@ // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan +use crate::extensions::triangulate::Triangulate; use crate::extensions::validation::InvalidPolygon; use algorithms::simplify_charshape::SimplifyCharshape; use algorithms::simplify_rdp::SimplifyRDP; @@ -26,9 +27,11 @@ use extensions::validation::Validate; use geo::{Polygon, Winding}; use pyo3::exceptions::PyValueError; use pyo3::prelude::*; +use spade::Triangulation; mod algorithms; mod extensions; +mod types; impl From for PyErr { fn from(err: InvalidPolygon) -> Self { @@ -124,6 +127,17 @@ fn is_valid(poly: Vec<[f64; 2]>) -> PyResult { Ok(poly.is_valid() && poly.exterior().is_cw()) } +#[pyfunction] +fn triangulate(poly: Vec<[f64; 2]>) -> PyResult> { + let poly = Polygon::new(poly.into(), vec![]); + let cdt = poly.triangulate(); + let edges = cdt + .undirected_edges() + .map(|edge| edge.positions().map(|p| p.into())) + .collect(); + Ok(edges) +} + #[pymodule] fn _polyshell(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(reduce_polygon_vw, m)?)?; diff --git a/src/types/mod.rs b/src/types/mod.rs new file mode 100644 index 0000000..5b347bc --- /dev/null +++ b/src/types/mod.rs @@ -0,0 +1,21 @@ +// Copyright 2025- European Centre for Medium-Range Weather Forecasts (ECMWF) + +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at + +// http://www.apache.org/licenses/LICENSE-2.0 + +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// In applying this licence, ECMWF does not waive the privileges and immunities +// granted to it by virtue of its status as an intergovernmental organisation nor +// does it submit to any jurisdiction. + +// Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan + +pub mod ord_triangle; diff --git a/src/extensions/ord_triangles.rs b/src/types/ord_triangle.rs similarity index 80% rename from src/extensions/ord_triangles.rs rename to src/types/ord_triangle.rs index 39237c8..76b59eb 100644 --- a/src/extensions/ord_triangles.rs +++ b/src/types/ord_triangle.rs @@ -19,13 +19,13 @@ // Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan use geo::algorithm::Area; -use geo::geometry::{Coord, LineString, Point, Triangle}; -use geo::{CoordFloat, CoordNum}; +use geo::geometry::{LineString, Triangle}; +use geo::{Coord, CoordFloat, CoordNum, Point}; // geo::Triangle sorts the vertices so they are always oriented ccw // this is not acceptable for our use case #[derive(Copy, Clone, Hash, Eq, PartialEq)] -pub struct OrdTriangle(pub Coord, pub Coord, pub Coord); +pub struct OrdTriangle(Coord, Coord, Coord); impl OrdTriangle { pub fn new(v1: Coord, v2: Coord, v3: Coord) -> Self { @@ -34,8 +34,8 @@ impl OrdTriangle { } impl From> for Triangle { - fn from(other: OrdTriangle) -> Self { - Self(other.0, other.1, other.2) + fn from(from: OrdTriangle) -> Self { + Triangle::new(from.0, from.1, from.2) } } @@ -55,12 +55,9 @@ pub trait OrdTriangles { impl OrdTriangles for LineString { fn ord_triangles(&'_ self) -> impl ExactSizeIterator> + '_ { - self.0.windows(3).map(|w| unsafe { - OrdTriangle::new( - *w.get_unchecked(0), - *w.get_unchecked(1), - *w.get_unchecked(2), - ) + self.0.windows(3).map(|w| { + let &[x, y, z] = w else { unreachable!() }; + OrdTriangle::new(x, y, z) }) } } diff --git a/tests/method_cases.py b/tests/method_cases.py index 895119c..89ad7fd 100644 --- a/tests/method_cases.py +++ b/tests/method_cases.py @@ -20,16 +20,13 @@ # Copyright 2025- Niall Oswald and Kenneth Martin and Jo Wayne Tan # -import pytest from polyshell import ReductionMethod -from pytest_cases import case def case_char() -> ReductionMethod: return ReductionMethod.CHARSHAPE -@case(marks=pytest.mark.xfail(reason="Known bug")) def case_rdp() -> ReductionMethod: return ReductionMethod.RDP