diff --git a/geo-bool-ops-benches/benches/boolean_ops.rs b/geo-bool-ops-benches/benches/boolean_ops.rs index c491a2e6f..865aa1aa5 100644 --- a/geo-bool-ops-benches/benches/boolean_ops.rs +++ b/geo-bool-ops-benches/benches/boolean_ops.rs @@ -1,7 +1,10 @@ use std::f64::consts::PI; use criterion::{measurement::Measurement, *}; -use geo::algorithm::{BooleanOps, Intersects, Rotate}; +use geo::{ + algorithm::{BooleanOps, Rotate}, + Relate, +}; use geo_booleanop::boolean::BooleanOp as OtherBooleanOp; use rand::{thread_rng, Rng}; @@ -77,10 +80,10 @@ fn run_complex(c: &mut Criterion) { }, ); - group.bench_with_input(BenchmarkId::new("geo::intersects", steps), &(), |b, _| { + group.bench_with_input(BenchmarkId::new("geo::relate", steps), &(), |b, _| { b.iter_batched( polys.sampler(), - |&(ref poly, ref poly2, _, _)| poly.intersects(poly2), + |&(ref poly, ref poly2, _, _)| poly.relate(poly2).is_intersects(), BatchSize::SmallInput, ); }); diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 8ea66dad9..0ba06b8c0 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -2,6 +2,8 @@ ## Unreleased +* Add `clip` boolean op. to clip a 1-D geometry with a 2-D geometry. + * * Add `Within` trait to determine if Geometry A is completely within Geometry B * * Add `Contains` impl for all remaining geometry types. diff --git a/geo/src/algorithm/bool_ops/assembly.rs b/geo/src/algorithm/bool_ops/assembly.rs index 544566985..53c1d42e3 100644 --- a/geo/src/algorithm/bool_ops/assembly.rs +++ b/geo/src/algorithm/bool_ops/assembly.rs @@ -1,6 +1,6 @@ use std::{ cell::Cell, - collections::{BTreeMap, HashMap}, + collections::{BTreeMap, HashMap, VecDeque}, }; use crate::{ @@ -20,11 +20,11 @@ use super::op::compare_crossings; /// describe a bounded region, do not intersect in their interior, and are not /// degenerate (not a point). #[derive(Debug)] -pub struct Assembly { +pub struct RegionAssembly { segments: Vec>, } -impl Default for Assembly { +impl Default for RegionAssembly { fn default() -> Self { Self { segments: Default::default(), @@ -32,13 +32,13 @@ impl Default for Assembly { } } -impl Assembly { +impl RegionAssembly { pub fn add_edge(&mut self, edge: LineOrPoint) { debug_assert!(edge.is_line()); self.segments.push(edge.into()); } pub fn finish(self) -> MultiPolygon { - let mut iter = CrossingsIter::from_iter(self.segments.iter()); + let mut iter = CrossingsIter::new_simple(self.segments.iter()); let mut snakes = vec![]; while let Some(pt) = iter.next() { @@ -175,6 +175,58 @@ impl Assembly { } } +#[derive(Debug)] +pub struct LineAssembly { + segments: Vec>>, + end_points: BTreeMap<(usize, SweepPoint), (usize, bool)>, +} + +impl LineAssembly { + pub fn add_edge(&mut self, geom: LineOrPoint, geom_idx: usize) { + // Try to find a line-string with either end-point + if let Some((seg_idx, at_front)) = self.end_points.remove(&(geom_idx, geom.left())) { + if at_front { + self.segments[seg_idx].push_front(geom.right()); + } else { + self.segments[seg_idx].push_back(geom.right()); + } + self.end_points + .insert((geom_idx, geom.right()), (seg_idx, at_front)); + } else if let Some((seg_idx, at_front)) = self.end_points.remove(&(geom_idx, geom.right())) + { + if at_front { + self.segments[seg_idx].push_front(geom.left()); + } else { + self.segments[seg_idx].push_back(geom.left()); + } + self.end_points + .insert((geom_idx, geom.left()), (seg_idx, at_front)); + } else { + let idx = self.segments.len(); + self.segments + .push(VecDeque::from_iter([geom.left(), geom.right()])); + self.end_points.insert((geom_idx, geom.left()), (idx, true)); + self.end_points + .insert((geom_idx, geom.right()), (idx, false)); + } + } + pub fn finish(self) -> Vec> { + self.segments + .into_iter() + .map(|pts| LineString::from_iter(pts.into_iter().map(|pt| *pt))) + .collect() + } +} + +impl Default for LineAssembly { + fn default() -> Self { + Self { + segments: Default::default(), + end_points: Default::default(), + } + } +} + #[derive(Debug, Clone)] struct Ring { ls: LineString, diff --git a/geo/src/algorithm/bool_ops/mod.rs b/geo/src/algorithm/bool_ops/mod.rs index e228a844c..8a9fdbab1 100644 --- a/geo/src/algorithm/bool_ops/mod.rs +++ b/geo/src/algorithm/bool_ops/mod.rs @@ -1,4 +1,4 @@ -use geo_types::MultiPolygon; +use geo_types::{MultiLineString, MultiPolygon}; use crate::{CoordsIter, GeoFloat, GeoNum, Polygon}; @@ -6,7 +6,8 @@ use crate::{CoordsIter, GeoFloat, GeoNum, Polygon}; /// /// Boolean operations are set operations on geometries considered as a subset /// of the 2-D plane. The operations supported are: intersection, union, xor or -/// symmetric difference, and set-difference. +/// symmetric difference, and set-difference on pairs of 2-D geometries and +/// clipping a 1-D geometry with self. /// /// These operations are implemented on [`Polygon`] and the [`MultiPolygon`] /// geometries. @@ -37,6 +38,16 @@ pub trait BooleanOps: Sized { fn difference(&self, other: &Self) -> MultiPolygon { self.boolean_op(other, OpType::Difference) } + + /// Clip a 1-D geometry with self. + /// + /// Returns the set-theoeretic intersection of `self` and `ls` if `invert` + /// is false, and the difference (`ls - self`) otherwise. + fn clip( + &self, + ls: &MultiLineString, + invert: bool, + ) -> MultiLineString; } #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)] @@ -51,9 +62,24 @@ impl BooleanOps for Polygon { type Scalar = T; fn boolean_op(&self, other: &Self, op: OpType) -> MultiPolygon { - let mut bop = Op::new(op, self.coords_count() + other.coords_count()); - bop.add_polygon(self, true); - bop.add_polygon(other, false); + let spec = BoolOp::from(op); + let mut bop = Proc::new(spec, self.coords_count() + other.coords_count()); + bop.add_polygon(self, 0); + bop.add_polygon(other, 1); + bop.sweep() + } + + fn clip( + &self, + ls: &MultiLineString, + invert: bool, + ) -> MultiLineString { + let spec = ClipOp::new(invert); + let mut bop = Proc::new(spec, self.coords_count() + ls.coords_count()); + bop.add_polygon(self, 0); + ls.0.iter().enumerate().for_each(|(idx, l)| { + bop.add_line_string(l, idx + 1); + }); bop.sweep() } } @@ -61,9 +87,24 @@ impl BooleanOps for MultiPolygon { type Scalar = T; fn boolean_op(&self, other: &Self, op: OpType) -> MultiPolygon { - let mut bop = Op::new(op, self.coords_count() + other.coords_count()); - bop.add_multi_polygon(self, true); - bop.add_multi_polygon(other, false); + let spec = BoolOp::from(op); + let mut bop = Proc::new(spec, self.coords_count() + other.coords_count()); + bop.add_multi_polygon(self, 0); + bop.add_multi_polygon(other, 1); + bop.sweep() + } + + fn clip( + &self, + ls: &MultiLineString, + invert: bool, + ) -> MultiLineString { + let spec = ClipOp::new(invert); + let mut bop = Proc::new(spec, self.coords_count() + ls.coords_count()); + bop.add_multi_polygon(self, 0); + ls.0.iter().enumerate().for_each(|(idx, l)| { + bop.add_line_string(l, idx + 1); + }); bop.sweep() } } @@ -71,6 +112,9 @@ impl BooleanOps for MultiPolygon { mod op; use op::*; mod assembly; +use assembly::*; +mod spec; +use spec::*; #[cfg(test)] mod tests; diff --git a/geo/src/algorithm/bool_ops/op.rs b/geo/src/algorithm/bool_ops/op.rs index c2007f084..38c495537 100644 --- a/geo/src/algorithm/bool_ops/op.rs +++ b/geo/src/algorithm/bool_ops/op.rs @@ -1,66 +1,69 @@ use std::{cell::Cell, cmp::Ordering, fmt::Debug}; -use super::{assembly::Assembly, *}; +use super::{MultiPolygon, Spec}; use crate::{ sweep::{Cross, Crossing, CrossingsIter, LineOrPoint}, CoordsIter, GeoFloat as Float, LineString, Polygon, }; #[derive(Debug, Clone)] -pub struct Op { - ty: OpType, - edges: Vec>, +pub struct Proc> { + spec: S, + edges: Vec>, } -impl Op { - pub fn new(ty: OpType, capacity: usize) -> Self { - Op { - ty, +impl> Proc { + pub fn new(spec: S, capacity: usize) -> Self { + Proc { + spec, edges: Vec::with_capacity(capacity), } } - // is_first -> whether it is from first input or second input - pub(crate) fn add_multi_polygon(&mut self, mp: &MultiPolygon, is_first: bool) { - mp.0.iter().for_each(|p| self.add_polygon(p, is_first)); + // idx: whether it is from first input or second input + pub(crate) fn add_multi_polygon(&mut self, mp: &MultiPolygon, idx: usize) { + mp.0.iter().for_each(|p| self.add_polygon(p, idx)); } - // is_first -> whether it is from first input or second input - pub(crate) fn add_polygon(&mut self, poly: &Polygon, is_first: bool) { - self.add_closed_ring(poly.exterior(), is_first, false); + // idx: whether it is from first input or second input + pub(crate) fn add_polygon(&mut self, poly: &Polygon, idx: usize) { + self.add_closed_ring(poly.exterior(), idx, false); for hole in poly.interiors() { - self.add_closed_ring(hole, is_first, true); + self.add_closed_ring(hole, idx, true); } } - // is_first -> whether it is from first input or second input - // _is_hole is not used rn; remove it once we fully handle fp issues - fn add_closed_ring(&mut self, ring: &LineString, is_first: bool, _is_hole: bool) { - assert!(ring.is_closed()); - if ring.coords_count() <= 3 { - return; - } - for line in ring.lines() { + pub(crate) fn add_line_string(&mut self, ls: &LineString, idx: usize) { + for line in ls.lines() { let lp: LineOrPoint<_> = line.into(); if !lp.is_line() { continue; } debug!("processing: {lp:?}"); - - let region = Region::infinity(self.ty); + let region = self.spec.infinity(); self.edges.push(Edge { geom: lp, - is_first, + idx, _region: region.into(), _region_2: region.into(), }); } } - pub fn sweep(&self) -> MultiPolygon { + // idx: whether it is from first input or second input + // _is_hole is not used rn; remove it once we fully handle fp issues + fn add_closed_ring(&mut self, ring: &LineString, idx: usize, _is_hole: bool) { + assert!(ring.is_closed()); + if ring.coords_count() <= 3 { + return; + } + + self.add_line_string(ring, idx); + } + + pub fn sweep(mut self) -> S::Output { let mut iter = CrossingsIter::from_iter(self.edges.iter()); - let mut output = Assembly::default(); while let Some(pt) = iter.next() { trace!( @@ -87,7 +90,7 @@ impl Op { geom = c.line, ); } - next_region.as_mut().unwrap().cross(cross.is_first); + next_region = Some(self.spec.cross(next_region.unwrap(), cross.idx)); let has_overlap = (idx + 1) < iter.intersections().len() && compare_crossings(c, &iter.intersections()[idx + 1]) == Ordering::Equal; if !has_overlap { @@ -97,11 +100,8 @@ impl Op { geom = c.line, next_region = next_region.unwrap() ); - let next_is_ty = next_region.unwrap().is_ty(self.ty); - if prev_region.is_ty(self.ty) ^ next_is_ty { - trace!("\tfull_geom: {geom:?}", geom = c.cross.geom); - output.add_edge(c.line) - } + self.spec + .output([prev_region, next_region.unwrap()], c.line, c.cross.idx); next_region = None; } idx += 1; @@ -127,14 +127,14 @@ impl Op { let mut region = prev .as_ref() .map(|(g, c)| c.get_region(*g)) - .unwrap_or_else(|| Region::infinity(self.ty)); + .unwrap_or_else(|| self.spec.infinity()); trace!("bot region: {region:?}"); while idx < iter.intersections().len() { let mut c = &iter.intersections()[idx]; let mut jdx = idx; loop { - region.cross(c.cross.is_first); + region = self.spec.cross(region, c.cross.idx); let has_overlap = (idx + 1) < iter.intersections().len() && compare_crossings(c, &iter.intersections()[idx + 1]) == Ordering::Equal; if !has_overlap { @@ -157,60 +157,20 @@ impl Op { idx += 1; } } - - output.finish() - } -} - -#[derive(Clone, Copy)] -struct Region { - is_first: bool, - is_second: bool, -} -impl Debug for Region { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!( - f, - "[{f}{s}]", - f = if self.is_first { "A" } else { "" }, - s = if self.is_second { "B" } else { "" }, - ) - } -} - -impl Region { - fn infinity(ty: OpType) -> Self { - Region { - is_first: false, - is_second: matches!(ty, OpType::Difference), - } - } - fn cross(&mut self, first: bool) { - if first { - self.is_first = !self.is_first; - } else { - self.is_second = !self.is_second; - } - } - fn is_ty(&self, ty: OpType) -> bool { - match ty { - OpType::Intersection | OpType::Difference => self.is_first && self.is_second, - OpType::Union => self.is_first || self.is_second, - OpType::Xor => self.is_first ^ self.is_second, - } + self.spec.finish() } } #[derive(Clone)] -struct Edge { +struct Edge> { geom: LineOrPoint, - is_first: bool, - _region: Cell, - _region_2: Cell, + idx: usize, + _region: Cell, + _region_2: Cell, } -impl Edge { - fn get_region(&self, piece: LineOrPoint) -> Region { +impl> Edge { + fn get_region(&self, piece: LineOrPoint) -> S::Region { // Note: This is related to the ordering of intersection // with respect to the complete geometry. Due to // finite-precision errors, intersection points might lie @@ -232,7 +192,7 @@ impl Edge { self._region_2.get() } } - fn set_region(&self, region: Region, piece: LineOrPoint) { + fn set_region(&self, region: S::Region, piece: LineOrPoint) { if piece.left() < self.geom.right() { self._region.set(region); } else { @@ -242,7 +202,7 @@ impl Edge { } } -impl std::fmt::Debug for Edge { +impl> std::fmt::Debug for Edge { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { let line = self.geom.line(); f.debug_struct("Edge") @@ -253,13 +213,13 @@ impl std::fmt::Debug for Edge { line.start.x, line.start.y, line.end.x, line.end.y ), ) - .field("is_first", &self.is_first) + .field("idx", &self.idx) .field("region", &self._region) .finish() } } -impl Cross for Edge { +impl> Cross for Edge { type Scalar = T; fn line(&self) -> LineOrPoint { diff --git a/geo/src/algorithm/bool_ops/spec.rs b/geo/src/algorithm/bool_ops/spec.rs new file mode 100644 index 000000000..839800094 --- /dev/null +++ b/geo/src/algorithm/bool_ops/spec.rs @@ -0,0 +1,130 @@ +use geo_types::{MultiLineString, MultiPolygon}; +use std::fmt::Debug; + +use super::*; +use crate::{sweep::LineOrPoint, GeoFloat, OpType}; + +pub trait Spec { + type Region: Copy + Debug; + type Output; + + fn infinity(&self) -> Self::Region; + fn cross(&self, prev_region: Self::Region, idx: usize) -> Self::Region; + fn output(&mut self, regions: [Self::Region; 2], geom: LineOrPoint, idx: usize); + fn finish(self) -> Self::Output; +} + +pub struct BoolOp { + ty: OpType, + assembly: RegionAssembly, +} +impl From for BoolOp { + fn from(ty: OpType) -> Self { + Self { + ty, + assembly: RegionAssembly::default(), + } + } +} +impl Spec for BoolOp { + type Region = Region; + type Output = MultiPolygon; + + fn infinity(&self) -> Self::Region { + Region { + is_first: false, + is_second: matches!(self.ty, OpType::Difference), + } + } + + fn cross(&self, mut prev_region: Self::Region, idx: usize) -> Self::Region { + prev_region.cross(idx == 0); + prev_region + } + + fn output(&mut self, regions: [Self::Region; 2], geom: LineOrPoint, _idx: usize) { + if regions[0].is_ty(self.ty) ^ regions[1].is_ty(self.ty) { + self.assembly.add_edge(geom) + } + } + + fn finish(self) -> Self::Output { + self.assembly.finish() + } +} + +pub struct ClipOp { + invert: bool, + assembly: LineAssembly, +} + +impl ClipOp { + pub fn new(invert: bool) -> Self { + Self { + invert, + assembly: Default::default(), + } + } +} + +impl Spec for ClipOp { + type Region = Region; + type Output = MultiLineString; + + fn infinity(&self) -> Self::Region { + Region { + is_first: false, + is_second: false, + } + } + + fn cross(&self, mut prev_region: Self::Region, idx: usize) -> Self::Region { + if idx == 0 { + prev_region.cross(true); + } + prev_region + } + + fn output(&mut self, regions: [Self::Region; 2], geom: LineOrPoint, idx: usize) { + if idx > 0 && (regions[0].is_first && regions[1].is_first) != self.invert { + self.assembly.add_edge(geom, idx); + } + } + + fn finish(self) -> Self::Output { + MultiLineString::new(self.assembly.finish()) + } +} + +#[derive(Clone, Copy)] +pub struct Region { + is_first: bool, + is_second: bool, +} +impl std::fmt::Debug for Region { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!( + f, + "[{f}{s}]", + f = if self.is_first { "A" } else { "" }, + s = if self.is_second { "B" } else { "" }, + ) + } +} + +impl Region { + fn cross(&mut self, first: bool) { + if first { + self.is_first = !self.is_first; + } else { + self.is_second = !self.is_second; + } + } + fn is_ty(&self, ty: OpType) -> bool { + match ty { + OpType::Intersection | OpType::Difference => self.is_first && self.is_second, + OpType::Union => self.is_first || self.is_second, + OpType::Xor => self.is_first ^ self.is_second, + } + } +} diff --git a/geo/src/algorithm/bool_ops/tests.rs b/geo/src/algorithm/bool_ops/tests.rs index ab8dd3308..a614104b2 100644 --- a/geo/src/algorithm/bool_ops/tests.rs +++ b/geo/src/algorithm/bool_ops/tests.rs @@ -1,5 +1,6 @@ use crate::{MultiPolygon, Polygon}; +use geo_types::LineString; use log::{error, info}; use std::{ @@ -27,9 +28,9 @@ fn check_sweep(wkt1: &str, wkt2: &str, ty: OpType) -> Result> let poly2 = MultiPolygon::try_from_wkt_str(wkt2) .or_else(|_| Polygon::::try_from_wkt_str(wkt2).map(MultiPolygon::from)) .unwrap(); - let mut bop = Op::new(ty, 0); - bop.add_multi_polygon(&poly1, true); - bop.add_multi_polygon(&poly2, false); + let mut bop = Proc::new(BoolOp::from(ty), 0); + bop.add_multi_polygon(&poly1, 0); + bop.add_multi_polygon(&poly2, 1); let geom = bop.sweep(); @@ -37,13 +38,6 @@ fn check_sweep(wkt1: &str, wkt2: &str, ty: OpType) -> Result> info!("{wkt}", wkt = geom.to_wkt()); Ok(geom) - - // let polygons = assemble(rings); - // info!("got {n} output polygons", n = polygons.len()); - // for p in polygons.iter() { - // info!("\t{wkt}", wkt = p.to_wkt()); - // } - // Ok(MultiPolygon::new(polygons)) } #[test] @@ -161,10 +155,17 @@ fn test_issue_865() -> Result<()> { } #[test] -fn test_jts_adhoc() -> Result<()> { +fn test_clip_adhoc() -> Result<()> { let wkt1 = "POLYGON ((20 0, 20 160, 200 160, 200 0, 20 0))"; - let wkt2 = "POLYGON ((220 80, 0 80, 0 240, 220 240, 220 80), (100 80, 120 120, 80 120, 100 80))"; + let wkt2 = "LINESTRING (0 0, 100 100, 240 100)"; - check_sweep(wkt1, wkt2, OpType::Intersection)?; + let poly1 = MultiPolygon::::try_from_wkt_str(wkt1) + .or_else(|_| Polygon::::try_from_wkt_str(wkt1).map(MultiPolygon::from)) + .unwrap(); + let mls = MultiLineString::try_from_wkt_str(wkt2) + .or_else(|_| LineString::::try_from_wkt_str(wkt2).map(MultiLineString::from)) + .unwrap(); + let output = poly1.clip(&mls, true); + eprintln!("{wkt}", wkt = output.to_wkt()); Ok(()) } diff --git a/geo/src/algorithm/sweep/iter.rs b/geo/src/algorithm/sweep/iter.rs index 026ae90e4..43a19b3a8 100644 --- a/geo/src/algorithm/sweep/iter.rs +++ b/geo/src/algorithm/sweep/iter.rs @@ -102,6 +102,12 @@ impl CrossingsIter where C: Cross + Clone, { + /// Faster sweep when input goemetries are known to not intersect except at + /// end-points. + pub fn new_simple>(iter: I) -> Self { + Self::new_ex(iter, true) + } + /// Returns the segments that intersect the last point yielded by /// the iterator. pub fn intersections_mut(&mut self) -> &mut [Crossing] { @@ -115,24 +121,28 @@ where pub(crate) fn prev_active(&self, c: &Crossing) -> Option<(LineOrPoint, &C)> { self.sweep.prev_active(c).map(|s| (s.geom, &s.cross)) } -} -impl FromIterator for CrossingsIter -where - C: Cross + Clone, -{ - fn from_iter>(iter: T) -> Self { + fn new_ex>(iter: T, is_simple: bool) -> Self { let iter = iter.into_iter(); let size = { let (min_size, max_size) = iter.size_hint(); max_size.unwrap_or(min_size) }; - let sweep = Sweep::new(iter); + let sweep = Sweep::new(iter, is_simple); let segments = Vec::with_capacity(4 * size); Self { sweep, segments } } } +impl FromIterator for CrossingsIter +where + C: Cross + Clone, +{ + fn from_iter>(iter: T) -> Self { + Self::new_ex(iter, false) + } +} + impl Iterator for CrossingsIter where C: Cross + Clone, diff --git a/geo/src/algorithm/sweep/proc.rs b/geo/src/algorithm/sweep/proc.rs index 0b4594aac..cfb4befa7 100644 --- a/geo/src/algorithm/sweep/proc.rs +++ b/geo/src/algorithm/sweep/proc.rs @@ -6,12 +6,13 @@ use std::{ use super::*; pub(crate) struct Sweep { + is_simple: bool, events: BinaryHeap>>, active_segments: BTreeSet>>, } impl Sweep { - pub(crate) fn new(iter: I) -> Self + pub(crate) fn new(iter: I, is_simple: bool) -> Self where I: IntoIterator, { @@ -24,6 +25,7 @@ impl Sweep { let mut sweep = Sweep { events: BinaryHeap::with_capacity(size), active_segments: Default::default(), + is_simple, }; for cr in iter { IMSegment::create_segment(cr, None, None, |ev| sweep.events.push(ev)); @@ -70,58 +72,61 @@ impl Sweep { match &event.ty { LineLeft => { let mut should_add = true; - for adj_segment in prev.into_iter().chain(next.into_iter()) { - if let Some(adj_intersection) = - segment.geom().intersect_line_ordered(&adj_segment.geom()) - { - trace!("Found intersection (LL):\n\tsegment1: {:?}\n\tsegment2: {:?}\n\tintersection: {:?}", segment, adj_segment, adj_intersection); - // 1. Split adj_segment, and extra splits to storage - let adj_overlap = adj_segment - .adjust_one_segment(adj_intersection, |e| self.events.push(e)); - - // A special case is if adj_segment was split, and the - // intersection is at the start of this segment. In this - // case, there is an right-end event in the heap, that - // needs to be handled before finishing up this event. - let handle_end_event = { - // Get first point of intersection - let int_pt = adj_intersection.left(); - // Check its not first point of the adjusted, but is - // first point of current segment - int_pt != adj_segment.geom().left() && int_pt == segment.geom().left() - }; - if handle_end_event { - let event = self.events.pop().unwrap(); - let done = self.handle_event(event, cb); - debug_assert!(done, "special right-end event handling failed") - } + if !self.is_simple { + for adj_segment in prev.into_iter().chain(next.into_iter()) { + if let Some(adj_intersection) = + segment.geom().intersect_line_ordered(&adj_segment.geom()) + { + trace!("Found intersection (LL):\n\tsegment1: {:?}\n\tsegment2: {:?}\n\tintersection: {:?}", segment, adj_segment, adj_intersection); + // 1. Split adj_segment, and extra splits to storage + let adj_overlap = adj_segment + .adjust_one_segment(adj_intersection, |e| self.events.push(e)); + + // A special case is if adj_segment was split, and the + // intersection is at the start of this segment. In this + // case, there is an right-end event in the heap, that + // needs to be handled before finishing up this event. + let handle_end_event = { + // Get first point of intersection + let int_pt = adj_intersection.left(); + // Check its not first point of the adjusted, but is + // first point of current segment + int_pt != adj_segment.geom().left() + && int_pt == segment.geom().left() + }; + if handle_end_event { + let event = self.events.pop().unwrap(); + let done = self.handle_event(event, cb); + debug_assert!(done, "special right-end event handling failed") + } - // 2. Split segment, adding extra segments as needed. - let seg_overlap_key = - segment.adjust_one_segment(adj_intersection, |e| self.events.push(e)); - - assert_eq!( - adj_overlap.is_some(), - seg_overlap_key.is_some(), - "one of the intersecting segments had an overlap, but not the other!" - ); - if let Some(adj_ovl) = adj_overlap { - let tgt = seg_overlap_key.unwrap(); - trace!("setting overlap: {adj_ovl:?} -> {tgt:?}"); - adj_ovl.chain_overlap(tgt.clone()); - - if tgt == segment { - // The whole event segment is now overlapping - // some other active segment. - // - // We do not need to continue iteration, but - // should callback if the left event of the - // now-parent has already been processed. - if Borrow::>::borrow(&adj_ovl).left_event_done { - should_add = false; - break; + // 2. Split segment, adding extra segments as needed. + let seg_overlap_key = segment + .adjust_one_segment(adj_intersection, |e| self.events.push(e)); + + assert_eq!( + adj_overlap.is_some(), + seg_overlap_key.is_some(), + "one of the intersecting segments had an overlap, but not the other!" + ); + if let Some(adj_ovl) = adj_overlap { + let tgt = seg_overlap_key.unwrap(); + trace!("setting overlap: {adj_ovl:?} -> {tgt:?}"); + adj_ovl.chain_overlap(tgt.clone()); + + if tgt == segment { + // The whole event segment is now overlapping + // some other active segment. + // + // We do not need to continue iteration, but + // should callback if the left event of the + // now-parent has already been processed. + if Borrow::>::borrow(&adj_ovl).left_event_done { + should_add = false; + break; + } + return true; } - return true; } } } @@ -152,35 +157,41 @@ impl Sweep { cb_seg = seg.overlapping().cloned(); } - if let (Some(prev), Some(next)) = (prev, next) { - let prev_geom = prev.geom(); - let next_geom = next.geom(); - if let Some(adj_intersection) = prev_geom.intersect_line_ordered(&next_geom) { - // 1. Split prev_segment, and extra splits to storage - let first = prev - .adjust_one_segment(adj_intersection, |e| self.events.push(e)) - .is_none(); - let second = next - .adjust_one_segment(adj_intersection, |e| self.events.push(e)) - .is_none(); - debug_assert!( - first && second, - "adjacent segments @ removal can't overlap!" - ); + if !self.is_simple { + if let (Some(prev), Some(next)) = (prev, next) { + let prev_geom = prev.geom(); + let next_geom = next.geom(); + if let Some(adj_intersection) = prev_geom.intersect_line_ordered(&next_geom) + { + // 1. Split prev_segment, and extra splits to storage + let first = prev + .adjust_one_segment(adj_intersection, |e| self.events.push(e)) + .is_none(); + let second = next + .adjust_one_segment(adj_intersection, |e| self.events.push(e)) + .is_none(); + debug_assert!( + first && second, + "adjacent segments @ removal can't overlap!" + ); + } } } } PointLeft => { - for adj_segment in prev.into_iter().chain(next.into_iter()) { - let geom = adj_segment.geom(); - if let Some(adj_intersection) = segment.geom().intersect_line_ordered(&geom) { - trace!("Found intersection:\n\tsegment1: {:?}\n\tsegment2: {:?}\n\tintersection: {:?}", segment, adj_segment, adj_intersection); - // 1. Split adj_segment, and extra splits to storage - let adj_overlap = adj_segment - .adjust_one_segment(adj_intersection, |e| self.events.push(e)); - - // Can't have overlap with a point - debug_assert!(adj_overlap.is_none()); + if !self.is_simple { + for adj_segment in prev.into_iter().chain(next.into_iter()) { + let geom = adj_segment.geom(); + if let Some(adj_intersection) = segment.geom().intersect_line_ordered(&geom) + { + trace!("Found intersection:\n\tsegment1: {:?}\n\tsegment2: {:?}\n\tintersection: {:?}", segment, adj_segment, adj_intersection); + // 1. Split adj_segment, and extra splits to storage + let adj_overlap = adj_segment + .adjust_one_segment(adj_intersection, |e| self.events.push(e)); + + // Can't have overlap with a point + debug_assert!(adj_overlap.is_none()); + } } } diff --git a/rfcs/2022-05-24-boolean-ops.md b/rfcs/2022-05-24-boolean-ops.md index ca8947b24..1711495fe 100644 --- a/rfcs/2022-05-24-boolean-ops.md +++ b/rfcs/2022-05-24-boolean-ops.md @@ -5,8 +5,10 @@ # Summary Boolean Ops. refer to constructive set-theoretic operations applied on -geometries; for eg., union, intersection. This RFC implements boolean-ops -on 2D geometries based on the algorithm of [Martinez-Rueda-Feito]. +geometries; for eg., union, intersection. The implementation supports +boolean-ops on 2D geometries based on the algorithm of +[Martinez-Rueda-Feito] as well as extensions to clip a 1D geometry by a 2D +geometry. # Implementation Break-down @@ -43,20 +45,46 @@ segments into segments that do not intersect in interiors. The sweep also provides a list of "active segments" which are ordered list of segments intersecting the current sweep line. This ordering allows us to infer the region above/below each segment, and thus the segments belong in the output -goemetry. Specifically, the outputs segments are whose either side have -different parity with the output region (i.e. one side of the segment -belongs to the output, and the other does not). These are collected and -used to assemble the final geometry. +goemetry. Pls. refer `geo/src/algorithm/bool_ops/op.rs`. +## Operation Specification + +The initial implementation focused on operation between two 2-d +geometries. However, the core logic to use the sweep to split segments, and +calculate regions based on segment arrangement lends itself to perform +boolean operations in multiple flavours. + +The general spec. of a boolean operation is extracted as the `Spec` trait. +This trait captures the calculation of regions from the segments and their +ordering, as well as the construction of the final output. + +Specifically, for operations between two 2D geometries, the outputs +segments are whose either side have different parity with the output region +(i.e. one side of the segment belongs to the output, and the other does +not). These are collected and used to assemble the final geometry. This is +captured in the `BoolOp` struct that implements the above trait. + +We also support clipping of a 1D geometry by a 2D geometry. This is the +same as intersection but between geometries of different dimensions. In +this case, we calculate region solely from the 2D geometry, and only output +segments from the 1D geomtries. These are then assembled to output a +`MultiLineString` in this case. This is captured in the `ClipOp` struct. + +Pls. ref `geo/src/algorithm/bool_ops/spec.rs`. + ## Output Construction Here, we construct a final geometry by assembling the segments obtained -from the previous section. These segments are guaranteed to represent a -bounded region, and thus can be decomposed into a set of cycles (the -eulerian graph condition). The only constraint is the ensure the output -satisfies the validity constraints of the OGC SFS. +from the previous section. For `BoolOp`, these segments are guaranteed to +represent a bounded region, and thus can be decomposed into a set of cycles +(the eulerian graph condition). The only constraint is to ensure the +output satisfies the validity constraints of the OGC SFS. + +For `ClipOp`, we simply assemble segments that belonged to the same +original `LineString` via a greedy algorithm to assemble a final list of 1D +geometries. Pls. ref `geo/src/algorithm/bool_ops/assembly.rs`.