diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index fd1a2fd..dff5cf1 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -12,6 +12,7 @@ Thanks for your interest in contributing! - Rust stable toolchain. - Run tests: `cargo test`. +- Run the original full-size randomized suites: `ITRIANGLE_FULL_RANDOM_TESTS=1 cargo test`. - Run formatter: `cargo fmt`. - Run lints (optional): `cargo clippy`. diff --git a/DebugApp/.gitignore b/DebugApp/.gitignore new file mode 100644 index 0000000..eb5a316 --- /dev/null +++ b/DebugApp/.gitignore @@ -0,0 +1 @@ +target diff --git a/DebugApp/README.md b/DebugApp/README.md new file mode 100644 index 0000000..3432ec8 --- /dev/null +++ b/DebugApp/README.md @@ -0,0 +1,19 @@ +# iTriangle debug applications + +`uniform_grid` is an interactive visual check for `UniformTriangulatable` and +`IntUniformGrid`. It shows the editable input boundary, the separately colored +boundary resampled by `SliceContour`, lattice candidates after point +containment, points remaining after the `edge_length / 3` edge-clearance +filter, and the resulting Delaunay mesh. Every resampled boundary edge is at +most `edge_length` long. Optional centroid-net relaxation can be enabled in the +sidebar; its iteration limit defaults to 40. Convex decomposition and centroid +net overlays can be toggled independently on top of the Delaunay mesh. + +Run it from the iTriangle repository root: + +```sh +cargo run --manifest-path DebugApp/uniform_grid/Cargo.toml +``` + +Drag a yellow vertex to edit the active contour. Use the mouse wheel to zoom, +and the middle or right mouse button to pan. diff --git a/DebugApp/debug_ui/Cargo.toml b/DebugApp/debug_ui/Cargo.toml new file mode 100644 index 0000000..b5807bc --- /dev/null +++ b/DebugApp/debug_ui/Cargo.toml @@ -0,0 +1,8 @@ +[package] +name = "debug_ui" +version = "0.1.0" +edition = "2024" +publish = false + +[dependencies] +eframe = "^0.34.0" diff --git a/DebugApp/debug_ui/src/camera.rs b/DebugApp/debug_ui/src/camera.rs new file mode 100644 index 0000000..d1d9bb3 --- /dev/null +++ b/DebugApp/debug_ui/src/camera.rs @@ -0,0 +1,62 @@ +use eframe::egui::{Pos2, Rect, Vec2}; + +#[derive(Clone, Copy, Debug)] +pub struct Camera { + pub center: Pos2, + pub zoom: f32, + pub min_zoom: f32, + pub max_zoom: f32, +} + +impl Default for Camera { + fn default() -> Self { + Self { + center: Pos2::ZERO, + zoom: 1.0, + min_zoom: 0.02, + max_zoom: 256.0, + } + } +} + +impl Camera { + pub fn screen_from_world(&self, rect: Rect, world: Pos2) -> Pos2 { + let origin = rect.center(); + + Pos2::new( + origin.x + (world.x - self.center.x) * self.zoom, + origin.y - (world.y - self.center.y) * self.zoom, + ) + } + + pub fn world_from_screen(&self, rect: Rect, screen: Pos2) -> Pos2 { + let origin = rect.center(); + + Pos2::new( + self.center.x + (screen.x - origin.x) / self.zoom, + self.center.y - (screen.y - origin.y) / self.zoom, + ) + } + + pub fn world_delta_from_screen_delta(&self, delta: Vec2) -> Vec2 { + Vec2::new(delta.x / self.zoom, -delta.y / self.zoom) + } + + pub fn pan_by_screen_delta(&mut self, delta: Vec2) { + self.center -= self.world_delta_from_screen_delta(delta); + } + + pub fn zoom_at_screen_pos(&mut self, rect: Rect, screen_pos: Pos2, factor: f32) { + let before = self.world_from_screen(rect, screen_pos); + self.zoom = (self.zoom * factor).clamp(self.min_zoom, self.max_zoom); + let after = self.world_from_screen(rect, screen_pos); + self.center += before - after; + } + + pub fn visible_world_rect(&self, rect: Rect) -> Rect { + Rect::from_two_pos( + self.world_from_screen(rect, rect.left_bottom()), + self.world_from_screen(rect, rect.right_top()), + ) + } +} diff --git a/DebugApp/debug_ui/src/grid.rs b/DebugApp/debug_ui/src/grid.rs new file mode 100644 index 0000000..e9a8e4a --- /dev/null +++ b/DebugApp/debug_ui/src/grid.rs @@ -0,0 +1,117 @@ +use crate::camera::Camera; +use eframe::egui::{ + Align2, Color32, FontId, Painter, PointerButton, Rect, Response, Stroke, Ui, Vec2, +}; + +#[derive(Clone, Copy, Debug)] +pub struct Grid { + pub base_step: f32, + pub min_screen_step: f32, + pub max_screen_step: f32, + pub minor_stroke: Stroke, + pub major_stroke: Stroke, + pub axis_stroke: Stroke, + pub background: Color32, +} + +impl Default for Grid { + fn default() -> Self { + Self { + base_step: 16.0, + min_screen_step: 12.0, + max_screen_step: 48.0, + minor_stroke: Stroke::new(1.0_f32, Color32::from_gray(36)), + major_stroke: Stroke::new(1.0_f32, Color32::from_gray(56)), + axis_stroke: Stroke::new(1.5_f32, Color32::from_rgb(88, 102, 124)), + background: Color32::from_rgb(18, 20, 24), + } + } +} + +impl Grid { + pub fn handle_input(&self, ui: &Ui, response: &Response, rect: Rect, camera: &mut Camera) { + if response.hovered() { + let scroll_y = ui.input(|input| input.smooth_scroll_delta.y); + + if scroll_y.abs() > f32::EPSILON { + let factor = (scroll_y * 0.0018).exp(); + + if let Some(pointer_pos) = ui.input(|input| input.pointer.hover_pos()) { + camera.zoom_at_screen_pos(rect, pointer_pos, factor); + } + } + } + + if response.dragged_by(PointerButton::Middle) + || response.dragged_by(PointerButton::Secondary) + { + camera.pan_by_screen_delta(response.drag_delta()); + } + } + + pub fn paint(&self, painter: &Painter, rect: Rect, camera: &Camera) { + painter.rect_filled(rect, 0.0, self.background); + + let world = camera.visible_world_rect(rect); + let step = self.step_for_zoom(camera.zoom); + let min_x_index = (world.left() / step).floor() as i32 - 1; + let max_x_index = (world.right() / step).ceil() as i32 + 1; + let min_y_index = (world.top() / step).floor() as i32 - 1; + let max_y_index = (world.bottom() / step).ceil() as i32 + 1; + + for index in min_x_index..=max_x_index { + let x = index as f32 * step; + let stroke = self.stroke_for_index(index); + let a = camera.screen_from_world(rect, eframe::egui::pos2(x, world.bottom())); + let b = camera.screen_from_world(rect, eframe::egui::pos2(x, world.top())); + painter.line_segment([a, b], stroke); + } + + for index in min_y_index..=max_y_index { + let y = index as f32 * step; + let stroke = self.stroke_for_index(index); + let a = camera.screen_from_world(rect, eframe::egui::pos2(world.left(), y)); + let b = camera.screen_from_world(rect, eframe::egui::pos2(world.right(), y)); + painter.line_segment([a, b], stroke); + } + } + + fn step_for_zoom(&self, zoom: f32) -> f32 { + let mut step = self.base_step; + + while step * zoom < self.min_screen_step { + step *= 2.0; + } + + while step * zoom > self.max_screen_step { + step *= 0.5; + } + + step + } + + fn stroke_for_index(&self, index: i32) -> Stroke { + if index == 0 { + self.axis_stroke + } else if index.rem_euclid(5) == 0 { + self.major_stroke + } else { + self.minor_stroke + } + } +} + +pub fn paint_camera_readout(painter: &Painter, rect: Rect, camera: &Camera) { + let text = format!( + "center ({:.1}, {:.1}) zoom {:.2}x", + camera.center.x, camera.center.y, camera.zoom + ); + + painter.text( + rect.left_top() + Vec2::new(12.0, 10.0), + Align2::LEFT_TOP, + text, + FontId::monospace(12.0), + Color32::from_rgb(196, 202, 214), + ); +} diff --git a/DebugApp/debug_ui/src/lib.rs b/DebugApp/debug_ui/src/lib.rs new file mode 100644 index 0000000..ab59c03 --- /dev/null +++ b/DebugApp/debug_ui/src/lib.rs @@ -0,0 +1,4 @@ +pub mod camera; +pub mod grid; + +pub use eframe::egui; diff --git a/DebugApp/uniform_grid/Cargo.toml b/DebugApp/uniform_grid/Cargo.toml new file mode 100644 index 0000000..be03551 --- /dev/null +++ b/DebugApp/uniform_grid/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "uniform_grid" +version = "0.1.0" +edition = "2024" +publish = false + +[dependencies] +eframe = "^0.34.0" +debug_ui = { path = "../debug_ui" } +i_triangle = { path = "../../iTriangle" } +i_overlay = { path = "../../../iOverlay/iOverlay" } diff --git a/DebugApp/uniform_grid/src/examples.rs b/DebugApp/uniform_grid/src/examples.rs new file mode 100644 index 0000000..96c536e --- /dev/null +++ b/DebugApp/uniform_grid/src/examples.rs @@ -0,0 +1,89 @@ +pub type Point = [f32; 2]; +pub type Shape = Vec>; + +#[derive(Clone)] +pub struct GridExample { + pub name: &'static str, + pub shape: Shape, + pub edge_length: f32, +} + +pub fn load_examples() -> Vec { + vec![ + GridExample { + name: "concave polygon", + shape: vec![vec![ + [-220.0, -135.0], + [35.0, -175.0], + [220.0, -70.0], + [80.0, 5.0], + [205.0, 155.0], + [-25.0, 105.0], + [-205.0, 175.0], + [-115.0, 10.0], + ]], + edge_length: 42.0, + }, + GridExample { + name: "shape with hole", + shape: vec![ + vec![ + [-225.0, -170.0], + [225.0, -170.0], + [225.0, 170.0], + [-225.0, 170.0], + ], + // Clockwise winding makes this contour a hole for NonZero fill. + vec![[-95.0, -65.0], [-95.0, 75.0], [105.0, 75.0], [105.0, -65.0]], + ], + edge_length: 38.0, + }, + GridExample { + name: "narrow contour", + shape: vec![vec![ + [-235.0, -28.0], + [235.0, -28.0], + [235.0, 28.0], + [-235.0, 28.0], + ]], + edge_length: 72.0, + }, + GridExample { + name: "narrow passage", + shape: vec![vec![ + [-230.0, -175.0], + [-25.0, -175.0], + [-25.0, 45.0], + [25.0, 45.0], + [25.0, -175.0], + [230.0, -175.0], + [230.0, 175.0], + [25.0, 175.0], + [25.0, 95.0], + [-25.0, 95.0], + [-25.0, 175.0], + [-230.0, 175.0], + ]], + edge_length: 50.0, + }, + GridExample { + name: "two holes", + shape: vec![ + vec![ + [-240.0, -175.0], + [240.0, -175.0], + [215.0, 175.0], + [-215.0, 175.0], + ], + vec![ + [-150.0, -55.0], + [-150.0, 70.0], + [-45.0, 70.0], + [-45.0, -55.0], + ], + vec![[45.0, -80.0], [45.0, 45.0], [160.0, 45.0], [160.0, -80.0]], + ], + edge_length: 34.0, + }, + ] +} diff --git a/DebugApp/uniform_grid/src/main.rs b/DebugApp/uniform_grid/src/main.rs new file mode 100644 index 0000000..2188b53 --- /dev/null +++ b/DebugApp/uniform_grid/src/main.rs @@ -0,0 +1,854 @@ +mod examples; + +use crate::examples::{GridExample, Point, Shape as PolygonShape, load_examples}; +use debug_ui::{ + camera::Camera, + egui::{ + self, Color32, CursorIcon, Id, Painter, Pos2, Rect, Sense, Shape, Stroke, Ui, Vec2, + epaint::PathShape, + }, + grid::{Grid, paint_camera_readout}, +}; +use i_overlay::{ + core::{ + fill_rule::FillRule, overlay::IntOverlayOptions, point_location::IntPointContainment, + simplify::Simplify, + }, + i_float::{ + adapter::FloatPointAdapter, + float::rect::FloatRect, + int::{point::IntPoint, rect::IntRect}, + }, + i_shape::{ + float::adapter::PathToInt, + int::shape::{IntShape, IntShapes}, + }, +}; +use i_triangle::{ + float::{ + relax::RelaxationOptions, triangulation::Triangulation, uniform::UniformTriangulatable, + }, + tessellation::{split::SliceContour, uniform::IntUniformGrid}, +}; + +const PANEL_WIDTH: f32 = 270.0; +const TRIANGLE_HEIGHT_NUMERATOR: i64 = 28_378; +const TRIANGLE_HEIGHT_SHIFT: u32 = 15; + +struct MeshResult { + mesh: Triangulation, + convex_polygons: PolygonShape, + centroid_net: PolygonShape, + resampled_boundary: Vec, + contained_candidates: Vec, + clearance_points: Vec, + max_boundary_edge: f32, + relaxation: Option, +} + +struct RelaxationStats { + iterations: usize, + converged: bool, +} + +struct UniformGridApp { + camera: Camera, + grid: Grid, + examples: Vec, + active_example: usize, + edge_length: f32, + relax_enabled: bool, + relax_iterations: usize, + show_fill: bool, + show_triangles: bool, + show_boundary: bool, + show_resampled_boundary: bool, + show_contained_candidates: bool, + show_clearance_points: bool, + show_vertices: bool, + show_convex_decomposition: bool, + show_centroid_net: bool, + result: Result, +} + +impl Default for UniformGridApp { + fn default() -> Self { + let examples = load_examples(); + let mut app = Self { + camera: Camera::default(), + grid: Grid::default(), + edge_length: examples[0].edge_length, + relax_enabled: false, + relax_iterations: 40, + examples, + active_example: 0, + show_fill: true, + show_triangles: true, + show_boundary: true, + show_resampled_boundary: true, + show_contained_candidates: true, + show_clearance_points: true, + show_vertices: false, + show_convex_decomposition: false, + show_centroid_net: false, + result: Err("not calculated".to_owned()), + }; + app.refresh_result(); + app.fit_active_example(); + app + } +} + +impl eframe::App for UniformGridApp { + fn ui(&mut self, ui: &mut Ui, _frame: &mut eframe::Frame) { + egui::Panel::left("uniform_grid_panel") + .resizable(false) + .default_size(PANEL_WIDTH) + .frame(egui::Frame::default().fill(Color32::from_rgb(24, 27, 32))) + .show_inside(ui, |ui| self.show_controls(ui)); + + egui::CentralPanel::default() + .frame(egui::Frame::default().fill(Color32::from_rgb(18, 20, 24))) + .show_inside(ui, |ui| self.show_canvas(ui)); + } +} + +impl UniformGridApp { + fn show_controls(&mut self, ui: &mut Ui) { + ui.spacing_mut().item_spacing = Vec2::new(6.0, 6.0); + ui.add_space(8.0); + + let mut selected = None; + ui.label("Contour"); + for (index, example) in self.examples.iter().enumerate() { + if ui + .selectable_label(index == self.active_example, example.name) + .clicked() + { + selected = Some(index); + } + } + + if let Some(index) = selected { + self.select_example(index); + } + + ui.add_space(8.0); + ui.separator(); + ui.label("Uniform grid"); + + let edge_changed = ui + .add( + egui::DragValue::new(&mut self.edge_length) + .prefix("edge_length ") + .range(1.0..=500.0) + .speed(1.0), + ) + .changed(); + if edge_changed { + self.refresh_result(); + } + + let relax_changed = ui.checkbox(&mut self.relax_enabled, "relax mesh").changed(); + let iterations_changed = ui + .add_enabled( + self.relax_enabled, + egui::DragValue::new(&mut self.relax_iterations) + .prefix("relax iterations ") + .range(0..=1_000) + .speed(1), + ) + .changed(); + if relax_changed || iterations_changed { + self.refresh_result(); + } + + ui.add_space(8.0); + ui.separator(); + ui.label("Layers"); + ui.checkbox(&mut self.show_fill, "triangle fill"); + ui.checkbox(&mut self.show_triangles, "Delaunay edges"); + ui.checkbox(&mut self.show_boundary, "input boundary"); + ui.checkbox(&mut self.show_resampled_boundary, "resampled boundary"); + ui.checkbox( + &mut self.show_contained_candidates, + "lattice after containment", + ); + ui.checkbox(&mut self.show_clearance_points, "after edge clearance"); + ui.checkbox(&mut self.show_vertices, "all mesh vertices"); + ui.checkbox(&mut self.show_convex_decomposition, "convex decomposition"); + ui.checkbox(&mut self.show_centroid_net, "centroid net"); + + ui.add_space(8.0); + ui.separator(); + match &self.result { + Ok(result) => { + ui.label(format!("Contours: {}", self.active_example().shape.len())); + ui.label(format!( + "After containment: {}", + result.contained_candidates.len() + )); + ui.label(format!( + "After clearance: {}", + result.clearance_points.len() + )); + ui.label(format!( + "Removed near edges: {}", + result.contained_candidates.len() - result.clearance_points.len() + )); + ui.label(format!( + "Boundary samples: {}", + result + .resampled_boundary + .iter() + .flatten() + .map(Vec::len) + .sum::() + )); + ui.colored_label( + Color32::from_rgb(128, 212, 156), + format!( + "Max boundary edge: {:.3} ≤ {:.3}", + result.max_boundary_edge, self.edge_length + ), + ); + ui.label(format!( + "Clearance: edge_length / 3 = {:.3}", + self.edge_length / 3.0 + )); + ui.label(format!("Mesh vertices: {}", result.mesh.points.len())); + ui.label(format!("Triangles: {}", result.mesh.indices.len() / 3)); + if self.show_convex_decomposition { + ui.label(format!("Convex polygons: {}", result.convex_polygons.len())); + } + if self.show_centroid_net { + ui.label(format!("Centroid cells: {}", result.centroid_net.len())); + } + if let Some(relaxation) = &result.relaxation { + ui.label(format!( + "Relax: {} iterations, converged: {}", + relaxation.iterations, relaxation.converged + )); + } + } + Err(error) => { + ui.colored_label(Color32::from_rgb(240, 118, 118), error); + } + } + + if ui.button("Fit view").clicked() { + self.fit_active_example(); + } + if ui.button("Reset example").clicked() { + let index = self.active_example; + self.examples[index] = load_examples().remove(index); + self.edge_length = self.examples[index].edge_length; + self.refresh_result(); + self.fit_active_example(); + } + + ui.add_space(8.0); + ui.separator(); + ui.small("Drag yellow boundary vertices to edit."); + ui.small("Wheel: zoom. Right/middle drag: pan."); + } + + fn show_canvas(&mut self, ui: &mut Ui) { + let available_size = ui.available_size(); + let (response, painter) = ui.allocate_painter(available_size, Sense::click_and_drag()); + let rect = response.rect; + + self.grid + .handle_input(ui, &response, rect, &mut self.camera); + self.grid.paint(&painter, rect, &self.camera); + + if let Ok(result) = &self.result { + paint_mesh( + &painter, + rect, + &self.camera, + result, + self.show_fill, + self.show_triangles, + self.show_vertices, + ); + + if self.show_contained_candidates { + paint_points( + &painter, + rect, + &self.camera, + &result.contained_candidates, + 3.5, + Color32::from_rgba_unmultiplied(190, 130, 255, 135), + ); + } + + if self.show_clearance_points { + paint_points( + &painter, + rect, + &self.camera, + &result.clearance_points, + 3.0, + Color32::from_rgb(233, 92, 132), + ); + } + + if self.show_convex_decomposition { + paint_contours( + &painter, + rect, + &self.camera, + result.convex_polygons.iter(), + Stroke::new(2.25_f32, Color32::from_rgb(255, 156, 72)), + ); + } + + if self.show_centroid_net { + paint_contours( + &painter, + rect, + &self.camera, + result.centroid_net.iter(), + Stroke::new(1.75_f32, Color32::from_rgb(115, 225, 150)), + ); + } + } + + let camera = self.camera; + let show_boundary = self.show_boundary; + let changed = edit_shape( + ui, + &painter, + rect, + &camera, + &mut self.examples[self.active_example].shape, + show_boundary, + ); + if changed { + self.refresh_result(); + } + + if self.show_resampled_boundary + && let Ok(result) = &self.result + { + paint_contours( + &painter, + rect, + &self.camera, + result.resampled_boundary.iter().flatten(), + Stroke::new(1.25_f32, Color32::from_rgb(80, 225, 220)), + ); + paint_resampled_boundary_points( + &painter, + rect, + &self.camera, + &result.resampled_boundary, + ); + } + + paint_camera_readout(&painter, rect, &self.camera); + } + + fn active_example(&self) -> &GridExample { + &self.examples[self.active_example] + } + + fn select_example(&mut self, index: usize) { + self.active_example = index; + self.edge_length = self.examples[index].edge_length; + self.refresh_result(); + self.fit_active_example(); + } + + fn refresh_result(&mut self) { + let shape = self.active_example().shape.clone(); + let edge_length = self.edge_length; + let relax_enabled = self.relax_enabled; + let relax_iterations = self.relax_iterations; + + self.result = match std::panic::catch_unwind(move || { + build_mesh_result(&shape, edge_length, relax_enabled, relax_iterations) + }) { + Ok(result) => result, + Err(payload) => Err(panic_message(payload)), + }; + } + + fn fit_active_example(&mut self) { + let Some(bounds) = shape_bounds(&self.active_example().shape) else { + return; + }; + + self.camera.center = Pos2::new( + 0.5 * (bounds.min_x + bounds.max_x), + 0.5 * (bounds.min_y + bounds.max_y), + ); + } +} + +fn build_mesh_result( + shape: &PolygonShape, + edge_length: f32, + relax_enabled: bool, + relax_iterations: usize, +) -> Result { + if !edge_length.is_finite() || edge_length <= 0.0 { + return Err("edge_length must be finite and positive".to_owned()); + } + // This is the public high-level API under test. + let mut delaunay = shape.uniform_triangulate(edge_length); + let relaxation = relax_enabled.then(|| { + let result = delaunay.relax_mut(RelaxationOptions::new(relax_iterations)); + RelaxationStats { + iterations: result.iterations, + converged: result.converged, + } + }); + let convex_polygons = delaunay.to_convex_polygons(); + let centroid_net = delaunay.to_centroid_net(0.0); + let mesh = delaunay.to_triangulation::(); + + // Reproduce the public float wrapper's single conversion into the integer pipeline. + let rect = FloatRect::with_iter(shape.iter().flatten()) + .ok_or_else(|| "input shape is empty".to_owned())?; + let adapter = FloatPointAdapter::::new(rect); + let int_edge_length = adapter.round_len_to_int(edge_length); + if int_edge_length <= 1 { + return Err("edge_length is below integer adapter precision".to_owned()); + } + + let int_shape: IntShape = shape.iter().map(|path| path.to_int(&adapter)).collect(); + + // These are the same integer stages used by IntUniformTriangulatable. + let split_boundary = int_shape.slice_contour(int_edge_length as u64); + let resampled_boundary = vec![int_shape_to_float(&split_boundary, &adapter)]; + let max_boundary_edge = max_contour_edge(&resampled_boundary); + let normalized = + split_boundary.simplify(FillRule::NonZero, IntOverlayOptions::keep_all_points()); + let contained_int = lattice_after_containment(&normalized, int_edge_length as u64); + let clearance_int = normalized.uniform_grid(int_edge_length as u64); + let contained_candidates = int_points_to_float(&contained_int, &adapter); + let clearance_points = int_points_to_float(&clearance_int, &adapter); + + Ok(MeshResult { + mesh, + convex_polygons, + centroid_net, + resampled_boundary, + contained_candidates, + clearance_points, + max_boundary_edge, + relaxation, + }) +} + +fn lattice_after_containment(shapes: &IntShapes, edge_length: u64) -> Vec> { + let Some(rect) = IntRect::with_iter(shapes.iter().flatten().flatten()) else { + return Vec::new(); + }; + let Ok(step) = i64::try_from(edge_length) else { + return Vec::new(); + }; + if step <= 1 { + return Vec::new(); + } + + // Keep this generator identical to IntUniformGrid's lattice stage. Only the + // subsequent edge-clearance filter is intentionally omitted here. + let row_step = (step * TRIANGLE_HEIGHT_NUMERATOR + (1_i64 << (TRIANGLE_HEIGHT_SHIFT - 1))) + >> TRIANGLE_HEIGHT_SHIFT; + if row_step <= 0 { + return Vec::new(); + } + + let half_step = step / 2; + let min_x = i64::from(rect.min_x); + let max_x = i64::from(rect.max_x); + let max_y = i64::from(rect.max_y); + let mut candidates = Vec::new(); + let mut row = 0usize; + let mut y = i64::from(rect.min_y) + row_step / 2; + + while y < max_y { + let row_offset = if row & 1 == 0 { half_step } else { step }; + let mut x = min_x + row_offset; + + while x < max_x { + candidates.push(IntPoint::new( + i32::try_from(x).expect("lattice x stays inside i32 bounds"), + i32::try_from(y).expect("lattice y stays inside i32 bounds"), + )); + x += step; + } + + row += 1; + y += row_step; + } + + let contains = shapes.contains_points(&candidates); + candidates + .into_iter() + .zip(contains) + .filter_map(|(point, is_inside)| is_inside.then_some(point)) + .collect() +} + +fn int_shape_to_float( + shape: &IntShape, + adapter: &FloatPointAdapter, +) -> PolygonShape { + shape + .iter() + .map(|contour| int_points_to_float(contour, adapter)) + .collect() +} + +fn int_points_to_float( + points: &[IntPoint], + adapter: &FloatPointAdapter, +) -> Vec { + points + .iter() + .map(|point| adapter.int_to_float(point)) + .collect() +} + +fn paint_mesh( + painter: &Painter, + rect: Rect, + camera: &Camera, + result: &MeshResult, + show_fill: bool, + show_edges: bool, + show_vertices: bool, +) { + let edge_stroke = Stroke::new(1.0_f32, Color32::from_rgba_unmultiplied(106, 185, 255, 190)); + let fill = Color32::from_rgba_unmultiplied(73, 170, 255, 30); + + for triangle in result.mesh.indices.chunks_exact(3) { + let points = [triangle[0], triangle[1], triangle[2]].map(|index| { + let point = result.mesh.points[index as usize]; + camera.screen_from_world(rect, point_to_pos(point)) + }); + + if show_fill { + painter.add(Shape::convex_polygon( + points.to_vec(), + fill, + Stroke::new(0.0_f32, Color32::TRANSPARENT), + )); + } + if show_edges { + painter.add(Shape::closed_line(points.to_vec(), edge_stroke)); + } + } + + if show_vertices { + paint_points( + painter, + rect, + camera, + &result.mesh.points, + 2.0, + Color32::from_rgb(128, 212, 156), + ); + } +} + +fn edit_shape( + ui: &mut Ui, + painter: &Painter, + rect: Rect, + camera: &Camera, + contours: &mut [Vec], + show_boundary: bool, +) -> bool { + let mut changed = false; + + for (contour_index, contour) in contours.iter_mut().enumerate() { + for (point_index, point) in contour.iter_mut().enumerate() { + let screen = camera.screen_from_world(rect, point_to_pos(*point)); + let hit_rect = Rect::from_center_size(screen, Vec2::splat(18.0)); + let response = ui + .interact( + hit_rect, + Id::new("boundary_point") + .with(contour_index) + .with(point_index), + Sense::drag(), + ) + .on_hover_cursor(CursorIcon::Grab); + + if response.dragged() + && let Some(screen_position) = ui.input(|input| input.pointer.interact_pos()) + { + let world = camera.world_from_screen(rect, screen_position); + *point = [world.x, world.y]; + changed = true; + } + + let fill = if response.dragged() || response.hovered() { + Color32::WHITE + } else { + Color32::from_rgb(255, 206, 102) + }; + painter.circle( + camera.screen_from_world(rect, point_to_pos(*point)), + 4.5, + fill, + Stroke::new(1.0_f32, Color32::from_rgb(18, 20, 24)), + ); + } + } + + if show_boundary { + paint_contours( + painter, + rect, + camera, + contours.iter(), + Stroke::new(2.5_f32, Color32::from_rgb(255, 206, 102)), + ); + } + + changed +} + +fn paint_contours<'a>( + painter: &Painter, + rect: Rect, + camera: &Camera, + contours: impl Iterator>, + stroke: Stroke, +) { + for contour in contours { + if contour.len() < 2 { + continue; + } + let screen_points = contour + .iter() + .map(|point| camera.screen_from_world(rect, point_to_pos(*point))) + .collect(); + painter.add(PathShape::closed_line(screen_points, stroke)); + } +} + +fn paint_points( + painter: &Painter, + rect: Rect, + camera: &Camera, + points: &[Point], + radius: f32, + color: Color32, +) { + for point in points { + painter.circle_filled( + camera.screen_from_world(rect, point_to_pos(*point)), + radius, + color, + ); + } +} + +fn paint_resampled_boundary_points( + painter: &Painter, + rect: Rect, + camera: &Camera, + shapes: &[PolygonShape], +) { + for contour in shapes.iter().flatten() { + paint_points( + painter, + rect, + camera, + contour, + 2.75, + Color32::from_rgb(80, 225, 220), + ); + } +} + +fn max_contour_edge(shapes: &[PolygonShape]) -> f32 { + shapes + .iter() + .flatten() + .filter(|contour| contour.len() > 1) + .flat_map(|contour| { + contour + .iter() + .zip(contour.iter().cycle().skip(1)) + .take(contour.len()) + .map(|(a, b)| { + let dx = b[0] - a[0]; + let dy = b[1] - a[1]; + (dx * dx + dy * dy).sqrt() + }) + }) + .fold(0.0_f32, f32::max) +} + +#[derive(Clone, Copy)] +struct Bounds { + min_x: f32, + max_x: f32, + min_y: f32, + max_y: f32, +} + +fn shape_bounds(shape: &PolygonShape) -> Option { + let mut points = shape.iter().flatten(); + let first = *points.next()?; + let mut bounds = Bounds { + min_x: first[0], + max_x: first[0], + min_y: first[1], + max_y: first[1], + }; + + for point in points { + bounds.min_x = bounds.min_x.min(point[0]); + bounds.max_x = bounds.max_x.max(point[0]); + bounds.min_y = bounds.min_y.min(point[1]); + bounds.max_y = bounds.max_y.max(point[1]); + } + + Some(bounds) +} + +fn point_to_pos(point: Point) -> Pos2 { + Pos2::new(point[0], point[1]) +} + +fn panic_message(payload: Box) -> String { + let message = payload + .downcast_ref::<&str>() + .copied() + .or_else(|| payload.downcast_ref::().map(String::as_str)) + .unwrap_or("unknown panic"); + + format!("triangulation panic: {message}") +} + +fn main() -> eframe::Result<()> { + let native_options = eframe::NativeOptions { + viewport: egui::ViewportBuilder::default() + .with_title("Uniform Triangular Grid") + .with_inner_size(Vec2::new(1100.0, 780.0)), + ..eframe::NativeOptions::default() + }; + + eframe::run_native( + "Uniform Triangular Grid", + native_options, + Box::new(|_cc| Ok(Box::new(UniformGridApp::default()))), + ) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn all_examples_build_meshes() { + for example in load_examples() { + let result = build_mesh_result(&example.shape, example.edge_length, false, 40) + .unwrap_or_else(|error| panic!("{}: {error}", example.name)); + assert!( + !result.mesh.indices.is_empty(), + "{} has no triangles", + example.name + ); + assert!( + !result.convex_polygons.is_empty(), + "{} has no convex polygons", + example.name + ); + assert!( + !result.centroid_net.is_empty(), + "{} has no centroid cells", + example.name + ); + assert!( + result.max_boundary_edge <= example.edge_length + 0.001, + "{} has a resampled boundary edge {} longer than edge_length {}", + example.name, + result.max_boundary_edge, + example.edge_length + ); + assert!( + result + .clearance_points + .iter() + .all(|point| result.contained_candidates.contains(point)) + ); + } + } + + #[test] + fn hole_case_has_no_lattice_points_or_triangles_in_hole() { + let example = load_examples() + .into_iter() + .find(|example| example.name == "shape with hole") + .expect("hole example"); + let result = build_mesh_result(&example.shape, example.edge_length, false, 40) + .expect("hole example mesh"); + + for points in [&result.contained_candidates, &result.clearance_points] { + assert!(points.iter().all(|point| { + point[0] <= -95.0 || point[0] >= 105.0 || point[1] <= -65.0 || point[1] >= 75.0 + })); + } + + for triangle in result.mesh.indices.chunks_exact(3) { + let a = result.mesh.points[triangle[0] as usize]; + let b = result.mesh.points[triangle[1] as usize]; + let c = result.mesh.points[triangle[2] as usize]; + let centroid = [(a[0] + b[0] + c[0]) / 3.0, (a[1] + b[1] + c[1]) / 3.0]; + assert!( + centroid[0] <= -95.0 + || centroid[0] >= 105.0 + || centroid[1] <= -65.0 + || centroid[1] >= 75.0, + "triangle centroid lies inside the hole: {centroid:?}" + ); + } + } + + #[test] + fn narrow_case_falls_back_to_boundary_mesh() { + let example = load_examples() + .into_iter() + .find(|example| example.name == "narrow contour") + .expect("narrow example"); + let result = build_mesh_result(&example.shape, example.edge_length, false, 40) + .expect("narrow example mesh"); + + assert!(result.clearance_points.len() <= result.contained_candidates.len()); + assert!(!result.mesh.indices.is_empty()); + } + + #[test] + fn clearance_stage_removes_near_boundary_candidates() { + let example = load_examples() + .into_iter() + .find(|example| example.name == "shape with hole") + .expect("hole example"); + let result = build_mesh_result(&example.shape, example.edge_length, false, 40) + .expect("hole example mesh"); + + assert!(result.contained_candidates.len() > result.clearance_points.len()); + } + + #[test] + fn optional_relax_uses_requested_iteration_limit() { + let example = load_examples().remove(0); + let result = build_mesh_result(&example.shape, example.edge_length, true, 40) + .expect("relaxed example mesh"); + let relaxation = result.relaxation.expect("relaxation stats"); + + assert!(relaxation.iterations <= 40); + assert!(!result.mesh.indices.is_empty()); + } +} diff --git a/editor/Cargo.toml b/editor/Cargo.toml index 4a586bc..26355ad 100644 --- a/editor/Cargo.toml +++ b/editor/Cargo.toml @@ -16,7 +16,7 @@ log = "0.4.22" console_log = "^1.0.0" console_error_panic_hook = "^0" -#i_mesh = { path = "../../iMesh/iMesh" } -#i_triangle = { path = "../iTriangle", default-features = true, features = ["serde"] } -i_mesh = "^0.4.0" -i_triangle = { version = "^0.37.0", features = ["serde"] } \ No newline at end of file +i_mesh = { path = "../../iMesh/iMesh" } +i_triangle = { path = "../iTriangle", default-features = true, features = ["serde"] } +#i_mesh = "^0.4.0" +#i_triangle = { version = "^0.37.0", features = ["serde"] } \ No newline at end of file diff --git a/iTriangle/Cargo.toml b/iTriangle/Cargo.toml index 65fb888..6038a7a 100644 --- a/iTriangle/Cargo.toml +++ b/iTriangle/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "i_triangle" -version = "0.46.0" +version = "0.47.0" edition = "2021" authors = ["Nail Sharipov "] description = "Polygon Triangulation Library: Efficient Delaunay Triangulation for Complex Shapes." @@ -20,7 +20,7 @@ serde = ["dep:serde", "i_overlay/serde"] [dependencies] serde = { version = "^1.0", default-features = false, features = ["derive"], optional = true } -i_overlay = { version = "^8.0.0"} +i_overlay = { version = "^8.1.0"} i_tree = "~0.19.0" i_key_sort = "~0.11.0" diff --git a/iTriangle/README.md b/iTriangle/README.md index cc14cb8..19b1390 100644 --- a/iTriangle/README.md +++ b/iTriangle/README.md @@ -7,7 +7,7 @@ [![codecov](https://codecov.io/gh/iShape-Rust/iTriangle/branch/main/graph/badge.svg)](https://codecov.io/gh/iShape-Rust/iTriangle) [![license](https://img.shields.io/crates/l/i_triangle.svg)](https://crates.io/crates/i_triangle) -iTriangle is a high-performance 2D polygon triangulation library for Rust. It turns real-world polygon input into triangle meshes, including shapes with holes, self-intersections, and mixed winding. The public API accepts `f32`/`f64` floating-point data and `i16`/`i32`/`i64` integer data, while the triangulation pipeline runs through a deterministic integer core for stable, reproducible output. +iTriangle provides reliable, deterministic 2D polygon triangulation for complex shapes with holes and self-intersections. *For detailed performance benchmarks, check out the* [Performance Comparison](https://ishape-rust.github.io/iShape-js/triangle/performance/performance.html) @@ -37,10 +37,12 @@ iTriangle is a high-performance 2D polygon triangulation library for Rust. It tu ## Features -- **Sweep-line Triangulation** - Fast and simple triangulation of polygons with or without holes. +- **Sweep-line Triangulation** - Fast polygon triangulation with `O(n log n)` complexity. - **Delaunay Triangulation** - Efficient and robust implementation for generating Delaunay triangulations. - **Self-Intersection Handling** – Fully supports self-intersecting polygons with automatic resolution. - **Adaptive Tessellation** - Refine Delaunay triangles using circumcenters for better shape quality. +- **Uniform Delaunay Triangulation** - Generate boundary-conforming meshes with a predictable target edge length. +- **Mesh Relaxation** - Improve mesh quality by moving interior vertices toward their centroid-net cell centroids. - **Convex Decomposition** - Convert triangulation into convex polygons. - **Centroidal Polygon Net**: Build per-vertex dual polygons using triangle centers and edge midpoints. - **Steiner Points**: Add custom inner points to influence triangulation. @@ -66,7 +68,7 @@ Add to your `Cargo.toml`: ```toml [dependencies] -i_triangle = "0.45" +i_triangle = "0.47" ``` Minimal example: @@ -194,6 +196,46 @@ println!("centroids: {:?}", centroids); > 💡 Output: Triangle indices and vertices, where all triangles oriented in a **counter-clockwise** direction. +### Uniform Mesh Relaxation and Centroid Net + +Use uniform triangulation when you need triangles with a predictable target edge +length. Relaxation moves only interior vertices toward their centroid-net cell +centers; boundary vertices stay fixed, and Delaunay edges are restored after each +iteration. + +```rust +use i_triangle::float::relax::RelaxationOptions; +use i_triangle::float::uniform::UniformTriangulatable; +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::float::single::SingleFloatOverlay; + +let contours = vec![ + vec![[0.0, 0.0], [12.0, 0.0], [12.0, 8.0], [0.0, 8.0]], + vec![[4.0, 2.0], [8.0, 2.0], [8.0, 6.0], [4.0, 6.0]], +]; +let empty: Vec> = Vec::new(); +let shapes = contours.overlay(&empty, OverlayRule::Union, FillRule::EvenOdd); +let shape = &shapes[0]; + +let mut delaunay = shape.uniform_triangulate(1.0); +let relaxation = delaunay.relax_mut(RelaxationOptions::new(24)); + +let triangles = delaunay.to_triangulation::(); +let centroid_net = delaunay.to_centroid_net(0.0); + +println!("triangles: {}", triangles.indices.len() / 3); +println!("centroid cells: {}", centroid_net.len()); +println!("relaxation: {relaxation:?}"); +``` + +| Relaxed Uniform Delaunay Mesh | Centroid Net | +| --- | --- | +| | | + +The complete reproducible renderer, including the eagle contours, is available in +[`examples/eagle_svg.rs`](examples/eagle_svg.rs). + ### Triangulating Multiple Shapes Efficiently If you need to triangulate many shapes, it is more efficient to use `Triangulator`. @@ -311,10 +353,6 @@ Benchmarks and interactive demos are available here: | --- | --- | --- | | | | | -| Tessellation | Centroid Net | | -| --- | --- | --- | -| | | | - ## Contributing See `CONTRIBUTING.md` for development setup, tests, and PR guidelines. diff --git a/iTriangle/examples/eagle_svg.rs b/iTriangle/examples/eagle_svg.rs new file mode 100644 index 0000000..c0da252 --- /dev/null +++ b/iTriangle/examples/eagle_svg.rs @@ -0,0 +1,372 @@ +use i_triangle::float::relax::{RelaxationOptions, RelaxationResult}; +use i_triangle::float::triangulation::Triangulation; +use i_triangle::float::uniform::UniformTriangulatable; +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::float::single::SingleFloatOverlay; +use std::collections::BTreeSet; +use std::fmt::Write as _; +use std::fs; +use std::io; +use std::path::Path; + +type Point = [f64; 2]; + +const EDGE_LENGTH: f64 = 3.5; +const RELAX_ITERATIONS: usize = 24; +const VIEW_BOX: &str = "0 0 132.29167 67.97879"; + +fn main() -> io::Result<()> { + let contours: Vec> = EAGLE_CONTOURS + .iter() + .map(|contour| contour.to_vec()) + .collect(); + let empty: Vec> = Vec::new(); + let shapes = contours.overlay(&empty, OverlayRule::Union, FillRule::EvenOdd); + assert_eq!(shapes.len(), 1, "the eagle must normalize to one shape"); + + let shape = &shapes[0]; + let mut delaunay = shape.uniform_triangulate(EDGE_LENGTH); + let relaxation = delaunay.relax_mut(RelaxationOptions::new(RELAX_ITERATIONS)); + let centroid_net = delaunay.to_centroid_net(0.0); + let triangulation = delaunay.to_triangulation::(); + + let root = Path::new(env!("CARGO_MANIFEST_DIR")).join("readme"); + fs::write( + root.join("eagle_tessellation.svg"), + tessellation_svg(shape, &triangulation, relaxation), + )?; + fs::write( + root.join("eagle_centroid.svg"), + centroid_svg(shape, ¢roid_net, relaxation), + )?; + + println!( + "generated: {} points, {} triangles, {} centroid cells, relax={:?}", + triangulation.points.len(), + triangulation.indices.len() / 3, + centroid_net.len(), + relaxation + ); + Ok(()) +} + +fn tessellation_svg( + shape: &[Vec], + triangulation: &Triangulation, + relaxation: RelaxationResult, +) -> String { + let mut svg = svg_header( + "Relaxed uniform Delaunay tessellation", + relaxation, + triangulation.points.len(), + ); + append_shape_fill(&mut svg, shape); + + let mut edges = BTreeSet::new(); + for triangle in triangulation.indices.chunks_exact(3) { + for (a, b) in [ + (triangle[0], triangle[1]), + (triangle[1], triangle[2]), + (triangle[2], triangle[0]), + ] { + edges.insert(if a < b { (a, b) } else { (b, a) }); + } + } + + svg.push_str(" \n", + ); + append_shape_outline(&mut svg, shape); + svg.push_str("\n"); + svg +} + +fn centroid_svg( + shape: &[Vec], + centroid_net: &[Vec], + relaxation: RelaxationResult, +) -> String { + let mut svg = svg_header( + "Centroid net of a relaxed uniform Delaunay mesh", + relaxation, + centroid_net.len(), + ); + append_shape_fill(&mut svg, shape); + + svg.push_str(" \n", + ); + append_shape_outline(&mut svg, shape); + svg.push_str("\n"); + svg +} + +fn svg_header(title: &str, relaxation: RelaxationResult, element_count: usize) -> String { + format!( + "\n\ + \n\ + \n\ + \n\ + {title}\n", + relaxation.iterations, relaxation.converged, element_count + ) +} + +fn append_shape_fill(svg: &mut String, shape: &[Vec]) { + svg.push_str(" \n"); +} + +fn append_shape_outline(svg: &mut String, shape: &[Vec]) { + svg.push_str(" \n", + ); +} + +fn append_contours(svg: &mut String, contours: &[Vec]) { + for contour in contours { + let Some(&first) = contour.first() else { + continue; + }; + append_move(svg, first); + for &point in &contour[1..] { + append_line(svg, point); + } + svg.push_str("Z\n "); + } +} + +fn append_segment(svg: &mut String, a: Point, b: Point) { + append_move(svg, a); + append_line(svg, b); + svg.push_str("\n "); +} + +fn append_move(svg: &mut String, point: Point) { + write!(svg, "M{:.5},{:.5} ", point[0], point[1]).unwrap(); +} + +fn append_line(svg: &mut String, point: Point) { + write!(svg, "L{:.5},{:.5}", point[0], point[1]).unwrap(); +} + +// Extracted once from the compound silhouette path in the historical SVG. +// The first contour is the eagle outline; the remaining contours are holes. +const EAGLE_CONTOURS: &[&[Point]] = &[ + &[ + [12.64125100, 0.21080764], + [12.64125100, 5.59553420], + [14.26289800, 9.36515110], + [16.42475200, 11.51914500], + [12.91135400, 10.71133300], + [9.39847060, 9.90352090], + [5.34538210, 8.28789720], + [1.29177810, 6.67278770], + [3.99435250, 9.63407880], + [8.31754400, 13.67313800], + [11.02011900, 14.21150700], + [6.69641170, 14.74987700], + [0.21085147, 14.74987700], + [2.91342590, 17.98112400], + [8.31754400, 19.05786400], + [1.83249920, 20.67348800], + [5.07527910, 22.82748000], + [10.47939800, 22.82748000], + [5.61548470, 24.98147400], + [9.39847060, 26.59658400], + [13.18197200, 26.59658400], + [10.47939800, 28.75057700], + [14.26289800, 29.28946000], + [18.04588400, 28.75057700], + [15.34382600, 30.90457100], + [20.20773800, 31.44345500], + [18.58609000, 33.05856400], + [24.26082600, 32.52019500], + [22.36959100, 36.28981200], + [27.77370900, 35.21255700], + [26.69329800, 37.90492000], + [32.09741500, 36.28981200], + [31.01649000, 38.98217400], + [34.79999200, 37.36655000], + [34.25926900, 40.59779800], + [38.04225400, 38.98217400], + [37.50205000, 41.67453700], + [41.28503500, 39.52054300], + [41.28503500, 42.75179100], + [44.52781500, 40.59779800], + [43.98761000, 43.82853000], + [47.77059500, 41.67453700], + [48.85152000, 39.52054300], + [48.85152000, 43.29016100], + [51.55358000, 41.13616700], + [52.09430100, 39.52054300], + [52.63450800, 42.75179100], + [54.25615600, 41.67453700], + [54.79636000, 40.05891300], + [55.87728800, 42.75179100], + [56.41749300, 46.52089500], + [52.63450800, 48.40595900], + [48.85152000, 50.29051200], + [42.90668200, 50.82888200], + [44.52781500, 54.06012700], + [47.22987300, 54.59849800], + [45.60874100, 55.67523700], + [47.77059500, 57.82923100], + [52.09430100, 58.36759900], + [48.85152000, 59.44485500], + [52.09430100, 62.13721700], + [55.87728800, 61.59884700], + [53.71543300, 63.75284100], + [57.49842100, 65.36846500], + [61.28192000, 64.82958200], + [59.66027300, 67.52245800], + [62.36284800, 67.52245800], + [66.14583300, 65.36846500], + [69.92881800, 67.52245800], + [72.63139300, 67.52245800], + [71.00974600, 64.82958200], + [74.79324600, 65.36846500], + [78.57623300, 63.75284100], + [76.41437800, 61.59884700], + [80.19736600, 62.13721700], + [83.44014600, 59.44485500], + [80.19736600, 58.36759900], + [84.52107200, 57.82923100], + [86.68292600, 55.67523700], + [85.06179300, 54.59849800], + [87.76385100, 54.06012700], + [89.38498500, 50.82888200], + [83.44014600, 50.29051200], + [79.65715900, 48.40595900], + [75.87365900, 46.52089500], + [76.41437800, 42.75179100], + [77.49530600, 40.05891300], + [78.03551100, 41.67453700], + [79.65715900, 42.75179100], + [80.19736600, 39.52054300], + [80.73808700, 41.13616700], + [83.44014600, 43.29016100], + [83.44014600, 39.52054300], + [84.52107200, 41.67453700], + [88.30405600, 43.82853000], + [87.76385100, 40.59779800], + [91.00663200, 42.75179100], + [91.00663200, 39.52054300], + [94.78961700, 41.67453700], + [94.24941200, 38.98217400], + [98.03239700, 40.59779800], + [97.49167500, 37.36655000], + [101.27518000, 38.98217400], + [100.19425000, 36.28981200], + [105.59836000, 37.90492000], + [104.51796000, 35.21255700], + [109.92207000, 36.28981200], + [108.03084000, 32.52019500], + [113.70558000, 33.05856400], + [112.08393000, 31.44345500], + [116.94784000, 30.90457100], + [114.24578000, 28.75057700], + [118.02876000, 29.28946000], + [121.81227000, 28.75057700], + [119.10970000, 26.59658400], + [122.89319000, 26.59658400], + [126.67619000, 24.98147400], + [121.81227000, 22.82748000], + [127.21639000, 22.82748000], + [130.45916000, 20.67348800], + [123.97412000, 19.05786400], + [129.37824000, 17.98112400], + [132.08082000, 14.74987700], + [125.59525000, 14.74987700], + [121.27155000, 14.21150700], + [123.97412000, 13.67313800], + [128.29731000, 9.63407880], + [130.99988000, 6.67278770], + [126.94628000, 8.28789720], + [122.89319000, 9.90352090], + [119.38031000, 10.71133300], + [115.86691000, 11.51914500], + [118.02876000, 9.36515110], + [119.65042000, 5.59553420], + [119.65042000, 0.21080764], + [116.94784000, 5.59553420], + [110.46279000, 8.28789720], + [103.97724000, 10.98026000], + [106.67929000, 8.28789720], + [106.67929000, 5.05716430], + [104.51796000, 7.74952740], + [99.65353000, 9.63407880], + [94.78961700, 11.51914500], + [91.36694000, 14.39096400], + [87.85405900, 17.17331300], + [84.43086400, 20.04513200], + [81.00818800, 22.82748000], + [77.49530600, 24.44259100], + [74.25252700, 24.44259100], + [71.55046800, 23.90422100], + [69.38861400, 22.28859700], + [66.14583300, 21.75022700], + [62.90305300, 22.28859700], + [60.74119900, 23.90422100], + [58.03914000, 24.44259100], + [54.79636000, 24.44259100], + [51.28347900, 22.82748000], + [47.86080200, 20.04513200], + [44.43760800, 17.17331300], + [40.92472600, 14.39096400], + [37.50205000, 11.51914500], + [32.63813700, 9.63407880], + [27.77370900, 7.74952740], + [25.61237100, 5.05716430], + [25.61237100, 8.28789720], + [28.31443100, 10.98026000], + [21.82887000, 8.28789720], + [15.34382600, 5.59553420], + ], + &[ + [62.90305300, 28.75057700], + [64.52470100, 30.36620100], + [63.98397900, 31.44345500], + [62.90305300, 30.36620100], + ], + &[ + [69.38861400, 28.75057700], + [69.38861400, 30.36620100], + [68.30768700, 31.44345500], + [67.76696600, 30.36620100], + ], + &[ + [66.14583300, 30.36620100], + [67.22675900, 30.90457100], + [67.76696600, 31.98182500], + [68.30768700, 32.52019500], + [70.46953900, 32.52019500], + [68.30768700, 33.59693300], + [67.22675900, 34.67418800], + [66.14583300, 37.36655000], + [65.06490800, 34.67418800], + [63.98397900, 33.59693300], + [61.82212700, 32.52019500], + [63.98397900, 32.52019500], + [64.52470100, 31.98182500], + [65.06490800, 30.90457100], + ], +]; diff --git a/iTriangle/readme/eagle_centroid.svg b/iTriangle/readme/eagle_centroid.svg index ac5bd65..8305649 100644 --- a/iTriangle/readme/eagle_centroid.svg +++ b/iTriangle/readme/eagle_centroid.svg @@ -1,1547 +1,577 @@ - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Centroid net of a relaxed uniform Delaunay mesh + + + diff --git a/iTriangle/readme/eagle_tessellation.svg b/iTriangle/readme/eagle_tessellation.svg index 71f594a..6c504ec 100644 --- a/iTriangle/readme/eagle_tessellation.svg +++ b/iTriangle/readme/eagle_tessellation.svg @@ -1,3849 +1,1403 @@ - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Relaxed uniform Delaunay tessellation + + + diff --git a/iTriangle/src/advanced/centroid.rs b/iTriangle/src/advanced/centroid.rs index bf46cac..1fa0815 100644 --- a/iTriangle/src/advanced/centroid.rs +++ b/iTriangle/src/advanced/centroid.rs @@ -110,7 +110,7 @@ impl IntTriangle { } #[inline] - fn left_neighbor_and_mid_edge(&self, vertex_index: usize) -> (usize, IntPoint) { + pub(crate) fn left_neighbor_and_mid_edge(&self, vertex_index: usize) -> (usize, IntPoint) { if self.vertices[0].index == vertex_index { let neighbor = self.neighbors[1]; let mid = middle(self.vertices[0].point, self.vertices[2].point); @@ -127,7 +127,7 @@ impl IntTriangle { } #[inline] - fn center(&self) -> IntPoint { + pub(crate) fn center(&self) -> IntPoint { let a = self.vertices[0].point; let b = self.vertices[1].point; let c = self.vertices[2].point; diff --git a/iTriangle/src/advanced/delaunay.rs b/iTriangle/src/advanced/delaunay.rs index 0ac4266..9d14cf4 100644 --- a/iTriangle/src/advanced/delaunay.rs +++ b/iTriangle/src/advanced/delaunay.rs @@ -183,7 +183,7 @@ impl DelaunayRefine for [IntTriangle] { } } -struct DelaunayCondition; +pub(crate) struct DelaunayCondition; impl DelaunayCondition { // if p is inside circumscribe circle of a, b, c return false @@ -191,7 +191,7 @@ impl DelaunayCondition { // return true if triangle satisfied condition and do not need flip triangles // more detail explanation and demo https://ishape-rust.github.io/iShape-js/triangle/delaunay.html #[inline] - fn is_flip_not_required( + pub(crate) fn is_flip_not_required( p: IntPoint, a: IntPoint, b: IntPoint, @@ -313,6 +313,7 @@ mod tests { use crate::geom::point::IndexPoint; use crate::geom::triangle::IntTriangle; use crate::int::triangulatable::IntTriangulatable; + use crate::test_util::random_cases; use alloc::vec; use i_overlay::core::fill_rule::FillRule; use i_overlay::core::overlay::IntOverlayOptions; @@ -455,7 +456,7 @@ mod tests { #[test] fn test_random_0() { - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { let shape = vec![random(8, 5)]; if let Some(first) = shape @@ -474,7 +475,7 @@ mod tests { #[test] fn test_random_1() { - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { let shape = vec![random(8, 12)]; if let Some(first) = shape @@ -493,7 +494,7 @@ mod tests { #[test] fn test_random_2() { - for _ in 0..2_000 { + for _ in 0..random_cases(2_000) { let main = random(50, 20); let mut shape = vec![main]; for _ in 0..10 { diff --git a/iTriangle/src/advanced/mod.rs b/iTriangle/src/advanced/mod.rs index 8a5c284..ee62ede 100644 --- a/iTriangle/src/advanced/mod.rs +++ b/iTriangle/src/advanced/mod.rs @@ -3,4 +3,5 @@ pub mod buffer; pub mod centroid; pub mod convex; pub mod delaunay; +pub mod relax; pub mod triangulation; diff --git a/iTriangle/src/advanced/relax.rs b/iTriangle/src/advanced/relax.rs new file mode 100644 index 0000000..e6f42fd --- /dev/null +++ b/iTriangle/src/advanced/relax.rs @@ -0,0 +1,685 @@ +use crate::advanced::buffer::DelaunayBuffer; +use crate::advanced::delaunay::{DelaunayRefine, IntDelaunay}; +use crate::geom::triangle::IntTriangle; +use alloc::vec; +use alloc::vec::Vec; +use i_overlay::i_float::int::number::int::IntNumber; +use i_overlay::i_float::int::number::product_uint::UIntProduct; +use i_overlay::i_float::int::number::signed_product::SignedProduct; +use i_overlay::i_float::int::number::uint::UIntNumber; +use i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_overlay::i_float::int::point::IntPoint; +use i_overlay::i_float::int::vector::IntVector; +use i_overlay::i_float::triangle::Triangle; + +/// Configuration for centroid-net relaxation of an integer Delaunay mesh. +/// +/// Boundary vertices, including vertices on hole boundaries, remain fixed. +/// Interior vertices move synchronously toward the area centroids of their +/// centroid-net cells. Every displacement is conservatively limited to less +/// than one quarter of the minimum altitude of each incident triangle before +/// the Delaunay property is restored with edge flips. +#[derive(Clone, Copy)] +pub struct RelaxationOptions { + /// Maximum number of relaxation iterations. + pub max_iterations: usize, + + /// Absolute convergence tolerance in input coordinate units. + /// + /// Relaxation stops before the next move when every unconstrained + /// centroid displacement is no greater than this value. Set it to zero to + /// stop only when integer rounding leaves no vertex to move. + pub tolerance: I::WideUInt, +} + +impl RelaxationOptions { + /// Creates options with zero convergence tolerance. + #[inline] + pub fn new(max_iterations: usize) -> Self { + Self { + max_iterations, + tolerance: I::WideUInt::ZERO, + } + } + + /// Sets the absolute convergence tolerance. + #[inline] + pub fn with_tolerance(mut self, tolerance: I::WideUInt) -> Self { + self.tolerance = tolerance; + self + } +} + +impl Default for RelaxationOptions { + #[inline] + fn default() -> Self { + Self::new(8) + } +} + +/// Outcome of an in-place relaxation run. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct RelaxationResult { + /// Number of completed vertex-movement iterations. + pub iterations: usize, + + /// `true` when relaxation stopped by tolerance or because no integer + /// vertex position changed; `false` when it reached the iteration limit. + pub converged: bool, +} + +struct RelaxationBuffer { + fixed: Vec, + starts: Vec, + proposals: Vec>, + order: Vec, + cell: Vec>, + ring: Vec, + delaunay: DelaunayBuffer, +} + +impl RelaxationBuffer { + fn new(delaunay: &IntDelaunay) -> Self { + let point_count = delaunay.points.len(); + let mut fixed = vec![false; point_count]; + + for triangle in &delaunay.triangles { + for edge in 0..3 { + if triangle.neighbors[edge] >= delaunay.triangles.len() { + fixed[triangle.vertices[(edge + 1) % 3].index] = true; + fixed[triangle.vertices[(edge + 2) % 3].index] = true; + } + } + } + + Self { + fixed, + starts: vec![usize::MAX; point_count], + proposals: delaunay.points.clone(), + order: (0..point_count).collect(), + cell: Vec::with_capacity(16), + ring: Vec::with_capacity(8), + delaunay: DelaunayBuffer::new(), + } + } + + fn rebuild_starts(&mut self, triangles: &[IntTriangle]) { + self.starts.fill(usize::MAX); + for (triangle_index, triangle) in triangles.iter().enumerate() { + for vertex in &triangle.vertices { + self.starts[vertex.index] = triangle_index; + } + } + } +} + +impl IntDelaunay { + /// Relaxes interior vertices toward their centroid-net area centroids. + /// + /// Boundary vertices are detected from triangle edges without a neighbor + /// and remain fixed. Vertex indices and the number of vertices are + /// preserved. The Delaunay property is restored with edge flips after each + /// synchronous move. + /// + /// Internal constrained edges are not retained because [`IntDelaunay`] + /// does not store which interior edges are constraints. + #[must_use] + #[inline] + pub fn relax(mut self, options: RelaxationOptions) -> Self { + self.relax_mut(options); + self + } + + /// In-place version of [`IntDelaunay::relax`]. + pub fn relax_mut(&mut self, options: RelaxationOptions) -> RelaxationResult { + if options.max_iterations == 0 { + return RelaxationResult { + iterations: 0, + converged: false, + }; + } + + debug_assert_positive_areas(&self.triangles); + + let mut buffer = RelaxationBuffer::new(self); + + for iteration in 0..options.max_iterations { + buffer.rebuild_starts(&self.triangles); + buffer.proposals.clone_from(&self.points); + + let mut within_tolerance = true; + + for vertex_index in 0..self.points.len() { + if buffer.fixed[vertex_index] { + continue; + } + + let start = buffer.starts[vertex_index]; + if start >= self.triangles.len() + || !self.collect_centroid_cell( + vertex_index, + start, + &mut buffer.cell, + &mut buffer.ring, + ) + { + continue; + } + + let point = self.points[vertex_index]; + let Some(centroid) = polygon_centroid(&buffer.cell, point) else { + continue; + }; + + let displacement = centroid - point; + let sqr_distance = displacement.sqr_length(); + if !is_within_tolerance::(sqr_distance, options.tolerance) { + within_tolerance = false; + } + + buffer.proposals[vertex_index] = + self.limit_displacement(point, displacement, &buffer.ring); + } + + if within_tolerance { + return RelaxationResult { + iterations: iteration, + converged: true, + }; + } + + resolve_collisions( + &self.points, + &mut buffer.proposals, + &buffer.fixed, + &mut buffer.order, + ); + + if buffer.proposals == self.points { + return RelaxationResult { + iterations: iteration, + converged: true, + }; + } + + self.points.clone_from(&buffer.proposals); + self.sync_triangle_points(); + + debug_assert_positive_areas(&self.triangles); + + self.triangles.build_with_buffer(&mut buffer.delaunay); + + debug_assert_positive_areas(&self.triangles); + } + + RelaxationResult { + iterations: options.max_iterations, + converged: false, + } + } + + fn collect_centroid_cell( + &self, + vertex_index: usize, + start: usize, + cell: &mut Vec>, + ring: &mut Vec, + ) -> bool { + cell.clear(); + ring.clear(); + + let mut triangle_index = start; + loop { + if triangle_index >= self.triangles.len() || ring.len() >= self.triangles.len() { + cell.clear(); + ring.clear(); + return false; + } + + let triangle = &self.triangles[triangle_index]; + debug_assert!(triangle + .vertices + .iter() + .any(|vertex| vertex.index == vertex_index)); + + let (next, mid) = triangle.left_neighbor_and_mid_edge(vertex_index); + ring.push(triangle_index); + cell.push(triangle.center()); + cell.push(mid); + + if next == start { + return true; + } + if next >= self.triangles.len() { + cell.clear(); + ring.clear(); + return false; + } + triangle_index = next; + } + } + + fn limit_displacement( + &self, + point: IntPoint, + displacement: IntVector, + ring: &[usize], + ) -> IntPoint { + let mut dx = displacement.x; + let mut dy = displacement.y; + + while dx != I::Wide::ZERO || dy != I::Wide::ZERO { + let candidate_displacement = IntVector::::new(dx, dy); + if ring.iter().all(|&triangle_index| { + displacement_is_safe(candidate_displacement, &self.triangles[triangle_index]) + }) { + break; + } + + dx = dx / I::Wide::TWO; + dy = dy / I::Wide::TWO; + } + + IntPoint::new( + I::from_wide(point.x.to_wide() + dx), + I::from_wide(point.y.to_wide() + dy), + ) + } + + fn sync_triangle_points(&mut self) { + for triangle in &mut self.triangles { + for vertex in &mut triangle.vertices { + debug_assert!(vertex.index < self.points.len()); + vertex.point = self.points[vertex.index]; + } + } + } +} + +fn polygon_centroid( + polygon: &[IntPoint], + origin: IntPoint, +) -> Option> { + if polygon.len() < 3 { + return None; + } + + let mut area_two = I::Wide::ZERO; + let mut x_numerator = SignedProduct::::multiply(I::Wide::ZERO, I::Wide::ZERO); + let mut y_numerator = SignedProduct::::multiply(I::Wide::ZERO, I::Wide::ZERO); + + for index in 0..polygon.len() { + let a = polygon[index] - origin; + let b = polygon[(index + 1) % polygon.len()] - origin; + let cross = a.cross_product(b); + + area_two = area_two + cross; + x_numerator = x_numerator.checked_add(SignedProduct::multiply(a.x + b.x, cross))?; + y_numerator = y_numerator.checked_add(SignedProduct::multiply(a.y + b.y, cross))?; + } + + if area_two <= I::Wide::ZERO { + return None; + } + + let three = I::WideUInt::from_u64(3); + let area = area_two.to_uint(); + if area > (I::WideUInt::LAST_BIT - I::WideUInt::ONE) / three { + return None; + } + let divisor = area * three; + + let dx = divide_signed_product(x_numerator, divisor)?; + let dy = divide_signed_product(y_numerator, divisor)?; + + Some(IntPoint::new( + I::from_wide(origin.x.to_wide() + dx), + I::from_wide(origin.y.to_wide() + dy), + )) +} + +fn divide_signed_product(value: SignedProduct, divisor: I::UInt) -> Option { + debug_assert!(divisor > I::UInt::ZERO); + debug_assert!(divisor < I::UInt::LAST_BIT); + + let magnitude = value.magnitude().divide_with_rounding(divisor); + if magnitude >= I::UInt::LAST_BIT { + return None; + } + + let result = I::from_uint(magnitude); + Some(if value.is_negative() { -result } else { result }) +} + +fn displacement_is_safe( + displacement: IntVector, + triangle: &IntTriangle, +) -> bool { + if displacement.x == I::Wide::ZERO && displacement.y == I::Wide::ZERO { + return true; + } + + let sqr_displacement = displacement.sqr_length(); + if sqr_displacement <= I::Wide::ZERO { + return false; + } + + let a = triangle.vertices[0].point; + let b = triangle.vertices[1].point; + let c = triangle.vertices[2].point; + + let area_two = Triangle::area_two(a, b, c); + debug_assert!(area_two > I::Wide::ZERO); + if area_two <= I::Wide::ZERO { + return false; + } + + let max_sqr_edge = a + .sqr_distance(b) + .max(b.sqr_distance(c)) + .max(c.sqr_distance(a)); + if max_sqr_edge <= I::Wide::ZERO { + return false; + } + + let left = ::Product::multiply( + sqr_displacement.to_uint(), + max_sqr_edge.to_uint(), + ); + + // |d| < h_min / 4, where h_min = area_two / longest_edge. + // Squaring and rearranging avoids both division and square roots: + // 16 * |d|^2 * longest_edge^2 < area_two^2. + let Some(left2) = left.checked_add(left) else { + return false; + }; + let Some(left4) = left2.checked_add(left2) else { + return false; + }; + let Some(left8) = left4.checked_add(left4) else { + return false; + }; + let Some(left16) = left8.checked_add(left8) else { + return false; + }; + + let area = area_two.to_uint(); + let right = ::Product::multiply(area, area); + + left16 < right +} + +fn is_within_tolerance(sqr_distance: I::Wide, tolerance: I::WideUInt) -> bool { + if sqr_distance < I::Wide::ZERO { + return false; + } + + let distance = ::Product::from_uint(sqr_distance.to_uint()); + let tolerance = ::Product::multiply(tolerance, tolerance); + distance <= tolerance +} + +fn resolve_collisions( + points: &[IntPoint], + proposals: &mut [IntPoint], + fixed: &[bool], + order: &mut [usize], +) { + loop { + // Topology normalization may represent regions touching at a former + // self-intersection with distinct vertex indices at the same point. + // Preserve those existing coincidences and reject only newly merged + // coordinates. + order.sort_unstable_by_key(|&index| (proposals[index], points[index])); + if !order.windows(2).any(|pair| { + proposals[pair[0]] == proposals[pair[1]] && points[pair[0]] != points[pair[1]] + }) { + return; + } + + // Back off every movable proposal together so the update remains + // synchronous and deterministic. + let mut changed = false; + for index in 0..proposals.len() { + if fixed[index] { + continue; + } + + let point = points[index]; + let displacement = proposals[index] - point; + let next = IntPoint::new( + I::from_wide(point.x.to_wide() + displacement.x / I::Wide::TWO), + I::from_wide(point.y.to_wide() + displacement.y / I::Wide::TWO), + ); + changed |= next != proposals[index]; + proposals[index] = next; + } + + if !changed { + proposals.clone_from_slice(points); + return; + } + } +} + +#[inline] +fn debug_assert_positive_areas(_triangles: &[IntTriangle]) { + #[cfg(debug_assertions)] + for triangle in _triangles { + debug_assert!( + Triangle::area_two( + triangle.vertices[0].point, + triangle.vertices[1].point, + triangle.vertices[2].point, + ) > I::Wide::ZERO + ); + } +} + +#[cfg(test)] +mod tests { + use super::{resolve_collisions, RelaxationOptions, RelaxationResult}; + use crate::advanced::delaunay::DelaunayCondition; + use crate::int::triangulatable::IntTriangulatable; + use crate::int::uniform::IntUniformTriangulatable; + use alloc::vec; + use alloc::vec::Vec; + use i_overlay::i_float::int::point::IntPoint; + + #[test] + fn zero_iterations_does_not_change_mesh() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + let mut delaunay = contour.uniform_triangulate(20u64); + let points = delaunay.points.clone(); + + let result = delaunay.relax_mut(RelaxationOptions::new(0)); + + assert_eq!( + result, + RelaxationResult { + iterations: 0, + converged: false, + } + ); + assert_eq!(delaunay.points, points); + } + + #[test] + fn large_tolerance_stops_before_moving() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + let mut delaunay = contour + .triangulate_with_steiner_points(&[ + IntPoint::new(25, 35), + IntPoint::new(72, 42), + IntPoint::new(43, 79), + ]) + .into_delaunay(); + let points = delaunay.points.clone(); + + let result = delaunay.relax_mut(RelaxationOptions::new(10).with_tolerance(1_000u64)); + + assert_eq!(result.iterations, 0); + assert!(result.converged); + assert_eq!(delaunay.points, points); + } + + #[test] + fn moves_interior_vertices_and_keeps_boundary_fixed() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + let mut delaunay = contour + .triangulate_with_steiner_points(&[ + IntPoint::new(25, 35), + IntPoint::new(72, 42), + IntPoint::new(43, 79), + ]) + .into_delaunay(); + let points = delaunay.points.clone(); + let boundary: Vec<_> = points + .iter() + .enumerate() + .filter(|(_, point)| { + (point.x == 0 || point.x == 100) && (point.y == 0 || point.y == 100) + }) + .map(|(index, _)| index) + .collect(); + + let result = delaunay.relax_mut(RelaxationOptions::new(2)); + + assert_eq!(result.iterations, 2); + assert!(!result.converged); + for &index in &boundary { + assert_eq!(delaunay.points[index], points[index]); + } + assert!(delaunay + .points + .iter() + .zip(&points) + .enumerate() + .any(|(index, (a, b))| !boundary.contains(&index) && a != b)); + assert_is_delaunay(&delaunay); + } + + #[test] + fn keeps_outer_and_hole_boundaries_fixed() { + let shape = vec![ + vec![ + IntPoint::new(0, 0), + IntPoint::new(120, 0), + IntPoint::new(120, 120), + IntPoint::new(0, 120), + ], + vec![ + IntPoint::new(45, 45), + IntPoint::new(45, 75), + IntPoint::new(75, 75), + IntPoint::new(75, 45), + ], + ]; + let mut delaunay = shape.uniform_triangulate(15u64); + let points = delaunay.points.clone(); + let boundary: Vec<_> = points + .iter() + .enumerate() + .filter(|(_, point)| { + point.x == 0 + || point.x == 120 + || point.y == 0 + || point.y == 120 + || ((point.x == 45 || point.x == 75) && (45..=75).contains(&point.y)) + || ((point.y == 45 || point.y == 75) && (45..=75).contains(&point.x)) + }) + .map(|(index, _)| index) + .collect(); + + delaunay.relax_mut(RelaxationOptions::new(3)); + + for index in boundary { + assert_eq!(delaunay.points[index], points[index]); + } + assert_is_delaunay(&delaunay); + } + + #[test] + fn consuming_relaxation_with_default_options_preserves_mesh() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + let delaunay = contour + .triangulate_with_steiner_points(&[ + IntPoint::new(25, 35), + IntPoint::new(72, 42), + IntPoint::new(43, 79), + ]) + .into_delaunay(); + let point_count = delaunay.points.len(); + let triangle_count = delaunay.triangles.len(); + + let relaxed = delaunay.relax(RelaxationOptions::default()); + + assert_eq!(relaxed.points.len(), point_count); + assert_eq!(relaxed.triangles.len(), triangle_count); + assert_is_delaunay(&relaxed); + } + + #[test] + fn collision_resolution_backs_off_movable_vertices() { + let points = [ + IntPoint::new(0i32, 0), + IntPoint::new(4, 0), + IntPoint::new(8, 0), + ]; + let mut proposals = [IntPoint::new(0, 0); 3]; + let fixed = [true, false, false]; + let mut order = [0, 1, 2]; + + resolve_collisions(&points, &mut proposals, &fixed, &mut order); + + assert_eq!( + proposals, + [ + IntPoint::new(0, 0), + IntPoint::new(2, 0), + IntPoint::new(4, 0), + ] + ); + } + + fn assert_is_delaunay( + delaunay: &crate::advanced::delaunay::IntDelaunay, + ) { + for (triangle_index, triangle) in delaunay.triangles.iter().enumerate() { + for &neighbor in &triangle.neighbors { + if neighbor <= triangle_index || neighbor >= delaunay.triangles.len() { + continue; + } + + let abc = triangle.abc_by_neighbor(neighbor); + let pcb = delaunay.triangles[neighbor].abc_by_neighbor(triangle_index); + assert!(DelaunayCondition::is_flip_not_required( + pcb.v0.vertex.point, + abc.v0.vertex.point, + abc.v1.vertex.point, + abc.v2.vertex.point, + )); + } + } + } +} diff --git a/iTriangle/src/float/mod.rs b/iTriangle/src/float/mod.rs index b601b5a..a7e8726 100644 --- a/iTriangle/src/float/mod.rs +++ b/iTriangle/src/float/mod.rs @@ -5,7 +5,9 @@ pub mod convex; pub mod custom; pub mod delaunay; pub mod locator; +pub mod relax; pub mod triangulatable; pub mod triangulation; pub mod triangulator; pub mod unchecked; +pub mod uniform; diff --git a/iTriangle/src/float/relax.rs b/iTriangle/src/float/relax.rs new file mode 100644 index 0000000..f910069 --- /dev/null +++ b/iTriangle/src/float/relax.rs @@ -0,0 +1,193 @@ +use crate::advanced::relax::RelaxationOptions as IntRelaxationOptions; +use crate::float::delaunay::Delaunay; +use i_overlay::i_float::float::compatible::FloatPointCompatible; +use i_overlay::i_float::float::number::FloatNumber; +use i_overlay::i_float::int::number::int::IntNumber; + +pub use crate::advanced::relax::RelaxationResult; + +/// Configuration for centroid-net relaxation in float coordinate units. +/// +/// This is a thin wrapper over the integer relaxation implementation. The +/// convergence tolerance is converted to the integer coordinate system by the +/// same adapter that was used to construct the mesh. +#[derive(Clone, Copy)] +pub struct RelaxationOptions { + /// Maximum number of relaxation iterations. + pub max_iterations: usize, + + /// Absolute convergence tolerance in input float coordinate units. + /// + /// Relaxation stops before the next move when every unconstrained + /// centroid displacement is no greater than this value. It must be finite + /// and non-negative. + pub tolerance: F, +} + +impl RelaxationOptions { + /// Creates options with zero convergence tolerance. + #[inline] + pub fn new(max_iterations: usize) -> Self { + Self { + max_iterations, + tolerance: F::ZERO, + } + } + + /// Sets the absolute convergence tolerance in float coordinate units. + #[inline] + pub fn with_tolerance(mut self, tolerance: F) -> Self { + self.tolerance = tolerance; + self + } +} + +impl Default for RelaxationOptions { + #[inline] + fn default() -> Self { + Self::new(8) + } +} + +impl Delaunay { + /// Relaxes interior vertices toward their centroid-net area centroids. + /// + /// Boundary vertices, including vertices on hole boundaries, remain fixed. + /// The operation is performed entirely by the integer mesh and then mapped + /// back to `P` when points are requested. + #[must_use] + #[inline] + pub fn relax(mut self, options: RelaxationOptions) -> Self { + self.relax_mut(options); + self + } + + /// In-place version of [`Delaunay::relax`]. + #[inline] + pub fn relax_mut(&mut self, options: RelaxationOptions) -> RelaxationResult { + assert!( + options.tolerance.is_finite() && options.tolerance >= P::Scalar::ZERO, + "tolerance must be finite and non-negative" + ); + + let tolerance = self.adapter.round_len_to_int(options.tolerance).to_uint(); + self.delaunay.relax_mut(IntRelaxationOptions { + max_iterations: options.max_iterations, + tolerance, + }) + } +} + +#[cfg(test)] +mod tests { + use super::RelaxationOptions; + use crate::float::triangulatable::Triangulatable; + use crate::float::uniform::UniformTriangulatable; + + #[test] + fn large_float_tolerance_stops_before_moving() { + let contour = [[0.0, 0.0], [100.0, 0.0], [100.0, 100.0], [0.0, 100.0]]; + let steiner = [[25.0, 35.0], [72.0, 42.0], [43.0, 79.0]]; + let mut delaunay = contour + .triangulate_with_steiner_points(&steiner) + .into_delaunay(); + let points = delaunay.points(); + + let result = delaunay.relax_mut(RelaxationOptions::new(10).with_tolerance(1_000.0)); + + assert_eq!(result.iterations, 0); + assert!(result.converged); + assert_eq!(delaunay.points(), points); + } + + #[test] + fn float_wrapper_moves_interior_and_keeps_boundary_fixed() { + let contour = [[0.0, 0.0], [100.0, 0.0], [100.0, 100.0], [0.0, 100.0]]; + let steiner = [[25.0, 35.0], [72.0, 42.0], [43.0, 79.0]]; + let mut delaunay = contour + .triangulate_with_steiner_points(&steiner) + .into_delaunay(); + let points = delaunay.points(); + let boundary: alloc::vec::Vec<_> = points + .iter() + .enumerate() + .filter(|(_, point)| { + (point[0] == 0.0 || point[0] == 100.0) && (point[1] == 0.0 || point[1] == 100.0) + }) + .map(|(index, _)| index) + .collect(); + + let result = delaunay.relax_mut(RelaxationOptions::new(2)); + let relaxed = delaunay.points(); + + assert_eq!(result.iterations, 2); + assert!(!result.converged); + for &index in &boundary { + assert_eq!(relaxed[index], points[index]); + } + assert!(relaxed + .iter() + .zip(&points) + .enumerate() + .any(|(index, (a, b))| !boundary.contains(&index) && a != b)); + delaunay + .to_triangulation::() + .validate(10_000.0, 0.000_001); + } + + #[test] + fn consuming_relaxation_with_default_options_returns_valid_mesh() { + let contour = [[0.0, 0.0], [100.0, 0.0], [100.0, 100.0], [0.0, 100.0]]; + let steiner = [[25.0, 35.0], [72.0, 42.0], [43.0, 79.0]]; + let delaunay = contour + .triangulate_with_steiner_points(&steiner) + .into_delaunay(); + + let relaxed = delaunay.relax(RelaxationOptions::default()); + + relaxed + .to_triangulation::() + .validate(10_000.0, 0.000_001); + } + + #[test] + #[should_panic(expected = "tolerance must be finite and non-negative")] + fn rejects_negative_tolerance() { + let contour = [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]]; + let mut delaunay = contour.triangulate().into_delaunay(); + + delaunay.relax_mut(RelaxationOptions::new(1).with_tolerance(-1.0)); + } + + #[test] + fn relaxes_uniform_mesh_from_self_intersecting_contour() { + struct RelaxationStats { + iterations: usize, + converged: bool, + } + + let shape = [[0.0, 0.0], [100.0, 100.0], [0.0, 100.0], [100.0, 0.0]]; + let edge_length = 10.0; + let relax_enabled = true; + let relax_iterations = 4; + + let mut delaunay = shape.uniform_triangulate(edge_length); + let relaxation = relax_enabled.then(|| { + let result = delaunay.relax_mut(RelaxationOptions::new(relax_iterations)); + RelaxationStats { + iterations: result.iterations, + converged: result.converged, + } + }); + + let relaxation = relaxation.unwrap(); + assert!(relaxation.iterations <= relax_iterations); + assert_eq!( + relaxation.converged, + relaxation.iterations < relax_iterations + ); + delaunay + .to_triangulation::() + .validate(5_000.0, 0.000_001); + } +} diff --git a/iTriangle/src/float/uniform.rs b/iTriangle/src/float/uniform.rs new file mode 100644 index 0000000..981b9a2 --- /dev/null +++ b/iTriangle/src/float/uniform.rs @@ -0,0 +1,118 @@ +use crate::float::delaunay::Delaunay; +use crate::int::uniform::IntUniformTriangulatable; +use i_overlay::core::integer::OverlayInt; +use i_overlay::i_float::adapter::FloatPointAdapter; +use i_overlay::i_float::float::compatible::FloatPointCompatible; +use i_overlay::i_float::float::number::FloatNumber; +use i_overlay::i_float::float::rect::FloatRect; +use i_overlay::i_shape::float::adapter::PathToInt; +use i_overlay::i_shape::int::shape::IntShape; +use i_overlay::i_shape::source::resource::ShapeResource; + +/// Float wrapper for the integer uniform triangulation pipeline. +/// +/// The input is converted once with a shared [`FloatPointAdapter`]. Boundary +/// splitting, topology normalization, lattice generation, boundary clearance, +/// and Delaunay triangulation are all performed by the integer implementation. +/// +/// # Example +/// +/// ``` +/// use i_triangle::float::uniform::UniformTriangulatable; +/// +/// let contour = [ +/// [0.0, 0.0], +/// [10.0, 0.0], +/// [10.0, 10.0], +/// [0.0, 10.0], +/// ]; +/// +/// let mesh = contour +/// .uniform_triangulate(2.0) +/// .to_triangulation::(); +/// +/// assert!(!mesh.indices.is_empty()); +/// ``` +pub trait UniformTriangulatable { + /// Triangulates with the default `i32` integer engine. + fn uniform_triangulate(&self, edge_length: P::Scalar) -> Delaunay

