From c565396c9c07b89f6806f49bc8d010028e404aff Mon Sep 17 00:00:00 2001 From: Nail Sharipov Date: Fri, 30 Oct 2026 16:21:47 +0300 Subject: [PATCH] minkowski ready --- iTriangle/Cargo.toml | 2 +- iTriangle/README.md | 71 +++ iTriangle/src/float/minkowski.rs | 514 +++++++++++++++++ iTriangle/src/float/mod.rs | 1 + iTriangle/src/int/minkowski.rs | 379 +++++++++++++ iTriangle/src/int/minkowski_sum.rs | 251 +++++++++ iTriangle/src/int/mod.rs | 2 + iTriangle/tests/float_minkowski_tests.rs | 524 ++++++++++++++++++ iTriangle/tests/minkowski_difference_tests.rs | 200 +++++++ iTriangle/tests/minkowski_operand_tests.rs | 90 +++ iTriangle/tests/minkowski_sum_tests.rs | 389 +++++++++++++ 11 files changed, 2422 insertions(+), 1 deletion(-) create mode 100644 iTriangle/src/float/minkowski.rs create mode 100644 iTriangle/src/int/minkowski.rs create mode 100644 iTriangle/src/int/minkowski_sum.rs create mode 100644 iTriangle/tests/float_minkowski_tests.rs create mode 100644 iTriangle/tests/minkowski_difference_tests.rs create mode 100644 iTriangle/tests/minkowski_operand_tests.rs create mode 100644 iTriangle/tests/minkowski_sum_tests.rs diff --git a/iTriangle/Cargo.toml b/iTriangle/Cargo.toml index b24c762..7b38852 100644 --- a/iTriangle/Cargo.toml +++ b/iTriangle/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "i_triangle" -version = "0.50.1" +version = "0.51.0" edition = "2021" authors = ["Nail Sharipov "] description = "Polygon Triangulation Library: Efficient Delaunay Triangulation for Complex Shapes." diff --git a/iTriangle/README.md b/iTriangle/README.md index cb996f8..5412824 100644 --- a/iTriangle/README.md +++ b/iTriangle/README.md @@ -23,6 +23,7 @@ iTriangle provides reliable, deterministic 2D polygon triangulation for complex - [Documentation](#documentation) - [Examples](#examples) - [Integer API](#integer-api) +- [Minkowski Operations](#minkowski-operations) - [Performance](#performance) - [Gallery](#gallery) - [Contributing](#contributing) @@ -44,6 +45,7 @@ iTriangle provides reliable, deterministic 2D polygon triangulation for complex - **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. +- **Minkowski Sum and Difference** - Prepare and reuse filled polygon regions, including holes and disconnected components. - **Centroidal Polygon Net**: Build per-vertex dual polygons using triangle centers and edge midpoints. - **Steiner Points**: Add custom inner points to influence triangulation. - **GPU-Friendly Layout**: Triangles and vertices are naturally ordered by X due to the sweep-line algorithm, improving cache locality for rendering. @@ -349,6 +351,75 @@ let triangulation = contour.uncheck_triangulate().into_triangulation::(); assert_eq!(triangulation.indices.len(), 6); ``` +## Minkowski Operations + +The Minkowski sum is `A + B = {a + b}`. Minkowski difference here means +`A + (-B) = {a - b}`: it reflects the second region through the origin. It is +not erosion or Boolean subtraction, and it does not undo a sum. + +For a single operation, the convenience functions use `i32` and automatically +choose a common domain containing both sources and the result. Both sources use +the same fill rule: + +```rust +use i_triangle::float::minkowski::{minkowski_difference, minkowski_sum}; +use i_triangle::i_overlay::core::fill_rule::FillRule; + +let square = [[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]]; +let sum = minkowski_sum(&square, &square, FillRule::NonZero).unwrap(); +let difference = minkowski_difference(&square, &square, FillRule::NonZero).unwrap(); +assert!(sum[0][0].contains(&[4.0, 4.0])); +assert!(difference[0][0].contains(&[-2.0, -2.0])); + +let empty: Vec<[f64; 2]> = Vec::new(); +assert!(minkowski_sum(&square, &empty, FillRule::NonZero).unwrap().is_empty()); +``` + +For repeated operations, select a fixed domain that contains the original +sources and every requested result. Each source is converted, normalized with +its fill rule, triangulated, and decomposed only during preparation: + +```rust +use i_triangle::float::minkowski::MinkowskiSolver; +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::i_float::float::rect::FloatRect; + +let domain = FloatRect::new(-16.0, 16.0, -16.0, 16.0).unwrap(); +let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain); +let source = [[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]; +let kernel = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; +let a = solver.prepare(&source, FillRule::NonZero).unwrap().unwrap(); +let b = solver.prepare(&kernel, FillRule::NonZero).unwrap().unwrap(); +let sum = solver.sum(&a, &b).unwrap(); +solver.clear(); +let difference = solver.difference(&a, &b).unwrap(); +assert!(sum[0][0].contains(&[5.0, 5.0])); +assert!(difference[0][0].contains(&[-1.0, -1.0])); +``` + +`prepare` returns `Ok(None)` for empty or unfilled regions, points, segments, +and geometry collapsed by the integer grid. A prepared float operand belongs +to its creating solver; another solver returns `MinkowskiError::ForeignOperand` +even with identical bounds. Moving or clearing the original solver preserves +its operands. Original source and result bounds outside the domain return +`OutsideDomain`; invalid coordinates or unsupported numeric result bounds return +`CoordinateOutOfRange`. Invalid constructor bounds panic. + +The adapter's grid stays fixed for the solver's lifetime. A larger domain can +reduce precision; `MinkowskiSolver::<[f64; 2], i64>` provides a larger integer +budget. Inspect `solver.adapter().inv_scale()` for the grid step. Rounding can +alter small features, and snapped output can lie slightly outside the domain. +`with_options(domain, capacity, true)` retains collinear vertices in the final +union; pairwise collinear edges are always merged, so original source vertices +are not all preserved. + +For already quantized geometry, use +`i_triangle::int::minkowski::{IntMinkowskiOperand, IntMinkowskiSolver}`. Integer +operands can be reused across integer solvers of the same coordinate type. +Input must satisfy the overlay integer coordinate range; conservative result +bounds are checked before pairwise operations. See the module rustdoc for +executable examples and the full range and precision contracts. + ## Performance Benchmarks and interactive demos are available here: diff --git a/iTriangle/src/float/minkowski.rs b/iTriangle/src/float/minkowski.rs new file mode 100644 index 0000000..eba775b --- /dev/null +++ b/iTriangle/src/float/minkowski.rs @@ -0,0 +1,514 @@ +//! Minkowski sums and differences on a shared, fixed integer grid. +//! +//! The sum is `A + B = {a + b | a in A, b in B}`. Here, **difference** means +//! `A + (-B) = {a - b | a in A, b in B}`. It is neither polygon subtraction +//! nor morphological erosion, and does not undo a sum. +//! +//! Use [`minkowski_sum`] or [`minkowski_difference`] for a single operation. +//! They choose a domain containing both sources and their result and use `i32` +//! coordinates. For repeated operations, choose a domain once and prepare each +//! source with [`MinkowskiSolver::prepare`]. Preparation applies the selected +//! [`FillRule`], resolves intersections, triangulates, and decomposes the filled +//! region into convex parts. Operations reuse these parts without converting +//! or triangulating the sources again. +//! +//! # Fixed domain and reusable operands +//! +//! The domain must contain the original source bounds and the bounds of every +//! requested result. Bounds checks use the original sources, including points +//! in contours later removed by the fill rule. An operand belongs to the solver +//! that prepared it; another solver rejects it even with an identical domain. +//! Moving or clearing a solver does not invalidate its operands. +//! +//! ``` +//! use i_triangle::float::minkowski::{MinkowskiError, MinkowskiSolver}; +//! use i_triangle::i_overlay::core::fill_rule::FillRule; +//! use i_triangle::i_overlay::i_float::float::rect::FloatRect; +//! +//! let domain = FloatRect::new(-16.0, 16.0, -16.0, 16.0).unwrap(); +//! let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain); +//! let source = [[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]; +//! let kernel = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; +//! let a = solver.prepare(&source, FillRule::NonZero)?.unwrap(); +//! let b = solver.prepare(&kernel, FillRule::NonZero)?.unwrap(); +//! let sum = solver.sum(&a, &b)?; +//! solver.clear(); +//! let difference = solver.difference(&a, &b)?; +//! assert_eq!(sum[0][0].len(), 4); +//! assert!(sum[0][0].contains(&[5.0, 5.0])); +//! assert!(difference[0][0].contains(&[-1.0, -1.0])); +//! # Ok::<(), MinkowskiError>(()) +//! ``` +//! +//! # Empty regions and errors +//! +//! Only filled, two-dimensional regions are represented. Empty input, an +//! unfilled source, or a region collapsed by quantization produces `Ok(None)` +//! during preparation. Points and segments are not operands. The one-shot +//! functions return an empty vector when either prepared region is empty. +//! Invalid coordinates are still validated; an empty source does not exempt +//! the other source from coordinate validation. +//! +//! ``` +//! use i_triangle::float::minkowski::{MinkowskiError, MinkowskiSolver}; +//! use i_triangle::i_overlay::core::fill_rule::FillRule; +//! use i_triangle::i_overlay::i_float::float::rect::FloatRect; +//! +//! let domain = FloatRect::new(-2.0, 2.0, -2.0, 2.0).unwrap(); +//! let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain); +//! let empty: Vec<[f64; 2]> = Vec::new(); +//! assert!(solver.prepare(&empty, FillRule::NonZero)?.is_none()); +//! let square = [[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0]]; +//! let a = solver.prepare(&square, FillRule::NonZero)?.unwrap(); +//! assert!(matches!(solver.sum(&a, &a), Err(MinkowskiError::OutsideDomain))); +//! let mut other = MinkowskiSolver::<[f64; 2]>::new(domain); +//! let b = other.prepare(&square, FillRule::NonZero)?.unwrap(); +//! assert!(matches!(solver.sum(&a, &b), Err(MinkowskiError::ForeignOperand))); +//! let invalid = [[f64::NAN, 0.0], [1.0, 0.0], [0.0, 1.0]]; +//! assert!(matches!(solver.prepare(&invalid, FillRule::NonZero), +//! Err(MinkowskiError::CoordinateOutOfRange))); +//! # Ok::<(), MinkowskiError>(()) +//! ``` +//! +//! # Precision and collinear vertices +//! +//! All operands share the solver's power-of-two scale and origin. A larger +//! domain generally means a coarser grid; a larger integer type can retain more +//! detail, subject to the precision of the floating-point scalar. The grid step +//! is [`FloatPointAdapter::inv_scale`], available through +//! [`MinkowskiSolver::adapter`]. Input rounding and integer intersection rounding +//! can alter small features or collapse them. Results approximate the operation +//! on the original floating-point regions; they are not exact real arithmetic. +//! Snapped output can extend slightly beyond the original domain bounds. +//! +//! [`MinkowskiSolver::with_options`] controls collinear vertices in the final +//! union. Enabling preservation retains vertices introduced by that union; it +//! does not retain every source vertex. Collinear edges of each convex pair +//! are merged before union, regardless of this setting. Output shapes contain +//! a counter-clockwise outer contour followed by clockwise hole contours. + +pub use crate::int::minkowski::MinkowskiError; +use crate::int::minkowski::{IntMinkowskiOperand, IntMinkowskiSolver}; +use alloc::sync::Arc; +use alloc::vec::Vec; +use i_overlay::core::fill_rule::FillRule; +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_float::int::rect::IntRect; +use i_overlay::i_shape::base::data::Shapes; +use i_overlay::i_shape::flat::buffer::FlatContoursBuffer; +use i_overlay::i_shape::int::shape::IntShapes; +use i_overlay::i_shape::source::float::resource::ShapeResource; + +/// An owned, prepared filled region on its creating solver's fixed grid. +/// +/// Create operands with [`MinkowskiSolver::prepare`]. The source need not remain +/// alive. Operands remain usable after further preparations, operations, solver +/// moves, and [`MinkowskiSolver::clear`], but only with their original solver. +pub struct MinkowskiOperand { + owner: Arc<()>, + bbox: FloatRect, + int_operand: IntMinkowskiOperand, +} + +impl MinkowskiOperand { + /// Borrows the prepared integer geometry in adapter coordinates. + /// + /// These coordinates have the adapter's origin subtracted and scale applied; + /// they are not coordinates in the original floating-point space. + pub fn int_operand(&self) -> &IntMinkowskiOperand { + &self.int_operand + } + + /// Returns the original source bounds, before quantization or fill filtering. + /// + /// Operations use these conservative bounds for domain and range checks. + pub fn bbox(&self) -> FloatRect { + self.bbox + } +} + +/// Reuses prepared regions and work buffers within one fixed coordinate domain. +/// +/// `P` may be a compatible `f32` or `f64` point type. `I` selects the integer +/// engine (`i16`, `i32`, or `i64`), with `i32` as the default. The adapter never +/// rescales after construction. See the [module documentation](self) for reuse, +/// ownership, empty-region, and precision contracts. +pub struct MinkowskiSolver { + owner: Arc<()>, + adapter: FloatPointAdapter, + int_solver: IntMinkowskiSolver, + source_buffer: FlatContoursBuffer, +} + +impl MinkowskiSolver { + /// Creates a solver with initial capacity 64 and collinear output simplified. + /// + /// `rect` must contain both sources and all requested result bounds. + /// + /// # Panics + /// Panics if `rect` has reversed bounds, non-finite coordinates, or absolute + /// coordinates above `2^60` for `f32` or `2^500` for `f64`. + pub fn new(rect: FloatRect) -> Self { + Self::with_options(rect, 64, false) + } + + /// Creates a solver with a fixed domain, capacity hint, and output policy. + /// + /// `capacity` is an initial allocation hint, not a limit. When + /// `preserve_output_collinear` is `true`, collinear vertices from the final + /// union are retained. Pairwise convex sums still merge collinear edges, + /// so this option does not preserve all original source vertices. + /// + /// The grid reserves an extra coordinate bit for adding or subtracting two + /// prepared operands: input magnitude is bounded by `2^(I::BITS - 4)`. + /// + /// # Panics + /// Panics for invalid `rect`, as in [`Self::new`]. + pub fn with_options( + rect: FloatRect, + capacity: usize, + preserve_output_collinear: bool, + ) -> Self { + Self { + owner: Arc::new(()), + adapter: FloatPointAdapter::with_coordinate_bits(rect, I::BITS - 4), + int_solver: IntMinkowskiSolver::with_options(capacity, preserve_output_collinear), + source_buffer: FlatContoursBuffer::with_capacity(capacity), + } + } + + /// Returns the fixed adapter, including its domain, origin, and grid scale. + /// + /// Operations restore sums with twice the adapter origin and differences + /// with zero origin. Directly applying `int_to_float` to an integer sum or + /// difference therefore does not in general reproduce the solver's output. + pub fn adapter(&self) -> &FloatPointAdapter { + &self.adapter + } + + /// Clears temporary geometry while retaining allocations and the fixed grid. + /// Previously prepared operands remain valid for this solver. + pub fn clear(&mut self) { + self.int_solver.clear(); + self.source_buffer.clear_and_reserve(0, 0); + } + + /// Quantizes, normalizes with `fill_rule`, and decomposes a source once. + /// + /// Accepts contours, shapes with holes, and collections of shapes through + /// [`ShapeResource`]. The returned operand owns its convex parts and can be + /// reused in any number of operations with this solver. Each source may use + /// a different fill rule. + /// + /// Returns `Ok(None)` when no filled area survives normalization on the + /// fixed grid, including empty, point, segment, and collapsed input. + /// + /// # Errors + /// Returns [`MinkowskiError::CoordinateOutOfRange`] for any non-finite source + /// coordinate or coordinate beyond the scalar's supported limit. Returns + /// [`MinkowskiError::OutsideDomain`] if the original source bounds extend + /// outside the solver's domain, even if those points would later be removed. + pub fn prepare( + &mut self, + source: &R, + fill_rule: FillRule, + ) -> Result>, MinkowskiError> + where + R: ShapeResource

+ ?Sized, + { + self.clear(); + let Some(bbox) = source_bbox(source)? else { + return Ok(None); + }; + self.check_rect(bbox)?; + let (points, paths) = source + .iter_paths() + .fold((0, 0), |(n, c), path| (n + path.len(), c + 1)); + self.source_buffer.clear_and_reserve(points, paths); + for path in source.iter_paths() { + let start = self.source_buffer.points.len(); + self.source_buffer + .points + .extend(path.iter().map(|p| self.adapter.float_to_int(p))); + self.source_buffer + .ranges + .push(start..self.source_buffer.points.len()); + } + Ok(self + .int_solver + .prepare(&self.source_buffer, fill_rule) + .map(|int_operand| MinkowskiOperand { + owner: Arc::clone(&self.owner), + bbox, + int_operand, + })) + } + + /// Returns `lhs + rhs = {a + b | a in lhs, b in rhs}`. + /// + /// Reuses the prepared parts and unions all pairwise sums. No source + /// conversion or triangulation is repeated. + /// + /// # Errors + /// Returns [`MinkowskiError::ForeignOperand`] if either operand was prepared + /// by another solver, [`MinkowskiError::OutsideDomain`] if the sum of the + /// original source bounds exceeds the fixed domain, or + /// [`MinkowskiError::CoordinateOutOfRange`] if numeric result bounds exceed + /// the supported coordinate range. The solver remains reusable on error. + pub fn sum( + &mut self, + lhs: &MinkowskiOperand, + rhs: &MinkowskiOperand, + ) -> Result, MinkowskiError> { + self.clear(); + self.check_operand(lhs)?; + self.check_operand(rhs)?; + self.check_rect(result_bbox(lhs.bbox, rhs.bbox, false)?)?; + let a = lhs.int_operand.bbox(); + let b = rhs.int_operand.bbox(); + let origin = self.adapter.offset(); + let offset = P::from_xy(origin.x() + origin.x(), origin.y() + origin.y()); + self.check_output_bounds( + IntRect::new( + a.min_x + b.min_x, + a.max_x + b.max_x, + a.min_y + b.min_y, + a.max_y + b.max_y, + ), + offset, + )?; + let shapes = self.int_solver.sum(&lhs.int_operand, &rhs.int_operand)?; + Ok(self.to_float(shapes, offset)) + } + + /// Returns `lhs + (-rhs) = {a - b | a in lhs, b in rhs}`. + /// + /// The right operand is reflected through the coordinate origin. This is + /// not erosion, Boolean subtraction, or an inverse of [`Self::sum`]. It + /// generally depends on operand order. Sources are not triangulated again. + /// + /// # Errors + /// Uses the same ownership, domain, and numeric checks as [`Self::sum`], + /// with result bounds computed by subtracting the original source bounds. + pub fn difference( + &mut self, + lhs: &MinkowskiOperand, + rhs: &MinkowskiOperand, + ) -> Result, MinkowskiError> { + self.clear(); + self.check_operand(lhs)?; + self.check_operand(rhs)?; + self.check_rect(result_bbox(lhs.bbox, rhs.bbox, true)?)?; + let a = lhs.int_operand.bbox(); + let b = rhs.int_operand.bbox(); + let offset = P::from_xy(P::Scalar::ZERO, P::Scalar::ZERO); + self.check_output_bounds( + IntRect::new( + a.min_x - b.max_x, + a.max_x - b.min_x, + a.min_y - b.max_y, + a.max_y - b.min_y, + ), + offset, + )?; + let shapes = self + .int_solver + .difference(&lhs.int_operand, &rhs.int_operand)?; + Ok(self.to_float(shapes, offset)) + } + + fn check_rect(&self, bbox: FloatRect) -> Result<(), MinkowskiError> { + let rect = self.adapter.rect(); + if bbox.min_x < rect.min_x + || bbox.max_x > rect.max_x + || bbox.min_y < rect.min_y + || bbox.max_y > rect.max_y + { + Err(MinkowskiError::OutsideDomain) + } else { + Ok(()) + } + } + + fn check_operand(&self, operand: &MinkowskiOperand) -> Result<(), MinkowskiError> { + if !Arc::ptr_eq(&self.owner, &operand.owner) { + return Err(MinkowskiError::ForeignOperand); + } + self.check_rect(operand.bbox) + } + + fn check_output_bounds(&self, bbox: IntRect, offset: P) -> Result<(), MinkowskiError> { + let inv = self.adapter.inv_scale(); + FloatRect::new( + P::Scalar::from_int(bbox.min_x) * inv + offset.x(), + P::Scalar::from_int(bbox.max_x) * inv + offset.x(), + P::Scalar::from_int(bbox.min_y) * inv + offset.y(), + P::Scalar::from_int(bbox.max_y) * inv + offset.y(), + ) + .map(|_| ()) + .map_err(|_| MinkowskiError::CoordinateOutOfRange) + } + + fn to_float(&self, shapes: IntShapes, offset: P) -> Shapes

{ + let inv = self.adapter.inv_scale(); + shapes + .into_iter() + .map(|shape| { + shape + .into_iter() + .map(|contour| { + contour + .into_iter() + .map(|p| { + P::from_xy( + P::Scalar::from_int(p.x) * inv + offset.x(), + P::Scalar::from_int(p.y) * inv + offset.y(), + ) + }) + .collect() + }) + .collect() + }) + .collect() + } +} + +/// Computes a one-shot Minkowski sum using an automatically selected `i32` grid. +/// +/// Both sources use `fill_rule`. The domain encloses their original bounds and +/// the sum bounds. Use [`MinkowskiSolver`] to choose a domain or integer type, +/// retain output collinear vertices, or reuse prepared operands. +/// +/// Returns an empty vector if either normalized source has no filled area, +/// including geometry collapsed by the common grid. +/// +/// # Errors +/// Returns [`MinkowskiError::CoordinateOutOfRange`] for invalid source +/// coordinates or nonempty result bounds beyond the supported numeric range. +/// +/// # Example +/// ``` +/// use i_triangle::float::minkowski::{minkowski_sum, MinkowskiError}; +/// use i_triangle::i_overlay::core::fill_rule::FillRule; +/// +/// let square = [[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]]; +/// let result = minkowski_sum(&square, &square, FillRule::NonZero)?; +/// assert!(result[0][0].contains(&[4.0, 4.0])); +/// let empty: Vec<[f64; 2]> = Vec::new(); +/// assert!(minkowski_sum(&square, &empty, FillRule::NonZero)?.is_empty()); +/// # Ok::<(), MinkowskiError>(()) +/// ``` +pub fn minkowski_sum( + lhs: &R0, + rhs: &R1, + fill_rule: FillRule, +) -> Result, MinkowskiError> +where + P: FloatPointCompatible, + R0: ShapeResource

