Improve robustness and performance of the boolean operation algorithm (#2191)

* Improve perf of path bool lib

* Use swap remove

* Use outer/inner bounding box for inclusion testing

* Reuse allocations for hit testing

* Use direct root finding for inclusion testing

* Reuse bounding box

* Use faster hash and specify capacities

* Use hashmap based approach for find vertices

* Unroll find_vertecies loop and use 32 bit positions

* Tune initial vec capacities

* Remove unused bounding boxes

* Use smallvec for storing outgoing edges

* Improve allocations for compute_minor

* Use approximate bounding box for edge finding

* Transition aabb to use glam vecs

* Make find vertecies use 64 bit again this is slower but less likely to cause issues

* Improve intersection candidate finding

* Remove loop check in bit vec iter

* Special case cubic line intersections

* Optimize grid rounding and add debug output

* Remove file write

* Remove faulty line intersection

* Fix grid rounding

* Improve robustness and cleanaup code

* Make elided lifetime explicit

* Fix tests

* Fix a boolean ops crash

* Add comment

---------

Co-authored-by: Keavon Chambers <keavon@keavon.com>

Cargo.lock

@@ -3859,7 +3859,10 @@ dependencies = [
  "lyon_geom",
  "regex",
  "resvg",
+ "roots",
+ "rustc-hash 2.1.1",
  "slotmap",
+ "smallvec",
  "svg",
 ]
 
@@ -4687,6 +4690,12 @@ dependencies = [
  "serde_derive",
 ]
 
+[[package]]
+name = "roots"
+version = "0.0.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "082f11ffa03bbef6c2c6ea6bea1acafaade2fd9050ae0234ab44a2153742b058"
+
 [[package]]
 name = "roxmltree"
 version = "0.20.0"

libraries/path-bool/Cargo.toml

@@ -28,6 +28,9 @@ glam = "0.29.0"
 regex = "1.10.6"
 slotmap = "1.0.7"
 lyon_geom = "1.0"
+roots = "0.0.8"
+rustc-hash = "2.0.0"
+smallvec = "1.13.2"
 
 [dev-dependencies]
 glob = "0.3"

libraries/path-bool/src/lib.rs

@@ -9,8 +9,8 @@ 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;
 }
 mod path;
 #[cfg(test)]
@@ -311,4 +311,23 @@ mod test {
 		assert_eq!(result.len(), 1);
 		assert!(!result[0].is_empty());
 	}
+	#[test]
+	fn real_01() {
+		let a = path_from_path_data(
+			"M 212.67152,105 A 64.171516,64.171516 0 0 1 148.5,169.17152 64.171516,64.171516 0 0 1 84.328484,105 64.171516,64.171516 0 0 1 148.5,40.828484 64.171516,64.171516 0 0 1 212.67152,105 Z",
+		)
+		.unwrap();
+		let b = path_from_path_data(
+			"m 83.755387,112.6962 h -7.62 v 37.0332 h 7.62 z m 22.097973,9.144 c -3.4544,0 -5.892801,1.4732 -7.620001,4.5212 v -4.064 H 91.12136 v 38.5064 h 7.111999 v -14.3256 c 1.7272,3.048 4.165601,4.4704 7.620001,4.4704 6.604,0 11.4808,-6.1976 11.4808,-14.5288 0,-8.0772 -4.3688,-14.5796 -11.4808,-14.5796 z m -1.6256,5.9436 c 3.6068,0 5.9944,3.5052 5.9944,8.7376 0,4.9784 -2.4892,8.4836 -5.9944,8.4836 -3.556,0 -5.994401,-3.4544 -5.994401,-8.5852 0,-5.1308 2.438401,-8.636 5.994401,-8.636 z m 23.62201,13.97 h -6.9596 c 0.2032,5.9944 4.6228,9.144 12.954,9.144 9.6012,0 11.9888,-5.4864 11.9888,-9.2964 0,-3.556 -1.778,-5.842 -5.3848,-6.9088 l -8.9916,-2.5908 c -1.9812,-0.6096 -2.4892,-1.016 -2.4892,-2.1336 0,-1.524 1.6256,-2.54 4.1148,-2.54 3.4036,0 5.08,1.2192 5.1308,3.7084 h 6.858 c -0.1016,-5.7912 -4.572,-9.2964 -11.938,-9.2964 -6.9596,0 -11.2776,3.5052 -11.2776,9.144 0,5.3848 4.4196,6.2484 5.9944,6.7564 l 8.4836,2.6416 c 1.778,0.5588 2.3876,1.1176 2.3876,2.2352 0,1.6764 -1.9812,2.6924 -5.2832,2.6924 -4.4704,0 -5.1816,-1.6764 -5.588,-3.556 z m 47.59959,7.9756 v -27.432 h -7.112 v 17.1704 c 0,3.2512 -2.286,5.3848 -5.7404,5.3848 -3.048,0 -4.572,-1.6256 -4.572,-4.9276 v -17.6276 h -7.112 v 19.1008 c 0,6.0452 3.3528,9.4996 9.1948,9.4996 3.7084,0 6.1976,-1.3716 8.2296,-4.4196 v 3.2512 z m 6.60404,-27.432 v 27.432 h 7.112 v -16.4592 c 0,-3.3528 1.8288,-5.3848 4.8768,-5.3848 2.3876,0 3.8608,1.3716 3.8608,3.556 v 18.288 h 7.112 v -16.4592 c 0,-3.3528 1.8288,-5.3848 4.8768,-5.3848 2.3876,0 3.8608,1.3716 3.8608,3.556 v 18.288 h 7.112 v -19.4056 c 0,-5.334 -3.2512,-8.4836 -8.7376,-8.4836 -3.4544,0 -5.8928,1.2192 -8.0264,4.064 -1.3208,-2.5908 -4.064,-4.064 -7.4676,-4.064 -3.1496,0 -5.1816,1.0668 -7.5184,3.8608 v -3.4036 z M 81.012192,49.196201 h -7.62 v 37.0332 h 25.3492 v -6.35 h -17.7292 z m 35.305978,9.144 c -8.382,0 -13.5128,5.5372 -13.5128,14.5288 0,9.0424 5.1308,14.5288 13.5636,14.5288 8.3312,0 13.5636,-5.5372 13.5636,-14.3256 0,-9.2964 -5.0292,-14.732 -13.6144,-14.732 z m 0.0508,5.7404 c 3.9116,0 6.4516,3.5052 6.4516,8.89 0,5.1308 -2.6416,8.6868 -6.4516,8.6868 -3.8608,0 -6.4516,-3.556 -6.4516,-8.7884 0,-5.2324 2.5908,-8.7884 6.4516,-8.7884 z m 18.89761,-5.2832 v 27.432 h 7.112 v -14.5796 c 0,-4.1656 2.0828,-6.2484 6.2484,-6.2484 0.762,0 1.27,0.0508 2.2352,0.2032 v -7.2136 c -0.4064,-0.0508 -0.6604,-0.0508 -0.8636,-0.0508 -3.2512,0 -6.096,2.1336 -7.62,5.842 v -5.3848 z m 31.34362,-0.4572 c -7.874,0 -12.7,5.6896 -12.7,14.8844 0,8.7884 4.7752,14.1732 12.5476,14.1732 6.14679,0 11.12519,-3.5052 12.69999,-8.89 h -7.0104 c -0.8636,2.1844 -2.8448,3.4544 -5.43559,3.4544 -4.7752,0 -5.5372,-3.5052 -5.6896,-7.2136 h 18.38959 c 0.0508,-0.6096 0.0508,-0.8636 0.0508,-1.2192 0,-12.1412 -7.366,-15.1892 -12.85239,-15.1892 z m 5.43559,11.684 H 161.1238 c 0.4572,-4.1656 2.2352,-6.2484 5.3848,-6.2484 3.30199,0 5.23239,2.2352 5.53719,6.2484 z m 12.75081,-11.2268 v 27.432 h 7.112 v -16.4592 c 0,-3.3528 1.8288,-5.3848 4.8768,-5.3848 2.3876,0 3.8608,1.3716 3.8608,3.556 v 18.288 h 7.112 v -16.4592 c 0,-3.3528 1.8288,-5.3848 4.8768,-5.3848 2.3876,0 3.8608,1.3716 3.8608,3.556 v 18.288 h 7.112 v -19.4056 c 0,-5.334 -3.2512,-8.4836 -8.7376,-8.4836 -3.4544,0 -5.8928,1.2192 -8.0264,4.064 -1.3208,-2.5908 -4.064,-4.064 -7.4676,-4.064 -3.1496,0 -5.1816,1.0668 -7.5184,3.8608 v -3.4036 z",
+		)
+		.unwrap();
+
+		let result = path_boolean(&a, FillRule::NonZero, &b, FillRule::NonZero, PathBooleanOperation::Union).unwrap();
+
+		// Add assertions here based on expected results
+		assert_eq!(result.len(), 1, "Expected 1 resulting path for Union operation");
+		dbg!(path_to_path_data(&result[0], 0.001));
+		// Add more specific assertions about the resulting path if needed
+		assert!(!result[0].is_empty());
+	}
 }

