From 83957c4bfcb20375bd0563aa0aa55c63939c7186 Mon Sep 17 00:00:00 2001 From: Nick P Date: Tue, 14 Jul 2026 14:58:40 -0600 Subject: [PATCH] feat: complete relate and buffer strategy parity --- crates/geometry-overlay/src/buffer.rs | 686 +++++++++++++++- crates/geometry-overlay/src/lib.rs | 1 + crates/geometry-overlay/src/relate.rs | 740 +++++++++++++++++- crates/geometry-strategy/src/buffer.rs | 96 +++ crates/geometry-strategy/src/lib.rs | 3 +- .../geometry/tests/buffer_strategy_parity.rs | 218 +++++- crates/geometry/tests/relate_pair_parity.rs | 154 +++- docs/03-overlay-engine.md | 26 +- docs/crates/geometry-overlay.md | 15 +- docs/crates/geometry-strategy.md | 2 +- 10 files changed, 1877 insertions(+), 64 deletions(-) diff --git a/crates/geometry-overlay/src/buffer.rs b/crates/geometry-overlay/src/buffer.rs index 6819817..e7b076a 100644 --- a/crates/geometry-overlay/src/buffer.rs +++ b/crates/geometry-overlay/src/buffer.rs @@ -5,7 +5,13 @@ //! the input outward by `distance`, rounding or mitering the corners, //! and unions the offset pieces into an output polygon. //! -//! Cartesian dispatch covers points, linestrings, and simple polygons. +//! Cartesian dispatch covers every static single and homogeneous multi kind. +//! Spherical and geographic inputs are projected into a local tangent plane, +//! buffered by the same Cartesian engine, and transformed back. The angular +//! path is intended for local buffers: unlike Boost's per-segment geodesic +//! offset formulas, its error grows with the geometry's angular extent and it +//! rejects projection centers at the poles. This deliberate approximation is +//! recorded in the project feature-parity map for later reassessment. //! Polygon offsets are signed, handle convex and reflex vertices, and move //! interior rings in the opposite topological direction from the exterior. //! @@ -32,13 +38,19 @@ use alloc::vec::Vec; use geometry_coords::{ CoordinateScalar, - math::{atan2, ceil, cos, hypot, mul_add, sin}, + math::{atan2, ceil, cos, hypot, mul_add, sin, sqrt}, +}; +use geometry_cs::{ + AngleUnit, Cartesian, CartesianFamily, CoordinateSystem, FromF64, Geographic, GeographicFamily, + Spherical, SphericalFamily, +}; +use geometry_model::{ + Box as ModelBox, Linestring, MultiLinestring, MultiPoint, MultiPolygon, Point2D, Polygon, Ring, }; -use geometry_cs::{CartesianFamily, CoordinateSystem, FromF64}; -use geometry_model::{Linestring, MultiPolygon, Polygon, Ring}; use geometry_strategy::buffer::{ BufferDistanceStrategy, BufferEndStrategy, BufferJoinStrategy, BufferPointStrategy, - BufferSettings, + BufferSettings, CartesianBuffer, DefaultBuffer, DefaultBufferStrategy, GeographicBuffer, + SphericalBuffer, }; use geometry_tag::{ BoxTag, LinestringTag, MultiLinestringTag, MultiPointTag, MultiPolygonTag, PointTag, @@ -100,11 +112,12 @@ pub enum PointStrategy { /// `boost::geometry::buffer` in /// `algorithms/detail/buffer/interface.hpp:246-273`. #[doc(hidden)] -pub trait BufferStrategy { +pub trait BufferStrategy { fn apply( &self, geometry: &G, settings: BufferSettings, + coordinate_strategy: &CoordinateStrategy, ) -> Result>, OverlayError>; } @@ -212,11 +225,11 @@ impl BufferStrategyForKind for MultiPolygonTag { /// Buffer a geometry using the public point and join strategies. /// /// Mirrors `boost::geometry::buffer` from -/// `boost/geometry/algorithms/detail/buffer/interface.hpp:246-273`. Cartesian -/// dispatch supports point, segment, linestring, ring, polygon, box, and all -/// three homogeneous multi-geometry kinds. Point inputs use `point`, linear -/// inputs use all five strategy roles, and areal inputs use signed distance -/// and join policies. +/// `boost/geometry/algorithms/detail/buffer/interface.hpp:246-273`. Cartesian, +/// spherical, and geographic dispatch supports point, segment, linestring, +/// ring, polygon, box, and all three homogeneous multi-geometry kinds. Point +/// inputs use `point`, linear inputs use all five strategy roles, and areal +/// inputs use signed distance and join policies. /// /// # Errors /// @@ -233,7 +246,9 @@ pub fn buffer( where G: Geometry, G::Kind: BufferStrategyForKind, - ::S: BufferStrategy, + <::Cs as CoordinateSystem>::Family: + DefaultBuffer<<::Cs as CoordinateSystem>::Family>, + ::S: BufferStrategy>, { let settings = BufferSettings { distance: BufferDistanceStrategy::Symmetric(distance), @@ -276,14 +291,56 @@ pub fn buffer_with( where G: Geometry, G::Kind: BufferStrategyForKind, - ::S: BufferStrategy, + <::Cs as CoordinateSystem>::Family: + DefaultBuffer<<::Cs as CoordinateSystem>::Family>, + ::S: BufferStrategy>, { - <::S as Default>::default().apply(geometry, settings) + buffer_with_strategy(geometry, settings, DefaultBufferStrategy::::default()) +} + +/// Buffer a geometry with explicit coordinate-system and five-role strategy +/// bundles. +/// +/// Mirrors the explicit strategy overload of `boost::geometry::buffer` from +/// `algorithms/detail/buffer/interface.hpp:246-273`, together with the +/// Cartesian, spherical, and geographic umbrella strategies under +/// `strategies/buffer/`. +/// +/// [`SphericalBuffer`] and [`GeographicBuffer`] use a geometry-centered local +/// tangent projection before invoking the Cartesian offset engine. This keeps +/// distance units explicit and `no_std` compatible, but is a local-extent +/// approximation rather than Boost's per-segment geodesic construction. +/// +/// # Errors +/// +/// Returns [`OverlayError::Unsupported`] for invalid strategy values, +/// non-finite/inapplicable distances, or degenerate linear input. +#[inline] +#[must_use = "buffering can fail and the generated geometry should be used"] +#[allow( + clippy::needless_pass_by_value, + reason = "Boost buffer coordinate strategies are small value objects passed explicitly" +)] +pub fn buffer_with_strategy( + geometry: &G, + settings: BufferSettings, + coordinate_strategy: CoordinateStrategy, +) -> Result>, OverlayError> +where + G: Geometry, + G::Kind: BufferStrategyForKind, + ::S: BufferStrategy, +{ + <::S as Default>::default().apply( + geometry, + settings, + &coordinate_strategy, + ) } /// Implements the point arm selected by `buffer_all` at /// `algorithms/detail/buffer/interface.hpp:269-273`. -impl BufferStrategy for PointBuffer +impl BufferStrategy for PointBuffer where G: Point + PointMut + Default + Copy, G::Scalar: CoordinateScalar + Into + FromF64, @@ -293,6 +350,7 @@ where &self, point_geometry: &G, settings: BufferSettings, + _coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let BufferDistanceStrategy::Symmetric(distance) = settings.distance else { return Err(OverlayError::Unsupported); @@ -316,7 +374,7 @@ where /// Implements the polygon arm selected by `buffer_all` at /// `algorithms/detail/buffer/interface.hpp:269-273`. -impl BufferStrategy for PolygonBuffer +impl BufferStrategy for PolygonBuffer where G: PolygonTrait, G::Point: PointMut + Default + Copy, @@ -327,6 +385,7 @@ where &self, polygon: &G, settings: BufferSettings, + _coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let BufferDistanceStrategy::Symmetric(distance) = settings.distance else { return Err(OverlayError::Unsupported); @@ -357,7 +416,7 @@ where } } -impl BufferStrategy for LinestringBuffer +impl BufferStrategy for LinestringBuffer where G: LinestringTrait, G::Point: PointMut + Default + Copy, @@ -368,6 +427,7 @@ where &self, line: &G, settings: BufferSettings, + _coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let (left, right) = match settings.distance { BufferDistanceStrategy::Symmetric(distance) => (distance, distance), @@ -384,7 +444,7 @@ where } } -impl BufferStrategy for SegmentBuffer +impl BufferStrategy for SegmentBuffer where G: SegmentTrait, G::Point: PointMut + Default + Copy, @@ -395,14 +455,15 @@ where &self, segment: &G, settings: BufferSettings, + coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let line: Linestring = Linestring::from_vec(alloc::vec![segment_start(segment), segment_end(segment)]); - LinestringBuffer.apply(&line, settings) + LinestringBuffer.apply(&line, settings, coordinate_strategy) } } -impl BufferStrategy for RingBuffer +impl BufferStrategy for RingBuffer where G: RingTrait, G::Point: PointMut + Default + Copy, @@ -413,6 +474,7 @@ where &self, ring: &G, settings: BufferSettings, + _coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let BufferDistanceStrategy::Symmetric(distance) = settings.distance else { return Err(OverlayError::Unsupported); @@ -427,7 +489,7 @@ where } } -impl BufferStrategy for BoxBuffer +impl BufferStrategy for BoxBuffer where G: BoxTrait, G::Point: PointMut + Default + Copy, @@ -438,6 +500,7 @@ where &self, bounds: &G, settings: BufferSettings, + coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let minimum = box_min(bounds); let maximum = box_max(bounds); @@ -452,11 +515,11 @@ where make_point(max_x, min_y), make_point(min_x, min_y), ]); - RingBuffer.apply(&ring, settings) + RingBuffer.apply(&ring, settings, coordinate_strategy) } } -impl BufferStrategy for MultiPointBuffer +impl BufferStrategy for MultiPointBuffer where G: MultiPointTrait::Point>, G::Point: PointMut + Default + Copy, @@ -467,16 +530,19 @@ where &self, points: &G, settings: BufferSettings, + coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let mut output = MultiPolygon::new(); for point in points.points() { - output.0.extend(PointBuffer.apply(point, settings)?.0); + output + .0 + .extend(PointBuffer.apply(point, settings, coordinate_strategy)?.0); } crate::merge::merge_polygons(output.0) } } -impl BufferStrategy for MultiLinestringBuffer +impl BufferStrategy for MultiLinestringBuffer where G: MultiLinestringTrait, G::Point: PointMut + Default + Copy, @@ -487,16 +553,21 @@ where &self, lines: &G, settings: BufferSettings, + coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let mut output = MultiPolygon::new(); for line in lines.linestrings() { - output.0.extend(LinestringBuffer.apply(line, settings)?.0); + output.0.extend( + LinestringBuffer + .apply(line, settings, coordinate_strategy)? + .0, + ); } crate::merge::merge_polygons(output.0) } } -impl BufferStrategy for MultiPolygonBuffer +impl BufferStrategy for MultiPolygonBuffer where G: MultiPolygonTrait, G::Point: PointMut + Default + Copy, @@ -507,15 +578,574 @@ where &self, polygons: &G, settings: BufferSettings, + coordinate_strategy: &CartesianBuffer, ) -> Result>, OverlayError> { let mut output = MultiPolygon::new(); for polygon in polygons.polygons() { - output.0.extend(PolygonBuffer.apply(polygon, settings)?.0); + output.0.extend( + PolygonBuffer + .apply(polygon, settings, coordinate_strategy)? + .0, + ); } crate::merge::merge_polygons(output.0) } } +trait AngularCoordinateSystem { + type Units: AngleUnit; +} + +impl AngularCoordinateSystem for Spherical { + type Units = Units; +} + +impl AngularCoordinateSystem for Geographic { + type Units = Units; +} + +#[derive(Debug, Clone, Copy)] +struct LocalProjection { + longitude: f64, + latitude: f64, + east_scale: f64, + north_scale: f64, +} + +impl LocalProjection { + fn project(self, longitude: f64, latitude: f64) -> (f64, f64) { + let mut delta_longitude = longitude - self.longitude; + if delta_longitude > core::f64::consts::PI { + delta_longitude -= 2.0 * core::f64::consts::PI; + } else if delta_longitude < -core::f64::consts::PI { + delta_longitude += 2.0 * core::f64::consts::PI; + } + ( + delta_longitude * self.east_scale, + (latitude - self.latitude) * self.north_scale, + ) + } + + fn unproject(self, x: f64, y: f64) -> (f64, f64) { + let mut longitude = self.longitude + x / self.east_scale; + if longitude > core::f64::consts::PI { + longitude -= 2.0 * core::f64::consts::PI; + } else if longitude < -core::f64::consts::PI { + longitude += 2.0 * core::f64::consts::PI; + } + (longitude, self.latitude + y / self.north_scale) + } +} + +trait AngularBufferProjection { + fn projection(&self, longitude: f64, latitude: f64) -> Result; +} + +impl AngularBufferProjection for SphericalBuffer { + fn projection(&self, longitude: f64, latitude: f64) -> Result { + if !self.radius.is_finite() || self.radius <= 0.0 { + return Err(OverlayError::Unsupported); + } + let east_scale = self.radius * cos(latitude); + if east_scale.abs() <= f64::EPSILON { + return Err(OverlayError::Unsupported); + } + Ok(LocalProjection { + longitude, + latitude, + east_scale, + north_scale: self.radius, + }) + } +} + +impl AngularBufferProjection for GeographicBuffer { + fn projection(&self, longitude: f64, latitude: f64) -> Result { + let spheroid = self.spheroid; + if !spheroid.equatorial_radius.is_finite() + || spheroid.equatorial_radius <= 0.0 + || !spheroid.flattening.is_finite() + || !(0.0..1.0).contains(&spheroid.flattening) + { + return Err(OverlayError::Unsupported); + } + + let eccentricity_squared = spheroid.eccentricity_squared(); + let sin_latitude = sin(latitude); + let denominator = sqrt(1.0 - eccentricity_squared * sin_latitude * sin_latitude); + let prime_vertical = spheroid.equatorial_radius / denominator; + let meridional = spheroid.equatorial_radius * (1.0 - eccentricity_squared) + / (denominator * denominator * denominator); + let east_scale = prime_vertical * cos(latitude); + if east_scale.abs() <= f64::EPSILON { + return Err(OverlayError::Unsupported); + } + Ok(LocalProjection { + longitude, + latitude, + east_scale, + north_scale: meridional, + }) + } +} + +fn angular_coordinates

