diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 9b603ecf0..ab0e24660 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -28,6 +28,8 @@ * * Add `wkt!` macro to define geometries at compile time. +* Add `TriangulateSpade` trait which provides (un)constrained Delaunay Triangulations for all `geo_types` via the `spade` crate + * ## 0.26.0 diff --git a/geo/Cargo.toml b/geo/Cargo.toml index f4cf7f648..e7f4a87dd 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -13,13 +13,14 @@ rust-version = "1.65" categories = ["science::geo"] [features] -default = ["earcutr"] +default = ["earcutr", "spade"] use-proj = ["proj"] proj-network = ["use-proj", "proj/network"] use-serde = ["serde", "geo-types/serde"] [dependencies] earcutr = { version = "0.4.2", optional = true } +spade = { version = "2.2.0", optional = true } float_next_after = "1.0.0" geo-types = { version = "0.7.9", features = ["approx", "use-rstar_0_11"] } geographiclib-rs = { version = "0.2.3", default-features = false } diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 386047f47..791ad98fe 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -252,6 +252,12 @@ pub mod triangulate_earcut; #[cfg(feature = "earcutr")] pub use triangulate_earcut::TriangulateEarcut; +/// Triangulate polygons using an (un)constrained [Delaunay Triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation) algorithm. +#[cfg(feature = "spade")] +pub mod triangulate_spade; +#[cfg(feature = "spade")] +pub use triangulate_spade::TriangulateSpade; + /// Vector Operations for 2D coordinates mod vector_ops; pub use vector_ops::Vector2DOps; diff --git a/geo/src/algorithm/triangulate_spade.rs b/geo/src/algorithm/triangulate_spade.rs new file mode 100644 index 000000000..fb38ff4f1 --- /dev/null +++ b/geo/src/algorithm/triangulate_spade.rs @@ -0,0 +1,782 @@ +use geo_types::{Coord, Line, Point, Triangle}; +use spade::{ + ConstrainedDelaunayTriangulation, DelaunayTriangulation, Point2, SpadeNum, Triangulation, +}; + +use crate::{ + line_intersection::line_intersection, CoordsIter, EuclideanDistance, GeoFloat, + LineIntersection, LinesIter, +}; +use crate::{Centroid, Contains}; + +// ======== Config ============ + +/// Collection of parameters that influence the precision of the algorithm in some sense (see +/// explanations on fields of this struct) +/// +/// This implements the `Default` trait and you can just use it most of the time +#[derive(Debug, Clone)] +pub struct SpadeTriangulationConfig { + /// Coordinates within this radius are snapped to the same position. For any two `Coords` there's + /// no real way to influence the decision when choosing the snapper and the snappee + snap_radius: T, +} + +impl Default for SpadeTriangulationConfig +where + T: SpadeTriangulationFloat, +{ + fn default() -> Self { + Self { + snap_radius: >::from(0.000_1), + } + } +} + +// ====== Error ======== + +#[derive(Debug)] +pub enum TriangulationError { + SpadeError(spade::InsertionError), + LoopTrap, + ConstraintFailure, +} + +impl std::fmt::Display for TriangulationError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{:?}", self) + } +} + +impl std::error::Error for TriangulationError {} + +pub type TriangulationResult = Result; + +// ======= Float trait ======== + +pub trait SpadeTriangulationFloat: GeoFloat + SpadeNum {} +impl SpadeTriangulationFloat for T {} + +// ======= Triangulation trait ========= + +pub type Triangles = Vec>; + +// seal the trait that needs to be implemented for TriangulateSpade to be implemented. This is done +// so that we don't leak these weird methods on the public interface. +mod private { + use super::*; + pub trait TriangulationRequirementTrait<'a, T> + where + T: SpadeTriangulationFloat, + { + /// collect all the lines that are relevant for triangulations from the geometric object that + /// should be triangulated. + /// + /// intersecting lines are allowed + fn lines(&'a self) -> Vec>; + /// collect all the coords that are relevant for triangulations from the geometric object that + /// should be triangulated + fn coords(&'a self) -> Vec>; + /// define a predicate that decides if a point is inside of the object (used for constrained triangulation) + fn contains_point(&'a self, p: Point) -> bool; + + // processing of the lines that prepare the lines for triangulation. + // + // `spade` has the general limitation that constraint lines cannot intersect or else it + // will panic. This is why we need to manually split up the lines into smaller parts at the + // intersection point + // + // there's also a preprocessing step which tries to minimize the risk of failure of the algo + // through edge cases (thin/flat triangles are prevented as much as possible & lines are deduped, ...) + fn cleanup_lines(lines: Vec>, snap_radius: T) -> TriangulationResult>> { + let (known_coords, lines) = preprocess_lines(lines, snap_radius); + prepare_intersection_contraint(lines, known_coords, snap_radius) + } + } +} + +/// Triangulate polygons using a [Delaunay +/// Triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation) +/// +/// This trait contains both constrained and unconstrained triangulation methods. To read more +/// about the differences of these methods also consult [this +/// page](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation) +pub trait TriangulateSpade<'a, T>: private::TriangulationRequirementTrait<'a, T> +where + T: SpadeTriangulationFloat, +{ + /// returns a triangulation that's solely based on the points of the geometric object + /// + /// The triangulation is guaranteed to be Delaunay + /// + /// Note that the lines of the triangulation don't necessarily follow the lines of the input + /// geometry. If you wish to achieve that take a look at the `constrained_triangulation` and the + /// `constrained_outer_triangulation` functions. + /// + /// ```rust + /// use geo::TriangulateSpade; + /// use geo::{Polygon, LineString, Coord}; + /// let u_shape = Polygon::new( + /// LineString::new(vec![ + /// Coord { x: 0.0, y: 0.0 }, + /// Coord { x: 1.0, y: 0.0 }, + /// Coord { x: 1.0, y: 1.0 }, + /// Coord { x: 2.0, y: 1.0 }, + /// Coord { x: 2.0, y: 0.0 }, + /// Coord { x: 3.0, y: 0.0 }, + /// Coord { x: 3.0, y: 3.0 }, + /// Coord { x: 0.0, y: 3.0 }, + /// ]), + /// vec![], + /// ); + /// let unconstrained_triangulation = u_shape.unconstrained_triangulation().unwrap(); + /// let num_triangles = unconstrained_triangulation.len(); + /// assert_eq!(num_triangles, 8); + /// ``` + /// + fn unconstrained_triangulation(&'a self) -> TriangulationResult> { + let points = self.coords(); + points + .into_iter() + .map(to_spade_point) + .try_fold(DelaunayTriangulation::>::new(), |mut tris, p| { + tris.insert(p).map_err(TriangulationError::SpadeError)?; + Ok(tris) + }) + .map(triangulation_to_triangles) + } + + /// returns triangulation that's based on the points of the geometric object and also + /// incorporates the lines of the input geometry + /// + /// The triangulation is not guaranteed to be Delaunay because of the constraint lines + /// + /// This outer triangulation also contains triangles that are not included in the input + /// geometry if it wasn't convex. Here's an example: + /// + /// ```text + /// ┌──────────────────┐ + /// │\ __/│ + /// │ \ __/ / │ + /// │ \ __/ / │ + /// │ \ __/ / │ + /// │ \/ / │ + /// │ ┌──────┐ │ + /// │ /│\:::::│\ │ + /// │ / │:\::::│ \ │ + /// │ / │::\:::│ \ │ + /// │ / │:::\::│ \ │ + /// │/ │::::\:│ \│ + /// └─────┘______└─────┘ + /// ``` + /// + /// ```rust + /// use geo::TriangulateSpade; + /// use geo::{Polygon, LineString, Coord}; + /// let u_shape = Polygon::new( + /// LineString::new(vec![ + /// Coord { x: 0.0, y: 0.0 }, + /// Coord { x: 1.0, y: 0.0 }, + /// Coord { x: 1.0, y: 1.0 }, + /// Coord { x: 2.0, y: 1.0 }, + /// Coord { x: 2.0, y: 0.0 }, + /// Coord { x: 3.0, y: 0.0 }, + /// Coord { x: 3.0, y: 3.0 }, + /// Coord { x: 0.0, y: 3.0 }, + /// ]), + /// vec![], + /// ); + /// // we use the default [`SpadeTriangulationConfig`] here + /// let constrained_outer_triangulation = + /// u_shape.constrained_outer_triangulation(Default::default()).unwrap(); + /// let num_triangles = constrained_outer_triangulation.len(); + /// assert_eq!(num_triangles, 8); + /// ``` + /// + /// The outer triangulation of the top down U-shape contains extra triangles marked + /// with ":". If you want to exclude those, take a look at `constrained_triangulation` + fn constrained_outer_triangulation( + &'a self, + config: SpadeTriangulationConfig, + ) -> TriangulationResult> { + let lines = self.lines(); + let lines = Self::cleanup_lines(lines, config.snap_radius)?; + lines + .into_iter() + .map(to_spade_line) + .try_fold( + ConstrainedDelaunayTriangulation::>::new(), + |mut cdt, [start, end]| { + let start = cdt.insert(start).map_err(TriangulationError::SpadeError)?; + let end = cdt.insert(end).map_err(TriangulationError::SpadeError)?; + // safety check (to prevent panic) whether we can add the line + if !cdt.can_add_constraint(start, end) { + return Err(TriangulationError::ConstraintFailure); + } + cdt.add_constraint(start, end); + Ok(cdt) + }, + ) + .map(triangulation_to_triangles) + } + + /// returns triangulation that's based on the points of the geometric object and also + /// incorporates the lines of the input geometry + /// + /// The triangulation is not guaranteed to be Delaunay because of the constraint lines + /// + /// This triangulation only contains triangles that are included in the input geometry. + /// Here's an example: + /// + /// ```text + /// ┌──────────────────┐ + /// │\ __/│ + /// │ \ __/ / │ + /// │ \ __/ / │ + /// │ \ __/ / │ + /// │ \/ / │ + /// │ ┌──────┐ │ + /// │ /│ │\ │ + /// │ / │ │ \ │ + /// │ / │ │ \ │ + /// │ / │ │ \ │ + /// │/ │ │ \│ + /// └─────┘ └─────┘ + /// ``` + /// + /// ```rust + /// use geo::TriangulateSpade; + /// use geo::{Polygon, LineString, Coord}; + /// let u_shape = Polygon::new( + /// LineString::new(vec![ + /// Coord { x: 0.0, y: 0.0 }, + /// Coord { x: 1.0, y: 0.0 }, + /// Coord { x: 1.0, y: 1.0 }, + /// Coord { x: 2.0, y: 1.0 }, + /// Coord { x: 2.0, y: 0.0 }, + /// Coord { x: 3.0, y: 0.0 }, + /// Coord { x: 3.0, y: 3.0 }, + /// Coord { x: 0.0, y: 3.0 }, + /// ]), + /// vec![], + /// ); + /// // we use the default [`SpadeTriangulationConfig`] here + /// let constrained_triangulation = u_shape.constrained_triangulation(Default::default()).unwrap(); + /// let num_triangles = constrained_triangulation.len(); + /// assert_eq!(num_triangles, 6); + /// ``` + /// + /// Compared to the `constrained_outer_triangulation` it only includes the triangles + /// inside of the input geometry + fn constrained_triangulation( + &'a self, + config: SpadeTriangulationConfig, + ) -> TriangulationResult> { + self.constrained_outer_triangulation(config) + .map(|triangles| { + triangles + .into_iter() + .filter(|triangle| { + let center = triangle.centroid(); + self.contains_point(center) + }) + .collect::>() + }) + } +} + +/// conversion from spade triangulation back to geo triangles +fn triangulation_to_triangles(triangulation: T) -> Triangles +where + T: Triangulation>, + F: SpadeTriangulationFloat, +{ + triangulation + .inner_faces() + .map(|face| face.positions()) + .map(|points| points.map(|p| Coord:: { x: p.x, y: p.y })) + .map(Triangle::from) + .collect::>() +} + +// ========== Triangulation trait impls ============ + +// everything that satisfies the requirement methods automatically implements the triangulation +impl<'a, T, G> TriangulateSpade<'a, T> for G +where + T: SpadeTriangulationFloat, + G: private::TriangulationRequirementTrait<'a, T>, +{ +} + +impl<'a, 'l, T, G> private::TriangulationRequirementTrait<'a, T> for G +where + 'a: 'l, + T: SpadeTriangulationFloat, + G: LinesIter<'l, Scalar = T> + CoordsIter + Contains>, +{ + fn coords(&'a self) -> Vec> { + self.coords_iter().collect::>() + } + + fn lines(&'a self) -> Vec> { + self.lines_iter().collect() + } + + fn contains_point(&'a self, p: Point) -> bool { + self.contains(&p) + } +} + +// it would be cool to impl the trait for GS: AsRef<[G]> but I wasn't able to get this to compile +// (yet) + +impl<'a, T, G> private::TriangulationRequirementTrait<'a, T> for Vec +where + T: SpadeTriangulationFloat, + G: TriangulateSpade<'a, T>, +{ + fn coords(&'a self) -> Vec> { + self.iter().flat_map(|g| g.coords()).collect::>() + } + + fn lines(&'a self) -> Vec> { + self.iter().flat_map(|g| g.lines()).collect::>() + } + + fn contains_point(&'a self, p: Point) -> bool { + self.iter().any(|g| g.contains_point(p)) + } +} + +impl<'a, T, G> private::TriangulationRequirementTrait<'a, T> for &[G] +where + T: SpadeTriangulationFloat, + G: TriangulateSpade<'a, T>, +{ + fn coords(&'a self) -> Vec> { + self.iter().flat_map(|g| g.coords()).collect::>() + } + + fn lines(&'a self) -> Vec> { + self.iter().flat_map(|g| g.lines()).collect::>() + } + + fn contains_point(&'a self, p: Point) -> bool { + self.iter().any(|g| g.contains_point(p)) + } +} + +// ========== Triangulation trait impl helpers ============ + +fn prepare_intersection_contraint( + mut lines: Vec>, + mut known_points: Vec>, + snap_radius: T, +) -> Result>, TriangulationError> { + // Rule 2 of "Power of 10" rules (NASA) + // safety net. We can't prove that the `while let` loop isn't going to run infinitely, so + // we abort after a fixed amount of iterations. In case that the iteration seems to loop + // indefinitely this check will return an Error indicating the infinite loop. + let mut loop_count = 1000; + let mut loop_check = || { + loop_count -= 1; + (loop_count != 0) + .then_some(()) + .ok_or(TriangulationError::LoopTrap) + }; + + while let Some((indices, intersection)) = { + let mut iter = iter_line_pairs(&lines); + iter.find_map(find_intersecting_lines_fn) + } { + loop_check()?; + let [l0, l1] = remove_lines_by_index(indices, &mut lines); + let new_lines = split_lines([l0, l1], intersection); + let new_lines = cleanup_filter_lines(new_lines, &lines, &mut known_points, snap_radius); + + lines.extend(new_lines); + } + + Ok(lines) +} + +/// iterates over all combinations (a,b) of lines in a vector where a != b +fn iter_line_pairs( + lines: &[Line], +) -> impl Iterator); 2]> { + lines.iter().enumerate().flat_map(|(idx0, line0)| { + lines + .iter() + .enumerate() + .filter(move |(idx1, line1)| *idx1 != idx0 && line0 != *line1) + .map(move |(idx1, line1)| [(idx0, line0), (idx1, line1)]) + }) +} + +/// checks whether two lines are intersecting and if so, checks the intersection to not be ill +/// formed +/// +/// returns +/// - [usize;2] : sorted indexes of lines, smaller one comes first +/// - intersection : type of intersection +fn find_intersecting_lines_fn( + [(idx0, line0), (idx1, line1)]: [(usize, &Line); 2], +) -> Option<([usize; 2], LineIntersection)> { + line_intersection(*line0, *line1) + .filter(|intersection| { + match intersection { + // intersection is not located in both lines + LineIntersection::SinglePoint { is_proper, .. } if !is_proper => false, + // collinear intersection is length zero line + LineIntersection::Collinear { intersection } + if intersection.start == intersection.end => + { + false + } + _ => true, + } + }) + .map(|intersection| ([idx0, idx1], intersection)) +} + +/// removes two lines by index in a safe way since the second index can be invalidated after +/// the first line was removed (remember `.remove(idx)` returns the element and shifts the tail +/// of the vector in direction of its start to close the gap) +fn remove_lines_by_index( + mut indices: [usize; 2], + lines: &mut Vec>, +) -> [Line; 2] { + indices.sort(); + let [idx0, idx1] = indices; + let l1 = lines.remove(idx1); + let l0 = lines.remove(idx0); + [l0, l1] +} + +/// split lines based on intersection kind: +/// +/// - intersection point: create 4 new lines from the existing line's endpoints to the intersection +/// point +/// - collinear: create 3 new lines (before overlap, overlap, after overlap) +fn split_lines( + [l0, l1]: [Line; 2], + intersection: LineIntersection, +) -> Vec> { + match intersection { + LineIntersection::SinglePoint { intersection, .. } => [ + (l0.start, intersection), + (l0.end, intersection), + (l1.start, intersection), + (l1.end, intersection), + ] + .map(|(a, b)| Line::new(a, b)) + .to_vec(), + LineIntersection::Collinear { .. } => { + let mut points = [l0.start, l0.end, l1.start, l1.end]; + // sort points by their coordinate values to resolve ambiguities + points.sort_by(|a, b| { + a.x.partial_cmp(&b.x) + .expect("sorting points by coordinate x failed") + .then_with(|| { + a.y.partial_cmp(&b.y) + .expect("sorting points by coordinate y failed") + }) + }); + // since all points are on one line we can just create new lines from consecutive + // points after sorting + points + .windows(2) + .map(|win| Line::new(win[0], win[1])) + .collect::>() + } + } +} + +/// new lines from the `split_lines` function may contain a variety of ill formed lines, this +/// function cleans all of these cases up +fn cleanup_filter_lines( + lines_need_check: Vec>, + existing_lines: &[Line], + known_points: &mut Vec>, + snap_radius: T, +) -> Vec> { + lines_need_check + .into_iter() + .map(|mut line| { + line.start = snap_or_register_point(line.start, known_points, snap_radius); + line.end = snap_or_register_point(line.end, known_points, snap_radius); + line + }) + .filter(|l| l.start != l.end) + .filter(|l| !existing_lines.contains(l)) + .filter(|l| !existing_lines.contains(&Line::new(l.end, l.start))) + .collect::>() +} + +/// snap point to the nearest existing point if it's close enough +/// +/// snap_radius can be configured via the third parameter of this function +fn snap_or_register_point( + point: Coord, + known_points: &mut Vec>, + snap_radius: T, +) -> Coord { + known_points + .iter() + // find closest + .min_by(|a, b| { + a.euclidean_distance(&point) + .partial_cmp(&b.euclidean_distance(&point)) + .expect("Couldn't compare coordinate distances") + }) + // only snap if closest is within epsilone range + .filter(|nearest_point| nearest_point.euclidean_distance(&point) < snap_radius) + .cloned() + // otherwise register and use input point + .unwrap_or_else(|| { + known_points.push(point); + point + }) +} + +/// preprocesses lines so that we're less likely to hit issues when using the spade triangulation +fn preprocess_lines( + lines: Vec>, + snap_radius: T, +) -> (Vec>, Vec>) { + let mut known_coords: Vec> = vec![]; + let capacity = lines.len(); + let lines = lines + .into_iter() + .fold(Vec::with_capacity(capacity), |mut lines, mut line| { + // deduplicating: + + // 1. snap coords of lines to existing coords + line.start = snap_or_register_point(line.start, &mut known_coords, snap_radius); + line.end = snap_or_register_point(line.end, &mut known_coords, snap_radius); + if + // 2. make sure line isn't degenerate (no length when start == end) + line.start != line.end + // 3. make sure line or flipped line wasn't already added + && !lines.contains(&line) + && !lines.contains(&Line::new(line.end, line.start)) + { + lines.push(line) + } + + lines + }); + (known_coords, lines) +} + +/// converts Line to something somewhat similar in the spade world +fn to_spade_line(line: Line) -> [Point2; 2] { + [to_spade_point(line.start), to_spade_point(line.end)] +} + +/// converts Coord to something somewhat similar in the spade world +fn to_spade_point(coord: Coord) -> Point2 { + Point2::new(coord.x, coord.y) +} + +#[cfg(test)] +mod spade_triangulation { + use super::*; + use geo_types::*; + + fn assert_num_triangles( + triangulation: &TriangulationResult>, + num: usize, + ) { + assert_eq!( + triangulation + .as_ref() + .map(|tris| tris.len()) + .expect("triangulation success"), + num + ) + } + + #[test] + fn basic_triangle_triangulates() { + let triangulation = Triangle::new( + Coord { x: 0.0, y: 0.0 }, + Coord { x: 1.0, y: 0.0 }, + Coord { x: 0.0, y: 1.0 }, + ) + .unconstrained_triangulation(); + + assert_num_triangles(&triangulation, 1); + } + + #[test] + fn basic_rectangle_triangulates() { + let triangulation = Rect::new(Coord { x: 0.0, y: 0.0 }, Coord { x: 1.0, y: 1.0 }) + .unconstrained_triangulation(); + + assert_num_triangles(&triangulation, 2); + } + + #[test] + fn basic_polygon_triangulates() { + let triangulation = Polygon::new( + LineString::new(vec![ + Coord { x: 0.0, y: 1.0 }, + Coord { x: -1.0, y: 0.0 }, + Coord { x: -0.5, y: -1.0 }, + Coord { x: 0.5, y: -1.0 }, + Coord { x: 1.0, y: 0.0 }, + ]), + vec![], + ) + .unconstrained_triangulation(); + + assert_num_triangles(&triangulation, 3); + } + + #[test] + fn overlapping_triangles_triangulate_unconstrained() { + let triangles = vec![ + Triangle::new( + Coord { x: 0.0, y: 0.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 0.0, y: 2.0 }, + ), + Triangle::new( + Coord { x: 1.0, y: 1.0 }, + Coord { x: -1.0, y: 1.0 }, + Coord { x: 1.0, y: -1.0 }, + ), + ]; + + let unconstrained_triangulation = triangles.unconstrained_triangulation(); + assert_num_triangles(&unconstrained_triangulation, 4); + } + + #[test] + fn overlapping_triangles_triangulate_constrained_outer() { + let triangles = vec![ + Triangle::new( + Coord { x: 0.0, y: 0.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 0.0, y: 2.0 }, + ), + Triangle::new( + Coord { x: 1.0, y: 1.0 }, + Coord { x: -1.0, y: 1.0 }, + Coord { x: 1.0, y: -1.0 }, + ), + ]; + + let constrained_outer_triangulation = + triangles.constrained_outer_triangulation(Default::default()); + assert_num_triangles(&constrained_outer_triangulation, 8); + } + + #[test] + fn overlapping_triangles_triangulate_constrained() { + let triangles = vec![ + Triangle::new( + Coord { x: 0.0, y: 0.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 0.0, y: 2.0 }, + ), + Triangle::new( + Coord { x: 1.0, y: 1.0 }, + Coord { x: -1.0, y: 1.0 }, + Coord { x: 1.0, y: -1.0 }, + ), + ]; + + let constrained_outer_triangulation = + triangles.constrained_triangulation(Default::default()); + assert_num_triangles(&constrained_outer_triangulation, 6); + } + + #[test] + fn u_shaped_polygon_triangulates_unconstrained() { + let u_shape = Polygon::new( + LineString::new(vec![ + Coord { x: 0.0, y: 0.0 }, + Coord { x: 1.0, y: 0.0 }, + Coord { x: 1.0, y: 1.0 }, + Coord { x: 2.0, y: 1.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 3.0, y: 0.0 }, + Coord { x: 3.0, y: 3.0 }, + Coord { x: 0.0, y: 3.0 }, + ]), + vec![], + ); + + let unconstrained_triangulation = u_shape.unconstrained_triangulation(); + assert_num_triangles(&unconstrained_triangulation, 8); + } + + #[test] + fn u_shaped_polygon_triangulates_constrained_outer() { + let u_shape = Polygon::new( + LineString::new(vec![ + Coord { x: 0.0, y: 0.0 }, + Coord { x: 1.0, y: 0.0 }, + Coord { x: 1.0, y: 1.0 }, + Coord { x: 2.0, y: 1.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 3.0, y: 0.0 }, + Coord { x: 3.0, y: 3.0 }, + Coord { x: 0.0, y: 3.0 }, + ]), + vec![], + ); + + let constrained_outer_triangulation = + u_shape.constrained_outer_triangulation(Default::default()); + assert_num_triangles(&constrained_outer_triangulation, 8); + } + + #[test] + fn u_shaped_polygon_triangulates_constrained_inner() { + let u_shape = Polygon::new( + LineString::new(vec![ + Coord { x: 0.0, y: 0.0 }, + Coord { x: 1.0, y: 0.0 }, + Coord { x: 1.0, y: 1.0 }, + Coord { x: 2.0, y: 1.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 3.0, y: 0.0 }, + Coord { x: 3.0, y: 3.0 }, + Coord { x: 0.0, y: 3.0 }, + ]), + vec![], + ); + + let constrained_triangulation = u_shape.constrained_triangulation(Default::default()); + assert_num_triangles(&constrained_triangulation, 6); + } + + #[test] + fn various_snap_radius_works() { + let u_shape = Polygon::new( + LineString::new(vec![ + Coord { x: 0.0, y: 0.0 }, + Coord { x: 1.0, y: 0.0 }, + Coord { x: 1.0, y: 1.0 }, + Coord { x: 2.0, y: 1.0 }, + Coord { x: 2.0, y: 0.0 }, + Coord { x: 3.0, y: 0.0 }, + Coord { x: 3.0, y: 3.0 }, + Coord { x: 0.0, y: 3.0 }, + ]), + vec![], + ); + + for snap_with in (1..6).map(|pow| 0.1_f64.powi(pow)) { + let constrained_triangulation = + u_shape.constrained_triangulation(SpadeTriangulationConfig { + snap_radius: snap_with, + }); + assert_num_triangles(&constrained_triangulation, 6); + } + } +}