libraries/path-bool/src/path/intersection_path_segment.rs

@@ -34,13 +34,13 @@ fn subdivide_intersection_segment(int_seg: &IntersectionSegment) -> [Intersectio
 			seg: seg0,
 			start_param: int_seg.start_param,
 			end_param: mid_param,
-			bounding_box: seg0.bounding_box(),
+			bounding_box: seg0.approx_bounding_box(),
 		},
 		IntersectionSegment {
 			seg: seg1,
 			start_param: mid_param,
 			end_param: int_seg.end_param,
-			bounding_box: seg1.bounding_box(),
+			bounding_box: seg1.approx_bounding_box(),
 		},
 	]
 }
@@ -116,8 +116,9 @@ pub fn path_segment_intersection(seg0: &PathSegment, seg1: &PathSegment, endpoin
 			return intersections;
 		}
 		_ => (),
-	}
+	};
 
+	// Fallback for quadratics and arc segments
 	// https://math.stackexchange.com/questions/20321/how-can-i-tell-when-two-cubic-b%C3%A9zier-curves-intersect
 
 	let mut pairs = vec![(
@@ -125,13 +126,13 @@ pub fn path_segment_intersection(seg0: &PathSegment, seg1: &PathSegment, endpoin
 			seg: *seg0,
 			start_param: 0.,
 			end_param: 1.,
-			bounding_box: seg0.bounding_box(),
+			bounding_box: seg0.approx_bounding_box(),
 		},
 		IntersectionSegment {
 			seg: *seg1,
 			start_param: 0.,
 			end_param: 1.,
-			bounding_box: seg1.bounding_box(),
+			bounding_box: seg1.approx_bounding_box(),
 		},
 	)];
 	let mut next_pairs = Vec::new();