(point: &P) -> (f64, f64) +where + P: Point, + P::Scalar: Into, + P::Cs: AngularCoordinateSystem, +{ + let longitude = ::Units::to_radians(point.get::<0>().into()); + let latitude = ::Units::to_radians(point.get::<1>().into()); + (longitude, latitude) +} + +fn angular_point

(longitude: f64, latitude: f64) -> P +where + P: PointMut + Default, + P::Scalar: FromF64, + P::Cs: AngularCoordinateSystem, +{ + let mut point = P::default(); + let longitude = ::Units::from_radians(longitude); + let latitude = ::Units::from_radians(latitude); + point.set::<0>(P::Scalar::from_f64(longitude)); + point.set::<1>(P::Scalar::from_f64(latitude)); + point +} + +fn projection_center(coordinates: &[(f64, f64)]) -> Result<(f64, f64), OverlayError> { + if coordinates.is_empty() { + return Err(OverlayError::Unsupported); + } + let mut longitude_sine = 0.0; + let mut longitude_cosine = 0.0; + let mut latitude = 0.0; + for &(longitude, point_latitude) in coordinates { + longitude_sine += sin(longitude); + longitude_cosine += cos(longitude); + latitude += point_latitude; + } + let count = coordinates.len() as f64; + Ok((atan2(longitude_sine, longitude_cosine), latitude / count)) +} + +type ProjectedPoint = Point2D; + +fn projected_point

(point: &P, projection: LocalProjection) -> ProjectedPoint +where + P: Point, + P::Scalar: Into, + P::Cs: AngularCoordinateSystem, +{ + let (longitude, latitude) = angular_coordinates(point); + let (x, y) = projection.project(longitude, latitude); + ProjectedPoint::new(x, y) +} + +fn projected_ring(ring: &R, projection: LocalProjection) -> Ring +where + R: RingTrait, + R::Point: Point, + ::Scalar: Into, + ::Cs: AngularCoordinateSystem, +{ + Ring::from_vec( + ring.points() + .map(|point| projected_point(point, projection)) + .collect(), + ) +} + +fn projected_polygon(polygon: &G, projection: LocalProjection) -> Polygon +where + G: PolygonTrait, + G::Point: Point, + ::Scalar: Into, + ::Cs: AngularCoordinateSystem, +{ + Polygon::with_inners( + projected_ring(polygon.exterior(), projection), + polygon + .interiors() + .map(|ring| projected_ring(ring, projection)) + .collect(), + ) +} + +fn unprojected_buffer

( + polygons: MultiPolygon>, + projection: LocalProjection, +) -> MultiPolygon> +where + P: PointMut + Default, + P::Scalar: FromF64, + P::Cs: AngularCoordinateSystem, +{ + MultiPolygon::from_vec( + polygons + .0 + .into_iter() + .map(|polygon| { + let outer = Ring::from_vec( + polygon + .outer + .0 + .into_iter() + .map(|point| { + let (longitude, latitude) = projection.unproject(point.x(), point.y()); + angular_point(longitude, latitude) + }) + .collect(), + ); + let inners = polygon + .inners + .into_iter() + .map(|ring| { + Ring::from_vec( + ring.0 + .into_iter() + .map(|point| { + let (longitude, latitude) = + projection.unproject(point.x(), point.y()); + angular_point(longitude, latitude) + }) + .collect(), + ) + }) + .collect(); + Polygon::with_inners(outer, inners) + }) + .collect(), + ) +} + +fn projection_for_points<'a, P>( + points: impl IntoIterator, + strategy: &impl AngularBufferProjection, +) -> Result +where + P: Point + 'a, + P::Scalar: Into, + P::Cs: AngularCoordinateSystem, +{ + let coordinates: Vec<_> = points.into_iter().map(angular_coordinates).collect(); + let (longitude, latitude) = projection_center(&coordinates)?; + strategy.projection(longitude, latitude) +} + +fn projected_point_apply

