Improve intersection candidate finding

libraries/path-bool/src/lib.rs

@@ -9,8 +9,10 @@ mod parsing {
 mod util {
 	pub(crate) mod aabb;
 	pub(crate) mod epsilons;
+	pub(crate) mod grid;
 	pub(crate) mod math;
 	pub(crate) mod quad_tree;
+	pub(crate) mod rtree;
 }
 mod path;
 #[cfg(test)]

libraries/path-bool/src/path_boolean.rs

@@ -65,6 +65,7 @@ new_key_type! {
 
 use crate::aabb::{bounding_box_around_point, bounding_box_max_extent, expand_bounding_box, extend_bounding_box, merge_bounding_boxes, Aabb};
 use crate::epsilons::Epsilons;
+use crate::grid::{BitVec, Grid};
 use crate::intersection_path_segment::{path_segment_intersection, segments_equal};
 use crate::path::Path;
 use crate::path_cubic_segment_self_intersection::path_cubic_segment_self_intersection;
@@ -431,29 +432,42 @@ fn split_at_self_intersections(edges: &mut Vec<MajorGraphEdgeStage1>) {
 /// A tuple containing:
 /// * A vector of split edges (MajorGraphEdgeStage2).
 /// * An optional overall bounding box (AaBb) for all edges.
-fn split_at_intersections(edges: &[MajorGraphEdgeStage1]) -> (Vec<MajorGraphEdgeStage1>, Option<Aabb>) {
-	if edges.is_empty() {
-		return (Vec::new(), None);
-	}
-
+fn split_at_intersections(edges: &[MajorGraphEdgeStage1]) -> Vec<MajorGraphEdgeStage1> {
 	// Step 1: Add bounding boxes to edges
 	let with_bounding_box: Vec<MajorGraphEdgeStage2> = edges.iter().map(|(seg, parent)| (*seg, *parent, seg.approx_bounding_box())).collect();
 	// Step 2: Calculate total bounding box
-	let total_bounding_box = with_bounding_box.iter().fold(Default::default(), |acc, (_, _, bb)| merge_bounding_boxes(&acc, &bb));
+	let total_bounding_box = with_bounding_box.iter().fold(Default::default(), |acc, (_, _, bb)| merge_bounding_boxes(&acc, bb));
+
+	let max_extent = bounding_box_max_extent(&total_bounding_box);
+	let cell_size = max_extent / (edges.len() as f64).sqrt();
+
+	// Step 3: Create grid for efficient intersection checks
+	let mut grid = Grid::new(cell_size, edges.len());
 
 	// Step 3: Create edge tree for efficient intersection checks
-	let mut edge_tree = QuadTree::new(total_bounding_box, INTERSECTION_TREE_DEPTH, 8);
+	// let mut edge_tree = QuadTree::new(total_bounding_box, INTERSECTION_TREE_DEPTH, 16);
+	// let mut rtree = crate::util::rtree::RTree::new(24);
 
-	let mut splits_per_edge: HashMap<usize, Vec<f64>> = HashMap::default();
+	let mut splits_per_edge: Vec<Vec<f64>> = vec![Vec::new(); edges.len()];
 
-	fn add_split(splits_per_edge: &mut HashMap<usize, Vec<f64>>, i: usize, t: f64) {
-		splits_per_edge.entry(i).or_default().push(t);
+	fn add_split(splits_per_edge: &mut [Vec<f64>], i: usize, t: f64) {
+		splits_per_edge[i].push(t);
 	}
+	// let mut candidates = Vec::with_capacity(8);
+	let mut candidates = BitVec::new(edges.len());
 
 	// Step 4: Find intersections and record split points
 	for (i, edge) in with_bounding_box.iter().enumerate() {
-		let candidates = edge_tree.find(&edge.2);
-		for &j in &candidates {
+		// let candidates = edge_tree.find(&edge.2);
+		// let mut quad_candidates: Vec<_> = quad_candidates.into_iter().collect();
+		// quad_candidates.sort_unstable();
+		// let mut candidates = rtree.query(&edge.2);
+		// candidates.sort_unstable();
+		// assert_eq!(candidates, quad_candidates);
+		candidates.clear();
+		grid.query(&edge.2, &mut candidates);
+
+		for j in candidates.iter_set_bits() {
 			let candidate: &(PathSegment, u8) = &edges[j];
 			let include_endpoints = edge.1 != candidate.1 || !(candidate.0.end().abs_diff_eq(edge.0.start(), EPS.point) || candidate.0.start().abs_diff_eq(edge.0.end(), EPS.point));
 			let intersection = path_segment_intersection(&edge.0, &candidate.0, include_endpoints, &EPS);
@@ -462,14 +476,16 @@ fn split_at_intersections(edges: &[MajorGraphEdgeStage1]) -> (Vec<MajorGraphEdge
 				add_split(&mut splits_per_edge, j, t1);
 			}
 		}
-		edge_tree.insert(edge.2, i);
+		grid.insert(&edge.2, i);
+		// edge_tree.insert(edge.2, i);
+		// rtree.insert(edge.2, i);
 	}
 
 	// Step 5: Apply splits to create new edges
 	let mut new_edges = Vec::new();
 
 	for (i, (seg, parent, _)) in with_bounding_box.into_iter().enumerate() {
-		if let Some(splits) = splits_per_edge.get(&i) {
+		if let Some(splits) = splits_per_edge.get(i) {
 			let mut splits = splits.clone();
 			splits.sort_by(|a, b| a.partial_cmp(b).unwrap());
 			let mut tmp_seg = seg;
@@ -496,7 +512,7 @@ fn split_at_intersections(edges: &[MajorGraphEdgeStage1]) -> (Vec<MajorGraphEdge
 		}
 	}
 
-	(new_edges, Some(total_bounding_box))
+	new_edges
 }
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
@@ -539,7 +555,7 @@ fn round_point(point: DVec2) -> I64Vec2 {
 }
 
 // TODO: Using 32bit values here might lead to incorrect results when the values collide. Even though this is very unlikely we should think about this case
-fn find_vertices(edges: &[MajorGraphEdgeStage1], total_bounding_box: Aabb) -> MajorGraph {
+fn find_vertices(edges: &[MajorGraphEdgeStage1]) -> MajorGraph {
 	let mut graph = MajorGraph {
 		edges: SlotMap::with_capacity_and_key(edges.len() * 2),
 		vertices: SlotMap::with_capacity_and_key(edges.len()),
@@ -1731,21 +1747,19 @@ impl Display for BooleanError {
 pub fn path_boolean(a: &Path, a_fill_rule: FillRule, b: &Path, b_fill_rule: FillRule, op: PathBooleanOperation) -> Result<Vec<Path>, BooleanError> {
 	let mut unsplit_edges: Vec<MajorGraphEdgeStage1> = a.iter().map(segment_to_edge(1)).chain(b.iter().map(segment_to_edge(2))).flatten().collect();
 
+	if unsplit_edges.is_empty() {
+		return Ok(Vec::new());
+	}
 	split_at_self_intersections(&mut unsplit_edges);
 
-	let (split_edges, total_bounding_box) = split_at_intersections(&unsplit_edges);
+	let split_edges = split_at_intersections(&unsplit_edges);
 
 	#[cfg(feature = "logging")]
 	for (edge, _) in split_edges.iter() {
 		eprintln!("{}", path_to_path_data(&vec![*edge], 0.001));
 	}
 
-	let total_bounding_box = match total_bounding_box {
-		Some(bb) => bb,
-		None => return Ok(Vec::new()), // Input geometry is empty
-	};
-
-	let major_graph = find_vertices(&split_edges, total_bounding_box);
+	let major_graph = find_vertices(&split_edges);
 
 	#[cfg(feature = "logging")]
 	eprintln!("Major graph:");
@@ -1832,10 +1846,7 @@ mod tests {
 	#[test]
 	fn test_split_at_intersections() {
 		let unsplit_edges = unsplit_edges();
-		let (split_edges, total_bounding_box) = split_at_intersections(&unsplit_edges);
-
-		// Check that we have a valid bounding box
-		assert!(total_bounding_box.is_some());
+		let split_edges = split_at_intersections(&unsplit_edges);
 
 		// Check that we have more edges after splitting (due to intersections)
 		assert!(split_edges.len() >= unsplit_edges.len());
@@ -1862,8 +1873,8 @@ mod tests {
 	fn test_compute_minor() {
 		// Set up the initial graph
 		let unsplit_edges = unsplit_edges();
-		let (split_edges, total_bounding_box) = split_at_intersections(&unsplit_edges);
-		let major_graph = find_vertices(&split_edges, total_bounding_box.unwrap());
+		let split_edges = split_at_intersections(&unsplit_edges);
+		let major_graph = find_vertices(&split_edges);
 
 		// Compute minor graph
 		let minor_graph = compute_minor(&major_graph);
@@ -1917,8 +1928,8 @@ mod tests {
 	fn test_sort_outgoing_edges_by_angle() {
 		// Set up the initial graph
 		let unsplit_edges = unsplit_edges();
-		let (split_edges, total_bounding_box) = split_at_intersections(&unsplit_edges);
-		let major_graph = find_vertices(&split_edges, total_bounding_box.unwrap());
+		let split_edges = split_at_intersections(&unsplit_edges);
+		let major_graph = find_vertices(&split_edges);
 		let mut minor_graph = compute_minor(&major_graph);
 
 		// Print initial state

libraries/path-bool/src/util/grid.rs

@@ -0,0 +1,129 @@
+use crate::aabb::Aabb;
+use glam::{DVec2, IVec2};
+use rustc_hash::FxHashMap;
+use smallvec::SmallVec;
+
+pub(crate) struct Grid {
+	cell_size: f64,
+	cells: FxHashMap<IVec2, SmallVec<[usize; 6]>>,
+}
+
+impl Grid {
+	pub(crate) fn new(cell_size: f64, edges: usize) -> Self {
+		Grid {
+			cell_size,
+			cells: FxHashMap::with_capacity_and_hasher(edges, Default::default()),
+		}
+	}
+
+	pub(crate) fn insert(&mut self, bbox: &Aabb, index: usize) {
+		let min_cell = self.point_to_cell(bbox.min());
+		let max_cell = self.point_to_cell(bbox.max());
+
+		for i in min_cell.x..=max_cell.x {
+			for j in min_cell.y..=max_cell.y {
+				self.cells.entry((i, j).into()).or_default().push(index);
+			}
+		}
+	}
+
+	pub(crate) fn query(&self, bbox: &Aabb, result: &mut BitVec) {
+		let min_cell = self.point_to_cell(bbox.min());
+		let max_cell = self.point_to_cell(bbox.max());
+
+		for i in min_cell.x..=max_cell.x {
+			for j in min_cell.y..=max_cell.y {
+				if let Some(indices) = self.cells.get(&(i, j).into()) {
+					for &index in indices {
+						result.set(index);
+					}
+				}
+			}
+		}
+		// result.sort_unstable();
+		// result.dedup();
+	}
+
+	fn point_to_cell(&self, point: DVec2) -> IVec2 {
+		(point / self.cell_size).as_ivec2()
+	}
+}
+
+use std::mem::size_of;
+
+pub struct BitVec {
+	data: Vec<u64>,
+	capacity: usize,
+}
+
+impl BitVec {
+	pub fn new(capacity: usize) -> Self {
+		let num_words = (capacity + 63) / 64;
+		BitVec { data: vec![0; num_words], capacity }
+	}
+
+	pub fn set(&mut self, index: usize) {
+		let word_index = index / 64;
+		let bit_index = index % 64;
+		self.data[word_index] |= 1u64 << bit_index;
+	}
+
+	pub fn clear(&mut self) {
+		self.data.fill(0);
+	}
+
+	pub fn iter_set_bits(&self) -> BitVecIterator {
+		BitVecIterator {
+			bit_vec: self,
+			current_word: self.data[0],
+			word_index: 0,
+		}
+	}
+}
+
+pub struct BitVecIterator<'a> {
+	bit_vec: &'a BitVec,
+	current_word: u64,
+	word_index: usize,
+}
+
+impl<'a> Iterator for BitVecIterator<'a> {
+	type Item = usize;
+
+	fn next(&mut self) -> Option<Self::Item> {
+		while self.word_index < self.bit_vec.data.len() {
+			if self.current_word == 0 {
+				self.word_index += 1;
+				if self.word_index == self.bit_vec.data.len() {
+					return None;
+				}
+				self.current_word = self.bit_vec.data[self.word_index];
+				continue;
+			}
+			let tz = self.current_word.trailing_zeros() as usize;
+			self.current_word ^= 1 << tz;
+
+			let result = self.word_index * 64 + tz;
+
+			return Some(result);
+		}
+		None
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+
+	#[test]
+	fn test_bitvec() {
+		let mut bv = BitVec::new(200);
+		bv.set(5);
+		bv.set(64);
+		bv.set(128);
+		bv.set(199);
+
+		let set_bits: Vec<usize> = bv.iter_set_bits().collect();
+		assert_eq!(set_bits, vec![5, 64, 128, 199]);
+	}
+}

libraries/path-bool/src/util/rtree.rs

@@ -0,0 +1,190 @@
+use glam::DVec2;
+use std::cmp::Ordering;
+
+use crate::aabb::{bounding_boxes_overlap, merge_bounding_boxes, Aabb};
+
+pub(crate) struct RTreeNode {
+	bbox: Aabb,
+	children: Vec<RTreeNode>,
+	entries: Vec<(Aabb, usize)>,
+	is_leaf: bool,
+}
+
+impl RTreeNode {
+	fn new(is_leaf: bool) -> Self {
+		RTreeNode {
+			bbox: Aabb::default(),
+			children: Vec::new(),
+			entries: Vec::new(),
+			is_leaf,
+		}
+	}
+
+	fn insert(&mut self, bbox: Aabb, index: usize, max_entries: usize) {
+		if self.is_leaf {
+			self.entries.push((bbox, index));
+			self.update_bbox(&bbox);
+			if self.entries.len() > max_entries {
+				self.split(max_entries);
+			}
+		} else {
+			let best_child = self.choose_subtree(&bbox);
+			self.children[best_child].insert(bbox, index, max_entries);
+			self.update_bbox(&self.children[best_child].bbox.clone());
+		}
+	}
+
+	fn choose_subtree(&self, bbox: &Aabb) -> usize {
+		self.children
+			.iter()
+			.enumerate()
+			.min_by(|(_, a), (_, b)| {
+				let area_increase_a = area(&merge_bounding_boxes(&a.bbox, bbox)) - area(&a.bbox);
+				let area_increase_b = area(&merge_bounding_boxes(&b.bbox, bbox)) - area(&b.bbox);
+				area_increase_a.partial_cmp(&area_increase_b).unwrap_or(Ordering::Equal)
+			})
+			.map(|(index, _)| index)
+			.unwrap_or(0)
+	}
+
+	fn update_bbox(&mut self, bbox: &Aabb) {
+		self.bbox = merge_bounding_boxes(&self.bbox, bbox);
+	}
+
+	fn split(&mut self, max_entries: usize) {
+		if !self.is_leaf {
+			return; // We only split leaf nodes in this implementation
+		}
+
+		let mut new_node = RTreeNode::new(true);
+		let (seed1, seed2) = self.pick_seeds();
+
+		new_node.entries.push(self.entries[seed2].clone());
+		new_node.update_bbox(&new_node.entries[0].0.clone());
+
+		self.bbox = self.entries[seed1].0;
+
+		let mut group1 = vec![self.entries[seed1].clone()];
+		let mut group2 = vec![self.entries[seed2].clone()];
+
+		let remaining_entries: Vec<_> = self.entries.iter().enumerate().filter(|&(i, _)| i != seed1 && i != seed2).map(|(_, e)| e.clone()).collect();
+
+		for entry in remaining_entries {
+			if group1.len() >= max_entries / 2 {
+				group2.push(entry);
+			} else if group2.len() >= max_entries / 2 {
+				group1.push(entry);
+			} else {
+				let increase1 = area(&merge_bounding_boxes(&self.bbox, &entry.0)) - area(&self.bbox);
+				let increase2 = area(&merge_bounding_boxes(&new_node.bbox, &entry.0)) - area(&new_node.bbox);
+				if increase1 < increase2 {
+					group1.push(entry);
+					self.update_bbox(&entry.0);
+				} else {
+					group2.push(entry);
+					new_node.update_bbox(&entry.0);
+				}
+			}
+		}
+
+		self.entries = group1;
+		new_node.entries = group2;
+
+		// Create a new parent if this was the root
+		if self.children.is_empty() {
+			let mut new_parent = RTreeNode::new(false);
+			new_parent.children.push(std::mem::replace(self, RTreeNode::new(true)));
+			new_parent.children.push(new_node);
+			new_parent.update_bbox(&new_parent.children[0].bbox.clone());
+			new_parent.update_bbox(&new_parent.children[1].bbox.clone());
+			*self = new_parent;
+		} else {
+			self.children.push(new_node);
+		}
+	}
+
+	fn pick_seeds(&self) -> (usize, usize) {
+		let mut max_waste = f64::NEG_INFINITY;
+		let mut seeds = (0, 1);
+
+		for i in 0..self.entries.len() {
+			for j in (i + 1)..self.entries.len() {
+				let combined_bbox = merge_bounding_boxes(&self.entries[i].0, &self.entries[j].0);
+				let waste = area(&combined_bbox) - area(&self.entries[i].0) - area(&self.entries[j].0);
+
+				if waste > max_waste {
+					max_waste = waste;
+					seeds = (i, j);
+				}
+			}
+		}
+
+		seeds
+	}
+
+	fn pick_next(&self, other: &RTreeNode) -> usize {
+		self.entries
+			.iter()
+			.enumerate()
+			.max_by(|(_, a), (_, b)| {
+				let diff_a = area(&merge_bounding_boxes(&self.bbox, &a.0)) - area(&merge_bounding_boxes(&other.bbox, &a.0));
+				let diff_b = area(&merge_bounding_boxes(&self.bbox, &b.0)) - area(&merge_bounding_boxes(&other.bbox, &b.0));
+				diff_a.partial_cmp(&diff_b).unwrap_or(Ordering::Equal)
+			})
+			.map(|(index, _)| index)
+			.unwrap_or(0)
+	}
+
+	fn assign_entry(&mut self, index: usize, other: &mut RTreeNode) {
+		let entry = self.entries.remove(index);
+		other.entries.push(entry);
+		other.update_bbox(&entry.0);
+	}
+
+	fn query(&self, search_bbox: &Aabb, results: &mut Vec<usize>) {
+		if !bounding_boxes_overlap(&self.bbox, search_bbox) {
+			return;
+		}
+
+		if self.is_leaf {
+			for (bbox, index) in &self.entries {
+				if bounding_boxes_overlap(bbox, search_bbox) {
+					results.push(*index);
+				}
+			}
+		} else {
+			for child in &self.children {
+				child.query(search_bbox, results);
+			}
+		}
+	}
+}
+
+pub(crate) struct RTree {
+	root: RTreeNode,
+	max_entries: usize,
+}
+
+impl RTree {
+	pub(crate) fn new(max_entries: usize) -> Self {
+		RTree {
+			root: RTreeNode::new(true),
+			max_entries,
+		}
+	}
+
+	pub(crate) fn insert(&mut self, bbox: Aabb, index: usize) {
+		self.root.insert(bbox, index, self.max_entries);
+	}
+
+	pub(crate) fn query(&self, search_bbox: &Aabb) -> Vec<usize> {
+		let mut results = Vec::new();
+		self.root.query(search_bbox, &mut results);
+		results
+	}
+}
+
+// Helper function to calculate the area of an Aabb
+fn area(bbox: &Aabb) -> f64 {
+	(bbox.right() - bbox.left()) * (bbox.bottom() - bbox.top())
+}