@@ -145,7 +146,7 @@ pub fn path_segment_intersection(seg0: &PathSegment, seg1: &PathSegment, endpoin
 	while !pairs.is_empty() {
 		next_pairs.clear();
 
-		if pairs.len() > 1000 {
+		if pairs.len() > 256 {
 			return calculate_overlap_intersections(seg0, seg1, eps);
 		}
 
@@ -191,10 +192,6 @@ pub fn path_segment_intersection(seg0: &PathSegment, seg1: &PathSegment, endpoin
 		std::mem::swap(&mut pairs, &mut next_pairs);
 	}
 
-	if !endpoints {
-		params.retain(|[s, t]| (s > &eps.param && s < &(1. - eps.param)) || (t > &eps.param && t < &(1. - eps.param)));
-	}
-
 	params
 }
 

libraries/path-bool/src/path/line_segment_aabb.rs

@@ -10,15 +10,15 @@ const TOP: u8 = 1 << 3;
 fn out_code(x: f64, y: f64, bounding_box: &Aabb) -> u8 {
 	let mut code = INSIDE;
 
-	if x < bounding_box.left {
+	if x < bounding_box.left() {
 		code |= LEFT;
-	} else if x > bounding_box.right {
+	} else if x > bounding_box.right() {
 		code |= RIGHT;
 	}
 
-	if y < bounding_box.top {
+	if y < bounding_box.top() {
 		code |= BOTTOM;
-	} else if y > bounding_box.bottom {
+	} else if y > bounding_box.bottom() {
 		code |= TOP;
 	}
 
@@ -57,20 +57,20 @@ pub(crate) fn line_segment_aabb_intersect(seg: LineSegment, bounding_box: &Aabb)
 			// outcode bit being tested guarantees the denominator is non-zero
 			if (outcode_out & TOP) != 0 {
 				// point is above the clip window
-				x = p0.x + (p1.x - p0.x) * (bounding_box.bottom - p0.y) / (p1.y - p0.y);
+				x = p0.x + (p1.x - p0.x) * (bounding_box.bottom() - p0.y) / (p1.y - p0.y);
-				y = bounding_box.bottom;
+				y = bounding_box.bottom();
 			} else if (outcode_out & BOTTOM) != 0 {
 				// point is below the clip window
-				x = p0.x + (p1.x - p0.x) * (bounding_box.top - p0.y) / (p1.y - p0.y);
+				x = p0.x + (p1.x - p0.x) * (bounding_box.top() - p0.y) / (p1.y - p0.y);
-				y = bounding_box.top;
+				y = bounding_box.top();
 			} else if (outcode_out & RIGHT) != 0 {
 				// point is to the right of clip window
-				y = p0.y + (p1.y - p0.y) * (bounding_box.right - p0.x) / (p1.x - p0.x);
+				y = p0.y + (p1.y - p0.y) * (bounding_box.right() - p0.x) / (p1.x - p0.x);
-				x = bounding_box.right;
+				x = bounding_box.right();
 			} else if (outcode_out & LEFT) != 0 {
 				// point is to the left of clip window
-				y = p0.y + (p1.y - p0.y) * (bounding_box.left - p0.x) / (p1.x - p0.x);
+				y = p0.y + (p1.y - p0.y) * (bounding_box.left() - p0.x) / (p1.x - p0.x);
-				x = bounding_box.left;
+				x = bounding_box.left();
 			}
 
 			// Now we move outside point to intersection point to clip

libraries/path-bool/src/path/path_segment.rs

@@ -588,21 +588,16 @@ impl PathSegment {
 	/// An [`Aabb`] representing the axis-aligned bounding box of the segment.
 	pub(crate) fn bounding_box(&self) -> Aabb {
 		match *self {
-			PathSegment::Line(start, end) => Aabb {
+			PathSegment::Line(start, end) => Aabb::new(start.x.min(end.x), start.y.min(end.y), start.x.max(end.x), start.y.max(end.y)),
-				top: start.y.min(end.y),
-				right: start.x.max(end.x),
-				bottom: start.y.max(end.y),
-				left: start.x.min(end.x),
-			},
 			PathSegment::Cubic(p1, p2, p3, p4) => {
 				let (left, right) = cubic_bounding_interval(p1.x, p2.x, p3.x, p4.x);
 				let (top, bottom) = cubic_bounding_interval(p1.y, p2.y, p3.y, p4.y);
-				Aabb { top, right, bottom, left }
+				Aabb::new(left, top, right, bottom)
 			}
 			PathSegment::Quadratic(p1, p2, p3) => {
 				let (left, right) = quadratic_bounding_interval(p1.x, p2.x, p3.x);
 				let (top, bottom) = quadratic_bounding_interval(p1.y, p2.y, p3.y);
-				Aabb { top, right, bottom, left }
+				Aabb::new(left, top, right, bottom)
 			}
 			PathSegment::Arc(start, rx, ry, phi, _, _, end) => {
 				if let Some(center_param) = self.arc_segment_to_center() {
@@ -627,11 +622,11 @@ impl PathSegment {
 					} else {
 						// TODO: Don't convert to cubics
 						let cubics = self.arc_segment_to_cubics(PI / 16.);
-						let mut bounding_box = None;
+						let mut bounding_box = bounding_box_around_point(start, 0.);
 						for cubic_seg in cubics {
-							bounding_box = Some(merge_bounding_boxes(bounding_box, &cubic_seg.bounding_box()));
+							bounding_box = merge_bounding_boxes(&bounding_box, &cubic_seg.bounding_box());
 						}
-						bounding_box.unwrap_or_else(|| bounding_box_around_point(start, 0.))
+						bounding_box
 					}
 				} else {
 					extend_bounding_box(Some(bounding_box_around_point(start, 0.)), end)
@@ -640,6 +635,30 @@ impl PathSegment {
 		}
 	}
 
+	/// Computes a loose bounding box that surrounds all anchors, but also the handles of cubic and quadratic segments.
+	/// This will usually be larger than the actual bounding box, but is faster to compute because it does not have to find where each curve reaches its maximum and minimum.
+	pub(crate) fn approx_bounding_box(&self) -> Aabb {
+		match *self {
+			PathSegment::Cubic(p1, p2, p3, p4) => {
+				// Use the control points to create a bounding box
+				let left = p1.x.min(p2.x).min(p3.x).min(p4.x);
+				let right = p1.x.max(p2.x).max(p3.x).max(p4.x);
+				let top = p1.y.min(p2.y).min(p3.y).min(p4.y);
+				let bottom = p1.y.max(p2.y).max(p3.y).max(p4.y);
+				Aabb::new(left, top, right, bottom)
+			}
+			PathSegment::Quadratic(p1, p2, p3) => {
+				// Use the control points to create a bounding box
+				let left = p1.x.min(p2.x).min(p3.x);
+				let right = p1.x.max(p2.x).max(p3.x);
+				let top = p1.y.min(p2.y).min(p3.y);
+				let bottom = p1.y.max(p2.y).max(p3.y);
+				Aabb::new(left, top, right, bottom)
+			}
+			seg => seg.bounding_box(),
+		}
+	}
+
 	/// Splits the path segment at a given parameter value.
 	///
 	/// # Arguments

libraries/path-bool/src/path_boolean.rs

@@ -63,21 +63,30 @@ new_key_type! {
 //
 // SPDX-License-Identifier: MIT
 
-use crate::aabb::{Aabb, bounding_box_around_point, bounding_box_max_extent, merge_bounding_boxes};
+use crate::aabb::{Aabb, bounding_box_max_extent, extend_bounding_box, merge_bounding_boxes};
 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;
 use crate::path_segment::PathSegment;
 #[cfg(feature = "logging")]
 use crate::path_to_path_data;
-use crate::quad_tree::QuadTree;
+
-use glam::DVec2;
+use glam::{BVec2, DVec2, I64Vec2};
+use roots::{Roots, find_roots_cubic};
+use rustc_hash::FxHashMap as HashMap;
+use rustc_hash::FxHashSet as HashSet;
 use slotmap::{SlotMap, new_key_type};
+use smallvec::SmallVec;
 use std::cmp::Ordering;
-use std::collections::{HashMap, HashSet, VecDeque};
+use std::collections::VecDeque;
 use std::fmt::Display;
 
+fn new_hash_map<K, V>(capacity: usize) -> HashMap<K, V> {
+	HashMap::with_capacity_and_hasher(capacity, Default::default())
+}
+
 /// Represents the types of boolean operations that can be performed on paths.
 #[derive(Debug, Clone, Copy)]
 pub enum PathBooleanOperation {
@@ -128,9 +137,6 @@ pub enum FillRule {
 	EvenOdd,
 }
 
-const INTERSECTION_TREE_DEPTH: usize = 8;
-const POINT_TREE_DEPTH: usize = 8;
-
 pub const EPS: Epsilons = Epsilons {
 	point: 1e-5,
 	linear: 1e-4,
@@ -153,7 +159,7 @@ pub struct MajorGraphEdge {
 pub struct MajorGraphVertex {
 	#[cfg_attr(not(feature = "logging"), expect(dead_code))]
 	pub point: DVec2,
-	outgoing_edges: Vec<MajorEdgeKey>,
+	outgoing_edges: SmallVec<[MajorEdgeKey; 4]>,
 }
 
 /// Represents the initial graph structure used in boolean operations.
@@ -167,7 +173,7 @@ struct MajorGraph {
 
 #[derive(Debug, Clone, PartialEq)]
 struct MinorGraphEdge {
-	segments: Vec<PathSegment>,
+	segments: SmallVec<[PathSegment; 4]>,
 	parent: u8,
 	incident_vertices: [MinorVertexKey; 2],
 	direction_flag: Direction,
@@ -215,7 +221,7 @@ impl PartialOrd for MinorGraphEdge {
 
 #[derive(Debug, Clone, Default)]
 struct MinorGraphVertex {
-	outgoing_edges: Vec<MinorEdgeKey>,
+	outgoing_edges: SmallVec<[MinorEdgeKey; 8]>,
 }
 
 #[derive(Debug, Clone)]
@@ -244,6 +250,23 @@ struct DualGraphHalfEdge {
 	twin: Option<DualEdgeKey>,
 }
 
+impl DualGraphHalfEdge {
+	fn start_segment(&self) -> PathSegment {
+		let segment = self.segments[0];
+		match self.direction_flag {
+			Direction::Forward => segment,
+			Direction::Backwards => segment.reverse(),
+		}
+	}
+
+	fn outer_boundnig_box(&self) -> Aabb {
+		self.segments
+			.iter()
+			.map(|seg| seg.approx_bounding_box())
+			.fold(Default::default(), |old, new| merge_bounding_boxes(&old, &new))
+	}
+}
+
 #[derive(Debug, Clone, PartialEq, Eq, Hash)]
 struct DualGraphVertex {
 	incident_edges: Vec<DualEdgeKey>,
@@ -258,6 +281,8 @@ struct DualGraphComponent {
 	edges: Vec<DualEdgeKey>,
 	vertices: Vec<DualVertexKey>,
 	outer_face: Option<DualVertexKey>,
+	inner_bb: Aabb,
+	outer_bb: Aabb,
 }
 
 /// Represents the dual graph of the MinorGraph.
@@ -334,7 +359,7 @@ fn dual_graph_to_dot(components: &[DualGraphComponent], edges: &SlotMap<DualEdge
 
 fn segment_to_edge(parent: u8) -> impl Fn(&PathSegment) -> Option<MajorGraphEdgeStage1> {
 	move |seg| {
-		if bounding_box_max_extent(&seg.bounding_box()) < EPS.point {
+		if bounding_box_max_extent(&seg.bounding_box()) < EPS.point / 2. {
 			return None;
 		}
 
@@ -404,31 +429,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<MajorGraphEdgeStage2>, Option<Aabb>) {
+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.bounding_box())).collect();
+	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(None, |acc, (_, _, bb)| Some(merge_bounding_boxes(acc, bb)));
+	let total_bounding_box = with_bounding_box.iter().fold(Default::default(), |acc, (_, _, bb)| merge_bounding_boxes(&acc, bb));
 
-	let total_bounding_box = match total_bounding_box {
+	let max_extent = bounding_box_max_extent(&total_bounding_box);
-		Some(bb) => bb,
+	let cell_size = max_extent / (edges.len() as f64).sqrt();
-		None => return (Vec::new(), None),
+
-	};
+	// 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::new();
+	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) {
+	fn add_split(splits_per_edge: &mut [Vec<f64>], i: usize, t: f64) {
-		splits_per_edge.entry(i).or_default().push(t);
+		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);
+		// let candidates = edge_tree.find(&edge.2);
-		for &j in &candidates {
+		// 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);
@@ -437,14 +473,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;
@@ -462,16 +500,16 @@ fn split_at_intersections(edges: &[MajorGraphEdgeStage1]) -> (Vec<MajorGraphEdge
 					continue;
 				}
 				let (seg1, seg2) = tmp_seg.split_at(tt);
-				new_edges.push((seg1, parent, seg1.bounding_box()));
+				new_edges.push((seg1, parent));
 				tmp_seg = seg2;
 			}
-			new_edges.push((tmp_seg, parent, tmp_seg.bounding_box()));
+			new_edges.push((tmp_seg, parent));
 		} else {
-			new_edges.push((seg, parent, seg.bounding_box()));
+			new_edges.push((seg, parent));
 		}
 	}
 
-	(new_edges, Some(total_bounding_box))
+	new_edges
 }
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
@@ -507,30 +545,41 @@ impl Direction {
 	}
 }
 
-// TODO:(@TrueDoctor) Optimize this by rounding each vertex up and down and then inserting them in a hashmap. This should remove the need for bbox calculations and the quad tree
+const ROUNDING_FACTOR: f64 = 1.0 / (2. * EPS.point);
-fn find_vertices(edges: &[MajorGraphEdgeStage2], bounding_box: Aabb) -> MajorGraph {
-	let mut vertex_tree = QuadTree::new(bounding_box, POINT_TREE_DEPTH, 8);
-	let mut graph = MajorGraph {
-		edges: SlotMap::with_key(),
-		vertices: SlotMap::with_key(),
-	};
 
-	let mut parents: HashMap<MajorEdgeKey, u8> = HashMap::new();
+fn round_point(point: DVec2) -> I64Vec2 {
-
+	(point * ROUNDING_FACTOR).round().as_i64vec2()
-	let mut vertex_pair_id_to_edges: HashMap<_, Vec<(MajorGraphEdgeStage2, MajorEdgeKey, MajorEdgeKey)>> = HashMap::new();
-
-	for (seg, parent, bounding_box) in edges {
-		let mut get_vertex = |point: DVec2| -> MajorVertexKey {
-			let box_around_point = bounding_box_around_point(point, EPS.point);
-			if let Some(&existing_vertex) = vertex_tree.find(&box_around_point).iter().next() {
-				existing_vertex
-			} else {
-				let vertex_key = graph.vertices.insert(MajorGraphVertex { point, outgoing_edges: Vec::new() });
-				vertex_tree.insert(box_around_point, vertex_key);
-				vertex_key
 }
+
+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()),
 	};
 
+	let mut vertex_pair_id_to_edges: HashMap<_, SmallVec<[(PathSegment, u8, MajorEdgeKey, MajorEdgeKey); 2]>> = new_hash_map(edges.len());
+	let mut vertex_hashmap: HashMap<I64Vec2, MajorVertexKey> = new_hash_map(edges.len() * 2);
+
+	for (seg, parent) in edges {
+		let mut get_vertex = |point: DVec2| -> MajorVertexKey {
+			let rounded = round_point(point);
+			for dx in -1..=1 {
+				for dy in -1..=1 {
+					let offset = I64Vec2::new(dx, dy);
+					if let Some(&vertex) = vertex_hashmap.get(&(rounded + offset)) {
+						return vertex;
+					}
+				}
+			}
+
+			let vertex_key = graph.vertices.insert(MajorGraphVertex {
+				point,
+				outgoing_edges: SmallVec::new(),
+			});
+			vertex_hashmap.insert(rounded, vertex_key);
+			vertex_key
+		};
+		// we should subtract the center instead here
 		let start_vertex = get_vertex(seg.start());
 		let end_vertex = get_vertex(seg.end());
 
@@ -556,10 +605,12 @@ fn find_vertices(edges: &[MajorGraphEdgeStage2], bounding_box: Aabb) -> MajorGra
 		if let Some(existing_edges) = vertex_pair_id_to_edges.get(&vertex_pair_id) {
 			if let Some(existing_edge) = existing_edges
 				.iter()
-				.find(|(other_seg, ..)| segments_equal(seg, &other_seg.0, EPS.point) || segments_equal(&seg.reverse(), &other_seg.0, EPS.point))
+				.find(|(other_seg, ..)| segments_equal(seg, other_seg, EPS.point) || segments_equal(&seg.reverse(), other_seg, EPS.point))
 			{
-				*parents.entry(existing_edge.1).or_default() |= parent;
+				if existing_edge.1 != *parent {
-				*parents.entry(existing_edge.2).or_default() |= parent;
+					graph.edges[existing_edge.2].parent = 0b11;
+					graph.edges[existing_edge.3].parent = 0b11;
+				}
 				continue;
 			}
 		}
@@ -585,13 +636,7 @@ fn find_vertices(edges: &[MajorGraphEdgeStage2], bounding_box: Aabb) -> MajorGra
 		graph.vertices[start_vertex].outgoing_edges.push(fwd_edge_key);
 		graph.vertices[end_vertex].outgoing_edges.push(bwd_edge_key);
 
-		vertex_pair_id_to_edges
+		vertex_pair_id_to_edges.entry(vertex_pair_id).or_default().push((*seg, *parent, fwd_edge_key, bwd_edge_key));
-			.entry(vertex_pair_id)
-			.or_default()
-			.push(((*seg, *parent, *bounding_box), fwd_edge_key, bwd_edge_key));
-	}
-	for (edge_key, parent) in parents {
-		graph.edges[edge_key].parent |= parent;
 	}
 
 	graph
@@ -624,11 +669,14 @@ fn get_order(vertex: &MajorGraphVertex) -> usize {
 ///
 /// A new MinorGraph representing the simplified structure.
 fn compute_minor(major_graph: &MajorGraph) -> MinorGraph {
-	let mut new_edges = SlotMap::with_key();
+	let vertex_count = major_graph.vertices.len();
-	let mut new_vertices = SlotMap::with_key();
+	let edge_count = major_graph.edges.len();
-	let mut to_minor_vertex = HashMap::new();
+	let mut new_edges = SlotMap::with_capacity_and_key(edge_count / 2);
-	let mut id_to_edge = HashMap::new();
+	let mut new_vertices = SlotMap::with_capacity_and_key(vertex_count.ilog2() as usize + 2);
-	let mut visited = HashSet::new();
+	let mut to_minor_vertex = new_hash_map(vertex_count);
+	let mut id_to_edge = new_hash_map(edge_count);
+	// merge with to_minor_vertex
+	let mut visited = HashSet::with_capacity_and_hasher(vertex_count, Default::default());
 
 	// Handle components that are not cycles
 	for (major_vertex_key, vertex) in &major_graph.vertices {
@@ -638,10 +686,10 @@ fn compute_minor(major_graph: &MajorGraph) -> MinorGraph {
 		}
 		let start_vertex = *to_minor_vertex
 			.entry(major_vertex_key)
-			.or_insert_with(|| new_vertices.insert(MinorGraphVertex { outgoing_edges: Vec::new() }));
+			.or_insert_with(|| new_vertices.insert(MinorGraphVertex { outgoing_edges: SmallVec::new() }));
 
 		for &start_edge_key in &vertex.outgoing_edges {
-			let mut segments = Vec::new();
+			let mut segments = SmallVec::new();
 			let mut edge_key = start_edge_key;
 			let mut edge = &major_graph.edges[edge_key];
 
@@ -660,7 +708,7 @@ fn compute_minor(major_graph: &MajorGraph) -> MinorGraph {
 
 			let end_vertex = *to_minor_vertex
 				.entry(edge.incident_vertices[1])
-				.or_insert_with(|| new_vertices.insert(MinorGraphVertex { outgoing_edges: Vec::new() }));
+				.or_insert_with(|| new_vertices.insert(MinorGraphVertex { outgoing_edges: SmallVec::new() }));
 			assert!(major_graph.edges[start_edge_key].twin.is_some());
 			assert!(edge.twin.is_some());
 
@@ -693,7 +741,7 @@ fn compute_minor(major_graph: &MajorGraph) -> MinorGraph {
 		let mut edge_key = vertex.outgoing_edges[0];
 		let mut edge = &major_graph.edges[edge_key];
 		let mut cycle = MinorGraphCycle {
-			segments: Vec::new(),
+			segments: Vec::with_capacity(4),
 			parent: edge.parent,
 			direction_flag: edge.direction_flag,
 		};
@@ -721,8 +769,10 @@ fn compute_minor(major_graph: &MajorGraph) -> MinorGraph {
 fn remove_dangling_edges(graph: &mut MinorGraph) {
 	// Basically DFS for each parent with BFS number
 	fn walk(parent: u8, graph: &MinorGraph) -> HashSet<MinorVertexKey> {
-		let mut kept_vertices = HashSet::new();
+		// merge
-		let mut vertex_to_level = HashMap::new();
+		let vertex_count = graph.vertices.len();
+		let mut kept_vertices = HashSet::with_capacity_and_hasher(vertex_count, Default::default());
+		let mut vertex_to_level = new_hash_map(vertex_count);
 
 		fn visit(
 			vertex: MinorVertexKey,
@@ -768,8 +818,8 @@ fn remove_dangling_edges(graph: &mut MinorGraph) {
 	graph.vertices.retain(|k, _| kept_vertices_a.contains(&k) || kept_vertices_b.contains(&k));
 
 	for vertex in graph.vertices.values_mut() {
-		vertex.outgoing_edges.retain(|&edge_key| {
+		vertex.outgoing_edges.retain(|edge_key| {
-			let edge = &graph.edges[edge_key];
+			let edge = &graph.edges[*edge_key];
 			(edge.parent & 1 == 1 && kept_vertices_a.contains(&edge.incident_vertices[0]) && kept_vertices_a.contains(&edge.incident_vertices[1]))
 				|| (edge.parent & 2 == 2 && kept_vertices_b.contains(&edge.incident_vertices[0]) && kept_vertices_b.contains(&edge.incident_vertices[1]))
 		});
@@ -786,7 +836,7 @@ fn sort_outgoing_edges_by_angle(graph: &mut MinorGraph) {
 		if vertex.outgoing_edges.len() > 2 {
 			vertex.outgoing_edges.sort_by(|&a, &b| graph.edges[a].partial_cmp(&graph.edges[b]).unwrap());
 			if cfg!(feature = "logging") {
-				eprintln!("Outgoing edges for {:?}:", vertex_key);
+				eprintln!("Outgoing edges for {vertex_key:?}:");
 				for &edge_key in &vertex.outgoing_edges {
 					let edge = &graph.edges[edge_key];
 					let angle = edge.start_segment().start_angle();
@@ -797,6 +847,7 @@ fn sort_outgoing_edges_by_angle(graph: &mut MinorGraph) {
 	}
 }
 
+#[cfg(feature = "logging")]
 fn face_to_polygon(face: &DualGraphVertex, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>) -> Vec<DVec2> {
 	const CNT: usize = 3;
 
@@ -829,6 +880,7 @@ fn line_segment_intersects_horizontal_ray(a: DVec2, b: DVec2, point: DVec2) -> b
 	x >= point.x
 }
 
+#[cfg(feature = "logging")]
 fn compute_point_winding(polygon: &[DVec2], tested_point: DVec2) -> i32 {
 	if polygon.len() <= 2 {
 		return 0;
@@ -844,6 +896,7 @@ fn compute_point_winding(polygon: &[DVec2], tested_point: DVec2) -> i32 {
 	winding
 }
 
+#[cfg(feature = "logging")]
 fn compute_winding(face: &DualGraphVertex, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>) -> Option<i32> {
 	let polygon = face_to_polygon(face, edges);
 
@@ -862,28 +915,49 @@ fn compute_winding(face: &DualGraphVertex, edges: &SlotMap<DualEdgeKey, DualGrap
 }
 
 fn compute_signed_area(face: &DualGraphVertex, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>) -> f64 {
-	let polygon = face_to_polygon(face, edges);
+	const CNT: usize = 3;
-	if polygon.len() <= 4 {
+	let mut area = 0.0;
-		return -1.;
+	let mut prev_point: Option<DVec2> = None;
+	let mut start_point = None;
+	let mut point_count = 0;
+
+	for &edge_key in &face.incident_edges {
+		let edge = &edges[edge_key];
+		for seg in &edge.segments {
+			for i in 0..CNT {
+				let t0 = i as f64 / CNT as f64;
+				let t = if edge.direction_flag.forward() { t0 } else { 1. - t0 };
+				let current_point = seg.sample_at(t);
+
+				if let Some(prev) = prev_point {
+					area += prev.x * current_point.y;
+					area -= current_point.x * prev.y;
+				} else {
+					start_point = Some(current_point);
+				}
+
+				prev_point = Some(current_point);
+				point_count += 1;
+			}
+		}
+	}
+
+	if point_count <= 4 {
+		return -f64::EPSILON;
+	}
+
+	// Close the polygon
+	if let (Some(start), Some(end)) = (start_point, prev_point) {
+		area += end.x * start.y;
+		area -= start.x * end.y;
 	}
 
 	#[cfg(feature = "logging")]
+	{
 		eprintln!("vertex: {:?}", face);
-	#[cfg(feature = "logging")]
-	for point in &polygon {
-		eprintln!("{}, {}", point.x, point.y);
-	}
-	let mut area = 0.;
-
-	for i in 0..polygon.len() {
-		let a = polygon[i];
-		let b = polygon[(i + 1) % polygon.len()];
-		area += a.x * b.y;
-		area -= b.x * a.y;
-	}
-
-	#[cfg(feature = "logging")]
 		eprintln!("winding: {}", area);
+	}
+
 	area
 }
 
@@ -912,8 +986,9 @@ fn compute_signed_area(face: &DualGraphVertex, edges: &SlotMap<DualEdgeKey, Dual
 /// operation cannot be completed successfully.
 fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 	let mut new_vertices: Vec<DualVertexKey> = Vec::new();
-	let mut minor_to_dual_edge: HashMap<MinorEdgeKey, DualEdgeKey> = HashMap::new();
+	let edge_count = minor_graph.edges.len();
-	let mut dual_edges = SlotMap::with_key();
+	let mut minor_to_dual_edge: HashMap<MinorEdgeKey, DualEdgeKey> = new_hash_map(edge_count);
+	let mut dual_edges = SlotMap::with_capacity_and_key(edge_count);
 	let mut dual_vertices = SlotMap::with_key();
 
 	for (start_edge_key, start_edge) in &minor_graph.edges {
@@ -935,7 +1010,7 @@ fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 			let twin_dual_key = minor_to_dual_edge.get(&twin).copied();
 
 			let new_edge_key = dual_edges.insert(DualGraphHalfEdge {
-				segments: edge.segments.clone(),
+				segments: edge.segments.to_vec(),
 				parent: edge.parent,
 				incident_vertex: face_key,
 				direction_flag: edge.direction_flag,
@@ -990,9 +1065,9 @@ fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 		new_vertices.push(outer_face_key);
 	}
 
-	let mut components = Vec::new();
+	let mut components = Vec::with_capacity(12);
-	let mut visited_vertices = HashSet::new();
+	let mut visited_vertices = HashSet::with_capacity_and_hasher(dual_vertices.len(), Default::default());
-	let mut visited_edges = HashSet::new();
+	let mut visited_edges = HashSet::with_capacity_and_hasher(edge_count, Default::default());
 
 	if cfg!(feature = "logging") {
 		eprintln!("faces: {}, dual-edges: {}, cycles: {}", new_vertices.len(), dual_edges.len(), minor_graph.cycles.len())
@@ -1046,14 +1121,17 @@ fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 		#[cfg(feature = "logging")]
 		eprintln!("component_vertices: {}", component_vertices.len());
 
+		#[cfg(feature = "logging")]
 		let windings: Option<Vec<_>> = component_vertices
 			.iter()
 			.map(|face_key| compute_winding(&dual_vertices[*face_key], &dual_edges).map(|w| (face_key, w)))
 			.collect();
-		let Some(windings) = windings else {
+		#[cfg(feature = "logging")]
+		let Some(_) = windings else {
 			return Err(BooleanError::NoEarInPolygon);
 		};
 
+		#[cfg(feature = "logging")]
 		let areas: Vec<_> = component_vertices
 			.iter()
 			.map(|face_key| (face_key, compute_signed_area(&dual_vertices[*face_key], &dual_edges)))
@@ -1070,31 +1148,38 @@ fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 						vertices: component_vertices.clone(),
 						edges: component_edges.clone(),
 						outer_face: None,
+						inner_bb: Default::default(),
+						outer_bb: Default::default(),
 					}],
 					&dual_edges,
 				)
 			);
 		}
 
-		let mut count = windings.iter().filter(|(_, winding)| winding < &0).count();
+		let areas: Vec<_> = component_vertices
-		let mut reverse_winding = false;
+			.iter()
-		// If the paths are reversed use positive winding as outer face
+			.map(|face_key| (face_key, (compute_signed_area(&dual_vertices[*face_key], &dual_edges) * 1000.) as i64))
-		if windings.len() > 2 && count == windings.len() - 1 {
+			.collect();
-			count = 1;
-			reverse_winding = true;
-		}
-		let outer_face_key = if count != 1 {
 		#[cfg(feature = "logging")]
 		eprintln!("Found multiple outer faces: {areas:?}, falling back to area calculation");
-			let (key, _) = *areas.iter().max_by_key(|(_, area)| (area.abs() * 1000.) as u64).unwrap();
+		let (&outer_face_key, _) = *areas
-			*key
-		} else {
-			*windings
 			.iter()
-				.find(|&&(&_, ref winding)| (winding < &0) ^ reverse_winding)
+			.max_by(|(_, a1), (_, a2)| match a1.abs().cmp(&a2.abs()) {
-				.expect("No outer face of a component found.")
+				Ordering::Equal => a1.cmp(a2),
-				.0
+				ord => ord,
-		};
+			})
+			.unwrap();
+		let inner_bb = dual_vertices[outer_face_key]
+			.incident_edges
+			.iter()
+			.map(|edge_key| dual_edges[*edge_key].start_segment().start())
+			.fold(Default::default(), |bbox, point| extend_bounding_box(Some(bbox), point));
+		let outer_bb = dual_vertices[outer_face_key]
+			.incident_edges
+			.iter()
+			.map(|edge_key| dual_edges[*edge_key].outer_boundnig_box())
+			.reduce(|acc, new| merge_bounding_boxes(&acc, &new))
+			.unwrap_or_default();
 		#[cfg(feature = "logging")]
 		dbg!(outer_face_key);
 
@@ -1102,6 +1187,8 @@ fn compute_dual(minor_graph: &MinorGraph) -> Result<DualGraph, BooleanError> {
 			vertices: component_vertices,
 			edges: component_edges,
 			outer_face: Some(outer_face_key),
+			inner_bb,
+			outer_bb,
 		});
 	}
 
@@ -1122,6 +1209,9 @@ fn get_next_edge(edge_key: MinorEdgeKey, graph: &MinorGraph) -> MinorEdgeKey {
 }
 
 fn test_inclusion(a: &DualGraphComponent, b: &DualGraphComponent, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>, vertices: &SlotMap<DualVertexKey, DualGraphVertex>) -> Option<DualVertexKey> {
+	if a.inner_bb.min().cmpge(b.outer_bb.min()) & a.inner_bb.max().cmple(b.outer_bb.max()) != BVec2::TRUE {
+		return None;
+	}
 	let tested_point = edges[a.edges[0]].segments[0].start();
 	for (face_key, face) in b.vertices.iter().map(|&key| (key, &vertices[key])) {
 		if Some(face_key) == b.outer_face {
@@ -1141,30 +1231,86 @@ fn test_inclusion(a: &DualGraphComponent, b: &DualGraphComponent, edges: &SlotMa
 	None
 }
 fn bounding_box_intersects_horizontal_ray(bounding_box: &Aabb, point: DVec2) -> bool {
-	interval_crosses_point(bounding_box.top, bounding_box.bottom, point[1]) && bounding_box.right >= point[0]
+	bounding_box.right() >= point[0] && (bounding_box.top()..bounding_box.bottom()).contains(&point[1])
 }
 
+#[derive(Copy, Clone)]
 struct IntersectionSegment {
 	bounding_box: Aabb,
 	seg: PathSegment,
 }
 
 pub fn path_segment_horizontal_ray_intersection_count(orig_seg: &PathSegment, point: DVec2) -> usize {
-	let total_bounding_box = orig_seg.bounding_box();
+	let total_bounding_box = orig_seg.approx_bounding_box();
-
 	if !bounding_box_intersects_horizontal_ray(&total_bounding_box, point) {
 		return 0;
 	}
 
-	let mut segments = vec![IntersectionSegment {
+	match orig_seg {
-		bounding_box: total_bounding_box,
+		PathSegment::Cubic(..) => cubic_bezier_horizontal_ray_intersection_count(orig_seg, point),
-		seg: *orig_seg,
+		_ => fallback_intersection_count(orig_seg, total_bounding_box, point),
-	}];
+	}
+}
+
+fn cubic_bezier_horizontal_ray_intersection_count(cubic: &PathSegment, point: DVec2) -> usize {
+	let y = point.y;
+	let PathSegment::Cubic(p0, p1, p2, p3) = cubic else { unreachable!() };
+
+	// Transform the curve so that the horizontal line is at y = 0
+	let a = -p0.y + 3.0 * p1.y - 3.0 * p2.y + p3.y;
+	let b = 3.0 * p0.y - 6.0 * p1.y + 3.0 * p2.y;
+	let c = -3.0 * p0.y + 3.0 * p1.y;
+	let d = p0.y - y;
+
+	let roots = find_roots_cubic(a, b, c, d);
+
 	let mut count = 0;
+	match roots {
+		Roots::Three(roots) => {
+			for &t in roots.iter() {
+				if (0.0..=1.0).contains(&t) {
+					let x = cubic.sample_at(t).x;
+					if x > point.x {
+						count += 1;
+					}
+				}
+			}
+		}
+		Roots::Two(roots) => {
+			for &t in roots.iter() {
+				if (0.0..=1.0).contains(&t) {
+					let x = cubic.sample_at(t).x;
+					if x > point.x {
+						count += 1;
+					}
+				}
+			}
+		}
+		Roots::One(roots) => {
+			for &t in roots.iter() {
+				if (0.0..=1.0).contains(&t) {
+					let x = cubic.sample_at(t).x;
+					if x > point.x {
+						count += 1;
+					}
+				}
+			}
+		}
+		_ => {}
+	}
+
+	count
+}
+
+fn fallback_intersection_count(orig_seg: &PathSegment, bounding_box: Aabb, point: DVec2) -> usize {
+	// Existing implementation for non-cubic segments
+	let mut segments = vec![IntersectionSegment { bounding_box, seg: *orig_seg }];
+	let mut count = 0;
+	let mut next_segments = Vec::new();
 
 	while !segments.is_empty() {
-		let mut next_segments = Vec::new();
+		next_segments.clear();
-		for segment in segments {
+		for segment in segments.iter() {
 			if bounding_box_max_extent(&segment.bounding_box) < EPS.linear {
 				if line_segment_intersects_horizontal_ray(segment.seg.start(), segment.seg.end(), point) {
 					count += 1;
@@ -1173,7 +1319,6 @@ pub fn path_segment_horizontal_ray_intersection_count(orig_seg: &PathSegment, po
 				let split = &segment.seg.split_at(0.5);
 				let bounding_box0 = split.0.bounding_box();
 				let bounding_box1 = split.1.bounding_box();
-
 				if bounding_box_intersects_horizontal_ray(&bounding_box0, point) {
 					next_segments.push(IntersectionSegment {
 						bounding_box: bounding_box0,
@@ -1188,9 +1333,8 @@ pub fn path_segment_horizontal_ray_intersection_count(orig_seg: &PathSegment, po
 				}
 			}
 		}
-		segments = next_segments;
+		std::mem::swap(&mut next_segments, &mut segments);
 	}
-
 	count
 }
 
@@ -1216,27 +1360,27 @@ pub fn path_segment_horizontal_ray_intersection_count(orig_seg: &PathSegment, po
 /// # Returns
 ///
 /// A vector of NestingTree structures representing the top-level components and their nested subcomponents.
-fn compute_nesting_tree(DualGraph { components, vertices, edges }: &DualGraph) -> Vec<NestingTree> {
+fn compute_nesting_tree(DualGraph { components, vertices, edges }: &mut DualGraph) -> Vec<NestingTree> {
-	let mut nesting_trees = Vec::new();
+	let mut nesting_trees = Vec::with_capacity(4);
 
-	for component in components {
+	for component in components.drain(..) {
 		insert_component(&mut nesting_trees, component, edges, vertices);
 	}
 
 	nesting_trees
 }
 
-fn insert_component(trees: &mut Vec<NestingTree>, component: &DualGraphComponent, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>, vertices: &SlotMap<DualVertexKey, DualGraphVertex>) {
+fn insert_component(trees: &mut Vec<NestingTree>, component: DualGraphComponent, edges: &SlotMap<DualEdgeKey, DualGraphHalfEdge>, vertices: &SlotMap<DualVertexKey, DualGraphVertex>) {
 	for tree in trees.iter_mut() {
-		if let Some(face_key) = test_inclusion(component, &tree.component, edges, vertices) {
+		if let Some(face_key) = test_inclusion(&component, &tree.component, edges, vertices) {
 			if let Some(children) = tree.outgoing_edges.get_mut(&face_key) {
 				insert_component(children, component, edges, vertices);
 			} else {
 				tree.outgoing_edges.insert(
 					face_key,
 					vec![NestingTree {
-						component: component.clone(),
+						component,
-						outgoing_edges: HashMap::new(),
+						outgoing_edges: HashMap::default(),
 					}],
 				);
 			}
@@ -1245,15 +1389,14 @@ fn insert_component(trees: &mut Vec<NestingTree>, component: &DualGraphComponent
 	}
 
 	let mut new_tree = NestingTree {
-		component: component.clone(),
+		component,
-		outgoing_edges: HashMap::new(),
+		outgoing_edges: HashMap::default(),
 	};
 
 	let mut i = 0;
 	while i < trees.len() {
 		if let Some(face_key) = test_inclusion(&trees[i].component, &new_tree.component, edges, vertices) {
-			// TODO: (@TrueDoctor) use swap remove
+			let tree = trees.swap_remove(i);
-			let tree = trees.remove(i);
 			new_tree.outgoing_edges.entry(face_key).or_default().push(tree);
 		} else {
 			i += 1;
@@ -1288,12 +1431,15 @@ fn flag_faces(
 	vertices: &SlotMap<DualVertexKey, DualGraphVertex>,
 	flags: &mut HashMap<DualVertexKey, u8>,
 ) {
+	let mut visited_faces = HashSet::default();
+	let mut face_stack = VecDeque::new();
 	for tree in nesting_trees.iter() {
 		let mut tree_stack = vec![(tree, 0, 0)];
 
 		while let Some((current_tree, a_running_count, b_running_count)) = tree_stack.pop() {
-			let mut visited_faces = HashSet::new();
+			// TODO: Test if clearing is faster
-			let mut face_stack = VecDeque::new();
+			visited_faces.clear();
+			face_stack.clear();
 
 			let outer_face_key = current_tree.component.outer_face.expect("Component doesn't have an outer face.");
 			face_stack.push_back((outer_face_key, a_running_count, b_running_count));
@@ -1346,7 +1492,7 @@ fn walk_faces<'a>(faces: &'a [DualVertexKey], edges: &SlotMap<DualEdgeKey, DualG
 	// TODO: Try using a binary search to avoid the hashset construction
 	let is_removed_edge = |edge: &DualGraphHalfEdge| face_set.contains(&edge.incident_vertex) == face_set.contains(&edges[edge.twin.unwrap()].incident_vertex);
 
-	let mut edge_to_next = HashMap::new();
+	let mut edge_to_next = HashMap::with_capacity_and_hasher(edges.len(), Default::default());
 	for face_key in faces {
 		let face = &vertices[*face_key];
 		let mut prev_edge = *face.incident_edges.last().unwrap();
@@ -1356,8 +1502,8 @@ fn walk_faces<'a>(faces: &'a [DualVertexKey], edges: &SlotMap<DualEdgeKey, DualG
 		}
 	}
 
-	let mut visited_edges = HashSet::new();
+	let mut visited_edges = HashSet::with_capacity_and_hasher(edges.len(), Default::default());
-	let mut result = Vec::new();
+	let mut result = Vec::with_capacity(edges.len());
 
 	for &face_key in faces {
 		let face = &vertices[face_key];
@@ -1544,24 +1690,23 @@ impl Display for BooleanError {
 /// - Input paths are invalid or cannot be processed.
 /// - The operation encounters an unsolvable geometric configuration.
 /// - Issues arise in determining the nesting structure of the paths.
+#[inline(never)]
 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() {
+	for (edge, _) in split_edges.iter() {
 		eprintln!("{}", path_to_path_data(&vec![*edge], 0.001));
 	}
 
-	let total_bounding_box = match total_bounding_box {
+	let major_graph = find_vertices(&split_edges);
-		Some(bb) => bb,
-		None => return Ok(Vec::new()), // Input geometry is empty
-	};
-
-	let major_graph = find_vertices(&split_edges, total_bounding_box);
 
 	#[cfg(feature = "logging")]
 	eprintln!("Major graph:");
@@ -1583,33 +1728,33 @@ pub fn path_boolean(a: &Path, a_fill_rule: FillRule, b: &Path, b_fill_rule: Fill
 
 	#[cfg(feature = "logging")]
 	for (key, edge) in minor_graph.edges.iter() {
-		eprintln!("{key:?}:\n{}", path_to_path_data(&edge.segments, 0.001));
+		eprintln!("{key:?}:\n{}", path_to_path_data(&edge.segments.to_vec(), 0.001));
 	}
 	#[cfg(feature = "logging")]
 	for vertex in minor_graph.vertices.values() {
-		eprintln!("{:?}", vertex);
+		eprintln!("{vertex:?}");
 	}
 	sort_outgoing_edges_by_angle(&mut minor_graph);
 	#[cfg(feature = "logging")]
 	for vertex in minor_graph.vertices.values() {
-		eprintln!("{:?}", vertex);
+		eprintln!("{vertex:?}");
 	}
 
 	for (edge_key, edge) in &minor_graph.edges {
-		assert!(minor_graph.vertices.contains_key(edge.incident_vertices[0]), "Edge {:?} has invalid start vertex", edge_key);
+		assert!(minor_graph.vertices.contains_key(edge.incident_vertices[0]), "Edge {edge_key:?} has invalid start vertex");
-		assert!(minor_graph.vertices.contains_key(edge.incident_vertices[1]), "Edge {:?} has invalid end vertex", edge_key);
+		assert!(minor_graph.vertices.contains_key(edge.incident_vertices[1]), "Edge {edge_key:?} has invalid end vertex");
-		assert!(edge.twin.is_some(), "Edge {:?} should have a twin", edge_key);
+		assert!(edge.twin.is_some(), "Edge {edge_key:?} should have a twin");
 		let twin = &minor_graph.edges[edge.twin.unwrap()];
-		assert_eq!(twin.twin.unwrap(), edge_key, "Twin relationship should be symmetrical for edge {:?}", edge_key);
+		assert_eq!(twin.twin.unwrap(), edge_key, "Twin relationship should be symmetrical for edge {edge_key:?}");
 	}
 
-	let dual_graph = compute_dual(&minor_graph)?;
+	let mut dual_graph = compute_dual(&minor_graph)?;
 
-	let nesting_trees = compute_nesting_tree(&dual_graph);
+	let nesting_trees = compute_nesting_tree(&mut dual_graph);
 
 	#[cfg(feature = "logging")]
 	for tree in &nesting_trees {
-		eprintln!("nesting_trees: {:?}", tree);
+		eprintln!("nesting_trees: {tree:?}");
 	}
 
 	let DualGraph { edges, vertices, .. } = &dual_graph;
@@ -1619,7 +1764,7 @@ pub fn path_boolean(a: &Path, a_fill_rule: FillRule, b: &Path, b_fill_rule: Fill
 	#[cfg(feature = "logging")]
 	eprintln!("{}", dual_graph_to_dot(&dual_graph.components, edges));
 
-	let mut flags = HashMap::new();
+	let mut flags = new_hash_map(vertices.len());
 	flag_faces(&nesting_trees, a_fill_rule, b_fill_rule, edges, vertices, &mut flags);
 
 	#[cfg(feature = "logging")]
@@ -1648,20 +1793,11 @@ mod tests {
 	#[test]
 	fn test_split_at_intersections() {
 		let unsplit_edges = unsplit_edges();
-		let (split_edges, total_bounding_box) = split_at_intersections(&unsplit_edges);
+		let split_edges = split_at_intersections(&unsplit_edges);
-
-		// Check that we have a valid bounding box
-		assert!(total_bounding_box.is_some());
 
 		// Check that we have more edges after splitting (due to intersections)
 		assert!(split_edges.len() >= unsplit_edges.len());
 
-		// Check that all edges have a valid bounding box
-		for (_, _, bb) in &split_edges {
-			assert!(bb.left <= bb.right);
-			assert!(bb.top <= bb.bottom);
-		}
-
 		// You might want to add more specific checks based on the expected behavior
 		// of your split_at_intersections function
 	}
@@ -1684,8 +1820,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 split_edges = split_at_intersections(&unsplit_edges);
-		let major_graph = find_vertices(&split_edges, total_bounding_box.unwrap());
+		let major_graph = find_vertices(&split_edges);
 
 		// Compute minor graph
 		let minor_graph = compute_minor(&major_graph);
@@ -1739,8 +1875,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 split_edges = split_at_intersections(&unsplit_edges);
-		let major_graph = find_vertices(&split_edges, total_bounding_box.unwrap());
+		let major_graph = find_vertices(&split_edges);
 		let mut minor_graph = compute_minor(&major_graph);
 
 		// Print initial state
@@ -1748,7 +1884,7 @@ mod tests {
 		print_minor_graph_state(&minor_graph);
 
 		// Store initial edge order
-		let initial_edge_order: HashMap<MinorVertexKey, Vec<MinorEdgeKey>> = minor_graph.vertices.iter().map(|(k, v)| (k, v.outgoing_edges.clone())).collect();
+		let initial_edge_order: HashMap<MinorVertexKey, SmallVec<[MinorEdgeKey; 8]>> = minor_graph.vertices.iter().map(|(k, v)| (k, v.outgoing_edges.clone())).collect();
 
 		// Apply sort_outgoing_edges_by_angle
 		sort_outgoing_edges_by_angle(&mut minor_graph);
@@ -1821,15 +1957,10 @@ mod tests {
 
 	#[test]
 	fn test_bounding_box_intersects_horizontal_ray() {
-		let bbox = Aabb {
+		let bbox = Aabb::new(20., 10., 40., 30.);
-			top: 10.,
-			right: 40.,
-			bottom: 30.,
-			left: 20.,
-		};
 
-		assert!(bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(0., 30.)));
+		assert!(bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(0., 29.)));
-		assert!(bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(20., 30.)));
+		assert!(bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(20., 29.)));
 		assert!(bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(10., 20.)));
 		assert!(!bounding_box_intersects_horizontal_ray(&bbox, DVec2::new(30., 40.)));
 	}

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

@@ -1,64 +1,92 @@
-use glam::DVec2;
+use glam::{BVec2, DVec2};
 
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub(crate) struct Aabb {
-	pub top: f64,
+	min: DVec2,
-	pub right: f64,
+	max: DVec2,
-	pub bottom: f64,
-	pub left: f64,
 }
 
+impl Default for Aabb {
+	fn default() -> Self {
+		Self {
+			min: DVec2::INFINITY,
+			max: DVec2::NEG_INFINITY,
+		}
+	}
+}
+
+impl Aabb {
+	#[inline]
+	pub(crate) fn min(&self) -> DVec2 {
+		self.min
+	}
+	#[inline]
+	pub(crate) fn max(&self) -> DVec2 {
+		self.max
+	}
+
+	pub(crate) const fn new(left: f64, top: f64, right: f64, bottom: f64) -> Self {
+		Aabb {
+			min: DVec2::new(left, top),
+			max: DVec2::new(right, bottom),
+		}
+	}
+	#[inline]
+	pub(crate) fn top(&self) -> f64 {
+		self.min.y
+	}
+	#[inline]
+	pub(crate) fn left(&self) -> f64 {
+		self.min.x
+	}
+	#[inline]
+	pub(crate) fn right(&self) -> f64 {
+		self.max.x
+	}
+	#[inline]
+	pub(crate) fn bottom(&self) -> f64 {
+		self.max.y
+	}
+}
+
+#[inline]
 pub(crate) fn bounding_boxes_overlap(a: &Aabb, b: &Aabb) -> bool {
-	a.left <= b.right && b.left <= a.right && a.top <= b.bottom && b.top <= a.bottom
+	(a.min.cmple(b.max) & b.min.cmple(a.max)) == BVec2::TRUE
 }
 
-pub(crate) fn merge_bounding_boxes(a: Option<Aabb>, b: &Aabb) -> Aabb {
+#[inline]
-	match a {
+pub(crate) fn merge_bounding_boxes(a: &Aabb, b: &Aabb) -> Aabb {
-		Some(a) => Aabb {
+	Aabb {
-			top: a.top.min(b.top),
+		min: a.min.min(b.min),
-			right: a.right.max(b.right),
+		max: a.max.max(b.max),
-			bottom: a.bottom.max(b.bottom),
-			left: a.left.min(b.left),
-		},
-		None => *b,
 	}
 }
 
+#[inline]
 pub(crate) fn extend_bounding_box(bounding_box: Option<Aabb>, point: DVec2) -> Aabb {
 	match bounding_box {
 		Some(bb) => Aabb {
-			top: bb.top.min(point.y),
+			min: bb.min.min(point),
-			right: bb.right.max(point.x),
+			max: bb.max.max(point),
-			bottom: bb.bottom.max(point.y),
-			left: bb.left.min(point.x),
-		},
-		None => Aabb {
-			top: point.y,
-			right: point.x,
-			bottom: point.y,
-			left: point.x,
 		},
+		None => Aabb { min: point, max: point },
 	}
 }
 
 pub(crate) fn bounding_box_max_extent(bounding_box: &Aabb) -> f64 {
-	(bounding_box.right - bounding_box.left).max(bounding_box.bottom - bounding_box.top)
+	(bounding_box.max - bounding_box.min).max_element()
 }
 
 pub(crate) fn bounding_box_around_point(point: DVec2, padding: f64) -> Aabb {
 	Aabb {
-		top: point.y - padding,
+		min: point - DVec2::splat(padding),
-		right: point.x + padding,
+		max: point + DVec2::splat(padding),
-		bottom: point.y + padding,
-		left: point.x - padding,
 	}
 }
 
 pub(crate) fn expand_bounding_box(bounding_box: &Aabb, padding: f64) -> Aabb {
 	Aabb {
-		top: bounding_box.top - padding,
+		min: bounding_box.min - DVec2::splat(padding),
-		right: bounding_box.right + padding,
+		max: bounding_box.max + DVec2::splat(padding),
-		bottom: bounding_box.bottom + padding,
-		left: bounding_box.left - padding,
 	}
 }

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

@@ -0,0 +1,128 @@
+use crate::aabb::Aabb;
+use glam::{DVec2, IVec2};
+use rustc_hash::FxHashMap;
+use smallvec::SmallVec;
+
+pub(crate) struct Grid {
+	cell_factor: f64,
+	cells: FxHashMap<IVec2, SmallVec<[usize; 6]>>,
+}
+
+impl Grid {
+	pub(crate) fn new(cell_size: f64, edges: usize) -> Self {
+		Grid {
+			cell_factor: cell_size.recip(),
+			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_floor(bbox.min());
+		let max_cell = self.point_to_cell_ceil(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_floor(bbox.min());
+		let max_cell = self.point_to_cell_ceil(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_ceil(&self, point: DVec2) -> IVec2 {
+		(point * self.cell_factor).ceil().as_ivec2()
+	}
+	fn point_to_cell_floor(&self, point: DVec2) -> IVec2 {
+		(point * self.cell_factor).floor().as_ivec2()
+	}
+}
+
+pub struct BitVec {
+	data: Vec<u64>,
+}
+
+impl BitVec {
+	pub fn new(capacity: usize) -> Self {
+		let num_words = capacity.div_ceil(64);
+		BitVec { data: vec![0; num_words] }
+	}
+
+	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> {
+		loop {
+			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);
+		}
+	}
+}
+
+#[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/quad_tree.rs

@@ -1,121 +0,0 @@
-use crate::aabb::Aabb;
-use std::collections::HashSet;
-
-pub struct QuadTree<T> {
-	bounding_box: Aabb,
-	depth: usize,
-	inner_node_capacity: usize,
-	subtrees: Option<Box<[QuadTree<T>; 4]>>,
-	pairs: Vec<(Aabb, T)>,
-}
-
-impl<T: Clone> QuadTree<T> {
-	pub fn new(bounding_box: Aabb, depth: usize, inner_node_capacity: usize) -> Self {
-		QuadTree {
-			bounding_box,
-			depth,
-			inner_node_capacity,
-			subtrees: None,
-			pairs: Vec::new(),
-		}
-	}
-
-	pub fn insert(&mut self, bounding_box: Aabb, value: T) -> bool {
-		if !crate::aabb::bounding_boxes_overlap(&bounding_box, &self.bounding_box) {
-			return false;
-		}
-
-		if self.depth > 0 && self.pairs.len() >= self.inner_node_capacity {
-			self.ensure_subtrees();
-			for tree in self.subtrees.as_mut().unwrap().iter_mut() {
-				tree.insert(bounding_box, value.clone());
-			}
-		} else {
-			self.pairs.push((bounding_box, value));
-		}
-
-		true
-	}
-
-	pub fn find(&self, bounding_box: &Aabb) -> HashSet<T>
-	where
-		T: Eq + std::hash::Hash + Clone,
-	{
-		let mut set = HashSet::new();
-		self.find_internal(bounding_box, &mut set);
-		set
-	}
-
-	fn find_internal(&self, bounding_box: &Aabb, set: &mut HashSet<T>)
-	where
-		T: Eq + std::hash::Hash + Clone,
-	{
-		if !crate::aabb::bounding_boxes_overlap(bounding_box, &self.bounding_box) {
-			return;
-		}
-
-		for (key, value) in &self.pairs {
-			if crate::aabb::bounding_boxes_overlap(bounding_box, key) {
-				set.insert(value.clone());
-			}
-		}
-
-		if let Some(subtrees) = &self.subtrees {
-			for tree in subtrees.iter() {
-				tree.find_internal(bounding_box, set);
-			}
-		}
-	}
-
-	fn ensure_subtrees(&mut self) {
-		if self.subtrees.is_some() {
-			return;
-		}
-
-		let mid_x = (self.bounding_box.left + self.bounding_box.right) / 2.;
-		let mid_y = (self.bounding_box.top + self.bounding_box.bottom) / 2.;
-
-		self.subtrees = Some(Box::new([
-			QuadTree::new(
-				Aabb {
-					top: self.bounding_box.top,
-					right: mid_x,
-					bottom: mid_y,
-					left: self.bounding_box.left,
-				},
-				self.depth - 1,
-				self.inner_node_capacity,
-			),
-			QuadTree::new(
-				Aabb {
-					top: self.bounding_box.top,
-					right: self.bounding_box.right,
-					bottom: mid_y,
-					left: mid_x,
-				},
-				self.depth - 1,
-				self.inner_node_capacity,
-			),
-			QuadTree::new(
-				Aabb {
-					top: mid_y,
-					right: mid_x,
-					bottom: self.bounding_box.bottom,
-					left: self.bounding_box.left,
-				},
-				self.depth - 1,
-				self.inner_node_capacity,
-			),
-			QuadTree::new(
-				Aabb {
-					top: mid_y,
-					right: self.bounding_box.right,
-					bottom: self.bounding_box.bottom,
-					left: mid_x,
-				},
-				self.depth - 1,
-				self.inner_node_capacity,
-			),
-		]));
-	}
-}

node-graph/gpath-bool/src/lib.rs

@@ -114,11 +114,16 @@ fn subtract<'a>(vector: impl Iterator<Item = TableRowRef<'a, Vector>>) -> Table<
 
 	let mut result_vector_table = Table::new_from_row(vector.next().map(|x| x.into_cloned()).unwrap_or_default());
 	let mut first_row = result_vector_table.iter_mut().next().expect("Expected the one row we just pushed");
+	let first_row_transform = if first_row.transform.matrix2.determinant() != 0. {
+		first_row.transform.inverse()
+	} else {
+		DAffine2::IDENTITY
+	};
 
 	let mut next_vector = vector.next();
 
 	while let Some(lower_vector) = next_vector {
-		let transform_of_lower_into_space_of_upper = first_row.transform.inverse() * *lower_vector.transform;
+		let transform_of_lower_into_space_of_upper = first_row_transform * *lower_vector.transform;
 
 		let result = &mut first_row.element;