( + point: &P, + settings: BufferSettings, + strategy: &impl AngularBufferProjection, +) -> Result>, OverlayError> +where + P: Point + PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into + FromF64, + P::Cs: AngularCoordinateSystem, +{ + let projection = projection_for_points(core::iter::once(point), strategy)?; + let point = projected_point(point, projection); + let output = PointBuffer.apply(&point, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) +} + +fn projected_linestring_apply( + line: &L, + settings: BufferSettings, + strategy: &impl AngularBufferProjection, +) -> Result>, OverlayError> +where + L: LinestringTrait, + L::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, +{ + let projection = projection_for_points(line.points(), strategy)?; + let projected = Linestring::from_vec( + line.points() + .map(|point| projected_point(point, projection)) + .collect(), + ); + let output = LinestringBuffer.apply(&projected, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) +} + +fn projected_ring_apply( + ring: &R, + settings: BufferSettings, + strategy: &impl AngularBufferProjection, +) -> Result>, OverlayError> +where + R: RingTrait, + R::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, +{ + let projection = projection_for_points(ring.points(), strategy)?; + let output = RingBuffer.apply( + &projected_ring(ring, projection), + settings, + &CartesianBuffer, + )?; + Ok(unprojected_buffer(output, projection)) +} + +fn projected_polygon_apply( + polygon: &G, + settings: BufferSettings, + strategy: &impl AngularBufferProjection, +) -> Result>, OverlayError> +where + G: PolygonTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, +{ + let mut coordinates = polygon + .exterior() + .points() + .map(angular_coordinates) + .collect::>(); + for ring in polygon.interiors() { + coordinates.extend(ring.points().map(angular_coordinates)); + } + let (longitude, latitude) = projection_center(&coordinates)?; + let projection = strategy.projection(longitude, latitude)?; + let output = PolygonBuffer.apply( + &projected_polygon(polygon, projection), + settings, + &CartesianBuffer, + )?; + Ok(unprojected_buffer(output, projection)) +} + +macro_rules! impl_angular_buffer_strategy { + ($strategy:ty, $family:ty) => { + impl BufferStrategy for PointBuffer + where + G: Point + PointMut + Default + Copy, + G::Scalar: CoordinateScalar + Into + FromF64, + G::Cs: AngularCoordinateSystem, + ::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + projected_point_apply(geometry, settings, coordinate_strategy) + } + } + + impl BufferStrategy for LinestringBuffer + where + G: LinestringTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + projected_linestring_apply(geometry, settings, coordinate_strategy) + } + } + + impl BufferStrategy for SegmentBuffer + where + G: SegmentTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + let line = Linestring::from_vec(alloc::vec![ + segment_start(geometry), + segment_end(geometry), + ]); + projected_linestring_apply(&line, settings, coordinate_strategy) + } + } + + impl BufferStrategy for RingBuffer + where + G: RingTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + projected_ring_apply(geometry, settings, coordinate_strategy) + } + } + + impl BufferStrategy for PolygonBuffer + where + G: PolygonTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + projected_polygon_apply(geometry, settings, coordinate_strategy) + } + } + + impl BufferStrategy for BoxBuffer + where + G: BoxTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + let minimum = box_min(geometry); + let maximum = box_max(geometry); + let projection = projection_for_points([&minimum, &maximum], coordinate_strategy)?; + let projected = ModelBox::from_corners( + projected_point(&minimum, projection), + projected_point(&maximum, projection), + ); + let output = BoxBuffer.apply(&projected, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) + } + } + + impl BufferStrategy for MultiPointBuffer + where + G: MultiPointTrait::Point>, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + let projection = projection_for_points(geometry.points(), coordinate_strategy)?; + let projected = MultiPoint::from_vec( + geometry + .points() + .map(|point| projected_point(point, projection)) + .collect(), + ); + let output = MultiPointBuffer.apply(&projected, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) + } + } + + impl BufferStrategy for MultiLinestringBuffer + where + G: MultiLinestringTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + let coordinates = geometry + .linestrings() + .flat_map(|line| line.points().map(angular_coordinates)) + .collect::>(); + let (longitude, latitude) = projection_center(&coordinates)?; + let projection = coordinate_strategy.projection(longitude, latitude)?; + let projected = MultiLinestring::from_vec( + geometry + .linestrings() + .map(|line| { + Linestring::from_vec( + line.points() + .map(|point| projected_point(point, projection)) + .collect(), + ) + }) + .collect(), + ); + let output = MultiLinestringBuffer.apply(&projected, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) + } + } + + impl BufferStrategy for MultiPolygonBuffer + where + G: MultiPolygonTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, + ::Cs: AngularCoordinateSystem, + <::Cs as CoordinateSystem>::Family: SameAs<$family>, + { + fn apply( + &self, + geometry: &G, + settings: BufferSettings, + coordinate_strategy: &$strategy, + ) -> Result>, OverlayError> { + let coordinates = geometry + .polygons() + .flat_map(|polygon| { + polygon + .exterior() + .points() + .chain(polygon.interiors().flat_map(RingTrait::points)) + .map(angular_coordinates) + }) + .collect::>(); + let (longitude, latitude) = projection_center(&coordinates)?; + let projection = coordinate_strategy.projection(longitude, latitude)?; + let projected = MultiPolygon::from_vec( + geometry + .polygons() + .map(|polygon| projected_polygon(polygon, projection)) + .collect(), + ); + let output = MultiPolygonBuffer.apply(&projected, settings, &CartesianBuffer)?; + Ok(unprojected_buffer(output, projection)) + } + } + }; +} + +impl_angular_buffer_strategy!(SphericalBuffer, SphericalFamily); +impl_angular_buffer_strategy!(GeographicBuffer, GeographicFamily); + /// Buffer a point by `distance`, producing the disc (or square) /// approximation. /// diff --git a/crates/geometry-overlay/src/lib.rs b/crates/geometry-overlay/src/lib.rs index 7d9fa06..f5f4bd7 100644 --- a/crates/geometry-overlay/src/lib.rs +++ b/crates/geometry-overlay/src/lib.rs @@ -44,6 +44,7 @@ pub mod validity; pub use buffer::{ JoinStrategy, PointStrategy, buffer, buffer_convex_polygon, buffer_point, buffer_with, + buffer_with_strategy, }; pub use merge::{merge_elements, merge_multipolygon, merge_polygons}; pub use operation::{OverlayError, difference, intersection, sym_difference, r#union, union_poly}; diff --git a/crates/geometry-overlay/src/relate.rs b/crates/geometry-overlay/src/relate.rs index 0b9cde7..456b903 100644 --- a/crates/geometry-overlay/src/relate.rs +++ b/crates/geometry-overlay/src/relate.rs @@ -9,26 +9,35 @@ //! (`algorithms/{crosses,overlaps,touches}.hpp`) are then thin tests on //! that matrix. //! -//! Cartesian point, linestring, and polygon pairs are dispatched through the -//! same matrix interface. Polygon relations include interior rings and exact -//! point/curve dimensions for colocated boundaries; the public mask-string -//! interface consumes the completed matrix. +//! Cartesian points, segments, linestrings, rings, boxes, polygons, +//! homogeneous multis, runtime geometries, and heterogeneous geometry +//! collections are dispatched through the same matrix interface. Union +//! topology follows OGC boundary rules, including the mod-2 boundary of a +//! multilinestring and absorption of members covered by areal interiors. The +//! public mask-string interface consumes the completed matrix. #![allow( clippy::float_cmp, reason = "exact equality identifies stored endpoint identity for DE-9IM boundary classification" )] +use alloc::vec::Vec; + use geometry_coords::{CoordinateScalar, precise_math}; -use geometry_cs::{CartesianFamily, CoordinateSystem}; -use geometry_tag::{LinestringTag, PointTag, PolygonTag, SameAs}; +use geometry_cs::{Cartesian, CartesianFamily, CoordinateSystem}; +use geometry_model::{DynGeometry, Point2D, Polygon, Ring}; +use geometry_tag::{ + BoxTag, DynamicGeometryTag, GeometryCollectionTag, LinestringTag, MultiLinestringTag, + MultiPointTag, MultiPolygonTag, PointTag, PolygonTag, RingTag, SameAs, SegmentTag, +}; use geometry_trait::{ - Geometry, Linestring as LinestringTrait, Point, PointMut, Polygon as PolygonTrait, - Ring as RingTrait, + Box as BoxTrait, Geometry, GeometryCollection, Linestring as LinestringTrait, MultiLinestring, + MultiPoint, MultiPolygon, Point, PointMut, Polygon as PolygonTrait, Ring as RingTrait, + Segment as SegmentTrait, corner, }; use crate::operation::OverlayError; -use crate::predicate::range_guard::polygon_in_range; +use crate::predicate::range_guard::{SAFE_ABS_MAX, polygon_in_range}; /// The dimension of an intersection cell in a [`De9im`] matrix. /// @@ -220,6 +229,9 @@ pub struct RelatePolygonLinestring; #[doc(hidden)] #[derive(Debug, Default, Clone, Copy)] pub struct RelatePolygonPolygon; +#[doc(hidden)] +#[derive(Debug, Default, Clone, Copy)] +pub struct RelateTopology; impl RelatePairStrategy for PointTag { type Strategy = RelatePointPoint; @@ -249,13 +261,97 @@ impl RelatePairStrategy for PolygonTag { type Strategy = RelatePolygonPolygon; } +trait TopologyKind {} + +impl TopologyKind for PointTag {} +impl TopologyKind for LinestringTag {} +impl TopologyKind for PolygonTag {} +impl TopologyKind for SegmentTag {} +impl TopologyKind for RingTag {} +impl TopologyKind for BoxTag {} +impl TopologyKind for MultiPointTag {} +impl TopologyKind for MultiLinestringTag {} +impl TopologyKind for MultiPolygonTag {} +impl TopologyKind for DynamicGeometryTag {} +impl TopologyKind for GeometryCollectionTag {} + +macro_rules! topology_pair_for_single { + ($single:ty, $($other:ty),+ $(,)?) => { + $( + impl RelatePairStrategy<$other> for $single { + type Strategy = RelateTopology; + } + )+ + }; +} + +topology_pair_for_single!( + PointTag, + SegmentTag, + RingTag, + BoxTag, + MultiPointTag, + MultiLinestringTag, + MultiPolygonTag, + DynamicGeometryTag, + GeometryCollectionTag, +); +topology_pair_for_single!( + LinestringTag, + SegmentTag, + RingTag, + BoxTag, + MultiPointTag, + MultiLinestringTag, + MultiPolygonTag, + DynamicGeometryTag, + GeometryCollectionTag, +); +topology_pair_for_single!( + PolygonTag, + SegmentTag, + RingTag, + BoxTag, + MultiPointTag, + MultiLinestringTag, + MultiPolygonTag, + DynamicGeometryTag, + GeometryCollectionTag, +); + +impl RelatePairStrategy for SegmentTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for RingTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for BoxTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for MultiPointTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for MultiLinestringTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for MultiPolygonTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for DynamicGeometryTag { + type Strategy = RelateTopology; +} +impl RelatePairStrategy for GeometryCollectionTag { + type Strategy = RelateTopology; +} + type PairStrategy = <::Kind as RelatePairStrategy<::Kind>>::Strategy; /// Compute the DE-9IM matrix for a supported pointlike, linear, or areal pair. /// /// Mirrors the pair dispatch in -/// `algorithms/detail/relate/interface.hpp:275-382`. +/// `algorithms/detail/relate/interface.hpp:275-382`, including the +/// geometry-collection path in `detail/relate/implementation_gc.hpp`. /// /// # Errors /// @@ -385,6 +481,286 @@ where } } +type TopologyPointModel = Point2D; + +#[derive(Debug, Default, Clone)] +struct Topology { + points: Vec<[f64; 2]>, + lines: Vec>, + polygons: Vec>, +} + +impl Topology { + fn in_range(&self) -> bool { + let in_range = |point: [f64; 2]| { + point[0].is_finite() + && point[1].is_finite() + && point[0].abs() <= SAFE_ABS_MAX + && point[1].abs() <= SAFE_ABS_MAX + }; + if !self + .points + .iter() + .chain(self.lines.iter().flatten()) + .copied() + .all(in_range) + { + return false; + } + self.polygons.iter().all(|polygon| { + polygon + .outer + .0 + .iter() + .chain(polygon.inners.iter().flat_map(|ring| ring.0.iter())) + .all(|point| in_range([point.x(), point.y()])) + }) + } +} + +trait TopologyBuilder { + fn append(&self, geometry: &G, topology: &mut Topology); +} + +trait TopologyBuilderForKind { + type Strategy: Default; +} + +#[derive(Debug, Default, Clone, Copy)] +struct TopologyPoint; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyLinestring; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyPolygon; +#[derive(Debug, Default, Clone, Copy)] +struct TopologySegment; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyRing; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyBox; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyMultiPoint; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyMultiLinestring; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyMultiPolygon; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyDynamic; +#[derive(Debug, Default, Clone, Copy)] +struct TopologyCollection; + +impl TopologyBuilderForKind for PointTag { + type Strategy = TopologyPoint; +} +impl TopologyBuilderForKind for LinestringTag { + type Strategy = TopologyLinestring; +} +impl TopologyBuilderForKind for PolygonTag { + type Strategy = TopologyPolygon; +} +impl TopologyBuilderForKind for SegmentTag { + type Strategy = TopologySegment; +} +impl TopologyBuilderForKind for RingTag { + type Strategy = TopologyRing; +} +impl TopologyBuilderForKind for BoxTag { + type Strategy = TopologyBox; +} +impl TopologyBuilderForKind for MultiPointTag { + type Strategy = TopologyMultiPoint; +} +impl TopologyBuilderForKind for MultiLinestringTag { + type Strategy = TopologyMultiLinestring; +} +impl TopologyBuilderForKind for MultiPolygonTag { + type Strategy = TopologyMultiPolygon; +} +impl TopologyBuilderForKind for DynamicGeometryTag { + type Strategy = TopologyDynamic; +} +impl TopologyBuilderForKind for GeometryCollectionTag { + type Strategy = TopologyCollection; +} + +type TopologyBuilderStrategy = <::Kind as TopologyBuilderForKind>::Strategy; + +impl TopologyBuilder for TopologyPoint +where + G: Point, + G::Scalar: Into, + ::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + topology.points.push(xy(geometry)); + } +} + +impl TopologyBuilder for TopologyLinestring +where + G: LinestringTrait, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + append_topology_line(geometry.points().map(xy).collect(), topology); + } +} + +impl TopologyBuilder for TopologyPolygon +where + G: PolygonTrait, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + append_topology_polygon(geometry, topology); + } +} + +impl TopologyBuilder for TopologySegment +where + G: SegmentTrait, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + append_topology_line( + alloc::vec![ + [ + geometry.get_indexed::<0, 0>().into(), + geometry.get_indexed::<0, 1>().into(), + ], + [ + geometry.get_indexed::<1, 0>().into(), + geometry.get_indexed::<1, 1>().into(), + ], + ], + topology, + ); + } +} + +impl TopologyBuilder for TopologyRing +where + G: RingTrait, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + topology + .polygons + .push(Polygon::new(topology_ring(geometry))); + } +} + +impl TopologyBuilder for TopologyBox +where + G: BoxTrait, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + let minimum = [ + geometry.get_indexed::<{ corner::MIN }, 0>().into(), + geometry.get_indexed::<{ corner::MIN }, 1>().into(), + ]; + let maximum = [ + geometry.get_indexed::<{ corner::MAX }, 0>().into(), + geometry.get_indexed::<{ corner::MAX }, 1>().into(), + ]; + topology + .polygons + .push(Polygon::new(Ring::from_vec(alloc::vec![ + topology_point(minimum), + topology_point([minimum[0], maximum[1]]), + topology_point(maximum), + topology_point([maximum[0], minimum[1]]), + topology_point(minimum), + ]))); + } +} + +impl TopologyBuilder for TopologyMultiPoint +where + G: MultiPoint, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + topology.points.extend(geometry.points().map(xy)); + } +} + +impl TopologyBuilder for TopologyMultiLinestring +where + G: MultiLinestring, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + for line in geometry.linestrings() { + append_topology_line(line.points().map(xy).collect(), topology); + } + } +} + +impl TopologyBuilder for TopologyMultiPolygon +where + G: MultiPolygon, + ::Scalar: Into, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + for polygon in geometry.polygons() { + append_topology_polygon(polygon, topology); + } + } +} + +impl TopologyBuilder> for TopologyDynamic +where + Scalar: CoordinateScalar + Into, + Cs: CoordinateSystem, + Cs::Family: SameAs, +{ + fn append(&self, geometry: &DynGeometry, topology: &mut Topology) { + append_dynamic_topology(geometry, topology); + } +} + +impl TopologyBuilder for TopologyCollection +where + G: GeometryCollection, + G::Item: Geometry, + ::Kind: TopologyBuilderForKind, + TopologyBuilderStrategy: TopologyBuilder + Default, +{ + fn append(&self, geometry: &G, topology: &mut Topology) { + for item in geometry.items() { + TopologyBuilderStrategy::::default().append(item, topology); + } + } +} + +impl RelateStrategy for RelateTopology +where + A: Geometry, + B: Geometry, + A::Kind: TopologyBuilderForKind, + B::Kind: TopologyBuilderForKind, + TopologyBuilderStrategy: TopologyBuilder + Default, + TopologyBuilderStrategy: TopologyBuilder + Default, +{ + fn relate(&self, first: &A, second: &B) -> Result { + let mut first_topology = Topology::default(); + TopologyBuilderStrategy::::default().append(first, &mut first_topology); + let mut second_topology = Topology::default(); + TopologyBuilderStrategy::::default().append(second, &mut second_topology); + relate_topologies(&first_topology, &second_topology) + } +} + #[derive(Debug, Clone, Copy, PartialEq, Eq)] enum Location { Interior, @@ -845,6 +1221,350 @@ fn interpolate(first: [f64; 2], second: [f64; 2], fraction: f64) -> [f64; 2] { ] } +fn topology_point(coordinates: [f64; 2]) -> TopologyPointModel { + TopologyPointModel::new(coordinates[0], coordinates[1]) +} + +fn topology_ring(ring: &R) -> Ring +where + R: RingTrait, + ::Scalar: Into, +{ + Ring::from_vec( + ring.points() + .map(|point| topology_point(xy(point))) + .collect(), + ) +} + +fn append_topology_line(mut points: Vec<[f64; 2]>, topology: &mut Topology) { + points.dedup_by(|first, second| xy_equal(*first, *second)); + if points.len() >= 2 { + topology.lines.push(points); + } else if let Some(point) = points.first() { + topology.points.push(*point); + } +} + +fn append_topology_polygon(polygon: &G, topology: &mut Topology) +where + G: PolygonTrait, + ::Scalar: Into, +{ + let outer = topology_ring(polygon.exterior()); + if outer.0.len() < 3 { + append_topology_line(outer.0.iter().map(xy).collect(), topology); + return; + } + topology.polygons.push(Polygon::with_inners( + outer, + polygon.interiors().map(topology_ring).collect(), + )); +} + +fn append_dynamic_topology(geometry: &DynGeometry, topology: &mut Topology) +where + Scalar: CoordinateScalar + Into, + Cs: CoordinateSystem, + Cs::Family: SameAs, +{ + match geometry { + DynGeometry::Point(point) => topology.points.push(xy(point)), + DynGeometry::LineString(line) => { + append_topology_line(line.points().map(xy).collect(), topology); + } + DynGeometry::Polygon(polygon) => append_topology_polygon(polygon, topology), + DynGeometry::MultiPoint(points) => topology.points.extend(points.points().map(xy)), + DynGeometry::MultiLineString(lines) => { + for line in lines.linestrings() { + append_topology_line(line.points().map(xy).collect(), topology); + } + } + DynGeometry::MultiPolygon(polygons) => { + for polygon in polygons.polygons() { + append_topology_polygon(polygon, topology); + } + } + DynGeometry::GeometryCollection(items) => { + for item in items { + append_dynamic_topology(item, topology); + } + } + } +} + +fn topology_segments(topology: &Topology) -> Vec<([f64; 2], [f64; 2])> { + let mut segments = Vec::new(); + for line in &topology.lines { + for points in line.windows(2) { + if !xy_equal(points[0], points[1]) { + segments.push((points[0], points[1])); + } + } + } + for polygon in &topology.polygons { + append_boundary_segments(&polygon.outer, &mut segments); + for ring in &polygon.inners { + append_boundary_segments(ring, &mut segments); + } + } + segments +} + +fn topology_location(topology: &Topology, point: [f64; 2]) -> Location { + let mut polygon_boundary = false; + for polygon in &topology.polygons { + match xy_location_polygon(point, polygon) { + Location::Interior => return Location::Interior, + Location::Boundary => polygon_boundary = true, + Location::Exterior => {} + } + } + + let mut on_line = false; + let mut endpoint_count = 0usize; + for line in &topology.lines { + for segment in line.windows(2) { + if point_on_segment(point, segment[0], segment[1]) { + on_line = true; + } + } + if let (Some(first), Some(last)) = (line.first(), line.last()) + && !xy_equal(*first, *last) + { + endpoint_count += usize::from(xy_equal(point, *first)); + endpoint_count += usize::from(xy_equal(point, *last)); + } + } + if on_line { + return if endpoint_count % 2 == 1 && !polygon_boundary { + Location::Boundary + } else { + Location::Interior + }; + } + if topology + .points + .iter() + .any(|candidate| xy_equal(*candidate, point)) + { + return Location::Interior; + } + if polygon_boundary { + Location::Boundary + } else { + Location::Exterior + } +} + +fn dimension_rank(dimension: Dimension) -> u8 { + match dimension { + Dimension::Empty => 0, + Dimension::Point => 1, + Dimension::Curve => 2, + Dimension::Area => 3, + } +} + +fn set_dimension(matrix: &mut De9im, row: Location, column: Location, dimension: Dimension) { + let cell = &mut matrix.m[row.index()][column.index()]; + if dimension_rank(dimension) > dimension_rank(*cell) { + *cell = dimension; + } +} + +fn segment_parameter(point: [f64; 2], start: [f64; 2], end: [f64; 2]) -> f64 { + let dx = end[0] - start[0]; + let dy = end[1] - start[1]; + if dx.abs() >= dy.abs() && dx != 0.0 { + (point[0] - start[0]) / dx + } else if dy != 0.0 { + (point[1] - start[1]) / dy + } else { + 0.0 + } +} + +fn segment_parameters( + segment: ([f64; 2], [f64; 2]), + all_segments: &[([f64; 2], [f64; 2])], +) -> Vec { + let mut parameters = alloc::vec![0.0, 1.0]; + for &(start, end) in all_segments { + match segment_relation(segment.0, segment.1, start, end) { + SegmentRelation::Point(point) => { + parameters.push(segment_parameter(point, segment.0, segment.1)); + } + SegmentRelation::Overlap => { + for point in [start, end] { + if point_on_segment(point, segment.0, segment.1) { + parameters.push(segment_parameter(point, segment.0, segment.1)); + } + } + } + SegmentRelation::Disjoint => {} + } + } + parameters.retain(|parameter| (-f64::EPSILON..=1.0 + f64::EPSILON).contains(parameter)); + parameters.sort_by(f64::total_cmp); + parameters.dedup_by(|first, second| (*first - *second).abs() <= f64::EPSILON); + parameters +} + +fn record_segment_cells( + matrix: &mut De9im, + first: &Topology, + second: &Topology, + segment: ([f64; 2], [f64; 2]), + all_segments: &[([f64; 2], [f64; 2])], +) { + for interval in segment_parameters(segment, all_segments).windows(2) { + if interval[1] - interval[0] <= f64::EPSILON { + continue; + } + let midpoint = interpolate(segment.0, segment.1, (interval[0] + interval[1]) * 0.5); + let first_location = topology_location(first, midpoint); + let second_location = topology_location(second, midpoint); + if first_location != Location::Exterior { + set_dimension(matrix, first_location, second_location, Dimension::Curve); + } + } +} + +fn append_topology_candidates( + topology: &Topology, + segments: &[([f64; 2], [f64; 2])], + output: &mut Vec<[f64; 2]>, +) { + output.extend(topology.points.iter().copied()); + for &(start, end) in segments { + output.push(start); + output.push(end); + } +} + +fn areas_intersect(first: &Topology, second: &Topology) -> Result { + for first_polygon in &first.polygons { + for second_polygon in &second.polygons { + if !crate::operation::intersection(first_polygon, second_polygon)? + .0 + .is_empty() + { + return Ok(true); + } + } + } + Ok(false) +} + +fn has_area_outside(first: &Topology, second: &Topology) -> Result { + for polygon in &first.polygons { + let mut pieces = alloc::vec![polygon.clone()]; + for clip in &second.polygons { + let mut remainder = Vec::new(); + for piece in pieces { + remainder.extend(crate::operation::difference(&piece, clip)?.0); + } + pieces = remainder; + if pieces.is_empty() { + break; + } + } + if !pieces.is_empty() { + return Ok(true); + } + } + Ok(false) +} + +fn relate_topologies(first: &Topology, second: &Topology) -> Result { + if !first.in_range() || !second.in_range() { + return Err(OverlayError::Unsupported); + } + + let first_segments = topology_segments(first); + let second_segments = topology_segments(second); + let all_segments = first_segments + .iter() + .chain(&second_segments) + .copied() + .collect::>(); + let mut matrix = empty_matrix(); + + if areas_intersect(first, second)? { + matrix.m[feature::INTERIOR][feature::INTERIOR] = Dimension::Area; + } + if has_area_outside(first, second)? { + matrix.m[feature::INTERIOR][feature::EXTERIOR] = Dimension::Area; + } + if has_area_outside(second, first)? { + matrix.m[feature::EXTERIOR][feature::INTERIOR] = Dimension::Area; + } + + for &segment in &first_segments { + record_segment_cells(&mut matrix, first, second, segment, &all_segments); + } + for &segment in &second_segments { + for interval in segment_parameters(segment, &all_segments).windows(2) { + if interval[1] - interval[0] <= f64::EPSILON { + continue; + } + let midpoint = interpolate(segment.0, segment.1, (interval[0] + interval[1]) * 0.5); + let first_location = topology_location(first, midpoint); + let second_location = topology_location(second, midpoint); + if second_location != Location::Exterior { + set_dimension( + &mut matrix, + first_location, + second_location, + Dimension::Curve, + ); + } + } + } + + let mut candidates = Vec::new(); + append_topology_candidates(first, &first_segments, &mut candidates); + append_topology_candidates(second, &second_segments, &mut candidates); + for &(first_start, first_end) in &first_segments { + for &(second_start, second_end) in &second_segments { + match segment_relation(first_start, first_end, second_start, second_end) { + SegmentRelation::Point(point) => candidates.push(point), + SegmentRelation::Overlap => { + for point in [first_start, first_end, second_start, second_end] { + if point_on_segment(point, first_start, first_end) + && point_on_segment(point, second_start, second_end) + { + candidates.push(point); + } + } + } + SegmentRelation::Disjoint => {} + } + } + } + candidates.sort_by(|first, second| { + first[0] + .total_cmp(&second[0]) + .then_with(|| first[1].total_cmp(&second[1])) + }); + candidates.dedup_by(|first, second| xy_equal(*first, *second)); + for point in candidates { + let first_location = topology_location(first, point); + let second_location = topology_location(second, point); + if first_location != Location::Exterior || second_location != Location::Exterior { + set_dimension( + &mut matrix, + first_location, + second_location, + Dimension::Point, + ); + } + } + + Ok(matrix) +} + /// Compute the DE-9IM matrix relating two polygons. /// /// Fills the matrix from Boolean interior regions and exact segment-pair diff --git a/crates/geometry-strategy/src/buffer.rs b/crates/geometry-strategy/src/buffer.rs index 0806a8a..939c16e 100644 --- a/crates/geometry-strategy/src/buffer.rs +++ b/crates/geometry-strategy/src/buffer.rs @@ -5,6 +5,102 @@ //! data-only and live below the overlay engine so algorithms can consume them //! without creating a dependency cycle. +use geometry_cs::{CartesianFamily, CoordinateSystem, GeographicFamily, SphericalFamily, Spheroid}; +use geometry_trait::Geometry; + +/// Cartesian coordinate strategy bundle for buffer construction. +/// +/// Mirrors `boost::geometry::strategies::buffer::cartesian<>` from +/// `strategies/buffer/cartesian.hpp:24-31`. +#[derive(Debug, Default, Clone, Copy, PartialEq, Eq)] +pub struct CartesianBuffer; + +/// Spherical coordinate strategy bundle for buffer construction. +/// +/// Mirrors `boost::geometry::strategies::buffer::spherical` from +/// `strategies/buffer/spherical.hpp:24-46`. Distances use the same unit as +/// `radius`. The overlay consumer applies this bundle through a local tangent +/// projection, an explicitly recorded approximation for local buffers. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SphericalBuffer { + /// Sphere radius in the buffer distance unit. + pub radius: f64, +} + +impl SphericalBuffer { + /// Unit-sphere strategy matching Boost's default radius. + pub const UNIT: Self = Self { radius: 1.0 }; + + /// Construct a spherical buffer strategy with an explicit radius. + #[must_use] + pub const fn new(radius: f64) -> Self { + Self { radius } + } +} + +impl Default for SphericalBuffer { + fn default() -> Self { + Self::UNIT + } +} + +/// Geographic coordinate strategy bundle for buffer construction. +/// +/// Mirrors `boost::geometry::strategies::buffer::geographic` from +/// `strategies/buffer/geographic.hpp:25-48`. The overlay consumer derives the +/// local meridional and prime-vertical scales from this spheroid before using +/// its Cartesian offset engine. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct GeographicBuffer { + /// Reference spheroid used to convert angular coordinates and metric + /// buffer distances. + pub spheroid: Spheroid, +} + +impl GeographicBuffer { + /// WGS84 geographic buffer strategy. + pub const WGS84: Self = Self { + spheroid: Spheroid::WGS84, + }; + + /// Construct a geographic buffer strategy with an explicit spheroid. + #[must_use] + pub const fn new(spheroid: Spheroid) -> Self { + Self { spheroid } + } +} + +impl Default for GeographicBuffer { + fn default() -> Self { + Self::WGS84 + } +} + +/// Select the default coordinate strategy bundle for a buffer input family. +/// +/// Mirrors `strategies::buffer::services::default_strategy` from +/// `strategies/buffer/services.hpp:24-40` and its Cartesian, spherical, and +/// geographic specializations. +pub trait DefaultBuffer { + /// Default coordinate strategy for this family. + type Strategy: Default; +} + +impl DefaultBuffer for CartesianFamily { + type Strategy = CartesianBuffer; +} + +impl DefaultBuffer for SphericalFamily { + type Strategy = SphericalBuffer; +} + +impl DefaultBuffer for GeographicFamily { + type Strategy = GeographicBuffer; +} + +/// Coordinate strategy selected by the point coordinate-system family of `G`. +pub type DefaultBufferStrategy = <<<::Point as geometry_trait::Point>::Cs as CoordinateSystem>::Family as DefaultBuffer<<<::Point as geometry_trait::Point>::Cs as CoordinateSystem>::Family>>::Strategy; + /// Signed offset distance policy. /// /// Mirrors `strategy::buffer::distance_symmetric` and `distance_asymmetric` diff --git a/crates/geometry-strategy/src/lib.rs b/crates/geometry-strategy/src/lib.rs index 434628c..72e522c 100644 --- a/crates/geometry-strategy/src/lib.rs +++ b/crates/geometry-strategy/src/lib.rs @@ -131,7 +131,8 @@ pub use area::{ pub use azimuth::{AzimuthStrategy, CartesianAzimuth, DefaultAzimuth, DefaultAzimuthStrategy}; pub use buffer::{ BufferDistanceStrategy, BufferEndStrategy, BufferJoinStrategy, BufferPointStrategy, - BufferSettings, BufferSideStrategy, + BufferSettings, BufferSideStrategy, CartesianBuffer, DefaultBuffer, DefaultBufferStrategy, + GeographicBuffer, SphericalBuffer, }; pub use cartesian::{ComparablePythagoras, PointToSegment, Pythagoras}; pub use centroid::{ diff --git a/crates/geometry/tests/buffer_strategy_parity.rs b/crates/geometry/tests/buffer_strategy_parity.rs index 1bff0a0..e8853e6 100644 --- a/crates/geometry/tests/buffer_strategy_parity.rs +++ b/crates/geometry/tests/buffer_strategy_parity.rs @@ -1,18 +1,20 @@ //! Public-facade tests for Boost's composable buffer strategy family. -use boost_geometry::cs::Cartesian; +use boost_geometry::cs::{Cartesian, Degree, Geographic, Spherical, Spheroid}; use boost_geometry::model::{ Box as ModelBox, Linestring, MultiLinestring, MultiPoint, MultiPolygon, Point2D, Polygon, Ring, Segment, polygon, }; use boost_geometry::overlay::{ JoinStrategy, OverlayError, PointStrategy, buffer, buffer_convex_polygon, buffer_with, + buffer_with_strategy, }; -use boost_geometry::prelude::area; +use boost_geometry::prelude::{area, distance_with}; use boost_geometry::strategy::buffer::{ BufferDistanceStrategy, BufferEndStrategy, BufferJoinStrategy, BufferPointStrategy, - BufferSettings, BufferSideStrategy, + BufferSettings, BufferSideStrategy, GeographicBuffer, SphericalBuffer, }; +use boost_geometry::strategy::{Haversine, Vincenty}; use boost_geometry::trait_::{MultiPolygon as _, Polygon as _, Ring as _}; type P = Point2D; @@ -463,3 +465,213 @@ fn areal_offset_handles_collinear_duplicate_and_collapsed_boundaries() { assert_eq!(collapsed_result.0.len(), 1); assert!((buffered_area(&collapsed_result) - 1.0).abs() < 1e-12); } + +/// `test/algorithms/buffer/buffer_point_geo.cpp:34-49` — the default +/// geographic coordinate strategy interprets buffer distance in metres and +/// constructs a geodesic point circle through the public facade. +#[test] +fn geographic_point_buffer_uses_wgs84_by_default() { + type GeographicPoint = Point2D>; + + let center = GeographicPoint::new(4.9, 52.0); + let result = buffer_with(¢er, BufferSettings::round(10.0, 360)).unwrap(); + let polygon = result.polygons().next().unwrap(); + assert_eq!(polygon.exterior().points().count(), 361); + for point in polygon.exterior().points().take(360) { + let distance = distance_with(¢er, point, Vincenty::WGS84); + assert!((distance - 10.0).abs() < 0.05); + } + let observed_area = area(polygon).abs(); + assert!((observed_area - 314.15).abs() < 314.15 * 0.005); +} + +/// `strategies/buffer/spherical.hpp:24-58` — an explicit sphere radius is +/// carried by the spherical strategy bundle. The great-circle distance of +/// each generated vertex is the self-contained oracle because Boost has no +/// spherical buffer-algorithm fixture. +#[test] +fn spherical_point_buffer_honors_the_explicit_radius_strategy() { + type SphericalPoint = Point2D>; + + let radius = 6_371_008.8; + let center = SphericalPoint::new(-113.49, 53.54); + let result = buffer_with_strategy( + ¢er, + BufferSettings::round(1_000.0, 72), + SphericalBuffer::new(radius), + ) + .unwrap(); + let polygon = result.polygons().next().unwrap(); + for point in polygon.exterior().points().take(72) { + let distance = distance_with(¢er, point, Haversine { radius }); + assert!((distance - 1_000.0).abs() < 0.5); + } +} + +/// `test/algorithms/buffer/buffer_geo_spheroid.cpp:107-121` — a caller can +/// replace WGS84 with the alternate spheroid used by Boost's oracle fixture. +#[test] +fn geographic_point_buffer_accepts_an_explicit_spheroid() { + type GeographicPoint = Point2D>; + + let spheroid = Spheroid { + equatorial_radius: 6_378_000.0, + flattening: (6_378_000.0 - 6_375_000.0) / 6_378_000.0, + }; + let center = GeographicPoint::new(10.393_775_9, 63.430_232_3); + let result = buffer_with_strategy( + ¢er, + BufferSettings::round(100.0, 360), + GeographicBuffer::new(spheroid), + ) + .unwrap(); + let polygon = result.polygons().next().unwrap(); + let distance_strategy = Vincenty { + spheroid, + max_iterations: 1_000, + tolerance: 1e-12, + }; + for point in polygon.exterior().points().take(360) { + let distance = distance_with(¢er, point, distance_strategy); + assert!((distance - 100.0).abs() < 0.5); + } + let observed_area = area(polygon).abs(); + assert!((observed_area - 31_414.33).abs() < 31_414.33 * 0.005); +} + +/// `test/algorithms/buffer/buffer_linestring_geo.cpp:15-64` and +/// `buffer_polygon_geo.cpp:15-55` — geographic linear and areal inputs use +/// the same five public strategy roles as Cartesian inputs. +#[test] +fn geographic_linear_and_areal_buffers_use_public_strategy_roles() { + type GeographicPoint = Point2D>; + + let line = Linestring::from_vec(vec![ + GeographicPoint::new(10.396_562_8, 63.427_678_6), + GeographicPoint::new(10.395_313_4, 63.429_963_4), + ]); + let line_settings = BufferSettings { + end: BufferEndStrategy::Flat, + ..BufferSettings::round(5.0, 360) + }; + let line_result = buffer_with(&line, line_settings).unwrap(); + let line_area: f64 = line_result + .polygons() + .map(|polygon| area(polygon).abs()) + .sum(); + assert!((line_area - 2_622.0).abs() < 35.0); + + let polygon: Polygon = Polygon::new(Ring::from_vec(vec![ + GeographicPoint::new(10.400_658_7, 63.437_798_2), + GeographicPoint::new(10.405_090_4, 63.439_599_3), + GeographicPoint::new(10.407_499_4, 63.438_252_7), + GeographicPoint::new(10.400_658_7, 63.437_798_2), + ])); + let polygon_result = buffer_with(&polygon, BufferSettings::round(5.0, 36)).unwrap(); + let polygon_area: f64 = polygon_result + .polygons() + .map(|polygon| area(polygon).abs()) + .sum(); + assert!((polygon_area - 32_940.0).abs() < 600.0); +} + +/// `test/algorithms/buffer/buffer_multi_linestring_geo.cpp:18-73` and +/// `buffer_multi_polygon_geo.cpp:59-122` — every static geometry-kind arm is +/// available with an angular coordinate strategy, not only point/polygon. +#[test] +fn angular_segment_ring_box_and_multi_dispatch_is_public() { + type GeographicPoint = Point2D>; + type SphericalPoint = Point2D>; + let spherical = SphericalBuffer::new(6_371_008.8); + let round = BufferSettings::round(100.0, 36); + + let segment = Segment::new( + SphericalPoint::new(-113.50, 53.54), + SphericalPoint::new(-113.49, 53.54), + ); + assert!( + !buffer_with_strategy(&segment, round, spherical) + .unwrap() + .0 + .is_empty() + ); + + let ring: Ring = Ring::from_vec(vec![ + SphericalPoint::new(-113.50, 53.53), + SphericalPoint::new(-113.50, 53.54), + SphericalPoint::new(-113.49, 53.54), + SphericalPoint::new(-113.49, 53.53), + SphericalPoint::new(-113.50, 53.53), + ]); + assert!( + !buffer_with_strategy(&ring, round, spherical) + .unwrap() + .0 + .is_empty() + ); + + let bounds = ModelBox::from_corners( + SphericalPoint::new(-113.50, 53.53), + SphericalPoint::new(-113.49, 53.54), + ); + assert!( + !buffer_with_strategy(&bounds, round, spherical) + .unwrap() + .0 + .is_empty() + ); + + let points = MultiPoint::from_vec(vec![ + SphericalPoint::new(-113.50, 53.54), + SphericalPoint::new(-113.48, 53.54), + ]); + assert_eq!( + buffer_with_strategy(&points, round, spherical) + .unwrap() + .polygons() + .count(), + 2 + ); + + let lines = MultiLinestring::from_vec(vec![ + Linestring::from_vec(vec![ + GeographicPoint::new(10.396, 63.427), + GeographicPoint::new(10.399, 63.428), + ]), + Linestring::from_vec(vec![ + GeographicPoint::new(10.406, 63.427), + GeographicPoint::new(10.409, 63.428), + ]), + ]); + assert_eq!( + buffer_with(&lines, BufferSettings::round(5.0, 36)) + .unwrap() + .polygons() + .count(), + 2 + ); + + let polygons: MultiPolygon> = MultiPolygon::from_vec(vec![ + Polygon::new(Ring::from_vec(vec![ + GeographicPoint::new(10.396, 63.427), + GeographicPoint::new(10.396, 63.428), + GeographicPoint::new(10.397, 63.428), + GeographicPoint::new(10.397, 63.427), + GeographicPoint::new(10.396, 63.427), + ])), + Polygon::new(Ring::from_vec(vec![ + GeographicPoint::new(10.406, 63.427), + GeographicPoint::new(10.406, 63.428), + GeographicPoint::new(10.407, 63.428), + GeographicPoint::new(10.407, 63.427), + GeographicPoint::new(10.406, 63.427), + ])), + ]); + assert_eq!( + buffer_with(&polygons, BufferSettings::round(5.0, 36)) + .unwrap() + .polygons() + .count(), + 2 + ); +} diff --git a/crates/geometry/tests/relate_pair_parity.rs b/crates/geometry/tests/relate_pair_parity.rs index 01e2ea5..0c00ac8 100644 --- a/crates/geometry/tests/relate_pair_parity.rs +++ b/crates/geometry/tests/relate_pair_parity.rs @@ -1,10 +1,14 @@ //! Public-facade DE-9IM tests across pointlike, linear, and areal pairs. use boost_geometry::cs::Cartesian; -use boost_geometry::model::{Linestring, Point2D, Polygon, Ring, polygon}; +use boost_geometry::model::{ + Box as ModelBox, DynGeometry, DynGeometryCollection, Linestring, MultiLinestring, MultiPoint, + MultiPolygon, Point2D, Polygon, Ring, Segment, polygon, +}; use boost_geometry::overlay::{ De9im, Dimension, OverlayError, RelateError, crosses, overlaps, relate, relation, touches, }; +use boost_geometry::trait_::Polygon as _; type P = Point2D; @@ -286,3 +290,151 @@ fn public_matrix_masks_cover_every_symbol_and_overlay_errors() { Err(RelateError::Overlay(OverlayError::Unsupported)) ); } + +/// `test/algorithms/overlaps/overlaps_box.cpp:21-35` and +/// `test/algorithms/relate/relate_pointlike_geometry.cpp:166-174` — boxes +/// participate in the same areal matrix dispatch as polygons. +#[test] +fn boxes_and_rings_use_public_areal_relate_dispatch() { + let first = ModelBox::from_corners(P::new(1.0, 1.0), P::new(3.0, 3.0)); + let second = ModelBox::from_corners(P::new(0.0, 0.0), P::new(2.0, 2.0)); + assert!( + relation(&first, &second) + .unwrap() + .matches("212101212") + .unwrap() + ); + assert!(overlaps(&first, &second).unwrap()); + + let ring = square().exterior().clone(); + assert!( + relation(&P::new(2.0, 2.0), &ring) + .unwrap() + .matches("0FFFFF212") + .unwrap() + ); + assert!(touches(&P::new(0.0, 2.0), &ring).unwrap()); +} + +/// `test/algorithms/relate/relate_pointlike_geometry.cpp:31-45,87-93` — +/// multi-point membership and the mod-2 boundary of a multi-linestring. +#[test] +fn pointlike_and_linear_multis_preserve_union_boundary_rules() { + let first = MultiPoint::from_vec(vec![P::new(0.0, 0.0), P::new(1.0, 1.0)]); + let second = MultiPoint::from_vec(vec![P::new(0.0, 0.0), P::new(1.0, 0.0)]); + assert!( + relation(&first, &second) + .unwrap() + .matches("0F0FFF0F2") + .unwrap() + ); + + let lines = MultiLinestring::from_vec(vec![ + Linestring::from_vec(vec![P::new(0.0, 0.0), P::new(2.0, 0.0), P::new(2.0, 2.0)]), + Linestring::from_vec(vec![P::new(0.0, 0.0), P::new(0.0, 2.0)]), + ]); + assert!( + relation(&P::new(0.0, 0.0), &lines) + .unwrap() + .matches("0FFFFF102") + .unwrap() + ); +} + +/// `test/algorithms/relate/relate_pointlike_geometry.cpp:181-219` — polygon +/// members are related as one multi-polygon point set, including a shared +/// boundary vertex. +#[test] +fn multipolygon_relate_uses_the_public_union_topology() { + let polygons = + MultiPolygon::from_vec(vec![box_at(0.0, 0.0, 5.0, 5.0), box_at(5.0, 5.0, 9.0, 9.0)]); + assert!( + relation(&P::new(5.0, 5.0), &polygons) + .unwrap() + .matches("F0FFFF212") + .unwrap() + ); + assert!( + relation(&P::new(6.0, 6.0), &polygons) + .unwrap() + .matches("0FFFFF212") + .unwrap() + ); +} + +/// `test/algorithms/relate/relate_gc.cpp:55-65` — heterogeneous collections +/// use the OGC union topology rather than treating members as independent +/// matrices. +#[test] +fn geometry_collections_relate_through_runtime_public_dispatch() { + let joined_lines = DynGeometryCollection(vec![ + DynGeometry::LineString(Linestring::from_vec(vec![ + P::new(0.0, 0.0), + P::new(1.0, 1.0), + ])), + DynGeometry::LineString(Linestring::from_vec(vec![ + P::new(1.0, 1.0), + P::new(2.0, 2.0), + ])), + ]); + let point = DynGeometryCollection(vec![DynGeometry::Point(P::new(1.0, 1.0))]); + assert!( + relation(&point, &joined_lines) + .unwrap() + .matches("0FFFFF102") + .unwrap() + ); + + let first = DynGeometryCollection(vec![ + DynGeometry::Polygon(box_at(0.0, 0.0, 5.0, 5.0)), + DynGeometry::LineString(Linestring::from_vec(vec![ + P::new(1.0, 1.0), + P::new(6.0, 6.0), + ])), + ]); + let second = DynGeometryCollection(vec![ + DynGeometry::Polygon(box_at(0.0, 0.0, 5.0, 5.0)), + DynGeometry::LineString(Linestring::from_vec(vec![ + P::new(5.0, 5.0), + P::new(6.0, 6.0), + ])), + ]); + let matrix = relation(&first, &second).unwrap(); + assert!(matrix.matches("2FFF1FFF2").unwrap()); +} + +/// `test/algorithms/relate/relate_linear_areal.cpp:44-64` and +/// `relate_gc.cpp:107-111` — segment, runtime-variant, and static-to-collection +/// ordered pairs all use the same public matrix contract. +#[test] +fn segment_dynamic_and_collection_reverse_pairs_are_public() { + let segment = Segment::new(P::new(-1.0, 2.0), P::new(5.0, 2.0)); + let bounds = ModelBox::from_corners(P::new(0.0, 0.0), P::new(4.0, 4.0)); + let segment_box = relation(&segment, &bounds).unwrap(); + assert!(segment_box.matches("101FF0212").unwrap()); + assert_eq!( + relation(&bounds, &segment).unwrap(), + segment_box.transposed() + ); + assert!(crosses(&segment, &bounds).unwrap()); + + let dynamic_point = DynGeometry::Point(P::new(2.0, 2.0)); + let dynamic_polygon = DynGeometry::Polygon(square()); + let dynamic_matrix = relation(&dynamic_point, &dynamic_polygon).unwrap(); + assert!(dynamic_matrix.matches("0FFFFF212").unwrap()); + assert_eq!( + relation(&dynamic_polygon, &dynamic_point).unwrap(), + dynamic_matrix.transposed() + ); + + let adjacent = DynGeometryCollection(vec![ + DynGeometry::Polygon(box_at(10.0, 0.0, 20.0, 10.0)), + DynGeometry::Point(P::new(15.0, 5.0)), + ]); + assert!( + relation(&box_at(0.0, 0.0, 10.0, 10.0), &adjacent) + .unwrap() + .matches("FF2F11212") + .unwrap() + ); +} diff --git a/docs/03-overlay-engine.md b/docs/03-overlay-engine.md index d50d30b..648be95 100644 --- a/docs/03-overlay-engine.md +++ b/docs/03-overlay-engine.md @@ -185,8 +185,10 @@ Thin orchestration over the pipeline above: from {Interior, Boundary, Exterior} of the two geometries, the *dimension* of their intersection (`Dimension::{Empty,Point,Curve,Area}`). **`relate`** tests that matrix against a DE-9IM mask; `touches`, `overlaps`, - and `crosses` are thin predicates over the same matrix. Built from the turn - graph plus interior-point sampling; areal × areal only in v1. + and `crosses` are thin predicates over the same matrix. Cartesian dispatch + covers static single kinds, homogeneous multis, runtime geometries, and + heterogeneous geometry collections. Collection topology uses OGC union + semantics, including mod-2 multiline boundaries. * **`is_valid`** tag-dispatches to the ring/polygon validators (and validates each multi-polygon member). The validators check the OGC simple-feature rules: finite in-range coordinates, enough points, closed boundary, no spikes, no @@ -204,17 +206,15 @@ exactly this failure mode being caught. ### OVL7 — Buffer (`buffer.rs`) -The public `buffer` entry tag-dispatches by geometry kind and grows a geometry -outward by a fixed distance. Its point and convex-polygon kernels are also -available directly as `buffer_point` and `buffer_convex_polygon`. - -**v1 scope:** positive-distance buffers of a **point** (→ circle, via -`PointStrategy::{Circle,Square}`) and a **convex polygon** (→ grown polygon -with `JoinStrategy::{Round,Miter}` corners). The general non-convex / -negative-distance buffer needs the full offset-and-self-union machinery and -is deferred. Note: `JoinStrategy::Miter` is currently **uncapped** — Boost's -`miter_limit` (default 5× distance) isn't implemented yet, so a near-180° -corner produces a proportionally long spike. +The public `buffer` entry tag-dispatches every static single and homogeneous +multi kind and grows or erodes it using explicit distance, side, join, end, +and point roles. Cartesian offsets are native and include holes, non-convex +polygons, signed distances, asymmetric linear widths, capped miters, and +round/flat ends. Spherical and geographic inputs use family-selected radius +or spheroid bundles, project into a local tangent plane, reuse the Cartesian +engine, and transform back. That angular path is an intentional local-extent +approximation; the feature-parity assumptions identify global/polar accuracy +as the revisit trigger. ## The recurring design principle: refuse, don't guess diff --git a/docs/crates/geometry-overlay.md b/docs/crates/geometry-overlay.md index 11e5a46..9d5dda2 100644 --- a/docs/crates/geometry-overlay.md +++ b/docs/crates/geometry-overlay.md @@ -22,10 +22,10 @@ functions (`intersection`, `union`, `difference`, `sym_difference`), plus | `traverse` (+ `enrich`, `state`) | OVL3 | Weiler–Atherton-style ring traversal — walk the turn graph, assemble output rings | | `assemble` | OVL4 | Nest traversed rings into `Polygon`/`MultiPolygon` by containment | | `operation` (+ `boolean`, `areal`) | OVL5 | split-edge arrangement; `intersection`, `r#union` (`union_poly` compatibility name), `difference`, `sym_difference`, `OverlayError` | -| `relate` | OVL6 | point/linestring/polygon `relation` matrix, `relate` mask, `De9im`, `crosses`/`overlaps`/`touches`, `Dimension` | +| `relate` | OVL6 | Cartesian static/multi/runtime/collection `relation` matrix, `relate` mask, `De9im`, `crosses`/`overlaps`/`touches`, `Dimension` | | `validity` | OVL6 | ring/polygon/multi-polygon `is_valid`, inter-ring/member topology, `ValidityFailure` | | `surface_point` | — | `point_on_surface` — a representative interior point, used by `assemble` and `relate` | -| `buffer` | OVL7 | Cartesian single/multi dispatch, `buffer`/`buffer_with`, signed areal offsets, linear ends, strategy compatibility enums | +| `buffer` | OVL7 | Cartesian/spherical/geographic single/multi dispatch, `buffer`/`buffer_with`/`buffer_with_strategy`, signed areal offsets, linear ends, strategy bundles | | `merge` | — | `merge_elements`, `merge_polygons`, `merge_multipolygon` | ## Robustness policy @@ -38,11 +38,12 @@ out-of-range coordinates are **refused**, not silently miscomputed. Boolean operations cover Cartesian polygon × polygon inputs with holes, containment, colocated vertices, shared edges, and collinear overlaps. Relate -covers point, linestring, and polygon pairs. Buffer covers Cartesian point, -segment, linestring, ring, polygon, box, and homogeneous multi kinds. Invalid -self-intersecting overlay inputs, linear/pointlike set-operation output, -geometry-collection relate, and non-Cartesian buffer specializations remain -outside this scope. +covers Cartesian static singles, homogeneous multis, runtime geometries, and +heterogeneous geometry collections. Buffer covers point, segment, linestring, +ring, polygon, box, and homogeneous multi kinds in all three coordinate +families. Angular buffer dispatch uses the recorded local-tangent +approximation. Invalid self-intersecting overlay inputs and linear/pointlike +set-operation output remain outside this scope. ## Why this is a separate crate, not part of `geometry-algorithm` diff --git a/docs/crates/geometry-strategy.md b/docs/crates/geometry-strategy.md index a5c786a..88d57d8 100644 --- a/docs/crates/geometry-strategy.md +++ b/docs/crates/geometry-strategy.md @@ -36,7 +36,7 @@ kind and coordinate-system family — combine). | `envelope` | `EnvelopeStrategy`, `EnvelopeStrategyForKind` | `Envelope{Point,Segment,Linestring,Ring,Polygon,Box,MultiPoint,MultiLinestring,MultiPolygon}` — see [tag-dispatch pattern](../02-tag-dispatch-pattern.md) | | `within` | `WithinStrategy`, `WithinStrategyForKind` | `WithinRing`, `WithinPoly`, `WithinBox` | | `intersects` | `IntersectsStrategy`, `IntersectsPairStrategy` | `CartesianIntersects` | -| `buffer` | data strategy bundle consumed by overlay | `BufferDistanceStrategy`, `BufferSideStrategy`, `BufferJoinStrategy`, `BufferEndStrategy`, `BufferPointStrategy`, `BufferSettings` | +| `buffer` | coordinate-family default plus five-role data bundle consumed by overlay | `CartesianBuffer`, `SphericalBuffer`, `GeographicBuffer`, `DefaultBuffer`; `BufferDistanceStrategy`, `BufferSideStrategy`, `BufferJoinStrategy`, `BufferEndStrategy`, `BufferPointStrategy`, `BufferSettings` | | `disjoint` | `DisjointStrategy` | `CartesianDisjoint` | | `equals` | `EqualsStrategy`, `EqualsPairStrategy` | `EqPointPoint`, `EqSegmentSegment`, `EqPolygonPolygon` | | `centroid` | `CentroidStrategy`, `CentroidStrategyForKind` | `Cartesian{Polygon,Ring,Linestring,Segment,Box,MultiPoint}Centroid` |