+ ?Sized, + R1: ShapeResource

+ ?Sized, +{ + calculate(lhs, rhs, fill_rule, false) +} + +/// Computes `{a - b | a in lhs, b in rhs}` on an automatically selected `i32` grid. +/// +/// This is the sum with `rhs` reflected through the origin, not erosion or +/// Boolean subtraction. It does not undo [`minkowski_sum`]. Both sources use +/// `fill_rule`; the domain encloses both original sources and their difference. +/// Empty-region behavior and numeric errors are the same as [`minkowski_sum`]. +/// +/// # Example +/// ``` +/// use i_triangle::float::minkowski::{minkowski_difference, MinkowskiError}; +/// use i_triangle::i_overlay::core::fill_rule::FillRule; +/// +/// let square = [[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]]; +/// let result = minkowski_difference(&square, &square, FillRule::NonZero)?; +/// assert!(result[0][0].contains(&[-2.0, -2.0])); +/// assert!(result[0][0].contains(&[2.0, 2.0])); +/// # Ok::<(), MinkowskiError>(()) +/// ``` +pub fn minkowski_difference( + lhs: &R0, + rhs: &R1, + fill_rule: FillRule, +) -> Result, MinkowskiError> +where + P: FloatPointCompatible, + R0: ShapeResource

+ ?Sized, + R1: ShapeResource