{ + self.uniform_triangulate_as::(edge_length) + } + + /// Triangulates with the requested integer engine. + fn uniform_triangulate_as(&self, edge_length: P::Scalar) -> Delaunay + where + I: OverlayInt; +} + +impl UniformTriangulatable

for S +where + S: ShapeResource

, + P: FloatPointCompatible, +{ + fn uniform_triangulate_as(&self, edge_length: P::Scalar) -> Delaunay + where + I: OverlayInt, + { + assert!( + edge_length.is_finite() && edge_length > P::Scalar::ZERO, + "edge_length must be finite and positive" + ); + + let rect = + FloatRect::with_iter(self.iter_paths().flatten()).unwrap_or_else(FloatRect::zero); + let adapter = FloatPointAdapter::::new(rect); + let int_edge_length = adapter.round_len_to_int(edge_length); + assert!( + int_edge_length > I::ONE, + "edge_length is below the precision of the selected integer engine" + ); + + let shape: IntShape = self + .iter_paths() + .map(|path| path.to_int(&adapter)) + .collect(); + let delaunay = shape.uniform_triangulate(int_edge_length.to_uint()); + + Delaunay { delaunay, adapter } + } +} + +#[cfg(test)] +mod tests { + use super::UniformTriangulatable; + use alloc::vec; + + #[test] + fn fills_square_with_boundary_and_grid_points() { + let contour = [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]]; + + let triangulation = contour.uniform_triangulate(2.0).to_triangulation::(); + + assert!(triangulation.points.len() > 20); + triangulation.validate(100.0, 0.000_001); + } + + #[test] + fn fills_shape_without_filling_hole() { + let shape = vec![ + vec![[0.0, 0.0], [20.0, 0.0], [20.0, 20.0], [0.0, 20.0]], + vec![[7.0, 7.0], [7.0, 13.0], [13.0, 13.0], [13.0, 7.0]], + ]; + + let triangulation = shape.uniform_triangulate(2.0).to_triangulation::(); + + assert!(triangulation.points.len() > 40); + triangulation.validate(364.0, 0.000_001); + } + + #[test] + fn narrow_shape_falls_back_to_split_boundary() { + let contour = [[0.0, 0.0], [10.0, 0.0], [10.0, 0.5], [0.0, 0.5]]; + + let triangulation = contour.uniform_triangulate(2.0).to_triangulation::(); + + assert!(!triangulation.indices.is_empty()); + triangulation.validate(5.0, 0.000_001); + } +} diff --git a/iTriangle/src/int/constraint.rs b/iTriangle/src/int/constraint.rs new file mode 100644 index 0000000..fa8b084 --- /dev/null +++ b/iTriangle/src/int/constraint.rs @@ -0,0 +1,283 @@ +use crate::advanced::delaunay::{DelaunayCondition, IntDelaunay}; +use crate::geom::triangle::IntTriangle; +use crate::int::triangulation::RawIntTriangulation; +use alloc::vec::Vec; +use i_overlay::i_float::int::number::int::IntNumber; +use i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_overlay::i_float::int::point::IntPoint; +use i_overlay::i_float::triangle::Triangle; +use i_overlay::string::line::IntLine; + +#[derive(Clone, Copy, PartialEq, Eq)] +struct ConstraintEdge { + a: usize, + b: usize, +} + +impl ConstraintEdge { + #[inline] + fn new(a: usize, b: usize) -> Self { + if a < b { + Self { a, b } + } else { + Self { a: b, b: a } + } + } +} + +pub(super) trait Constrain { + fn into_constrained_delaunay(self, constraints: &[IntLine]) -> IntDelaunay; +} + +pub(super) fn constraint_points(constraints: &[IntLine]) -> Vec> { + let mut points = Vec::with_capacity(2 * constraints.len()); + for &line in constraints { + if line[0] != line[1] { + points.extend_from_slice(&line); + } + } + points.sort_unstable(); + points.dedup(); + points +} + +impl Constrain for RawIntTriangulation { + fn into_constrained_delaunay(mut self, constraints: &[IntLine]) -> IntDelaunay { + let mut locked = Vec::new(); + let mut vertices = Vec::new(); + + for &line in constraints { + if line[0] == line[1] { + continue; + } + + self.vertices_on_line(line, &mut vertices); + assert!( + vertices.len() >= 2 + && self.points[vertices[0]] == line[0].min(line[1]) + && self.points[*vertices.last().unwrap()] == line[0].max(line[1]), + "constraint endpoints must lie strictly inside the triangulated geometry" + ); + + for pair in vertices.windows(2) { + let edge = ConstraintEdge::new(pair[0], pair[1]); + self.recover_edge(edge, &locked); + if !locked.contains(&edge) { + locked.push(edge); + } + } + } + + self.refine_delaunay(&locked); + + IntDelaunay { + triangles: self.triangles, + points: self.points, + } + } +} + +impl RawIntTriangulation { + fn vertices_on_line(&self, line: IntLine, result: &mut Vec) { + result.clear(); + let min = line[0].min(line[1]); + let max = line[0].max(line[1]); + + for (index, &point) in self.points.iter().enumerate() { + if min <= point + && point <= max + && Triangle::area_two(line[0], line[1], point) == I::Wide::ZERO + { + result.push(index); + } + } + + result.sort_unstable_by_key(|&index| self.points[index]); + } + + fn recover_edge(&mut self, constraint: ConstraintEdge, locked: &[ConstraintEdge]) { + if self.has_edge(constraint) { + return; + } + + let flip_limit = self.triangles.len().saturating_mul(self.triangles.len()) + 1; + for _ in 0..flip_limit { + let mut candidate = None; + + 'scan: for triangle_index in 0..self.triangles.len() { + let neighbors = self.triangles[triangle_index].neighbors; + for neighbor_index in neighbors { + if neighbor_index <= triangle_index || neighbor_index >= self.triangles.len() { + continue; + } + + let abc = self.triangles[triangle_index].abc_by_neighbor(neighbor_index); + let edge = ConstraintEdge::new(abc.v1.vertex.index, abc.v2.vertex.index); + if locked.contains(&edge) { + continue; + } + + let p0 = self.points[constraint.a]; + let p1 = self.points[constraint.b]; + if proper_intersection(p0, p1, abc.v1.vertex.point, abc.v2.vertex.point) + && can_flip(&self.triangles, triangle_index, neighbor_index) + { + let pcb = self.triangles[neighbor_index].abc_by_neighbor(triangle_index); + if !proper_intersection(p0, p1, abc.v0.vertex.point, pcb.v0.vertex.point) { + candidate = Some((triangle_index, neighbor_index)); + break 'scan; + } + } + } + } + + if let Some((triangle_index, neighbor_index)) = candidate { + flip(&mut self.triangles, triangle_index, neighbor_index); + if self.has_edge(constraint) { + return; + } + } else { + break; + } + } + + panic!("unable to recover constraint edge; constraints must not intersect"); + } + + #[inline] + fn has_edge(&self, edge: ConstraintEdge) -> bool { + self.triangles.iter().any(|triangle| { + let mut has_a = false; + let mut has_b = false; + for vertex in triangle.vertices { + has_a |= vertex.index == edge.a; + has_b |= vertex.index == edge.b; + } + has_a && has_b + }) + } + + fn refine_delaunay(&mut self, locked: &[ConstraintEdge]) { + loop { + let mut candidate = None; + + 'scan: for triangle_index in 0..self.triangles.len() { + let neighbors = self.triangles[triangle_index].neighbors; + for neighbor_index in neighbors { + if neighbor_index <= triangle_index || neighbor_index >= self.triangles.len() { + continue; + } + + let abc = self.triangles[triangle_index].abc_by_neighbor(neighbor_index); + let edge = ConstraintEdge::new(abc.v1.vertex.index, abc.v2.vertex.index); + if locked.contains(&edge) + || !can_flip(&self.triangles, triangle_index, neighbor_index) + { + continue; + } + + let pcb = self.triangles[neighbor_index].abc_by_neighbor(triangle_index); + if !DelaunayCondition::is_flip_not_required( + pcb.v0.vertex.point, + abc.v0.vertex.point, + abc.v1.vertex.point, + abc.v2.vertex.point, + ) { + candidate = Some((triangle_index, neighbor_index)); + break 'scan; + } + } + } + + if let Some((triangle_index, neighbor_index)) = candidate { + flip(&mut self.triangles, triangle_index, neighbor_index); + } else { + return; + } + } + } +} + +#[inline] +fn proper_intersection( + a: IntPoint, + b: IntPoint, + c: IntPoint, + d: IntPoint, +) -> bool { + opposite_signs(Triangle::area_two(a, b, c), Triangle::area_two(a, b, d)) + && opposite_signs(Triangle::area_two(c, d, a), Triangle::area_two(c, d, b)) +} + +#[inline] +fn opposite_signs(a: W, b: W) -> bool { + a < W::ZERO && b > W::ZERO || a > W::ZERO && b < W::ZERO +} + +#[inline] +fn can_flip( + triangles: &[IntTriangle], + triangle_index: usize, + neighbor_index: usize, +) -> bool { + let abc = triangles[triangle_index].abc_by_neighbor(neighbor_index); + let pcb = triangles[neighbor_index].abc_by_neighbor(triangle_index); + opposite_signs( + Triangle::area_two( + abc.v0.vertex.point, + pcb.v0.vertex.point, + abc.v1.vertex.point, + ), + Triangle::area_two( + abc.v0.vertex.point, + pcb.v0.vertex.point, + abc.v2.vertex.point, + ), + ) +} + +fn flip( + triangles: &mut [IntTriangle], + triangle_index: usize, + neighbor_index: usize, +) { + let abc = triangles[triangle_index].abc_by_neighbor(neighbor_index); + let pcb = triangles[neighbor_index].abc_by_neighbor(triangle_index); + + update_neighbor(triangles, abc.v1.neighbor, triangle_index, neighbor_index); + update_neighbor(triangles, pcb.v1.neighbor, neighbor_index, triangle_index); + + let abp = &mut triangles[triangle_index]; + abp.neighbors[abc.v0.position] = pcb.v1.neighbor; + abp.neighbors[abc.v1.position] = neighbor_index; + abp.neighbors[abc.v2.position] = abc.v2.neighbor; + abp.vertices[abc.v2.position] = pcb.v0.vertex; + + let pca = &mut triangles[neighbor_index]; + pca.neighbors[pcb.v0.position] = abc.v1.neighbor; + pca.neighbors[pcb.v1.position] = triangle_index; + pca.neighbors[pcb.v2.position] = pcb.v2.neighbor; + pca.vertices[pcb.v2.position] = abc.v0.vertex; +} + +#[inline] +fn update_neighbor( + triangles: &mut [IntTriangle], + neighbor_index: usize, + old_index: usize, + new_index: usize, +) { + if neighbor_index >= triangles.len() { + return; + } + + let triangle = &mut triangles[neighbor_index]; + if triangle.neighbors[0] == old_index { + triangle.neighbors[0] = new_index; + } else if triangle.neighbors[1] == old_index { + triangle.neighbors[1] = new_index; + } else { + debug_assert_eq!(triangle.neighbors[2], old_index); + triangle.neighbors[2] = new_index; + } +} diff --git a/iTriangle/src/int/earcut/earcut_64.rs b/iTriangle/src/int/earcut/earcut_64.rs index 4f818b7..1d833d0 100644 --- a/iTriangle/src/int/earcut/earcut_64.rs +++ b/iTriangle/src/int/earcut/earcut_64.rs @@ -553,6 +553,7 @@ mod tests { use crate::int::earcut::earcut_64::{Bit, Ear, Earcut64, EarcutSolver}; use crate::int::earcut::flat::FlatEarcutStore; use crate::int::triangulation::{IntTriangulation, RawIntTriangulation}; + use crate::test_util::random_cases; use alloc::vec; use alloc::vec::Vec; use i_overlay::core::fill_rule::FillRule; @@ -1558,7 +1559,7 @@ mod tests { #[test] fn test_random_0() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { if let Some(first) = random(8, 5) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1574,7 +1575,7 @@ mod tests { #[test] fn test_random_1() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { if let Some(first) = random(8, 7) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1590,7 +1591,7 @@ mod tests { #[test] fn test_random_2() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { if let Some(first) = random(8, 10) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1606,7 +1607,7 @@ mod tests { #[test] fn test_random_3() { - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { if let Some(first) = random(8, 12) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1622,7 +1623,7 @@ mod tests { #[test] fn test_random_4() { - for _ in 0..10_000 { + for _ in 0..random_cases(10_000) { if let Some(first) = random(16, 32) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1639,7 +1640,7 @@ mod tests { #[test] fn test_random_5() { - for _ in 0..5_000 { + for _ in 0..random_cases(5_000) { if let Some(first) = random(16, 48) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() @@ -1656,7 +1657,7 @@ mod tests { #[test] fn test_random_6() { - for _ in 0..2_000 { + for _ in 0..random_cases(2_000) { if let Some(first) = random(16, 64) .simplify(FillRule::NonZero, IntOverlayOptions::keep_output_points()) .first() diff --git a/iTriangle/src/int/mod.rs b/iTriangle/src/int/mod.rs index fd64741..438445f 100644 --- a/iTriangle/src/int/mod.rs +++ b/iTriangle/src/int/mod.rs @@ -1,4 +1,5 @@ mod binder; +mod constraint; pub mod custom; pub mod earcut; pub mod locator; @@ -9,4 +10,5 @@ pub mod triangulatable; pub mod triangulation; pub mod triangulator; pub mod unchecked; +pub mod uniform; pub mod validation; diff --git a/iTriangle/src/int/monotone/flat/triangulator.rs b/iTriangle/src/int/monotone/flat/triangulator.rs index 894c804..3476a08 100644 --- a/iTriangle/src/int/monotone/flat/triangulator.rs +++ b/iTriangle/src/int/monotone/flat/triangulator.rs @@ -387,6 +387,7 @@ mod tests { use crate::int::monotone::triangulator::MonotoneTriangulator; use crate::int::triangulation::IntTriangulation; + use crate::test_util::random_cases; use alloc::vec; use alloc::vec::Vec; use i_overlay::core::fill_rule::FillRule; @@ -820,7 +821,7 @@ mod tests { #[test] fn test_random_0() { let mut raw = IntTriangulation::::default(); - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { let path = random(8, 5); let shape = vec![path]; if let Some(first) = shape @@ -839,7 +840,7 @@ mod tests { #[test] fn test_random_1() { let mut raw = IntTriangulation::::default(); - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { let path = random(10, 6); let shape = vec![path]; if let Some(first) = shape @@ -858,7 +859,7 @@ mod tests { #[test] fn test_random_2() { let mut raw = IntTriangulation::::default(); - for _ in 0..100_000 { + for _ in 0..random_cases(100_000) { let path = random(10, 12); let shape = vec![path]; if let Some(first) = shape @@ -877,7 +878,7 @@ mod tests { #[test] fn test_random_3() { let mut raw = IntTriangulation::::default(); - for _ in 0..50_000 { + for _ in 0..random_cases(50_000) { let path = random(20, 20); let shape = vec![path]; if let Some(first) = shape @@ -896,7 +897,7 @@ mod tests { #[test] fn test_random_4() { let mut raw = IntTriangulation::::default(); - for _ in 0..5_000 { + for _ in 0..random_cases(5_000) { let path = random(30, 50); let shape = vec![path]; if let Some(first) = shape @@ -915,7 +916,7 @@ mod tests { #[test] fn test_random_5() { let mut raw = IntTriangulation::::default(); - for _ in 0..2_000 { + for _ in 0..random_cases(2_000) { let main = random(50, 20); let mut shape = vec![main]; for _ in 0..10 { diff --git a/iTriangle/src/int/monotone/net/triangulator.rs b/iTriangle/src/int/monotone/net/triangulator.rs index b783dc4..a63fa43 100644 --- a/iTriangle/src/int/monotone/net/triangulator.rs +++ b/iTriangle/src/int/monotone/net/triangulator.rs @@ -629,6 +629,7 @@ mod tests { use crate::int::binder::SteinerInference; use crate::int::monotone::triangulator::MonotoneTriangulator; use crate::int::triangulation::RawIntTriangulation; + use crate::test_util::random_cases; use alloc::vec; use alloc::vec::Vec; use i_overlay::core::fill_rule::FillRule; @@ -1291,7 +1292,7 @@ mod tests { #[test] fn test_random_0() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let path = random(8, 5); let shape = vec![path]; if let Some(first) = shape @@ -1312,7 +1313,7 @@ mod tests { #[test] fn test_random_1() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let path = random(10, 6); let shape = vec![path]; if let Some(first) = shape @@ -1333,7 +1334,7 @@ mod tests { #[test] fn test_random_2() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let path = random(10, 12); let shape = vec![path]; if let Some(first) = shape @@ -1354,7 +1355,7 @@ mod tests { #[test] fn test_random_3() { - for _ in 0..10_000 { + for _ in 0..random_cases(10_000) { let path = random(20, 20); let shape = vec![path]; if let Some(first) = shape @@ -1375,7 +1376,7 @@ mod tests { #[test] fn test_random_4() { - for _ in 0..2_000 { + for _ in 0..random_cases(2_000) { let path = random(30, 50); let shape = vec![path]; if let Some(first) = shape @@ -1396,7 +1397,7 @@ mod tests { #[test] fn test_random_5() { - for _ in 0..1_000 { + for _ in 0..random_cases(1_000) { let main = random(50, 20); let mut shape = vec![main]; for _ in 0..10 { @@ -1423,7 +1424,7 @@ mod tests { fn test_random_6() { let shape = vec![path(&[[-10, 0], [0, -10], [10, 0], [0, 10]])]; let shape_area = shape.area_two(); - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let points = random_points(5, 10); let mut raw = RawIntTriangulation::default(); @@ -1442,7 +1443,7 @@ mod tests { fn test_random_7() { let shapes = vec![vec![path(&[[-5, 0], [0, -5], [5, 0], [0, 5]])]]; let shape_area = shapes.area_two(); - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let points = random_points(8, 2); let group = shapes.group_by_shapes(&points); @@ -1460,7 +1461,7 @@ mod tests { #[test] fn test_random_8() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let points = random_points(15, 1); let shape = random(10, 4); @@ -1488,7 +1489,7 @@ mod tests { #[test] fn test_random_9() { - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let points = random_points(10, 2); let shape = random(10, 4); @@ -1516,7 +1517,7 @@ mod tests { #[test] fn test_random_10() { - for _ in 0..5_000 { + for _ in 0..random_cases(5_000) { let points = random_points(10, 8); let shape = random(10, 8); @@ -1544,7 +1545,7 @@ mod tests { #[test] fn test_random_11() { - for _ in 0..2_000 { + for _ in 0..random_cases(2_000) { let main = random(50, 20); let mut shape = vec![main]; for _ in 0..10 { diff --git a/iTriangle/src/int/triangulatable.rs b/iTriangle/src/int/triangulatable.rs index c29c602..fb67f3e 100644 --- a/iTriangle/src/int/triangulatable.rs +++ b/iTriangle/src/int/triangulatable.rs @@ -1,9 +1,12 @@ +use crate::advanced::delaunay::IntDelaunay; +use crate::int::constraint::{constraint_points, Constrain}; use crate::int::solver::{ContourSolver, ShapeSolver, ShapesSolver}; use crate::int::triangulation::RawIntTriangulation; use i_overlay::core::integer::OverlayInt; use i_overlay::i_float::int::number::int::IntNumber; use i_overlay::i_float::int::point::IntPoint; use i_overlay::i_shape::int::shape::{IntContour, IntShape, IntShapes}; +use i_overlay::string::line::IntLine; /// A trait for performing triangulation with default validation settings. /// /// Provides a simplified interface for converting shapes or contours into triangle meshes. @@ -69,3 +72,117 @@ impl IntTriangulatable for IntShapes { ShapesSolver::triangulate_with_steiner_points(Default::default(), self, points) } } + +/// A trait for building a constrained Delaunay triangulation. +/// +/// Unlike holes, constraints do not remove any area. They force the requested internal +/// line segments to appear as edges in the resulting triangle mesh. +pub trait IntConstrainedTriangulatable { + /// Builds a constrained Delaunay triangulation containing every requested internal edge. + /// + /// Constraint endpoints must lie strictly inside the geometry. Constraints may share + /// endpoints, but must not cross each other or the shape boundary. + fn triangulate_with_constraints(&self, constraints: &[IntLine]) -> IntDelaunay; +} + +impl IntConstrainedTriangulatable for IntContour { + #[inline] + fn triangulate_with_constraints(&self, constraints: &[IntLine]) -> IntDelaunay { + let points = constraint_points(constraints); + ContourSolver::triangulate_with_steiner_points(Default::default(), self, &points) + .into_constrained_delaunay(constraints) + } +} + +impl IntConstrainedTriangulatable for IntShape { + #[inline] + fn triangulate_with_constraints(&self, constraints: &[IntLine]) -> IntDelaunay { + let points = constraint_points(constraints); + ShapeSolver::triangulate_with_steiner_points(Default::default(), self, &points) + .into_constrained_delaunay(constraints) + } +} + +impl IntConstrainedTriangulatable for IntShapes { + #[inline] + fn triangulate_with_constraints(&self, constraints: &[IntLine]) -> IntDelaunay { + let points = constraint_points(constraints); + ShapesSolver::triangulate_with_steiner_points(Default::default(), self, &points) + .into_constrained_delaunay(constraints) + } +} + +#[cfg(test)] +mod tests { + extern crate std; + + use super::{IntConstrainedTriangulatable, IntTriangulatable}; + use crate::int::triangulation::IntTriangulation; + use i_overlay::i_float::int::point::IntPoint; + use i_overlay::i_shape::int::shape::IntShapes; + use i_overlay::i_shape::int_shapes; + + fn has_edge(triangulation: &IntTriangulation, edge: [IntPoint; 2]) -> bool { + let a = triangulation + .points + .iter() + .position(|&point| point == edge[0]) + .unwrap() as u16; + let b = triangulation + .points + .iter() + .position(|&point| point == edge[1]) + .unwrap() as u16; + triangulation + .indices + .chunks_exact(3) + .any(|triangle| triangle.contains(&a) && triangle.contains(&b)) + } + + fn assert_has_edge(triangulation: &IntTriangulation, edge: [IntPoint; 2]) { + assert!( + has_edge(triangulation, edge), + "missing constraint edge {edge:?}" + ); + } + + #[test] + fn test_0() { + let shapes: IntShapes = int_shapes![[[[-5, -5], [5, -5], [5, 5], [-5, 5]],],]; + let constraints = [[IntPoint::new(-2, 0), IntPoint::new(2, 0)]]; + + let triangulation: IntTriangulation<_, u16> = shapes + .triangulate_with_constraints(&constraints) + .into_triangulation(); + + assert_has_edge(&triangulation, constraints[0]); + + std::println!("points: {:#?}", triangulation.points); + std::println!("triangles: {:#?}", triangulation.indices); + } + + #[test] + fn multiple_constraints_are_preserved() { + let shapes: IntShapes = int_shapes![[[[-10, -10], [10, -10], [10, 10], [-10, 10]],],]; + let constraints = [ + [IntPoint::new(-8, 0), IntPoint::new(8, 0)], + [IntPoint::new(8, 0), IntPoint::new(0, -4)], + [IntPoint::new(0, 1), IntPoint::new(0, 4)], + ]; + + let points: alloc::vec::Vec<_> = constraints.iter().flatten().copied().collect(); + let ordinary: IntTriangulation<_, u16> = shapes + .triangulate_with_steiner_points(&points) + .into_delaunay() + .into_triangulation(); + assert!(!has_edge(&ordinary, constraints[0])); + + let triangulation: IntTriangulation<_, u16> = shapes + .triangulate_with_constraints(&constraints) + .into_triangulation(); + + for edge in constraints { + assert_has_edge(&triangulation, edge); + } + } +} diff --git a/iTriangle/src/int/uniform.rs b/iTriangle/src/int/uniform.rs new file mode 100644 index 0000000..07fb77c --- /dev/null +++ b/iTriangle/src/int/uniform.rs @@ -0,0 +1,167 @@ +use crate::advanced::delaunay::IntDelaunay; +use crate::int::unchecked::IntUncheckedTriangulatable; +use crate::tessellation::split::SliceContour; +use crate::tessellation::uniform::IntUniformGrid; +use i_overlay::core::fill_rule::FillRule; +use i_overlay::core::integer::OverlayInt; +use i_overlay::core::overlay::IntOverlayOptions; +use i_overlay::core::simplify::Simplify; +use i_overlay::i_float::int::number::uint::UIntNumber; +use i_overlay::i_shape::int::shape::{IntContour, IntShape, IntShapes}; + +/// Builds a boundary-conforming Delaunay mesh from a uniform triangular lattice. +/// +/// Boundary edges are split first. The split geometry is then normalized while +/// preserving every inserted collinear point. Interior lattice points that are too +/// close to any boundary edge are discarded before triangulation. +pub trait IntUniformTriangulatable { + /// Triangulates using `edge_length` both as the maximum boundary segment length + /// and as the horizontal spacing of the interior lattice. + fn uniform_triangulate(&self, edge_length: I::WideUInt) -> IntDelaunay; +} + +impl IntUniformTriangulatable for IntContour { + #[inline] + fn uniform_triangulate(&self, edge_length: I::WideUInt) -> IntDelaunay { + validate_edge_length::(edge_length); + let sliced = self.slice_contour(edge_length); + build_uniform( + sliced.simplify(FillRule::NonZero, IntOverlayOptions::keep_all_points()), + edge_length, + ) + } +} + +impl IntUniformTriangulatable for IntShape { + #[inline] + fn uniform_triangulate(&self, edge_length: I::WideUInt) -> IntDelaunay { + validate_edge_length::(edge_length); + let sliced = self.slice_contour(edge_length); + build_uniform( + sliced.simplify(FillRule::NonZero, IntOverlayOptions::keep_all_points()), + edge_length, + ) + } +} + +impl IntUniformTriangulatable for IntShapes { + #[inline] + fn uniform_triangulate(&self, edge_length: I::WideUInt) -> IntDelaunay { + validate_edge_length::(edge_length); + let sliced = self.slice_contour(edge_length); + build_uniform( + sliced.simplify(FillRule::NonZero, IntOverlayOptions::keep_all_points()), + edge_length, + ) + } +} + +#[inline] +fn validate_edge_length(edge_length: I::WideUInt) { + assert!( + edge_length > I::WideUInt::ONE && edge_length <= I::WideUInt::HALF_MASK, + "edge_length must be greater than one and fit the integer coordinate budget" + ); +} + +#[inline] +fn build_uniform(shapes: IntShapes, edge_length: I::WideUInt) -> IntDelaunay { + let steiner_points = shapes.uniform_grid(edge_length); + shapes + .uncheck_triangulate_with_steiner_points(&steiner_points) + .into_delaunay() +} + +#[cfg(test)] +mod tests { + use super::IntUniformTriangulatable; + use alloc::vec; + use i_overlay::i_float::int::point::IntPoint; + use i_overlay::i_shape::int::shape::IntShapes; + + #[test] + fn preserves_split_boundary_points() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + + let delaunay = contour.uniform_triangulate(20u64); + + for x in [20, 40, 60, 80] { + assert!(delaunay.points.contains(&IntPoint::new(x, 0))); + } + assert!(!delaunay.triangles.is_empty()); + } + + #[test] + fn triangulates_shape_with_hole() { + let shape = vec![ + vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ], + vec![ + IntPoint::new(40, 40), + IntPoint::new(40, 60), + IntPoint::new(60, 60), + IntPoint::new(60, 40), + ], + ]; + + let delaunay = shape.uniform_triangulate(10u64); + + assert!(!delaunay.triangles.is_empty()); + assert!(delaunay + .points + .iter() + .all(|p| p.x <= 40 || 60 <= p.x || p.y <= 40 || 60 <= p.y)); + } + + #[test] + fn triangulates_multiple_disjoint_shapes() { + let shapes: IntShapes = vec![ + vec![vec![ + IntPoint::new(0, 0), + IntPoint::new(20, 0), + IntPoint::new(20, 20), + IntPoint::new(0, 20), + ]], + vec![vec![ + IntPoint::new(40, 0), + IntPoint::new(60, 0), + IntPoint::new(60, 20), + IntPoint::new(40, 20), + ]], + ]; + + let delaunay = shapes.uniform_triangulate(5u64); + + assert!(delaunay.points.iter().any(|point| point.x <= 20)); + assert!(delaunay.points.iter().any(|point| point.x >= 40)); + assert!(delaunay + .points + .iter() + .all(|point| point.x <= 20 || point.x >= 40)); + assert!(!delaunay.triangles.is_empty()); + } + + #[test] + #[should_panic( + expected = "edge_length must be greater than one and fit the integer coordinate budget" + )] + fn rejects_edge_length_of_one() { + let contour = vec![ + IntPoint::new(0i32, 0), + IntPoint::new(10, 0), + IntPoint::new(10, 10), + IntPoint::new(0, 10), + ]; + + contour.uniform_triangulate(1u64); + } +} diff --git a/iTriangle/src/lib.rs b/iTriangle/src/lib.rs index 8bb0e90..28a5ccf 100644 --- a/iTriangle/src/lib.rs +++ b/iTriangle/src/lib.rs @@ -1,6 +1,23 @@ #![no_std] extern crate alloc; +#[cfg(test)] +extern crate std; + +#[cfg(test)] +pub(crate) mod test_util { + const QUICK_RANDOM_DIVISOR: usize = 20; + const FULL_RANDOM_TESTS_ENV: &str = "ITRIANGLE_FULL_RANDOM_TESTS"; + + pub(crate) fn random_cases(full_count: usize) -> usize { + if std::env::var_os(FULL_RANDOM_TESTS_ENV).is_some() { + full_count + } else { + full_count.div_ceil(QUICK_RANDOM_DIVISOR) + } + } +} + pub mod advanced; pub mod float; pub mod geom; diff --git a/iTriangle/src/tessellation/mod.rs b/iTriangle/src/tessellation/mod.rs index 59ab6a7..619d373 100644 --- a/iTriangle/src/tessellation/mod.rs +++ b/iTriangle/src/tessellation/mod.rs @@ -1,2 +1,3 @@ pub mod circumcenter; pub mod split; +pub mod uniform; diff --git a/iTriangle/src/tessellation/split.rs b/iTriangle/src/tessellation/split.rs index 2e5e4f4..aa35dcd 100644 --- a/iTriangle/src/tessellation/split.rs +++ b/iTriangle/src/tessellation/split.rs @@ -6,6 +6,8 @@ use i_overlay::i_float::int::point::IntPoint; use i_overlay::i_shape::int::shape::{IntContour, IntShape, IntShapes}; pub trait SliceContour { + /// Splits every contour edge so that no resulting segment is longer than + /// `max_edge_length`. fn slice_contour(&self, max_edge_length: I::WideUInt) -> Self; } @@ -19,7 +21,7 @@ impl SliceContour for IntContour { }; let radius = max_edge_length; - if radius > I::WideUInt::HALF_MASK { + if radius == I::WideUInt::ZERO || radius > I::WideUInt::HALF_MASK { return self.clone(); } @@ -76,13 +78,8 @@ fn extract( contour.push(b); return; } - let len = i_overlay::i_float::float::number::FloatNumber::sqrt(sqr_len.to_f64()); - let n = ((len + 0.5 * radius.to_f64()) / radius.to_f64()) as usize; - if n <= 1 { - contour.push(b); - return; - } + let n = ((sqr_len - I::WideUInt::ONE).isqrt() / radius + I::WideUInt::ONE).to_usize(); if n == 2 { let x = I::from_wide((a.x.to_wide() + b.x.to_wide()) / I::Wide::TWO); @@ -120,18 +117,63 @@ mod tests { ]; let s0 = contour.slice_contour(8u64); - assert_eq!(s0.len(), 4); + assert_eq!(s0.len(), 8); + assert_max_edge_length(&s0, 8); let s1 = contour.slice_contour(7u64); - assert_eq!(s1.len(), 4); + assert_eq!(s1.len(), 8); + assert_max_edge_length(&s1, 7); let s2 = contour.slice_contour(6u64); assert_eq!(s2.len(), 8); + assert_max_edge_length(&s2, 6); let s3 = contour.slice_contour(5u64); assert_eq!(s3.len(), 8); + assert_max_edge_length(&s3, 5); let s4 = contour.slice_contour(3u64); - assert_eq!(s4.len(), 12); + assert_eq!(s4.len(), 16); + assert_max_edge_length(&s4, 3); + } + + #[test] + fn uses_integer_sqrt_segment_count() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(25, 0), + IntPoint::new(25, 10), + IntPoint::new(0, 10), + ]; + + let sliced = contour.slice_contour(6u64); + + // 25 / 6 produces five parts and 10 / 6 produces two parts. + assert_eq!(sliced.len(), 14); + } + + #[test] + fn exact_multiple_does_not_add_an_extra_segment() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(4, 0), + IntPoint::new(4, 2), + IntPoint::new(0, 2), + ]; + + let sliced = contour.slice_contour(2u64); + + assert_eq!(sliced.len(), 6); + } + + fn assert_max_edge_length(contour: &[IntPoint], max_edge_length: i32) { + let mut a = *contour.last().unwrap(); + let sqr_max = max_edge_length * max_edge_length; + for &b in contour { + let dx = b.x - a.x; + let dy = b.y - a.y; + assert!(dx * dx + dy * dy <= sqr_max); + a = b; + } } } diff --git a/iTriangle/src/tessellation/uniform.rs b/iTriangle/src/tessellation/uniform.rs new file mode 100644 index 0000000..005ee00 --- /dev/null +++ b/iTriangle/src/tessellation/uniform.rs @@ -0,0 +1,309 @@ +use alloc::vec::Vec; +use i_key_sort::sort::two_keys::TwoKeysSort; +use i_overlay::core::integer::OverlayInt; +use i_overlay::core::point_location::IntPointContainment; +use i_overlay::i_float::int::number::product_uint::UIntProduct; +use i_overlay::i_float::int::number::uint::UIntNumber; +use i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_overlay::i_float::int::point::IntPoint; +use i_overlay::i_float::int::rect::IntRect; +use i_overlay::i_shape::int::shape::{IntContour, IntShape}; + +const TRIANGLE_HEIGHT_NUMERATOR: u32 = 28_378; +const TRIANGLE_HEIGHT_SHIFT: u32 = 15; + +/// Generates vertices of an equilateral triangular lattice inside integer geometry. +/// +/// The input geometry must have resolved topology. A shape may contain holes and a +/// collection of shapes is treated as their union. Candidate points are tested in one +/// batch with [`IntPointContainment`], then points close to any boundary edge are removed. +pub trait IntUniformGrid { + /// Returns lattice points contained by the geometry. + /// + /// `edge_length` is the horizontal lattice spacing. Consecutive rows are separated by + /// `round(sqrt(3) / 2 * edge_length)` and shifted by half an edge. + fn uniform_grid(&self, edge_length: I::WideUInt) -> Vec>; +} + +impl IntUniformGrid for [IntPoint] { + #[inline] + fn uniform_grid(&self, edge_length: I::WideUInt) -> Vec> { + let mut edges = Vec::with_capacity(self.len()); + append_edges(self, &mut edges); + build_grid(self, self.iter(), &edges, edge_length) + } +} + +impl IntUniformGrid for [IntContour] { + #[inline] + fn uniform_grid(&self, edge_length: I::WideUInt) -> Vec> { + let mut edges = Vec::new(); + for contour in self { + append_edges(contour, &mut edges); + } + build_grid(self, self.iter().flatten(), &edges, edge_length) + } +} + +impl IntUniformGrid for [IntShape] { + #[inline] + fn uniform_grid(&self, edge_length: I::WideUInt) -> Vec> { + let mut edges = Vec::new(); + for contour in self.iter().flatten() { + append_edges(contour, &mut edges); + } + build_grid(self, self.iter().flatten().flatten(), &edges, edge_length) + } +} + +#[derive(Clone, Copy)] +struct Edge { + a: IntPoint, + b: IntPoint, +} + +fn append_edges(contour: &[IntPoint], edges: &mut Vec>) { + let Some(&mut_a) = contour.last() else { + return; + }; + let mut a = mut_a; + for &b in contour { + edges.push(Edge { a, b }); + a = b; + } +} + +fn build_grid<'a, I, G, It>( + geometry: &G, + points: It, + edges: &[Edge], + edge_length: I::WideUInt, +) -> Vec> +where + I: OverlayInt + 'a, + G: IntPointContainment + ?Sized, + It: Iterator>, +{ + let Some(rect) = IntRect::with_iter(points) else { + return Vec::new(); + }; + + let step = I::Wide::from_uint(edge_length); + if step <= I::Wide::ONE { + return Vec::new(); + } + + // Fixed-point approximation of sqrt(3) / 2, + // the height-to-edge ratio of an equilateral triangle. + let row_step = (step * I::Wide::from_u32(TRIANGLE_HEIGHT_NUMERATOR) + + I::Wide::from_u32(1 << (TRIANGLE_HEIGHT_SHIFT - 1))) + >> TRIANGLE_HEIGHT_SHIFT; + if row_step <= I::Wide::ZERO { + return Vec::new(); + } + + let half_step = step / I::Wide::TWO; + let min_x = rect.min_x.to_wide(); + let max_x = rect.max_x.to_wide(); + let max_y = rect.max_y.to_wide(); + + let mut candidates = Vec::new(); + let mut row = 0usize; + let mut y = rect.min_y.to_wide() + row_step / I::Wide::TWO; + + while y < max_y { + let row_offset = if row & 1 == 0 { half_step } else { step }; + let mut x = min_x + row_offset; + + while x < max_x { + candidates.push(IntPoint::new(I::from_wide(x), I::from_wide(y))); + x = x + step; + } + + row += 1; + y = y + row_step; + } + + let contains = geometry.contains_points(&candidates); + let contained = candidates + .into_iter() + .zip(contains) + .filter_map(|(point, is_inside)| is_inside.then_some(point)) + .collect(); + + let third = edge_length / I::WideUInt::from_u64(3); + let clearance = if third == I::WideUInt::ZERO { + I::WideUInt::ONE + } else { + third + }; + + filter_near_edges(contained, edges, rect, edge_length, clearance) +} + +fn filter_near_edges( + points: Vec>, + edges: &[Edge], + rect: IntRect, + cell_size: I::WideUInt, + clearance: I::WideUInt, +) -> Vec> { + if points.is_empty() || edges.is_empty() { + return points; + } + + let origin_x = rect.min_x.to_wide(); + let origin_y = rect.min_y.to_wide(); + let limit_x = rect.max_x.to_wide(); + let limit_y = rect.max_y.to_wide(); + let cell_size_wide = I::Wide::from_uint(cell_size); + let clearance_wide = I::Wide::from_uint(clearance); + + // (cell_y, cell_x, edge_index). Key sorting groups edge references by cell + // without relying on hashing in this no_std crate. + let mut cell_edges = Vec::new(); + for (edge_index, edge) in edges.iter().enumerate() { + let min_x = (edge.a.x.min(edge.b.x).to_wide() - clearance_wide).max(origin_x); + let max_x = (edge.a.x.max(edge.b.x).to_wide() + clearance_wide).min(limit_x); + let min_y = (edge.a.y.min(edge.b.y).to_wide() - clearance_wide).max(origin_y); + let max_y = (edge.a.y.max(edge.b.y).to_wide() + clearance_wide).min(limit_y); + + let min_cell_x = ((min_x - origin_x) / cell_size_wide).to_usize(); + let max_cell_x = ((max_x - origin_x) / cell_size_wide).to_usize(); + let min_cell_y = ((min_y - origin_y) / cell_size_wide).to_usize(); + let max_cell_y = ((max_y - origin_y) / cell_size_wide).to_usize(); + + for cell_y in min_cell_y..=max_cell_y { + for cell_x in min_cell_x..=max_cell_x { + cell_edges.push((cell_y, cell_x, edge_index)); + } + } + } + cell_edges.sort_by_two_keys(false, |entry| entry.0, |entry| entry.1); + + points + .into_iter() + .filter(|point| { + let cell_x = ((point.x.to_wide() - origin_x) / cell_size_wide).to_usize(); + let cell_y = ((point.y.to_wide() - origin_y) / cell_size_wide).to_usize(); + let cell = (cell_y, cell_x); + let start = cell_edges.partition_point(|entry| (entry.0, entry.1) < cell); + let end = cell_edges.partition_point(|entry| (entry.0, entry.1) <= cell); + + !cell_edges[start..end] + .iter() + .any(|entry| is_close_to_edge(*point, edges[entry.2], clearance)) + }) + .collect() +} + +#[inline] +fn is_close_to_edge( + point: IntPoint, + edge: Edge, + clearance: I::WideUInt, +) -> bool { + let ab = edge.b - edge.a; + let ap = point - edge.a; + let length_sqr = ab.sqr_length(); + let clearance_sqr = clearance * clearance; + + if length_sqr <= I::Wide::ZERO { + return ap.sqr_length().to_uint() <= clearance_sqr; + } + + let projection = ap.dot_product(ab); + if projection <= I::Wide::ZERO { + return ap.sqr_length().to_uint() <= clearance_sqr; + } + if projection >= length_sqr { + return (point - edge.b).sqr_length().to_uint() <= clearance_sqr; + } + + let cross = ab.cross_product(ap).unsigned_abs(); + let distance_product = ::Product::multiply(cross, cross); + let limit_product = + ::Product::multiply(clearance_sqr, length_sqr.to_uint()); + + distance_product <= limit_product +} + +#[cfg(test)] +mod tests { + use super::{filter_near_edges, Edge, IntUniformGrid}; + use alloc::vec; + use i_overlay::i_float::int::point::IntPoint; + use i_overlay::i_float::int::rect::IntRect; + + #[test] + fn square_grid_has_staggered_rows() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + + let points = contour.uniform_grid(20u64); + + assert!(!points.is_empty()); + assert!(points + .iter() + .all(|p| 0 < p.x && p.x < 100 && 0 < p.y && p.y < 100)); + assert!(points.iter().any(|p| p.x == 10)); + assert!(points.iter().any(|p| p.x == 20)); + } + + #[test] + fn shape_grid_excludes_hole() { + let shape = vec![ + vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ], + vec![ + IntPoint::new(40, 40), + IntPoint::new(40, 60), + IntPoint::new(60, 60), + IntPoint::new(60, 40), + ], + ]; + + let points = shape.uniform_grid(10u64); + + assert!(!points.is_empty()); + assert!(points + .iter() + .all(|p| p.x <= 40 || 60 <= p.x || p.y <= 40 || 60 <= p.y)); + } + + #[test] + fn grid_removes_points_in_edge_influence_across_cell_boundary() { + let contour = vec![ + IntPoint::new(0, 0), + IntPoint::new(100, 0), + IntPoint::new(100, 100), + IntPoint::new(0, 100), + ]; + + let points = contour.uniform_grid(20u64); + + assert!(points.iter().all(|p| 6 < p.x && p.x < 94)); + assert!(points.iter().all(|p| 6 < p.y && p.y < 94)); + } + + #[test] + fn edge_influences_a_cell_it_does_not_cross() { + let points = vec![IntPoint::new(9, 5), IntPoint::new(7, 5)]; + let edges = [Edge { + a: IntPoint::new(10, 0), + b: IntPoint::new(10, 10), + }]; + + let filtered = filter_near_edges(points, &edges, IntRect::new(0, 20, 0, 20), 10u64, 2u64); + + assert_eq!(filtered, vec![IntPoint::new(7, 5)]); + } +} diff --git a/iTriangle/tests/doc_tests.rs b/iTriangle/tests/doc_tests.rs index 727c9f8..f163f94 100644 --- a/iTriangle/tests/doc_tests.rs +++ b/iTriangle/tests/doc_tests.rs @@ -1,9 +1,14 @@ #[cfg(test)] mod tests { use i_overlay::i_shape::base::data::Contour; + use i_triangle::float::relax::RelaxationOptions; use i_triangle::float::triangulatable::Triangulatable; use i_triangle::float::triangulation::Triangulation; use i_triangle::float::triangulator::Triangulator; + use i_triangle::float::uniform::UniformTriangulatable; + use i_triangle::i_overlay::core::fill_rule::FillRule; + use i_triangle::i_overlay::core::overlay_rule::OverlayRule; + use i_triangle::i_overlay::float::single::SingleFloatOverlay; use rand::RngExt; #[test] @@ -99,6 +104,26 @@ mod tests { } } + #[test] + fn uniform_relaxation_and_centroid_net() { + let contours = vec![ + vec![[0.0, 0.0], [12.0, 0.0], [12.0, 8.0], [0.0, 8.0]], + vec![[4.0, 2.0], [8.0, 2.0], [8.0, 6.0], [4.0, 6.0]], + ]; + let empty: Vec> = Vec::new(); + let shapes = contours.overlay(&empty, OverlayRule::Union, FillRule::EvenOdd); + let shape = &shapes[0]; + + let mut delaunay = shape.uniform_triangulate(1.0); + let relaxation = delaunay.relax_mut(RelaxationOptions::new(24)); + let triangles = delaunay.to_triangulation::(); + let centroid_net = delaunay.to_centroid_net(0.0); + + assert!(!triangles.indices.is_empty()); + assert!(!centroid_net.is_empty()); + assert!(relaxation.iterations <= 24); + } + fn random_contours(count: usize) -> Vec> { let mut contours = Vec::with_capacity(count); for _ in 0..count { diff --git a/iTriangle/tests/float_tests.rs b/iTriangle/tests/float_tests.rs index 806f2de..1159cd7 100644 --- a/iTriangle/tests/float_tests.rs +++ b/iTriangle/tests/float_tests.rs @@ -14,6 +14,17 @@ mod tests { impl TestInt for I {} + const QUICK_RANDOM_DIVISOR: usize = 20; + const FULL_RANDOM_TESTS_ENV: &str = "ITRIANGLE_FULL_RANDOM_TESTS"; + + fn random_cases(full_count: usize) -> usize { + if std::env::var_os(FULL_RANDOM_TESTS_ENV).is_some() { + full_count + } else { + full_count.div_ceil(QUICK_RANDOM_DIVISOR) + } + } + #[test] fn test_0() { test_0_as::(); @@ -132,7 +143,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let contour = random(8, 5); let area = contour.simplify_shape(FillRule::NonZero).area(); @@ -157,7 +168,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let contour = random(10, 6); let area = contour.simplify_shape(FillRule::NonZero).area(); @@ -182,7 +193,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..20_000 { + for _ in 0..random_cases(20_000) { let contour = random(10, 12); let area = contour.simplify_shape(FillRule::NonZero).area(); @@ -207,7 +218,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..10_000 { + for _ in 0..random_cases(10_000) { let contour = random(20, 20); let area = contour.simplify_shape(FillRule::NonZero).area(); @@ -232,7 +243,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..1_000 { + for _ in 0..random_cases(1_000) { let contour = random(30, 50); let area = contour.simplify_shape(FillRule::NonZero).area(); @@ -257,7 +268,7 @@ mod tests { let mut triangulator = Triangulator::::default(); let mut t = Triangulation::with_capacity(8); - for _ in 0..500 { + for _ in 0..random_cases(500) { let main = random(50, 20); let mut shape = vec![main]; for _ in 0..10 { diff --git a/iTriangle/tests/relax_stress.rs b/iTriangle/tests/relax_stress.rs new file mode 100644 index 0000000..9f1830b --- /dev/null +++ b/iTriangle/tests/relax_stress.rs @@ -0,0 +1,165 @@ +use i_triangle::float::relax::RelaxationOptions; +use i_triangle::float::uniform::UniformTriangulatable; +use i_triangle::i_overlay::core::integer::OverlayInt; + +const DEFAULT_CASE_COUNT: usize = 3_000; +const CASE_COUNT_ENV: &str = "ITRIANGLE_RELAX_STRESS_CASES"; + +#[test] +#[ignore = "deterministic stress test; run with `cargo test --test relax_stress -- --ignored`"] +fn self_intersecting_uniform_relax_stress() { + let case_count = std::env::var(CASE_COUNT_ENV) + .ok() + .and_then(|value| value.parse().ok()) + .unwrap_or(DEFAULT_CASE_COUNT); + + for case in 0..case_count { + match case % 3 { + 0 => run_case::(case), + 1 => run_case::(case), + _ => run_case::(case), + } + } +} + +fn run_case(case: usize) { + let mut rng = Rng::new(case as u64 + 1); + let vertex_count = [5, 7, 9, 11, 13][case % 5]; + let radius = 30.0 + 70.0 * rng.unit(); + let center = [200.0 * rng.signed_unit(), 200.0 * rng.signed_unit()]; + let angle_step = core::f64::consts::TAU / vertex_count as f64; + + let mut ring = Vec::with_capacity(vertex_count); + for index in 0..vertex_count { + let angle = index as f64 * angle_step + 0.16 * angle_step * rng.signed_unit(); + let local_radius = radius * (0.82 + 0.36 * rng.unit()); + ring.push([ + center[0] + local_radius * angle.cos(), + center[1] + local_radius * angle.sin(), + ]); + } + + // Visiting every second point of an odd ring creates a single + // self-intersecting star contour. + let mut shape = Vec::with_capacity(vertex_count); + let mut index = 0; + for _ in 0..vertex_count { + shape.push(ring[index]); + index = (index + 2) % vertex_count; + } + + let edge_length = radius / (4 + case % 6) as f64; + let relax_iterations = 1 + case % 8; + let mut delaunay = shape.uniform_triangulate_as::(edge_length); + + let before = delaunay.points(); + let indices_before = delaunay.triangle_indices::(); + let neighbors_before = delaunay.triangle_neighbors(); + assert!(!indices_before.is_empty(), "case {case}: empty mesh"); + + let boundary = boundary_vertices(before.len(), &indices_before, &neighbors_before); + let area_before = mesh_area(case, &before, &indices_before); + + let result = delaunay.relax_mut(RelaxationOptions::new(relax_iterations)); + + assert!( + result.iterations <= relax_iterations, + "case {case}: iteration limit exceeded" + ); + assert_eq!( + result.converged, + result.iterations < relax_iterations, + "case {case}: inconsistent convergence result" + ); + + let after = delaunay.points(); + let indices_after = delaunay.triangle_indices::(); + assert_eq!( + after.len(), + before.len(), + "case {case}: point count changed" + ); + assert_eq!( + indices_after.len(), + indices_before.len(), + "case {case}: triangle count changed" + ); + + for (index, is_boundary) in boundary.into_iter().enumerate() { + if is_boundary { + assert_eq!(after[index], before[index], "case {case}: boundary moved"); + } + } + + assert_no_new_collisions(case, &before, &after); + + let area_after = mesh_area(case, &after, &indices_after); + let area_tolerance = 1.0e-10 * area_before.max(1.0); + assert!( + (area_after - area_before).abs() <= area_tolerance, + "case {case}: area changed from {area_before} to {area_after}" + ); +} + +fn boundary_vertices(point_count: usize, indices: &[u32], neighbors: &[[usize; 3]]) -> Vec { + let triangle_count = indices.len() / 3; + assert_eq!(neighbors.len(), triangle_count); + + let mut result = vec![false; point_count]; + for (triangle_index, triangle) in indices.chunks_exact(3).enumerate() { + for edge in 0..3 { + if neighbors[triangle_index][edge] >= triangle_count { + result[triangle[(edge + 1) % 3] as usize] = true; + result[triangle[(edge + 2) % 3] as usize] = true; + } + } + } + result +} + +fn mesh_area(case: usize, points: &[[f64; 2]], indices: &[u32]) -> f64 { + let mut area_two = 0.0; + for triangle in indices.chunks_exact(3) { + let a = points[triangle[0] as usize]; + let b = points[triangle[1] as usize]; + let c = points[triangle[2] as usize]; + let triangle_area_two = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]); + assert!( + triangle_area_two > 0.0, + "case {case}: non-positive triangle area" + ); + area_two += triangle_area_two; + } + 0.5 * area_two +} + +fn assert_no_new_collisions(case: usize, before: &[[f64; 2]], after: &[[f64; 2]]) { + for left in 0..after.len() { + for right in left + 1..after.len() { + assert!( + after[left] != after[right] || before[left] == before[right], + "case {case}: vertices {left} and {right} collided" + ); + } + } +} + +struct Rng(u64); + +impl Rng { + fn new(seed: u64) -> Self { + Self(seed) + } + + fn unit(&mut self) -> f64 { + self.0 = self + .0 + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + ((self.0 >> 11) as f64) * (1.0 / ((1u64 << 53) as f64)) + } + + fn signed_unit(&mut self) -> f64 { + 2.0 * self.unit() - 1.0 + } +} diff --git a/readme/eagle_centroid.svg b/readme/eagle_centroid.svg index ac5bd65..8305649 100644 --- a/readme/eagle_centroid.svg +++ b/readme/eagle_centroid.svg @@ -1,1547 +1,577 @@ - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Centroid net of a relaxed uniform Delaunay mesh + + + diff --git a/readme/eagle_tessellation.svg b/readme/eagle_tessellation.svg index 71f594a..6c504ec 100644 --- a/readme/eagle_tessellation.svg +++ b/readme/eagle_tessellation.svg @@ -1,3849 +1,1403 @@ - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Relaxed uniform Delaunay tessellation + + +