From c49e41c8d1e0aca259a85503ade9e003be2c8cad Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sun, 13 Feb 2022 17:24:13 +0100 Subject: [PATCH] Adds documentation for `refine`, `with_angle_limit` and `exclude_outer_faces` Various bug fixes & improvements to the refinement --- examples/svg_renderer/main.rs | 19 +- examples/svg_renderer/scenario_list.rs | 95 ++- images/angle_limit_00.svg | 74 +++ images/angle_limit_20.svg | 120 ++++ images/angle_limit_30.svg | 232 ++++++++ images/angle_limit_34.svg | 472 +++++++++++++++ images/exclude_outer_faces_refined.svg | 274 +++++++++ images/exclude_outer_faces_unrefined.svg | 190 ++++++ images/refined.svg | 706 +++++++++++++++++++++++ images/unrefined.svg | 190 ++++++ src/delaunay_core/dcel_operations.rs | 11 +- src/delaunay_core/mod.rs | 2 +- src/delaunay_core/refinement.rs | 214 ++++++- src/delaunay_core/triangulation_ext.rs | 20 +- src/lib.rs | 2 +- 15 files changed, 2560 insertions(+), 61 deletions(-) create mode 100644 images/angle_limit_00.svg create mode 100644 images/angle_limit_20.svg create mode 100644 images/angle_limit_30.svg create mode 100644 images/angle_limit_34.svg create mode 100644 images/exclude_outer_faces_refined.svg create mode 100644 images/exclude_outer_faces_unrefined.svg create mode 100644 images/refined.svg create mode 100644 images/unrefined.svg diff --git a/examples/svg_renderer/main.rs b/examples/svg_renderer/main.rs index d15c078..f59d146 100644 --- a/examples/svg_renderer/main.rs +++ b/examples/svg_renderer/main.rs @@ -38,7 +38,24 @@ fn main() -> Result { scenario_list::dual_edge_example().save_to_svg("dual_edges", "images/dual_edges.svg")?; scenario_list::project_point_scenario() .save_to_svg("project_point", "images/project_point.svg")?; - scenario_list::refinement_scenario().save_to_svg("refinement", "images/refinement.svg")?; + scenario_list::refinement_scenario(true).save_to_svg("refined", "images/refined.svg")?; + scenario_list::refinement_scenario(false).save_to_svg("unrefined", "images/unrefined.svg")?; + scenario_list::angle_limit_scenario(0.0) + .save_to_svg("angle_limit_0", "images/angle_limit_00.svg")?; + scenario_list::angle_limit_scenario(20.0) + .save_to_svg("angle_limit_20", "images/angle_limit_20.svg")?; + scenario_list::angle_limit_scenario(30.0) + .save_to_svg("angle_limit_30", "images/angle_limit_30.svg")?; + scenario_list::angle_limit_scenario(34.0) + .save_to_svg("angle_limit_34", "images/angle_limit_34.svg")?; + + scenario_list::exclude_outer_faces_scenario(false).save_to_svg( + "exclude_unrefined", + "images/exclude_outer_faces_unrefined.svg", + )?; + scenario_list::exclude_outer_faces_scenario(true) + .save_to_svg("exclude_refined", "images/exclude_outer_faces_refined.svg")?; + Ok(()) } diff --git a/examples/svg_renderer/scenario_list.rs b/examples/svg_renderer/scenario_list.rs index b5c558f..46c34ca 100644 --- a/examples/svg_renderer/scenario_list.rs +++ b/examples/svg_renderer/scenario_list.rs @@ -2,13 +2,14 @@ use super::quicksketch::{ ArrowType, HorizontalAlignment, Point, Sketch, SketchColor, SketchElement, SketchFill, StrokeStyle, Vector, }; + use cgmath::{Angle, Bounded, Deg, EuclideanSpace, InnerSpace, Point2, Vector2}; use spade::{ handles::{ FixedDirectedEdgeHandle, VoronoiVertex::{self, Inner, Outer}, }, - InsertionError, RefinementParameters, Triangulation as _, + AngleLimit, InsertionError, RefinementParameters, Triangulation as _, }; use crate::{ @@ -901,7 +902,52 @@ pub fn project_point_scenario() -> Sketch { sketch } -pub fn refinement_scenario() -> Sketch { +pub fn refinement_scenario(do_refine: bool) -> Sketch { + let mut cdt = create_refinement_cdt(); + + let parameters = RefinementParameters::default(); + + if do_refine { + cdt.refine(parameters); + } + + convert_refinement_cdt(&mut cdt) +} + +pub fn exclude_outer_faces_scenario(do_refine: bool) -> Sketch { + let mut cdt = create_refinement_cdt(); + + let num_additional_vertices = if do_refine { 500 } else { 0 }; + + let parameters = RefinementParameters::::default() + .exclude_outer_faces(&cdt) + .with_max_additional_vertices(num_additional_vertices); + + let result = cdt.refine(parameters); + for face in &result.excluded_faces { + cdt.face_data_mut(*face).fill = SketchFill::solid(SketchColor::TAN); + } + + convert_refinement_cdt(&mut cdt) +} + +fn convert_refinement_cdt(cdt: &mut Cdt) -> Sketch { + for vertex in cdt.fixed_vertices() { + cdt.vertex_data_mut(vertex).radius = 0.4; + } + for edge in cdt.fixed_undirected_edges() { + if cdt.is_constraint_edge(edge) { + cdt.undirected_edge_data_mut(edge).data_mut().color = SketchColor::DARK_RED; + } else { + cdt.undirected_edge_data_mut(edge).data_mut().color = SketchColor::DARK_GRAY; + } + } + let mut sketch = convert_triangulation(cdt, &ConversionOptions::default()); + sketch.set_width(360); + sketch +} + +fn create_refinement_cdt() -> Cdt { let mut cdt = Cdt::new(); let v0 = VertexType::new(0.0, 0.0); let v1 = VertexType::new(0.0, 100.0); @@ -968,22 +1014,45 @@ pub fn refinement_scenario() -> Sketch { for window in inner_vertices.windows(2) { cdt.add_constraint_edge(window[0], window[1]).unwrap(); } + cdt +} - let parameters = RefinementParameters::default(); +pub fn angle_limit_scenario(angle_limit_degrees: f64) -> Sketch { + let mut cdt = create_angle_limit_cdt(); - cdt.refine(parameters); + cdt.refine( + RefinementParameters::default().with_angle_limit(AngleLimit::from_deg(angle_limit_degrees)), + ); + let mut result = convert_refinement_cdt(&mut cdt); + result.set_width(200); + result +} - for vertex in cdt.fixed_vertices() { - cdt.vertex_data_mut(vertex).radius = 0.4; - } +fn create_angle_limit_cdt() -> Cdt { + let mut cdt = Cdt::new(); - for edge in cdt.fixed_undirected_edges() { - if cdt.is_constraint_edge(edge) { - cdt.undirected_edge_data_mut(edge).data_mut().color = SketchColor::DARK_RED; - } else { - cdt.undirected_edge_data_mut(edge).data_mut().color = SketchColor::DARK_GRAY; + let num_slices = 22; + + cdt.insert(VertexType::new(0.0, 0.0)).unwrap(); + for index in 0..num_slices { + if index == 2 || index == 5 || index == 6 || index == 15 || index == 16 || index == 17 { + // Add some arbitrary irregularities to make the result look more interesting + continue; } + + let angle = std::f64::consts::PI * 0.9 * index as f64 / num_slices as f64; + let distance = 50.0; + let (sin, cos) = angle.sin_cos(); + cdt.insert(VertexType::new(sin * distance, cos * distance)) + .unwrap(); + } + + let mut handles = cdt.fixed_vertices().collect::>(); + + handles.push(handles[0]); + for vertices in handles.windows(2) { + cdt.add_constraint(vertices[0], vertices[1]); } - convert_triangulation(&cdt, &ConversionOptions::default()) + cdt } diff --git a/images/angle_limit_00.svg b/images/angle_limit_00.svg new file mode 100644 index 0000000..3dc60f8 --- /dev/null +++ b/images/angle_limit_00.svg @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/angle_limit_20.svg b/images/angle_limit_20.svg new file mode 100644 index 0000000..1dee8e0 --- /dev/null +++ b/images/angle_limit_20.svg @@ -0,0 +1,120 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/angle_limit_30.svg b/images/angle_limit_30.svg new file mode 100644 index 0000000..0203e78 --- /dev/null +++ b/images/angle_limit_30.svg @@ -0,0 +1,232 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/angle_limit_34.svg b/images/angle_limit_34.svg new file mode 100644 index 0000000..194bcfe --- /dev/null +++ b/images/angle_limit_34.svg @@ -0,0 +1,472 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/exclude_outer_faces_refined.svg b/images/exclude_outer_faces_refined.svg new file mode 100644 index 0000000..194bed8 --- /dev/null +++ b/images/exclude_outer_faces_refined.svg @@ -0,0 +1,274 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/exclude_outer_faces_unrefined.svg b/images/exclude_outer_faces_unrefined.svg new file mode 100644 index 0000000..2924c94 --- /dev/null +++ b/images/exclude_outer_faces_unrefined.svg @@ -0,0 +1,190 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/refined.svg b/images/refined.svg new file mode 100644 index 0000000..4798866 --- /dev/null +++ b/images/refined.svg @@ -0,0 +1,706 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/images/unrefined.svg b/images/unrefined.svg new file mode 100644 index 0000000..4e1b06b --- /dev/null +++ b/images/unrefined.svg @@ -0,0 +1,190 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/src/delaunay_core/dcel_operations.rs b/src/delaunay_core/dcel_operations.rs index f710203..aac6fef 100644 --- a/src/delaunay_core/dcel_operations.rs +++ b/src/delaunay_core/dcel_operations.rs @@ -532,7 +532,10 @@ where /// Splits `edge_handle` only one side. Used to split edges on the convex hull. /// /// Returns the newly inserted vertex and the two resulting parts of the split edge. -/// The returned edges will point away from the inserted vertex +/// +/// The returned edges will point in the same direction as the input edge. The first +/// returned edge has the same origin as the input edge, the second returned edge has the +/// same destination as the input edge. pub fn split_half_edge( dcel: &mut Dcel, edge_handle: FixedDirectedEdgeHandle, @@ -543,7 +546,7 @@ where UE: Default, F: Default, { - // Original quad: + // Original face: // // to // + @@ -583,6 +586,8 @@ where // nf = new face // e1, e2 = new edges // nv = new vertex + // + // This would return [e, e2] let edge = dcel.directed_edge(edge_handle); let v = edge.prev().from().fix(); let to = edge.to().fix(); @@ -660,7 +665,7 @@ where dcel.vertices[to.index()].out_edge = optional::some(e2.rev()); dcel.faces[f1.index()].adjacent_edge = optional::some(edge_handle); - (nv, [e2, edge_handle]) + (nv, [edge_handle, e2]) } /// Splits `edge_handle`, introducing 6 new half edges, two new faces and one diff --git a/src/delaunay_core/mod.rs b/src/delaunay_core/mod.rs index dcb377c..9dbb5fa 100644 --- a/src/delaunay_core/mod.rs +++ b/src/delaunay_core/mod.rs @@ -27,7 +27,7 @@ pub use hint_generator::{ LastUsedVertexHintGenerator, }; -pub use refinement::{AngleLimit, RefinementParameters}; +pub use refinement::{AngleLimit, RefinementParameters, RefinementResult}; pub use line_side_info::LineSideInfo; diff --git a/src/delaunay_core/refinement.rs b/src/delaunay_core/refinement.rs index b21be3b..216ccf9 100644 --- a/src/delaunay_core/refinement.rs +++ b/src/delaunay_core/refinement.rs @@ -12,12 +12,17 @@ use super::{ TriangulationExt, UndirectedEdgeHandle, }; +pub struct RefinementResult { + pub excluded_faces: HashSet>, + pub refinement_complete: bool, +} + /// Specifies the minimum allowed angle that should be kept after a refinement procedure. /// /// The refinement algorithm will attempt to keep the *minimum angle in the triangulation* greater than /// an angle limit specified with this struct. /// -/// *See method [refine](crate::ConstrainedDelaunayTriangulation::refine) * +/// *See [ConstrainedDelaunayTriangulation::refine], [RefinementParameters::with_angle_limit]* #[derive(Copy, Clone, PartialEq, PartialOrd)] pub struct AngleLimit { radius_to_shortest_edge_limit: f64, @@ -148,17 +153,38 @@ impl RefinementParameters { Self::default() } + /// Specifies the smallest allowed inner angle in a refined triangulation. + /// + /// The refinement algorithm will attempt to insert additional points (called steiner points) until the + /// minimum angle is larger than the angle bound specified by this function. + /// + /// Defaults to 30 degrees. + /// + /// Note that angle limits much larger than 30 degrees may not always terminate successfully - check + /// [RefinementResult::refinement_complete] to make sure that the angle limit did not lead to an overly strong + /// refined triangulation. + /// + /// # Examples of different angle limits + /// + /// + ///
0° (no angle refinement)20°30°34°
+ #[doc = concat!(include_str!("../../images/angle_limit_00.svg"), "")] + #[doc = concat!(include_str!("../../images/angle_limit_20.svg"), "")] + #[doc = concat!(include_str!("../../images/angle_limit_30.svg"), "")] + #[doc = concat!(include_str!("../../images/angle_limit_34.svg"), "
")] + /// + /// *See also [ConstrainedDelaunayTriangulation::refine]* pub fn with_angle_limit(mut self, angle_limit: AngleLimit) -> Self { self.angle_limit = angle_limit; self } - pub fn with_min_area(mut self, min_area: S) -> Self { + pub fn with_min_required_area(mut self, min_area: S) -> Self { self.min_area = Some(min_area); self } - pub fn with_max_area(mut self, max_area: S) -> Self { + pub fn with_max_allowed_area(mut self, max_area: S) -> Self { self.max_area = Some(max_area); self } @@ -168,11 +194,42 @@ impl RefinementParameters { self } + /// Prevents constraint edges from being split during refinement. + /// + /// By default, constraint edges may be split in order to restore the triangulation's Delaunay property. + /// The resulting two new edges will *become* new constraint edges, hence the original shape outlined by + /// constraint edges remains the same - no "gaps" or deviations are introduced. + /// + /// Enabling this option will, in general, reduce the quality of the resulting mesh - it is not necessarily + /// Delaunay anymore and faces adjacent to long constraint edges may violate the configured [AngleLimit]. pub fn keep_constraint_edges(mut self) -> Self { self.keep_constraint_edges = true; self } + /// Allows to exclude outer faces from the refinement process. + /// + /// This is useful if the constraint edges form a *closed shape* with a clearly defined inner and outer part. + /// Spade will determine inner and outer faces by identifying which faces can be reached from the outer face + /// without "crossing" a constraint edge, similar to a flood fill algorithm. + /// + /// Any holes in the triangulation will also be excluded. More specifically, any point with an odd winding number + /// is considered to be inner (see e.g. [Wikipedia](https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm)). + /// + /// Note that excluded faces may still be subdivided if a neighboring edge needs to be split. However, they will never be the + /// *cause* for a subdivision - their angle and area is ignored. + /// + /// The resulting outer faces of the triangulation are returned by the call to [refine](ConstrainedDelaunayTriangulation::refine), + /// see [RefinementResult::excluded_faces]. + /// + /// # Example + /// + /// + ///
UnrefinedRefined
+ #[doc = concat!(include_str!("../../images/exclude_outer_faces_unrefined.svg"), "",include_str!("../../images/exclude_outer_faces_refined.svg"), "
")] + /// + /// *A refinement operation configured to exclude outer faces. All colored faces are considered outer faces and are + /// ignored during refinement. Note that the inner part of the "A" shape forms a hole and is also excluded.* pub fn exclude_outer_faces( mut self, triangulation: &ConstrainedDelaunayTriangulation, @@ -249,7 +306,6 @@ impl RefinementParameters { let (_, radius2) = face.circumcircle(); let ratio2 = radius2 / length2; - let angle_limit = self.angle_limit.radius_to_shortest_edge_limit; if ratio2.into() > angle_limit * angle_limit { RefinementHint::ShouldRefine @@ -268,10 +324,94 @@ where L: HintGenerator<::Scalar>, ::Scalar: Float, { - pub fn refine( - &mut self, - mut parameters: RefinementParameters, - ) -> HashSet> { + /// Refines a triangulation by inserting additional points to improve the quality of its mesh. + /// + /// *Mesh quality*, when applied to constrained delaunay triangulations (CDT), usually refers to how skewed its + /// triangles are. A skewed triangle is a triangle with very large or very small (acute) inner angles. + /// Some applications (e.g. interpolation and finite element methods) perform poorly in the presence of skewed triangles. + /// + /// Refining by inserting additional points (called "steiner points") may increase the minimum angle. The given + /// [RefinementParameters] should be used to fine-tune the refinement behavior. + /// + /// # General usage + /// + /// The vertex type must implement `From>` - otherwise, Spade cannot construct new steiner points at a + /// certain location. The refinement itself happens *in place* and will result in a valid CDT. + /// + /// ``` + /// use spade::{ConstrainedDelaunayTriangulation, RefinementParameters, Point2, InsertionError, Triangulation}; + /// + /// fn get_refined_triangulation(vertices: Vec>) -> + /// Result>, InsertionError> + /// { + /// let mut cdt = ConstrainedDelaunayTriangulation::bulk_load(vertices)?; + /// let result = cdt.refine(RefinementParameters::default()); + /// if !result.refinement_complete { + /// panic!("Refinement failed - I should consider using different parameters.") + /// } + /// + /// return Ok(cdt) + /// } + /// ``` + /// + /// # Example image + /// + /// + /// + ///
UnrefinedRefined
+ #[doc = concat!(include_str!("../../images/unrefined.svg"), "",include_str!("../../images/refined.svg"), "
")] + /// + /// *A refinement example. The CDT on the left has some acute angles and skewed triangles. + /// The refined CDT on the right contains several additional points that prevents such triangles from appearing while keeping + /// all input vertices and constraint edges.* + /// + /// # Quality guarantees + /// + /// Refinement will ensure that the resulting triangulation fulfills a few properties: + /// - The triangulation's minimum angle will be larger than the angle specified by + /// [with_angle_limit](crate::RefinementParameters::with_angle_limit).
+ /// *Exception*: Acute input angles (small angles between initial constraint edges) cannot be resolved as the constraint edges + /// must be kept intact. The algorithm will, for the most part, leave those places unchanged.
+ /// *Exception*: The refinement will often not be able to increase the minimal angle much above 30 degrees as every newly + /// inserted steiner point may create additional skewed triangles. + /// - The refinement will fullfil the *Delaunay Property*: Every triangle's circumcenter will not contain any other vertex.
+ /// *Exception*: Refining with [keep_constraint_edges](crate::RefinementParameters::keep_constraint_edges) cannot restore + /// the Delaunay property if doing so would require splitting a constraint edge.
+ /// *Exception*: Refining with [exclude_outer_faces](crate::RefinementParameters::exclude_outer_faces) will not + /// ensure the Delaunay property of any outer face. + /// - Spade allows to specify a [maximum allowed triangle area](crate::RefinementParameters::with_max_allowed_area). + /// The algorithm will attempt to subdivide any triangle with an area larger than this, independent of its smallest angle. + /// - Spade allows to specify a [minimum required triangle area](crate::RefinementParameters::with_min_required_area). + /// The refinement will attempt to ignore any triangle with an area smaller than this parameter. This can prevent the + /// refinement algorithm from over-refining in some cases. + /// + /// # General limitations + /// + /// The algorithm may fail to terminate in some cases for a minimum angle limit larger than 30 degrees. Such a limit can + /// result in an endless loop: Every additionally inserted point creates more triangles that need to be refined. + /// + /// To prevent this, spade limits the number of additionally inserted steiner points + /// (see [RefinementParameters::with_max_additional_vertices]). However, this may leave the refinement in an incomplete state - + /// some areas of the input mesh may not have been triangulated at all, some will be overly refined. + /// Use [RefinementResult::refinement_complete] to identify if a refinement operation has succeeded without running out of + /// vertices. + /// + /// To prevent this from happening, consider either lowering the minimum angle limit + /// (see [RefinementParameters::with_angle_limit]) or introduce a + /// [minimum required area](RefinementParameters::with_min_required_area). + /// + /// # References + /// + /// This is an adaption of the classical refinement algorithms introduced by Jim Ruppert and Paul Chew. + /// + /// For a good introduction to the topic, refer to the slides from a short course at the thirteenth and fourteenth + /// International Meshing Roundtables (2005) by Jonathan Richard Shewchuk: + /// + /// + /// Wikipedia: + /// + /// + pub fn refine(&mut self, mut parameters: RefinementParameters) -> RefinementResult { use PositionInTriangulation::*; let mut excluded_faces = std::mem::take(&mut parameters.excluded_faces); @@ -294,18 +434,21 @@ where .map(|edge| edge.fix()), ); - let mut encroached_face_candidates: Vec<_> = self.fixed_inner_faces().collect(); + let mut encroached_face_candidates: VecDeque<_> = self.fixed_inner_faces().collect(); let num_initial_vertices: usize = self.num_vertices(); let num_additional_vertices = parameters .max_additional_vertices - .unwrap_or_else(|| num_initial_vertices * 4); + .unwrap_or_else(|| num_initial_vertices * 16); let max_allowed_vertices = num_initial_vertices + num_additional_vertices; - let is_original_vertex = |vertex: FixedVertexHandle| vertex.index() <= num_initial_vertices; + let is_original_vertex = |vertex: FixedVertexHandle| vertex.index() < num_initial_vertices; + + let mut refinement_complete = true; 'outer: loop { if self.num_vertices() >= max_allowed_vertices { + refinement_complete = false; break; } @@ -358,7 +501,7 @@ where } // Take the next triangle candidate - if let Some(face) = encroached_face_candidates.pop() { + if let Some(face) = encroached_face_candidates.pop_front() { if excluded_faces.contains(&face) { continue; } @@ -474,13 +617,16 @@ where ); } else if !forcibly_splitted_segments_buffer.is_empty() { // Revisit this face later - encroached_face_candidates.push(face.fix()); + encroached_face_candidates.push_back(face.fix()); } } else { break; } } - excluded_faces + RefinementResult { + excluded_faces, + refinement_complete, + } } fn is_fixed_edge(edge: UndirectedEdgeHandle, F>) -> bool { @@ -490,12 +636,13 @@ where fn resolve_encroachment( &mut self, encroached_segments_buffer: &mut VecDeque, - encroached_faces_buffer: &mut Vec>, + encroached_faces_buffer: &mut VecDeque>, encroached_edge: FixedUndirectedEdgeHandle, is_original_vertex: impl Fn(FixedVertexHandle) -> bool, excluded_faces: &mut HashSet>, ) { let segment = self.directed_edge(encroached_edge.as_directed()); + let [v0, v1] = segment.vertices(); let half = Into::::into(0.5f32); @@ -550,16 +697,19 @@ where } let (h1, h2) = (self.directed_edge(e1), self.directed_edge(e2)); + if is_left_side_protected == Some(true) { - excluded_faces.insert(h1.face().as_inner().unwrap().fix()); - excluded_faces.insert(h2.face().as_inner().unwrap().fix()); + excluded_faces.insert(h1.face().fix().as_inner().unwrap()); + excluded_faces.insert(h2.face().fix().as_inner().unwrap()); } if is_right_side_protected == Some(true) { - excluded_faces.insert(h1.rev().face().as_inner().unwrap().fix()); - excluded_faces.insert(h2.rev().face().as_inner().unwrap().fix()); + excluded_faces.insert(h1.rev().face().fix().as_inner().unwrap()); + excluded_faces.insert(h2.rev().face().fix().as_inner().unwrap()); } + self.legalize_vertex(new_vertex); + encroached_faces_buffer.extend( self.vertex(new_vertex) .out_edges() @@ -570,7 +720,7 @@ where encroached_segments_buffer.push_back(e1.as_undirected()); encroached_segments_buffer.push_back(e2.as_undirected()); - // Neighboring edges may have become encroached. Check if the need to be added to the encroached segment buffer. + // Neighboring edges may have become encroached. Check if they need to be added to the encroached segment buffer. encroached_segments_buffer.extend( [e1.rev(), e2] .into_iter() @@ -593,8 +743,8 @@ fn is_encroaching_edge( query_point.distance_2(edge_center) < radius_2 } -fn nearest_power_of_two(input: S) -> S { - input.log2().floor().exp2() +fn nearest_power_of_two(input: S) -> S { + input.log2().round().exp2() } #[cfg(test)] @@ -654,16 +804,18 @@ mod test { for i in 1..400u32 { let log = (i as f64).log2() as u32; - assert_eq!((1 << log) as f64, nearest_power_of_two(i as f64)) + let nearest = nearest_power_of_two(i as f64); + + assert!((1 << log) as f64 == nearest || (1 << (log + 1)) as f64 == nearest); } assert_eq!(0.25, nearest_power_of_two(0.25)); assert_eq!(0.25, nearest_power_of_two(0.27)); assert_eq!(0.5, nearest_power_of_two(0.5)); - assert_eq!(0.5, nearest_power_of_two(0.75)); - assert_eq!(1.0, nearest_power_of_two(1.5)); + assert_eq!(1.0, nearest_power_of_two(0.75)); + assert_eq!(2.0, nearest_power_of_two(1.5)); assert_eq!(2.0, nearest_power_of_two(2.5)); - assert_eq!(2.0, nearest_power_of_two(3.231)); + assert_eq!(4.0, nearest_power_of_two(3.231)); assert_eq!(4.0, nearest_power_of_two(4.0)); } @@ -695,14 +847,12 @@ mod test { cdt.add_constraint_edge(Point2::new(0.0, 0.0), Point2::new(4.0, 90.0))?; let refinement_parameters = RefinementParameters::new() - .with_max_additional_vertices(20) - .with_angle_limit(AngleLimit::from_radius_to_shortest_edge_ratio(0.7)); - - cdt.refine(refinement_parameters); + .with_max_additional_vertices(10) + .with_angle_limit(AngleLimit::from_radius_to_shortest_edge_ratio(1.0)); - assert_eq!(cdt.num_vertices(), 5); - assert_eq!(cdt.num_inner_faces(), 3); + let result = cdt.refine(refinement_parameters); + assert!(result.refinement_complete); cdt.cdt_sanity_check(); Ok(()) } diff --git a/src/delaunay_core/triangulation_ext.rs b/src/delaunay_core/triangulation_ext.rs index 57ca01c..2348e6c 100644 --- a/src/delaunay_core/triangulation_ext.rs +++ b/src/delaunay_core/triangulation_ext.rs @@ -69,18 +69,18 @@ pub trait TriangulationExt: Triangulation { InsertionResult::NewlyInserted(self.insert_into_face(face, t)) } OnEdge(edge) => { - if self.is_defined_legal(edge.as_undirected()) { - let (new_handle, split_parts) = self.insert_on_edge(edge, t); + let (new_handle, split_parts) = self.insert_on_edge(edge, t); + if self.is_defined_legal(edge.as_undirected()) { // If the edge is defined as legal the resulting edges must // be redefined as legal self.handle_legal_edge_split( split_parts.map(|edge| edge.as_undirected()), ); - InsertionResult::NewlyInserted(new_handle) - } else { - InsertionResult::NewlyInserted(self.insert_on_edge(edge, t).0) } + self.legalize_vertex(new_handle); + + InsertionResult::NewlyInserted(new_handle) } OnVertex(vertex) => { self.s_mut().update_vertex(vertex, t); @@ -311,15 +311,15 @@ pub trait TriangulationExt: Triangulation { new_vertex: Self::Vertex, ) -> (FixedVertexHandle, [FixedDirectedEdgeHandle; 2]) { let edge_handle = self.directed_edge(edge); - let result = if edge_handle.is_outer_edge() { - dcel_operations::split_half_edge(self.s_mut(), edge.rev(), new_vertex) + if edge_handle.is_outer_edge() { + let (new_vertex, [e0, e1]) = + dcel_operations::split_half_edge(self.s_mut(), edge.rev(), new_vertex); + (new_vertex, [e1.rev(), e0.rev()]) } else if edge_handle.rev().is_outer_edge() { dcel_operations::split_half_edge(self.s_mut(), edge, new_vertex) } else { dcel_operations::split_edge(self.s_mut(), edge, new_vertex) - }; - self.legalize_vertex(result.0); - result + } } fn legalize_vertex(&mut self, new_handle: FixedVertexHandle) { diff --git a/src/lib.rs b/src/lib.rs index adee610..6944f8f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,7 +33,7 @@ pub use crate::delaunay_core::math::{ pub use delaunay_core::{ AngleLimit, HierarchyHintGenerator, HierarchyHintGeneratorWithBranchFactor, HintGenerator, - LastUsedVertexHintGenerator, RefinementParameters, + LastUsedVertexHintGenerator, RefinementParameters, RefinementResult, }; pub use delaunay_core::LineSideInfo;