+ ?Sized, +{ + calculate(lhs, rhs, fill_rule, true) +} + +fn source_bbox(source: &R) -> Result>, MinkowskiError> +where + P: FloatPointCompatible, + R: ShapeResource

+ ?Sized, +{ + FloatRect::with_iter(source.iter_paths().flatten()) + .map_err(|_| MinkowskiError::CoordinateOutOfRange) +} + +fn result_bbox( + a: FloatRect, + b: FloatRect, + difference: bool, +) -> Result, MinkowskiError> { + let rect = if difference { + FloatRect::new( + a.min_x - b.max_x, + a.max_x - b.min_x, + a.min_y - b.max_y, + a.max_y - b.min_y, + ) + } else { + FloatRect::new( + a.min_x + b.min_x, + a.max_x + b.max_x, + a.min_y + b.min_y, + a.max_y + b.max_y, + ) + }; + rect.map_err(|_| MinkowskiError::CoordinateOutOfRange) +} + +fn calculate( + lhs: &R0, + rhs: &R1, + fill_rule: FillRule, + difference: bool, +) -> Result, MinkowskiError> +where + P: FloatPointCompatible, + R0: ShapeResource

+ ?Sized, + R1: ShapeResource

+ ?Sized, +{ + let a = source_bbox(lhs)?; + let b = source_bbox(rhs)?; + let (Some(a), Some(b)) = (a, b) else { + return Ok(Vec::new()); + }; + let inputs = FloatRect::with_rects(a, b).map_err(|_| MinkowskiError::CoordinateOutOfRange)?; + let rect = match result_bbox(a, b, difference) { + Ok(output) => FloatRect::with_rects(inputs, output) + .map_err(|_| MinkowskiError::CoordinateOutOfRange)?, + Err(_) => inputs, + }; + let mut solver = MinkowskiSolver::

::new(rect); + let a = solver.prepare(lhs, fill_rule)?; + let b = solver.prepare(rhs, fill_rule)?; + let (Some(a), Some(b)) = (a, b) else { + return Ok(Vec::new()); + }; + if difference { + solver.difference(&a, &b) + } else { + solver.sum(&a, &b) + } +} diff --git a/iTriangle/src/float/mod.rs b/iTriangle/src/float/mod.rs index 51e07f0..dd121d4 100644 --- a/iTriangle/src/float/mod.rs +++ b/iTriangle/src/float/mod.rs @@ -15,6 +15,7 @@ pub mod convex; pub mod custom; pub mod delaunay; pub mod locator; +pub mod minkowski; pub mod relax; pub mod triangulatable; pub mod triangulation; diff --git a/iTriangle/src/int/minkowski.rs b/iTriangle/src/int/minkowski.rs new file mode 100644 index 0000000..fe35061 --- /dev/null +++ b/iTriangle/src/int/minkowski.rs @@ -0,0 +1,379 @@ +//! Minkowski operations on prepared, filled integer regions. +//! +//! Preparation resolves the selected [`FillRule`], triangulates the normalized +//! region, and groups triangles into convex parts. Sums and differences merge +//! each pair's edge directions and union the resulting polygons. Concavities, +//! holes, and disconnected components are handled through this decomposition. +//! +//! Difference means `A + (-B) = {a - b | a in A, b in B}`, not erosion or +//! Boolean subtraction. It is not the inverse of a Minkowski sum. +//! +//! # Reuse and empty regions +//! +//! ``` +//! use i_triangle::int::minkowski::{IntMinkowskiSolver, MinkowskiError}; +//! use i_triangle::i_overlay::core::fill_rule::FillRule; +//! use i_triangle::i_overlay::i_float::int::point::IntPoint; +//! use i_triangle::i_overlay::i_shape::int_path; +//! +//! let source = int_path![[0, 0], [4, 0], [4, 4], [0, 4]]; +//! let kernel = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; +//! let mut solver = IntMinkowskiSolver::::default(); +//! let a = solver.prepare(&source, FillRule::NonZero).unwrap(); +//! let b = solver.prepare(&kernel, FillRule::NonZero).unwrap(); +//! let sum = solver.sum(&a, &b)?; +//! solver.clear(); +//! let difference = solver.difference(&a, &b)?; +//! assert!(sum[0][0].contains(&(source[2] + kernel[2]))); +//! assert!(difference[0][0].contains(&IntPoint::new(-1, -1))); +//! let line = int_path![[0, 0], [1, 0], [2, 0]]; +//! assert!(solver.prepare(&line, FillRule::NonZero).is_none()); +//! # Ok::<(), MinkowskiError>(()) +//! ``` +//! +//! Only two-dimensional filled regions are represented: points, segments, +//! empty sources, and sources with no filled area produce `None` in preparation. +//! Integer operands own their parts and can be used with any integer solver of +//! the same coordinate type. The floating-point API adds solver ownership to +//! protect its fixed coordinate transform. +//! +//! # Coordinate range and rounding +//! +//! Input coordinates must satisfy the [overlay integer range](i_overlay::core::integer). +//! Preparation does not validate this arithmetic precondition. Operations check +//! conservative result bounds before generating pairs, requiring every result +//! coordinate to be strictly between `-2^(I::BITS - 2)` and `2^(I::BITS - 2)`. +//! Bounds come from the original source, including contours removed by filling. +//! Intersections are rounded to the integer grid, so small features can change +//! during normalization and union. Increase coordinate precision before +//! preparation if one grid unit is too coarse for the application. + +use crate::advanced::delaunay::IntDelaunay; +use crate::int::minkowski_sum::ConvexSum; +use crate::int::triangulator::IntTriangulator; +use alloc::vec::Vec; +use i_overlay::core::fill_rule::FillRule; +use i_overlay::core::integer::OverlayInt; +use i_overlay::core::overlay::{IntOverlayOptions, ShapeType}; +use i_overlay::core::overlay_rule::OverlayRule; +use i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_overlay::i_float::int::rect::IntRect; +use i_overlay::i_shape::int::shape::{IntContour, IntShapes}; +use i_overlay::i_shape::source::int::resource::IntShapeResource; + +/// A numeric, domain, or operand-ownership failure in a Minkowski operation. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum MinkowskiError { + /// A coordinate or conservative result bound exceeds the supported range. + /// Floating-point input must also be finite. + CoordinateOutOfRange, + /// Original source or result bounds extend outside a float solver's domain. + OutsideDomain, + /// A floating-point operand was prepared by a different solver instance. + ForeignOperand, +} + +/// A reusable integer triangulator and union workspace for Minkowski operations. +/// +/// Prepared operands are independent of the workspace and survive clearing it. +/// `I` can be `i16`, `i32`, or `i64`; callers must respect the module's input +/// [coordinate-range contract](self#coordinate-range-and-rounding). +pub struct IntMinkowskiSolver { + triangulator: IntTriangulator, + output_options: IntOverlayOptions, +} + +/// An owned convex decomposition of a normalized, nonempty filled region. +/// +/// Prepare with [`Self::with_source`] or [`IntMinkowskiSolver::prepare`]. The +/// original source may be dropped, and the operand may be reused with any +/// [`IntMinkowskiSolver`] of the same integer type without further triangulation. +pub struct IntMinkowskiOperand { + convex_polygons: Vec>, + bbox: IntRect, +} + +impl IntMinkowskiOperand { + /// Prepares a source using a temporary default solver and `fill_rule`. + /// + /// Returns `None` for an empty or unfilled region, including points and + /// segments. Input must satisfy the module's coordinate-range contract. + /// For repeated preparations, reuse [`IntMinkowskiSolver::prepare`]. + pub fn with_source(source: &R, fill_rule: FillRule) -> Option + where + R: IntShapeResource + ?Sized, + { + IntMinkowskiSolver::default().prepare(source, fill_rule) + } + + /// Borrows the prepared convex parts in counter-clockwise order. + /// + /// These describe the filled region after normalization and integer rounding, + /// rather than the original contour grouping or vertex sequence. + pub fn convex_polygons(&self) -> &[IntContour] { + &self.convex_polygons + } + + /// Returns the original source bounds, before fill filtering. + /// These conservative bounds are used for result-range checks. + pub fn bbox(&self) -> IntRect { + self.bbox + } +} + +impl IntMinkowskiSolver { + /// Creates a solver with an initial capacity hint and simplified collinear output. + /// Capacity is not a limit on source or result size. + pub fn new(capacity: usize) -> Self { + Self::with_options(capacity, false) + } + + /// Creates a solver with an initial capacity hint and a final-output policy. + /// + /// `preserve_output_collinear` retains collinear vertices in the final union. + /// Collinear edges of each convex pair are always merged first, so enabling + /// this option does not preserve every original source vertex. Outer output + /// contours are counter-clockwise; hole contours are clockwise. + pub fn with_options(capacity: usize, preserve_output_collinear: bool) -> Self { + Self { + output_options: IntOverlayOptions { + preserve_output_collinear, + ..Default::default() + }, + triangulator: IntTriangulator::new(capacity, Default::default(), Default::default()), + } + } + + /// Clears temporary overlay geometry while retaining its allocations. + /// Prepared operands and the configured output policy remain usable. + pub fn clear(&mut self) { + self.triangulator.overlay.clear(); + } + + /// Normalizes a source with `fill_rule` and prepares its convex parts once. + /// + /// Accepts contours, shapes with holes, and collections through + /// [`IntShapeResource`]. Returns `None` if no filled area remains, including + /// empty, point, and segment input. Existing operands remain valid. + /// + /// Input coordinates must satisfy the [overlay integer range](i_overlay::core::integer); + /// this method does not perform numeric range validation. + pub fn prepare(&mut self, source: &R, fill_rule: FillRule) -> Option> + where + R: IntShapeResource + ?Sized, + { + self.clear(); + let bbox = IntRect::with_iter(source.iter_paths().flatten())?; + self.triangulator.overlay.options = IntOverlayOptions::keep_output_points(); + let shapes = self.triangulator.overlay.simplify_source(source, fill_rule); + if shapes.is_empty() { + return None; + } + + let mut delaunay = IntDelaunay { + points: Vec::new(), + triangles: Vec::new(), + }; + self.triangulator + .uncheck_delaunay_shapes_into(&shapes, &mut delaunay); + let convex_polygons = delaunay.to_convex_polygons(); + debug_assert!(!convex_polygons.is_empty()); + + Some(IntMinkowskiOperand { + convex_polygons, + bbox, + }) + } + + /// Returns `{a + b | a in lhs, b in rhs}` using the prepared convex parts. + /// + /// # Errors + /// Returns [`MinkowskiError::CoordinateOutOfRange`] when the sum of the + /// original source bounds reaches or exceeds the strict coordinate limits + /// documented in this module. No pairwise points are generated on that error, + /// and the solver remains reusable. + pub fn sum( + &mut self, + lhs: &IntMinkowskiOperand, + rhs: &IntMinkowskiOperand, + ) -> Result, MinkowskiError> { + self.sum_parts( + &lhs.convex_polygons, + lhs.bbox, + &rhs.convex_polygons, + rhs.bbox, + ) + } + + pub(crate) fn sum_parts( + &mut self, + lhs: &R0, + lhs_bbox: IntRect, + rhs: &R1, + rhs_bbox: IntRect, + ) -> Result, MinkowskiError> + where + R0: IntShapeResource + ?Sized, + R1: IntShapeResource + ?Sized, + { + self.clear(); + self.triangulator.overlay.options = self.output_options; + self.triangulator.overlay.options.preserve_input_collinear = true; + + Self::validate_bounds([ + lhs_bbox.min_x.to_wide() + rhs_bbox.min_x.to_wide(), + lhs_bbox.max_x.to_wide() + rhs_bbox.max_x.to_wide(), + lhs_bbox.min_y.to_wide() + rhs_bbox.min_y.to_wide(), + lhs_bbox.max_y.to_wide() + rhs_bbox.max_y.to_wide(), + ])?; + + for a in lhs.iter_paths() { + for b in rhs.iter_paths() { + if let Some(path) = ConvexSum::new(a, b) { + self.triangulator + .overlay + .add_path_iter(path, ShapeType::Subject); + } + } + } + + Ok(self + .triangulator + .overlay + .overlay(OverlayRule::Subject, FillRule::NonZero)) + } + + /// Returns `{a - b | a in lhs, b in rhs}`, reflecting `rhs` through the origin. + /// + /// This is `lhs + (-rhs)`, not erosion, Boolean subtraction, or an inverse of + /// [`Self::sum`]. No repeated triangulation is performed. + /// + /// # Errors + /// Returns [`MinkowskiError::CoordinateOutOfRange`] if the difference of the + /// original source bounds reaches or exceeds the module's strict coordinate + /// limits. The solver remains reusable after an error. + pub fn difference( + &mut self, + lhs: &IntMinkowskiOperand, + rhs: &IntMinkowskiOperand, + ) -> Result, MinkowskiError> { + self.difference_parts( + &lhs.convex_polygons, + lhs.bbox, + &rhs.convex_polygons, + rhs.bbox, + ) + } + + pub(crate) fn difference_parts( + &mut self, + lhs: &R0, + lhs_bbox: IntRect, + rhs: &R1, + rhs_bbox: IntRect, + ) -> Result, MinkowskiError> + where + R0: IntShapeResource + ?Sized, + R1: IntShapeResource + ?Sized, + { + self.clear(); + self.triangulator.overlay.options = self.output_options; + self.triangulator.overlay.options.preserve_input_collinear = true; + + Self::validate_bounds([ + lhs_bbox.min_x.to_wide() - rhs_bbox.max_x.to_wide(), + lhs_bbox.max_x.to_wide() - rhs_bbox.min_x.to_wide(), + lhs_bbox.min_y.to_wide() - rhs_bbox.max_y.to_wide(), + lhs_bbox.max_y.to_wide() - rhs_bbox.min_y.to_wide(), + ])?; + + for a in lhs.iter_paths() { + for b in rhs.iter_paths() { + if let Some(path) = ConvexSum::difference(a, b) { + self.triangulator + .overlay + .add_path_iter(path, ShapeType::Subject); + } + } + } + + Ok(self + .triangulator + .overlay + .overlay(OverlayRule::Subject, FillRule::NonZero)) + } + + pub(crate) fn validate_bounds(bounds: [I::Wide; 4]) -> Result<(), MinkowskiError> { + let limit = I::Wide::ONE << (I::BITS - 2); + if bounds + .iter() + .any(|&value| value <= -limit || value >= limit) + { + Err(MinkowskiError::CoordinateOutOfRange) + } else { + Ok(()) + } + } +} + +impl Default for IntMinkowskiSolver { + fn default() -> Self { + Self::new(64) + } +} + +#[cfg(test)] +mod tests { + use super::IntMinkowskiSolver; + use i_overlay::core::fill_rule::FillRule; + use i_overlay::i_float::int::point::IntPoint; + use i_overlay::i_float::int::rect::IntRect; + use i_overlay::i_shape::int::shape::IntContour; + use i_overlay::i_shape::int_path; + + #[test] + fn repeated_operations_preserve_results() { + let contour: IntContour = (-32..=32).map(|x| IntPoint::new(x, x * x)).collect(); + let mut solver = IntMinkowskiSolver::new(0); + let operand = solver.prepare(&contour, FillRule::NonZero).unwrap(); + assert_eq!(operand.convex_polygons.len(), 1); + let expected = solver.sum(&operand, &operand).unwrap(); + for _ in 0..3 { + let operand = solver.prepare(&contour, FillRule::NonZero).unwrap(); + assert_eq!(solver.sum(&operand, &operand).unwrap(), expected); + } + } + + #[test] + fn degenerate_parts_are_combined_before_discarding_zero_area_results() { + let point = int_path![[2, 3], [2, 3], [2, 3]]; + let horizontal = int_path![[0, 0], [2, 0], [2, 0]]; + let vertical = int_path![[0, 0], [0, 2], [0, 2]]; + let point_bbox = IntRect::with_points(&point).unwrap(); + let h_bbox = IntRect::with_points(&horizontal).unwrap(); + let v_bbox = IntRect::with_points(&vertical).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert!(solver + .sum_parts(&point, point_bbox, &point, point_bbox) + .unwrap() + .is_empty()); + assert!(solver + .difference_parts(&horizontal, h_bbox, &horizontal, h_bbox) + .unwrap() + .is_empty()); + let sum = solver + .sum_parts(&horizontal, h_bbox, &vertical, v_bbox) + .unwrap(); + let difference = solver + .difference_parts(&horizontal, h_bbox, &vertical, v_bbox) + .unwrap(); + let expected_sum = int_path![[0, 0], [2, 0], [2, 2], [0, 2]]; + let expected_difference = int_path![[0, -2], [2, -2], [2, 0], [0, 0]]; + assert_eq!(sum[0][0].len(), 4); + assert!(sum[0][0].iter().all(|p| expected_sum.contains(p))); + assert_eq!(difference[0][0].len(), 4); + assert!(difference[0][0] + .iter() + .all(|p| expected_difference.contains(p))); + } +} diff --git a/iTriangle/src/int/minkowski_sum.rs b/iTriangle/src/int/minkowski_sum.rs new file mode 100644 index 0000000..d7805b1 --- /dev/null +++ b/iTriangle/src/int/minkowski_sum.rs @@ -0,0 +1,251 @@ +use core::cmp::Ordering; +use i_overlay::core::integer::OverlayInt; +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; + +struct Walk<'a, I: OverlayInt> { + points: &'a [IntPoint], + index: usize, + reverse: bool, + remaining: usize, + reflected: bool, + segment: Option<(usize, usize)>, +} + +fn angle_cmp(a: IntVector, b: IntVector) -> Ordering { + let lower = + |v: IntVector| v.y < I::Wide::ZERO || (v.y == I::Wide::ZERO && v.x < I::Wide::ZERO); + lower(a) + .cmp(&lower(b)) + .then_with(|| I::Wide::ZERO.cmp(&a.cross_product(b))) +} + +impl<'a, I: OverlayInt> Walk<'a, I> { + fn new(points: &'a [IntPoint], reflected: bool) -> Option { + let key = |&i: &usize| (points[i].y, points[i].x); + let index = if reflected { + (0..points.len()).max_by_key(key)? + } else { + (0..points.len()).min_by_key(key)? + }; + let mut walk = Self { + points, + index, + reverse: false, + remaining: points.len(), + reflected, + segment: None, + }; + let origin = points[index]; + let mut winding = I::Wide::ZERO; + for _ in 0..points.len() { + let a = points[walk.index] - origin; + let b = points[walk.next_index()] - origin; + winding = winding + a.cross_product(b); + walk.advance(); + } + if winding == I::Wide::ZERO { + let end = if reflected { + (0..points.len()).min_by_key(key)? + } else { + (0..points.len()).max_by_key(key)? + }; + walk.index = index; + walk.segment = Some((index, end)); + walk.remaining = if points[index] == points[end] { 0 } else { 2 }; + return Some(walk); + } + walk.index = index; + walk.remaining = points.len(); + walk.reverse = winding < I::Wide::ZERO; + Some(walk) + } + + #[inline] + fn next_index(&self) -> usize { + if let Some((start, end)) = self.segment { + return if self.index == start { end } else { start }; + } + if self.reverse { + if self.index == 0 { + self.points.len() - 1 + } else { + self.index - 1 + } + } else if self.index + 1 == self.points.len() { + 0 + } else { + self.index + 1 + } + } + + #[inline] + fn edge(&self) -> IntVector { + let a = self.points[self.index]; + let b = self.points[self.next_index()]; + if self.reflected { + a - b + } else { + b - a + } + } + + #[inline] + fn point(&self) -> IntPoint { + let p = self.points[self.index]; + if self.reflected { + IntPoint::new(-p.x, -p.y) + } else { + p + } + } + + #[inline] + fn advance(&mut self) { + self.index = self.next_index(); + self.remaining -= 1; + } +} + +pub(super) struct ConvexSum<'a, I: OverlayInt> { + a: Walk<'a, I>, + b: Walk<'a, I>, +} + +impl<'a, I: OverlayInt> ConvexSum<'a, I> { + pub(super) fn new(a: &'a [IntPoint], b: &'a [IntPoint]) -> Option { + Some(Self { + a: Walk::new(a, false)?, + b: Walk::new(b, false)?, + }) + } + + pub(super) fn difference(a: &'a [IntPoint], b: &'a [IntPoint]) -> Option { + Some(Self { + a: Walk::new(a, false)?, + b: Walk::new(b, true)?, + }) + } + + fn next_step(&self) -> Option<(bool, IntVector)> { + if self.a.remaining == 0 && self.b.remaining == 0 { + return None; + } + if self.a.remaining == 0 { + return Some((false, self.b.edge())); + } + if self.b.remaining == 0 { + return Some((true, self.a.edge())); + } + let a = self.a.edge(); + let b = self.b.edge(); + let zero = IntVector::::new(I::Wide::ZERO, I::Wide::ZERO); + if a == zero || (b != zero && angle_cmp(a, b) != Ordering::Greater) { + Some((true, a)) + } else { + Some((false, b)) + } + } +} + +impl Iterator for ConvexSum<'_, I> { + type Item = IntPoint; + + fn next(&mut self) -> Option { + let point = self.a.point() + self.b.point(); + let zero = IntVector::::new(I::Wide::ZERO, I::Wide::ZERO); + let mut direction = None; + while let Some((advance_a, edge)) = self.next_step() { + if edge != zero { + if let Some(previous) = direction { + if angle_cmp(previous, edge) != Ordering::Equal { + break; + } + } else { + direction = Some(edge); + } + } + if advance_a { + self.a.advance(); + } else { + self.b.advance(); + } + } + direction.map(|_| point) + } +} + +#[cfg(test)] +mod tests { + use super::ConvexSum; + use alloc::vec::Vec; + use i_overlay::i_float::int::point::IntPoint; + use i_overlay::i_shape::int_path; + + #[test] + fn collinear_edges_and_repeated_vertices_are_merged() { + let a = int_path![[0, 0], [2, 0], [2, 0], [4, 0], [4, 4], [0, 4]]; + let b = int_path![[0, 0], [2, 0], [2, 2], [0, 2]]; + let sum: Vec<_> = ConvexSum::new(&a, &b).unwrap().collect(); + assert_eq!(sum, int_path![[0, 0], [6, 0], [6, 6], [0, 6]]); + } + + #[test] + fn merges_closing_edges_for_all_starts_and_windings() { + let contour = int_path![ + [0, 0], + [2, 0], + [4, 0], + [4, 2], + [4, 4], + [2, 4], + [0, 4], + [0, 2], + [0, 0] + ]; + let expected = int_path![[0, 0], [8, 0], [8, 8], [0, 8]]; + for reverse_a in [false, true] { + for reverse_b in [false, true] { + for start_a in 0..contour.len() { + for start_b in 0..contour.len() { + let mut a = contour.clone(); + let mut b = contour.clone(); + if reverse_a { + a.reverse(); + } + if reverse_b { + b.reverse(); + } + a.rotate_left(start_a); + b.rotate_left(start_b); + let mut sum = ConvexSum::new(&a, &b).unwrap(); + assert_eq!(sum.by_ref().collect::>(), expected); + assert_eq!(sum.next(), None); + assert_eq!(sum.next(), None); + } + } + } + } + } + + #[test] + fn slightly_nonconvex_parts_are_accepted() { + let a = int_path![[0, 0], [3000, 0], [2999, 1], [3000, 3000], [0, 3000]]; + let point = int_path![[2, 3], [2, 3], [2, 3]]; + let sum: Vec<_> = ConvexSum::new(&a, &point).unwrap().collect(); + let difference: Vec<_> = ConvexSum::difference(&a, &point).unwrap().collect(); + assert_eq!( + sum, + a.iter() + .map(|p| IntPoint::new(p.x + 2, p.y + 3)) + .collect::>() + ); + assert_eq!( + difference, + a.iter() + .map(|p| IntPoint::new(p.x - 2, p.y - 3)) + .collect::>() + ); + } +} diff --git a/iTriangle/src/int/mod.rs b/iTriangle/src/int/mod.rs index 438445f..bb8aca0 100644 --- a/iTriangle/src/int/mod.rs +++ b/iTriangle/src/int/mod.rs @@ -4,6 +4,8 @@ pub mod custom; pub mod earcut; pub mod locator; mod meta; +pub mod minkowski; +mod minkowski_sum; pub(crate) mod monotone; mod solver; pub mod triangulatable; diff --git a/iTriangle/tests/float_minkowski_tests.rs b/iTriangle/tests/float_minkowski_tests.rs new file mode 100644 index 0000000..de2bb70 --- /dev/null +++ b/iTriangle/tests/float_minkowski_tests.rs @@ -0,0 +1,524 @@ +use i_triangle::float::minkowski::{ + minkowski_difference, minkowski_sum, MinkowskiError, MinkowskiSolver, +}; +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::integer::OverlayInt; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::float::overlay::FloatOverlay; +use i_triangle::i_overlay::i_float::float::compatible::FloatPointCompatible; +use i_triangle::i_overlay::i_float::float::number::FloatNumber; +use i_triangle::i_overlay::i_float::float::rect::FloatRect; +use i_triangle::i_overlay::i_shape::base::data::Shapes; + +fn domain(min: f64, max: f64) -> FloatRect { + FloatRect::new( + F::from_float(min), + F::from_float(max), + F::from_float(min), + F::from_float(max), + ) + .unwrap() +} + +fn rectangle(min_x: f64, max_x: f64, min_y: f64, max_y: f64) -> Vec<[F; 2]> { + [ + [min_x, min_y], + [max_x, min_y], + [max_x, max_y], + [min_x, max_y], + ] + .map(|[x, y]| [F::from_float(x), F::from_float(y)]) + .into() +} + +fn assert_rectangle(result: &Shapes<[F; 2]>, expected: &[[F; 2]], tolerance: f64) { + assert_eq!(result.len(), 1); + assert_eq!(result[0].len(), 1); + assert_eq!(result[0][0].len(), 4); + for p in &result[0][0] { + assert!( + expected + .iter() + .any(|q| (p[0].to_f64() - q[0].to_f64()).abs() <= tolerance + && (p[1].to_f64() - q[1].to_f64()).abs() <= tolerance), + "unexpected point ({}, {})", + p[0], + p[1] + ); + } +} + +#[test] +fn fixed_grid_operands_work_for_all_coordinate_types() { + fn check() { + let a = rectangle::(-10.5, -4.5, 3.25, 8.25); + let b = rectangle::(2.0, 6.0, -5.0, -2.0); + let mut solver = MinkowskiSolver::<[F; 2], I>::new(domain(-32.0, 32.0)); + let lhs = solver + .prepare(a.as_slice(), FillRule::NonZero) + .unwrap() + .unwrap(); + let rhs = solver.prepare(&b, FillRule::NonZero).unwrap().unwrap(); + for _ in 0..3 { + assert_rectangle( + &solver.sum(&lhs, &rhs).unwrap(), + &rectangle::(-8.5, 1.5, -1.75, 6.25), + 1e-6, + ); + assert_rectangle( + &solver.difference(&lhs, &rhs).unwrap(), + &rectangle::(-16.5, -6.5, 5.25, 13.25), + 1e-6, + ); + assert_rectangle( + &solver.difference(&rhs, &lhs).unwrap(), + &rectangle::(6.5, 16.5, -13.25, -5.25), + 1e-6, + ); + } + } + check::(); + check::(); + check::(); + check::(); + check::(); + check::(); +} + +#[test] +fn fractional_rectangles_use_the_fixed_coordinate_system() { + fn check() { + let mut solver = MinkowskiSolver::<[F; 2], I>::new(domain(-32.0, 32.0)); + for i in 0..64 { + let x = i as f64 * 0.137 - 4.1; + let y = i as f64 * -0.083 + 2.3; + let a = rectangle::(x, x + 1.31, y, y + 0.79); + let b = rectangle::(-0.23, 0.48, 0.17, 0.88); + let lhs = solver.prepare(&a, FillRule::NonZero).unwrap().unwrap(); + let rhs = solver.prepare(&b, FillRule::NonZero).unwrap().unwrap(); + assert_rectangle( + &solver.sum(&lhs, &rhs).unwrap_or_else(|e| { + panic!( + "{e:?}, F={}, I={}, i={i}", + core::any::type_name::(), + core::any::type_name::() + ) + }), + &rectangle::(x - 0.23, x + 1.79, y + 0.17, y + 1.67), + 0.01, + ); + assert_rectangle( + &solver.difference(&lhs, &rhs).unwrap(), + &rectangle::(x - 0.48, x + 1.54, y - 0.88, y + 0.62), + 0.01, + ); + } + } + check::(); + check::(); + check::(); + check::(); + check::(); + check::(); +} + +#[test] +fn holes_concavities_and_collinear_options_work_on_fixed_grid() { + let concave = vec![ + [0.0, 0.0], + [4.0, 0.0], + [4.0, 1.0], + [1.0, 1.0], + [1.0, 4.0], + [0.0, 4.0], + ]; + let source = vec![ + rectangle::(0.0, 10.0, 0.0, 10.0), + vec![[3.0, 3.0], [3.0, 7.0], [7.0, 7.0], [7.0, 3.0]], + ]; + let kernel = rectangle::(0.0, 1.0, 0.0, 1.0); + let expected = vec![ + rectangle::(0.0, 11.0, 0.0, 11.0), + vec![[4.0, 4.0], [4.0, 7.0], [7.0, 7.0], [7.0, 4.0]], + ]; + let expected_concave = vec![ + [-1.0, -1.0], + [4.0, -1.0], + [4.0, 1.0], + [1.0, 1.0], + [1.0, 4.0], + [-1.0, 4.0], + ]; + for keep_collinear in [false, true] { + let mut solver = + MinkowskiSolver::<[f64; 2]>::with_options(domain(-16.0, 16.0), 0, keep_collinear); + let rhs = solver.prepare(&kernel, FillRule::NonZero).unwrap().unwrap(); + for _ in 0..3 { + let lhs = solver.prepare(&source, FillRule::NonZero).unwrap().unwrap(); + let result = solver.sum(&lhs, &rhs).unwrap(); + assert!( + FloatOverlay::<[f64; 2]>::from_subj_and_clip(&result, &expected) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty() + ); + let lhs = solver + .prepare(&concave, FillRule::NonZero) + .unwrap() + .unwrap(); + let result = solver.difference(&lhs, &rhs).unwrap(); + assert!( + FloatOverlay::<[f64; 2]>::from_subj_and_clip(&result, &expected_concave) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty() + ); + assert_eq!(result[0][0].len() > expected_concave.len(), keep_collinear); + } + } +} + +#[test] +fn empty_and_unfilled_sources_have_no_operand() { + let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain(-32.0, 32.0)); + assert!(solver + .prepare(&Vec::<[f64; 2]>::new(), FillRule::NonZero) + .unwrap() + .is_none()); + assert!(solver + .prepare(&[[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]], FillRule::NonZero) + .unwrap() + .is_none()); + let square = rectangle::(0.0, 1.0, 0.0, 1.0); + assert!(solver + .prepare(&square, FillRule::Negative) + .unwrap() + .is_none()); + assert!(solver + .prepare(&[square.clone(), square.clone()], FillRule::EvenOdd) + .unwrap() + .is_none()); + assert!(solver + .prepare(&square, FillRule::NonZero) + .unwrap() + .is_some()); +} + +#[test] +fn out_of_range_result_is_reported_without_poisoning_the_solver() { + fn check() { + let limit = F::MAX_COORDINATE.to_f64(); + let a = rectangle::(limit * 0.75, limit, 0.0, limit * 0.25); + let mut solver = MinkowskiSolver::<[F; 2]>::new(domain(-limit, limit)); + let operand = solver.prepare(&a, FillRule::NonZero).unwrap().unwrap(); + assert!(matches!( + solver.sum(&operand, &operand), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert_rectangle( + &solver.difference(&operand, &operand).unwrap(), + &rectangle::(-limit * 0.25, limit * 0.25, -limit * 0.25, limit * 0.25), + limit * 1e-6, + ); + } + check::(); + check::(); +} + +#[derive(Clone, Copy)] +struct Point { + x: f64, + y: f64, +} + +impl FloatPointCompatible for Point { + type Scalar = f64; + fn from_xy(x: f64, y: f64) -> Self { + Self { x, y } + } + fn x(&self) -> f64 { + self.x + } + fn y(&self) -> f64 { + self.y + } +} + +#[test] +fn supports_custom_float_points() { + let points = rectangle::(0.0, 2.0, 0.0, 2.0) + .into_iter() + .map(|[x, y]| Point { x, y }) + .collect::>(); + let mut solver = MinkowskiSolver::::new(domain(-32.0, 32.0)); + let operand = solver.prepare(&points, FillRule::NonZero).unwrap().unwrap(); + let result = solver.sum(&operand, &operand).unwrap(); + assert_eq!(result[0][0].len(), 4); + assert!(result[0][0] + .iter() + .all(|p| (p.x == 0.0 || p.x == 4.0) && (p.y == 0.0 || p.y == 4.0))); +} + +#[test] +fn nonzero_grid_origin_is_accounted_for_in_both_operations() { + let mut solver = MinkowskiSolver::<[f64; 2], i16>::new(domain(-4.0, 24.0)); + let a = solver + .prepare(&rectangle::(8.0, 10.0, 8.0, 10.0), FillRule::NonZero) + .unwrap() + .unwrap(); + let b = solver + .prepare(&rectangle::(10.0, 12.0, 10.0, 12.0), FillRule::NonZero) + .unwrap() + .unwrap(); + let original = a.int_operand().convex_polygons().to_vec(); + for _ in 0..3 { + assert_rectangle( + &solver.sum(&a, &b).unwrap(), + &rectangle::(18.0, 22.0, 18.0, 22.0), + 1e-8, + ); + assert_rectangle( + &solver.difference(&a, &b).unwrap(), + &rectangle::(-4.0, 0.0, -4.0, 0.0), + 1e-8, + ); + assert_eq!(original, a.int_operand().convex_polygons()); + solver.clear(); + } +} + +#[test] +fn rejects_source_and_result_bounds_outside_the_domain() { + let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain(-8.0, 8.0)); + for source in [ + rectangle::(-9.0, 0.0, 0.0, 1.0), + rectangle::(0.0, 9.0, 0.0, 1.0), + rectangle::(0.0, 1.0, -9.0, 0.0), + rectangle::(0.0, 1.0, 0.0, 9.0), + ] { + assert!(matches!( + solver.prepare(&source, FillRule::NonZero), + Err(MinkowskiError::OutsideDomain) + )); + } + let a = solver + .prepare(&rectangle::(4.0, 8.0, 4.0, 8.0), FillRule::NonZero) + .unwrap() + .unwrap(); + let b = solver + .prepare(&rectangle::(-8.0, -4.0, -8.0, -4.0), FillRule::NonZero) + .unwrap() + .unwrap(); + assert!(matches!( + solver.sum(&a, &a), + Err(MinkowskiError::OutsideDomain) + )); + assert!(matches!( + solver.difference(&a, &b), + Err(MinkowskiError::OutsideDomain) + )); + assert_rectangle( + &solver.sum(&a, &b).unwrap(), + &rectangle::(-4.0, 4.0, -4.0, 4.0), + 1e-8, + ); +} + +#[test] +fn foreign_operands_are_rejected_even_with_identical_grids() { + let source = rectangle::(0.0, 1.0, 0.0, 1.0); + let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain(-16.0, 16.0)); + let a = solver.prepare(&source, FillRule::NonZero).unwrap().unwrap(); + for rect in [ + domain(-16.0, 16.0), + domain(-15.0, 15.0), + domain(-8.0, 8.0), + domain(-15.0, 17.0), + ] { + let mut foreign = MinkowskiSolver::<[f64; 2]>::new(rect); + let b = foreign + .prepare(&source, FillRule::NonZero) + .unwrap() + .unwrap(); + for (lhs, rhs) in [(&a, &b), (&b, &a)] { + assert!(matches!( + solver.sum(lhs, rhs), + Err(MinkowskiError::ForeignOperand) + )); + assert!(matches!( + solver.difference(lhs, rhs), + Err(MinkowskiError::ForeignOperand) + )); + } + } + let mut foreign = MinkowskiSolver::<[f64; 2]>::new(domain(100.0, 120.0)); + let b = foreign + .prepare( + &rectangle::(101.0, 102.0, 101.0, 102.0), + FillRule::NonZero, + ) + .unwrap() + .unwrap(); + assert!(matches!( + solver.sum(&a, &b), + Err(MinkowskiError::ForeignOperand) + )); + assert_rectangle( + &solver.sum(&a, &a).unwrap(), + &rectangle::(0.0, 2.0, 0.0, 2.0), + 1e-8, + ); +} + +#[test] +fn operand_owner_survives_move_clear_and_further_preparation() { + let rect = domain(-8.0, 8.0); + let source = rectangle::(0.0, 1.0, 0.0, 1.0); + let mut solver = MinkowskiSolver::<[f64; 2]>::new(rect); + let a = solver.prepare(&source, FillRule::NonZero).unwrap().unwrap(); + let mut moved = core::mem::replace(&mut solver, MinkowskiSolver::new(rect)); + moved.clear(); + let b = moved.prepare(&source, FillRule::NonZero).unwrap().unwrap(); + moved.clear(); + assert_rectangle( + &moved.sum(&a, &b).unwrap(), + &rectangle::(0.0, 2.0, 0.0, 2.0), + 1e-8, + ); + assert_rectangle( + &moved.difference(&a, &b).unwrap(), + &rectangle::(-1.0, 1.0, -1.0, 1.0), + 1e-8, + ); + assert!(matches!( + solver.sum(&a, &b), + Err(MinkowskiError::ForeignOperand) + )); + drop(moved); + let mut replacement = MinkowskiSolver::<[f64; 2]>::new(rect); + assert!(matches!( + replacement.difference(&a, &b), + Err(MinkowskiError::ForeignOperand) + )); +} + +#[test] +fn sugar_selects_domain_for_sources_and_each_operation() { + for (a, b) in [ + ( + rectangle::(8.0, 10.0, 3.0, 5.0), + rectangle::(10.0, 12.0, 6.0, 8.0), + ), + ( + rectangle::(1e6, 1e6 + 4.0, 1e6, 1e6 + 4.0), + rectangle::(-1e6, -1e6 + 0.25, -1e6, -1e6 + 0.25), + ), + ] { + assert_rectangle( + &minkowski_sum(a.as_slice(), &b, FillRule::NonZero).unwrap(), + &rectangle::( + a[0][0] + b[0][0], + a[2][0] + b[2][0], + a[0][1] + b[0][1], + a[2][1] + b[2][1], + ), + 1e-8, + ); + assert_rectangle( + &minkowski_difference(&a, &b, FillRule::NonZero).unwrap(), + &rectangle::( + a[0][0] - b[2][0], + a[2][0] - b[0][0], + a[0][1] - b[2][1], + a[2][1] - b[0][1], + ), + 1e-8, + ); + } + let a = rectangle::(-3.0, 1.0, 2.0, 4.0); + assert_rectangle( + &minkowski_sum(&a, &a, FillRule::NonZero).unwrap(), + &rectangle::(-6.0, 2.0, 4.0, 8.0), + 1e-6, + ); +} + +#[test] +fn sugar_returns_empty_for_empty_unfilled_and_collapsed_sources() { + let square = rectangle::(0.0, 1.0, 0.0, 1.0); + let empty: Vec<[f64; 2]> = Vec::new(); + let line = vec![[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]; + for source in [&empty, &line] { + assert!(minkowski_sum(source, &square, FillRule::NonZero) + .unwrap() + .is_empty()); + assert!(minkowski_difference(&square, source, FillRule::NonZero) + .unwrap() + .is_empty()); + } + assert!(minkowski_sum(&square, &square, FillRule::Negative) + .unwrap() + .is_empty()); + assert!(minkowski_difference( + &[square.clone(), square.clone()], + &square, + FillRule::EvenOdd + ) + .unwrap() + .is_empty()); + let tiny = rectangle::(0.0, 1e-12, 0.0, 1e-12); + assert!(minkowski_sum(&tiny, &square, FillRule::NonZero) + .unwrap() + .is_empty()); +} + +#[test] +fn sugar_reports_numeric_errors_and_prepare_rejects_invalid_points() { + let limit = f64::MAX_COORDINATE; + let large = rectangle::(limit * 0.75, limit, 0.0, limit * 0.25); + assert!(matches!( + minkowski_sum(&large, &large, FillRule::NonZero), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert!(minkowski_sum(&large, &large, FillRule::Negative) + .unwrap() + .is_empty()); + let negative = rectangle::(-limit, -limit * 0.75, 0.0, limit * 0.25); + assert!(matches!( + minkowski_difference(&large, &negative, FillRule::NonZero), + Err(MinkowskiError::CoordinateOutOfRange) + )); + let mut solver = MinkowskiSolver::<[f64; 2]>::new(domain(-8.0, 8.0)); + for x in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY, limit * 2.0] { + let invalid = [[0.0, 0.0], [x, 0.0], [0.0, 1.0]]; + assert!(matches!( + solver.prepare(&invalid, FillRule::NonZero), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert!(matches!( + minkowski_sum(&invalid, &invalid, FillRule::NonZero), + Err(MinkowskiError::CoordinateOutOfRange) + )); + } +} + +#[test] +fn exact_grid_extremes_do_not_overflow_integer_sum() { + let mut solver = MinkowskiSolver::<[f64; 2], i16>::new(domain(-1.0, 1.0)); + let a = solver + .prepare(&rectangle::(0.0, 0.5, 0.0, 0.5), FillRule::NonZero) + .unwrap() + .unwrap(); + assert_rectangle( + &solver.sum(&a, &a).unwrap(), + &rectangle::(0.0, 1.0, 0.0, 1.0), + 0.0, + ); + let a = solver + .prepare(&rectangle::(-1.0, 0.0, -1.0, 0.0), FillRule::NonZero) + .unwrap() + .unwrap(); + assert_rectangle( + &solver.difference(&a, &a).unwrap(), + &rectangle::(-1.0, 1.0, -1.0, 1.0), + 0.0, + ); +} diff --git a/iTriangle/tests/minkowski_difference_tests.rs b/iTriangle/tests/minkowski_difference_tests.rs new file mode 100644 index 0000000..8fbcb4c --- /dev/null +++ b/iTriangle/tests/minkowski_difference_tests.rs @@ -0,0 +1,200 @@ +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::integer::OverlayInt; +use i_triangle::i_overlay::core::overlay::Overlay; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_triangle::i_overlay::i_float::int::point::IntPoint; +use i_triangle::i_overlay::i_shape::int::shape::{IntContour, IntShapes}; +use i_triangle::i_overlay::i_shape::source::int::resource::IntShapeResource; +use i_triangle::i_overlay::i_shape::{int_path, int_shapes}; +use i_triangle::int::minkowski::{IntMinkowskiSolver, MinkowskiError}; + +fn assert_geometry + ?Sized>( + actual: &IntShapes, + expected: &R, +) { + assert!(Overlay::from_subj_and_clip(actual, expected) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty()); +} + +#[test] +fn subtracts_translated_rectangles_for_all_integer_types() { + fn check>() { + let contour = |points: [(i16, i16); 4]| -> IntContour { + points + .map(|(x, y)| IntPoint::new(I::from(x), I::from(y))) + .into() + }; + let a = contour([(-10, 3), (-4, 3), (-4, 8), (-10, 8)]); + let b = contour([(2, -5), (6, -5), (6, -2), (2, -2)]); + let expected = contour([(-16, 5), (-6, 5), (-6, 13), (-16, 13)]); + let reversed = contour([(6, -13), (16, -13), (16, -5), (6, -5)]); + let mut solver = IntMinkowskiSolver::default(); + let lhs = solver.prepare(&a, FillRule::NonZero).unwrap(); + let rhs = solver.prepare(&b, FillRule::NonZero).unwrap(); + assert_geometry(&solver.difference(&lhs, &rhs).unwrap(), &expected); + assert_geometry(&solver.difference(&rhs, &lhs).unwrap(), &reversed); + } + + check::(); + check::(); + check::(); +} + +#[test] +fn triangle_difference_matches_sum_with_reflection() { + let a = int_path![[0, 0], [4, 0], [0, 4]]; + let b = int_path![[0, 0], [2, 0], [0, 2]]; + let reflected = int_path![[0, 0], [-2, 0], [0, -2]]; + let expected = int_path![[-2, 0], [0, -2], [4, -2], [4, 0], [0, 4], [-2, 4]]; + let mut solver = IntMinkowskiSolver::default(); + let lhs = solver.prepare(&a, FillRule::NonZero).unwrap(); + let rhs = solver.prepare(&b, FillRule::NonZero).unwrap(); + let reflected = solver.prepare(&reflected, FillRule::NonZero).unwrap(); + let difference = solver.difference(&lhs, &rhs).unwrap(); + assert_geometry(&difference, &expected); + assert_geometry(&difference, &solver.sum(&lhs, &reflected).unwrap()); +} + +#[test] +fn difference_is_not_the_inverse_of_sum() { + let a = int_path![[0, 0], [10, 0], [10, 10], [0, 10]]; + let b = int_path![[0, 0], [2, 0], [2, 2], [0, 2]]; + let mut solver = IntMinkowskiSolver::default(); + let lhs = solver.prepare(&a, FillRule::NonZero).unwrap(); + let rhs = solver.prepare(&b, FillRule::NonZero).unwrap(); + let sum = solver.sum(&lhs, &rhs).unwrap(); + assert_geometry(&sum, &int_path![[0, 0], [12, 0], [12, 12], [0, 12]]); + let sum = solver.prepare(&sum, FillRule::NonZero).unwrap(); + let difference = solver.difference(&sum, &rhs).unwrap(); + assert_geometry( + &difference, + &int_path![[-2, -2], [12, -2], [12, 12], [-2, 12]], + ); + assert!(!Overlay::from_subj_and_clip(&difference, &a) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty()); +} + +#[test] +fn preserves_concavity_holes_and_output_options_across_operations() { + let concave = int_path![[0, 0], [4, 0], [4, 1], [1, 1], [1, 4], [0, 4]]; + let holed = int_shapes![[ + [[0, 0], [10, 0], [10, 10], [0, 10]], + [[3, 3], [3, 7], [7, 7], [7, 3]] + ]]; + let square = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let expected_concave = int_path![[-1, -1], [4, -1], [4, 1], [1, 1], [1, 4], [-1, 4]]; + let expected_holed = int_shapes![[ + [[-1, -1], [10, -1], [10, 10], [-1, 10]], + [[3, 3], [3, 6], [6, 6], [6, 3]] + ]]; + for preserve_output_collinear in [false, true] { + let mut solver = IntMinkowskiSolver::with_options(0, preserve_output_collinear); + let rhs = solver.prepare(&square, FillRule::NonZero).unwrap(); + for _ in 0..3 { + let lhs = solver.prepare(&concave, FillRule::NonZero).unwrap(); + let result = solver.difference(&lhs, &rhs).unwrap(); + assert_geometry(&result, &expected_concave); + if preserve_output_collinear { + assert!(result[0][0].len() > expected_concave.len()); + } else { + assert_eq!(result[0][0].len(), expected_concave.len()); + } + let lhs = solver.prepare(&holed, FillRule::NonZero).unwrap(); + assert_geometry(&solver.difference(&lhs, &rhs).unwrap(), &expected_holed); + assert_geometry( + &solver.sum(&rhs, &rhs).unwrap(), + &int_path![[0, 0], [2, 0], [2, 2], [0, 2]], + ); + } + } +} + +#[test] +fn subtraction_checks_all_bbox_bounds_and_recovers_after_errors() { + fn check() { + let limit = I::Wide::ONE << (I::BITS - 2); + let one = I::Wide::ONE; + let two = I::Wide::TWO; + let four = I::Wide::FOUR; + for swap_axes in [false, true] { + let rectangle = |min: I::Wide, max: I::Wide| -> IntContour { + [ + (min, I::Wide::ZERO), + (max, I::Wide::ZERO), + (max, two), + (min, two), + ] + .map(|(x, y)| { + let (x, y) = if swap_axes { (y, x) } else { (x, y) }; + IntPoint::new(I::from_wide(x), I::from_wide(y)) + }) + .into() + }; + let mut solver = IntMinkowskiSolver::default(); + for positive in [false, true] { + let a = if positive { + rectangle(limit - four, limit - two) + } else { + rectangle(-limit + two, -limit + four) + }; + let valid = if positive { + rectangle(-one, I::Wide::ZERO) + } else { + rectangle(I::Wide::ZERO, one) + }; + let invalid = if positive { + rectangle(-two, -one) + } else { + rectangle(one, two) + }; + let lhs = solver.prepare(&a, FillRule::NonZero).unwrap(); + let valid = solver.prepare(&valid, FillRule::NonZero).unwrap(); + let invalid = solver.prepare(&invalid, FillRule::NonZero).unwrap(); + let before = solver.difference(&lhs, &valid).unwrap(); + assert!(matches!( + solver.difference(&lhs, &invalid), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert_geometry(&solver.difference(&lhs, &valid).unwrap(), &before); + let self_difference = solver.difference(&lhs, &lhs).unwrap(); + let expected: IntContour = [(-two, -two), (two, -two), (two, two), (-two, two)] + .map(|(x, y)| IntPoint::new(I::from_wide(x), I::from_wide(y))) + .into(); + assert_geometry(&self_difference, &expected); + } + } + } + + check::(); + check::(); + check::(); +} + +#[test] +fn rejects_original_bbox_even_when_retained_geometry_would_fit() { + let limit = 1_i32 << 30; + let square = int_path![[0, 0], [2, 0], [2, 2], [0, 2]]; + let remote = int_path![ + [limit - 4, 0], + [limit - 2, 0], + [limit - 2, 2], + [limit - 4, 2] + ]; + let source = vec![square.clone(), remote.clone(), remote]; + let kernel = int_path![[-5, -2], [-3, -2], [-3, 0], [-5, 0]]; + let mut solver = IntMinkowskiSolver::default(); + let lhs = solver.prepare(&source, FillRule::EvenOdd).unwrap(); + let rhs = solver.prepare(&kernel, FillRule::NonZero).unwrap(); + assert!(matches!( + solver.difference(&lhs, &rhs), + Err(MinkowskiError::CoordinateOutOfRange) + )); + let retained = solver.prepare(&square, FillRule::NonZero).unwrap(); + assert_geometry( + &solver.difference(&retained, &rhs).unwrap(), + &int_path![[3, 0], [7, 0], [7, 4], [3, 4]], + ); +} diff --git a/iTriangle/tests/minkowski_operand_tests.rs b/iTriangle/tests/minkowski_operand_tests.rs new file mode 100644 index 0000000..e6d98c0 --- /dev/null +++ b/iTriangle/tests/minkowski_operand_tests.rs @@ -0,0 +1,90 @@ +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::integer::OverlayInt; +use i_triangle::i_overlay::core::overlay::Overlay; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::i_float::int::point::IntPoint; +use i_triangle::i_overlay::i_float::int::rect::IntRect; +use i_triangle::i_overlay::i_shape::int::path::ContourExtension; +use i_triangle::i_overlay::i_shape::int::shape::IntContour; +use i_triangle::i_overlay::i_shape::{int_path, int_shapes}; +use i_triangle::int::minkowski::IntMinkowskiOperand; + +#[test] +fn prepares_borrowed_contours_for_all_integer_types() { + fn check>() { + let contour: IntContour = [(0, 0), (10, 0), (10, 10), (0, 10)] + .map(|(x, y)| IntPoint::new(I::from(x), I::from(y))) + .into(); + let operand = + IntMinkowskiOperand::with_source(contour.as_slice(), FillRule::NonZero).unwrap(); + + assert!(operand.bbox() == IntRect::new(I::from(0), I::from(10), I::from(0), I::from(10))); + assert_eq!(operand.convex_polygons().len(), 1); + let polygon = &operand.convex_polygons()[0]; + assert_eq!(polygon.len(), contour.len()); + assert!(polygon.iter().all(|p| contour.contains(p))); + assert!(polygon.is_convex()); + assert!(!polygon.is_clockwise_ordered()); + } + + check::(); + check::(); + check::(); +} + +#[test] +fn empty_and_degenerate_sources_have_no_operand() { + let contours: Vec> = vec![ + vec![], + int_path![[0, 0]], + int_path![[0, 0], [10, 0]], + int_path![[0, 0], [5, 0], [10, 0]], + ]; + for contour in &contours { + assert!(IntMinkowskiOperand::with_source(contour, FillRule::NonZero).is_none()); + } + assert!(IntMinkowskiOperand::with_source(&contours, FillRule::NonZero).is_none()); +} + +#[test] +fn decomposition_preserves_concavity_holes_and_separate_components() { + let shapes = int_shapes![ + [ + [[0, 0], [20, 0], [20, 10], [10, 10], [10, 20], [0, 20]], + [[2, 2], [2, 6], [6, 6], [6, 2]] + ], + [[[30, 0], [40, 0], [40, 10], [30, 10]]] + ]; + let operand = IntMinkowskiOperand::with_source(&shapes, FillRule::NonZero).unwrap(); + assert_eq!(operand.bbox(), IntRect::new(0, 40, 0, 20)); + for polygon in operand.convex_polygons() { + assert!(polygon.len() >= 3); + assert!(polygon.is_convex()); + assert!(!polygon.is_clockwise_ordered()); + } + + let difference = Overlay::from_subj_and_clip(&shapes, operand.convex_polygons()) + .overlay(OverlayRule::Xor, FillRule::NonZero); + assert!(difference.is_empty()); +} + +#[test] +fn fill_rule_does_not_shrink_the_source_bbox() { + let square = int_path![[0, 0], [10, 0], [10, 10], [0, 10]]; + let remote = int_path![[100, 100], [110, 100], [110, 110], [100, 110]]; + let paths = vec![square.clone(), remote.clone(), remote.clone()]; + let operand = IntMinkowskiOperand::with_source(&paths, FillRule::EvenOdd).unwrap(); + + assert_eq!(operand.bbox(), IntRect::new(0, 110, 0, 110)); + assert!( + Overlay::from_subj_and_clip(&square, operand.convex_polygons()) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty() + ); + + assert!( + IntMinkowskiOperand::with_source(&[remote.clone(), remote], FillRule::EvenOdd).is_none() + ); + assert!(IntMinkowskiOperand::with_source(&square, FillRule::Positive).is_some()); + assert!(IntMinkowskiOperand::with_source(&square, FillRule::Negative).is_none()); +} diff --git a/iTriangle/tests/minkowski_sum_tests.rs b/iTriangle/tests/minkowski_sum_tests.rs new file mode 100644 index 0000000..27f9986 --- /dev/null +++ b/iTriangle/tests/minkowski_sum_tests.rs @@ -0,0 +1,389 @@ +use i_triangle::i_overlay::core::fill_rule::FillRule; +use i_triangle::i_overlay::core::integer::OverlayInt; +use i_triangle::i_overlay::core::overlay::Overlay; +use i_triangle::i_overlay::core::overlay_rule::OverlayRule; +use i_triangle::i_overlay::i_float::int::number::wide_int::WideIntNumber; +use i_triangle::i_overlay::i_float::int::point::IntPoint; +use i_triangle::i_overlay::i_shape::int::path::ContourExtension; +use i_triangle::i_overlay::i_shape::int::shape::{IntContour, IntShapes}; +use i_triangle::i_overlay::i_shape::source::int::resource::IntShapeResource; +use i_triangle::i_overlay::i_shape::{int_path, int_shapes}; +use i_triangle::int::minkowski::{IntMinkowskiOperand, IntMinkowskiSolver, MinkowskiError}; +use rand::rngs::StdRng; +use rand::{RngExt, SeedableRng}; +use std::collections::BTreeSet; + +fn assert_geometry + ?Sized>( + actual: &IntShapes, + expected: &R, +) { + assert!(Overlay::from_subj_and_clip(actual, expected) + .overlay(OverlayRule::Xor, FillRule::NonZero) + .is_empty()); +} + +#[test] +fn adds_translated_rectangles_for_all_integer_types() { + fn check>() { + let contour = |points: [(i16, i16); 4]| -> IntContour { + points + .map(|(x, y)| IntPoint::new(I::from(x), I::from(y))) + .into() + }; + let a = contour([(-10, 3), (-4, 3), (-4, 8), (-10, 8)]); + let b = contour([(2, -5), (6, -5), (6, -2), (2, -2)]); + let expected = contour([(-8, -2), (2, -2), (2, 6), (-8, 6)]); + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let rhs = IntMinkowskiOperand::with_source(&b, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert_geometry(&solver.sum(&lhs, &rhs).unwrap(), &expected); + assert_geometry(&solver.sum(&rhs, &lhs).unwrap(), &expected); + } + + check::(); + check::(); + check::(); +} + +#[test] +fn builds_a_hexagon_from_triangles_with_opposite_directions() { + let a = int_path![[0, 0], [4, 0], [0, 4]]; + let b = int_path![[0, 0], [-2, 0], [0, -2]]; + let expected = int_path![[-2, 0], [0, -2], [4, -2], [4, 0], [0, 4], [-2, 4]]; + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let rhs = IntMinkowskiOperand::with_source(&b, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert_geometry(&solver.sum(&lhs, &rhs).unwrap(), &expected); + assert_geometry(&solver.sum(&rhs, &lhs).unwrap(), &expected); +} + +#[test] +fn merges_partial_sums_without_filling_concavities() { + let a = int_path![[0, 0], [4, 0], [4, 1], [1, 1], [1, 4], [0, 4]]; + let b = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let expected = int_path![[0, 0], [5, 0], [5, 2], [2, 2], [2, 5], [0, 5]]; + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let rhs = IntMinkowskiOperand::with_source(&b, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert_geometry(&solver.sum(&lhs, &rhs).unwrap(), &expected); + assert_geometry(&solver.sum(&rhs, &lhs).unwrap(), &expected); +} + +#[test] +fn preserves_and_closes_holes() { + let a = int_shapes![[ + [[0, 0], [10, 0], [10, 10], [0, 10]], + [[3, 3], [3, 7], [7, 7], [7, 3]] + ]]; + let small = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let large = int_path![[0, 0], [4, 0], [4, 4], [0, 4]]; + let expected = int_shapes![[ + [[0, 0], [11, 0], [11, 11], [0, 11]], + [[4, 4], [4, 7], [7, 7], [7, 4]] + ]]; + let filled = int_path![[0, 0], [14, 0], [14, 14], [0, 14]]; + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let small = IntMinkowskiOperand::with_source(&small, FillRule::NonZero).unwrap(); + let large = IntMinkowskiOperand::with_source(&large, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert_geometry(&solver.sum(&lhs, &small).unwrap(), &expected); + assert_geometry(&solver.sum(&lhs, &large).unwrap(), &filled); +} + +#[test] +fn reuses_solver_for_disconnected_and_merged_results() { + let a = int_shapes![ + [[[0, 0], [1, 0], [1, 1], [0, 1]]], + [[[4, 0], [5, 0], [5, 1], [4, 1]]] + ]; + let small = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let large = int_path![[0, 0], [4, 0], [4, 1], [0, 1]]; + let separate = int_shapes![ + [[[0, 0], [2, 0], [2, 2], [0, 2]]], + [[[4, 0], [6, 0], [6, 2], [4, 2]]] + ]; + let merged = int_path![[0, 0], [9, 0], [9, 2], [0, 2]]; + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let small = IntMinkowskiOperand::with_source(&small, FillRule::NonZero).unwrap(); + let large = IntMinkowskiOperand::with_source(&large, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::new(0); + for _ in 0..3 { + let result = solver.sum(&lhs, &small).unwrap(); + assert_eq!(result.len(), 2); + assert_geometry(&result, &separate); + let result = solver.sum(&lhs, &large).unwrap(); + assert_eq!(result.len(), 1); + assert_geometry(&result, &merged); + solver.clear(); + } +} + +#[test] +fn checks_all_bbox_bounds_and_recovers_after_errors() { + fn check() { + let limit = I::Wide::ONE << (I::BITS - 2); + let one = I::Wide::ONE; + let two = I::Wide::TWO; + let four = I::Wide::FOUR; + for swap_axes in [false, true] { + let rectangle = |min: I::Wide, max: I::Wide| -> IntContour { + [ + (min, I::Wide::ZERO), + (max, I::Wide::ZERO), + (max, two), + (min, two), + ] + .map(|(x, y)| { + let (x, y) = if swap_axes { (y, x) } else { (x, y) }; + IntPoint::new(I::from_wide(x), I::from_wide(y)) + }) + .into() + }; + let mut solver = IntMinkowskiSolver::default(); + for positive in [false, true] { + let a = if positive { + rectangle(limit - four, limit - two) + } else { + rectangle(-limit + two, -limit + four) + }; + let valid = if positive { + rectangle(I::Wide::ZERO, one) + } else { + rectangle(-one, I::Wide::ZERO) + }; + let invalid = if positive { + rectangle(one, two) + } else { + rectangle(-two, -one) + }; + let lhs = IntMinkowskiOperand::with_source(&a, FillRule::NonZero).unwrap(); + let valid = IntMinkowskiOperand::with_source(&valid, FillRule::NonZero).unwrap(); + let invalid = + IntMinkowskiOperand::with_source(&invalid, FillRule::NonZero).unwrap(); + let before = solver.sum(&lhs, &valid).unwrap(); + assert!(matches!( + solver.sum(&lhs, &invalid), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert!(matches!( + solver.sum(&lhs, &lhs), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert_geometry(&solver.sum(&lhs, &valid).unwrap(), &before); + } + } + } + + check::(); + check::(); + check::(); +} + +#[test] +fn rejects_source_bbox_even_when_normalized_sum_would_fit() { + let limit = 1_i32 << 30; + let square = int_path![[0, 0], [2, 0], [2, 2], [0, 2]]; + let remote = int_path![ + [limit - 4, 0], + [limit - 2, 0], + [limit - 2, 2], + [limit - 4, 2] + ]; + let source = vec![square.clone(), remote.clone(), remote]; + let kernel = int_path![[3, 0], [5, 0], [5, 2], [3, 2]]; + let lhs = IntMinkowskiOperand::with_source(&source, FillRule::EvenOdd).unwrap(); + let rhs = IntMinkowskiOperand::with_source(&kernel, FillRule::NonZero).unwrap(); + let retained = IntMinkowskiOperand::with_source(&square, FillRule::NonZero).unwrap(); + let mut solver = IntMinkowskiSolver::default(); + assert!(matches!( + solver.sum(&lhs, &rhs), + Err(MinkowskiError::CoordinateOutOfRange) + )); + assert_geometry( + &solver.sum(&retained, &rhs).unwrap(), + &int_path![[3, 0], [7, 0], [7, 4], [3, 4]], + ); +} + +#[test] +fn preparation_and_sum_share_overlay_without_leaking_fill_rules_or_geometry() { + let paths = vec![ + int_path![[0, 0], [10, 0], [10, 10], [0, 10]], + int_path![[3, 3], [7, 3], [7, 7], [3, 7]], + ]; + let square = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let filled = int_path![[0, 0], [11, 0], [11, 11], [0, 11]]; + let hole = int_shapes![[ + [[0, 0], [11, 0], [11, 11], [0, 11]], + [[4, 4], [4, 7], [7, 7], [7, 4]] + ]]; + let empty: IntContour = Vec::new(); + let mut solver = IntMinkowskiSolver::default(); + let kernel = solver + .prepare(square.as_slice(), FillRule::NonZero) + .unwrap(); + let original_kernel = kernel.convex_polygons().to_vec(); + + for _ in 0..3 { + for fill_rule in [FillRule::EvenOdd, FillRule::NonZero, FillRule::Positive] { + let operand = solver.prepare(&paths, fill_rule).unwrap(); + assert!(solver.prepare(&empty, FillRule::NonZero).is_none()); + assert!(solver.prepare(&paths, FillRule::Negative).is_none()); + let result = solver.sum(&operand, &kernel).unwrap(); + if fill_rule == FillRule::EvenOdd { + assert_geometry(&result, &hole); + } else { + assert_geometry(&result, &filled); + } + } + let reversed: Vec<_> = paths + .iter() + .map(|path| path.iter().rev().copied().collect::>()) + .collect(); + let operand = solver.prepare(&reversed, FillRule::Negative).unwrap(); + assert_geometry(&solver.sum(&operand, &kernel).unwrap(), &filled); + assert_eq!(kernel.convex_polygons(), original_kernel); + } +} + +#[test] +fn output_collinear_option_survives_preparation() { + let source = int_path![[0, 0], [4, 0], [4, 1], [1, 1], [1, 4], [0, 4]]; + let square = int_path![[0, 0], [1, 0], [1, 1], [0, 1]]; + let expected = int_path![[0, 0], [5, 0], [5, 2], [2, 2], [2, 5], [0, 5]]; + let empty: IntContour = Vec::new(); + for preserve_output_collinear in [false, true] { + let mut solver = IntMinkowskiSolver::with_options(0, preserve_output_collinear); + for _ in 0..3 { + let operand = solver.prepare(&source, FillRule::NonZero).unwrap(); + let kernel = solver.prepare(&square, FillRule::NonZero).unwrap(); + assert!(operand + .convex_polygons() + .iter() + .all(|p| !p.is_clockwise_ordered())); + assert!(solver.prepare(&empty, FillRule::NonZero).is_none()); + assert!(solver.prepare(&square, FillRule::Negative).is_none()); + + let result = solver.sum(&operand, &kernel).unwrap(); + assert_eq!(result.len(), 1); + assert!(!result[0][0].is_clockwise_ordered()); + assert_geometry(&result, &expected); + if preserve_output_collinear { + assert!(result[0][0].len() > expected.len()); + } else { + assert_eq!(result[0][0].len(), expected.len()); + } + solver.clear(); + } + } +} + +fn cells_to_contours(cells: &BTreeSet<(i32, i32)>) -> Vec> { + cells + .iter() + .map(|&(x, y)| int_path![[x, y], [x + 1, y], [x + 1, y + 1], [x, y + 1]]) + .collect() +} + +fn contains_cell_center(shapes: &IntShapes, x: i32, y: i32) -> bool { + let px = 2 * x + 1; + let py = 2 * y + 1; + let mut winding = 0; + for contour in shapes.iter().flatten() { + for (&a, &b) in contour.iter().zip(contour.iter().cycle().skip(1)) { + let cross = (b.x - a.x) * (py - 2 * a.y) - (b.y - a.y) * (px - 2 * a.x); + if 2 * a.y <= py && 2 * b.y > py && cross > 0 { + winding += 1; + } else if 2 * a.y > py && 2 * b.y <= py && cross < 0 { + winding -= 1; + } + } + } + winding != 0 +} + +#[test] +fn random_cell_shapes_match_independent_discrete_sum_and_difference() { + let mut rng = StdRng::seed_from_u64(0x4d_49_4e_4b); + let mut solver = IntMinkowskiSolver::new(0); + for case in 0..128 { + let mut make_cells = || { + let mut cells = BTreeSet::new(); + for x in -2..=2 { + for y in -2..=2 { + if rng.random_bool(0.3) { + cells.insert((x, y)); + } + } + } + if cells.is_empty() { + cells.insert((0, 0)); + } + cells + }; + let a = make_cells(); + let b = make_cells(); + let mut expected_sum = BTreeSet::new(); + let mut expected_difference = BTreeSet::new(); + for &(ax, ay) in &a { + for &(bx, by) in &b { + for dx in 0..=1 { + for dy in 0..=1 { + expected_sum.insert((ax + bx + dx, ay + by + dy)); + expected_difference.insert((ax - bx - 1 + dx, ay - by - 1 + dy)); + } + } + } + } + + let oblique = case % 2 != 0; + let transform = |cells: &BTreeSet<(i32, i32)>| { + let mut contours = cells_to_contours(cells); + if oblique { + for p in contours.iter_mut().flatten() { + *p = IntPoint::new(2 * p.x + p.y, p.x + p.y); + } + } + contours + }; + let lhs = solver.prepare(&transform(&a), FillRule::NonZero).unwrap(); + let rhs = solver.prepare(&transform(&b), FillRule::NonZero).unwrap(); + for (lhs, rhs, difference) in [(&lhs, &rhs, false), (&rhs, &lhs, false), (&lhs, &rhs, true)] + { + let (mut result, expected) = if difference { + (solver.difference(lhs, rhs).unwrap(), &expected_difference) + } else { + (solver.sum(lhs, rhs).unwrap(), &expected_sum) + }; + if oblique { + for p in result.iter_mut().flatten().flatten() { + *p = IntPoint::new(p.x - p.y, 2 * p.y - p.x); + } + } + let mut area_two = 0; + for contour in result.iter().flatten() { + for (&p, &q) in contour.iter().zip(contour.iter().cycle().skip(1)) { + assert!( + p.x == q.x || p.y == q.y, + "case {case}: non-axis-aligned boundary" + ); + area_two += p.x * q.y - p.y * q.x; + } + } + assert_eq!( + area_two, + 2 * expected.len() as i32, + "case {case}: area mismatch" + ); + for x in -5..=6 { + for y in -5..=6 { + assert_eq!( + contains_cell_center(&result, x, y), + expected.contains(&(x, y)), + "case {case}, cell ({x}, {y}), a={a:?}, b={b:?}" + ); + } + } + } + } +}