Special case cubic line intersections

Author: TrueDoctor

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

@@ -6,6 +6,7 @@ use crate::math::lerp;
 use crate::path_segment::PathSegment;
 
 use glam::DVec2;
+use roots::{find_roots_cubic, Roots};
 
 #[derive(Clone)]
 struct IntersectionSegment {
@@ -88,32 +89,41 @@ pub fn segments_equal(seg0: &PathSegment, seg1: &PathSegment, point_epsilon: f64
 }
 
 pub fn path_segment_intersection(seg0: &PathSegment, seg1: &PathSegment, endpoints: bool, eps: &Epsilons) -> Vec<[f64; 2]> {
-	if let (PathSegment::Line(start0, end0), PathSegment::Line(start1, end1)) = (seg0, seg1) {
-		if let Some(st) = line_segment_intersection([*start0, *end0], [*start1, *end1], eps.param) {
-			if !endpoints && (st.0 < eps.param || st.0 > 1. - eps.param) && (st.1 < eps.param || st.1 > 1. - eps.param) {
-				return vec![];
+	match (seg0, seg1) {
+		(PathSegment::Line(start0, end0), PathSegment::Line(start1, end1)) => {
+			if let Some(st) = line_segment_intersection([*start0, *end0], [*start1, *end1], eps.param) {
+				if !endpoints && (st.0 < eps.param || st.0 > 1. - eps.param) && (st.1 < eps.param || st.1 > 1. - eps.param) {
+					return vec![];
+				}
+				return vec![st.into()];
 			}
-			return vec![st.into()];
 		}
-	}
-	if let (PathSegment::Cubic(s1, c11, c21, e1), PathSegment::Cubic(s2, c12, c22, e2)) = (seg0, seg1) {
-		let path1 = lyon_geom::CubicBezierSegment {
-			from: (s1.x, s1.y).into(),
-			ctrl1: (c11.x, c11.y).into(),
-			ctrl2: (c21.x, c21.y).into(),
-			to: (e1.x, e1.y).into(),
-		};
-		let path2 = lyon_geom::CubicBezierSegment {
-			from: (s2.x, s2.y).into(),
-			ctrl1: (c12.x, c12.y).into(),
-			ctrl2: (c22.x, c22.y).into(),
-			to: (e2.x, e2.y).into(),
-		};
-
-		let intersections = path1.cubic_intersections_t(&path2);
-		let intersections: Vec<_> = intersections.into_iter().map(|(s, t)| [s, t]).collect();
-		return intersections;
-	}
+		(PathSegment::Cubic(s1, c11, c21, e1), PathSegment::Cubic(s2, c12, c22, e2)) => {
+			let path1 = lyon_geom::CubicBezierSegment {
+				from: (s1.x, s1.y).into(),
+				ctrl1: (c11.x, c11.y).into(),
+				ctrl2: (c21.x, c21.y).into(),
+				to: (e1.x, e1.y).into(),
+			};
+			let path2 = lyon_geom::CubicBezierSegment {
+				from: (s2.x, s2.y).into(),
+				ctrl1: (c12.x, c12.y).into(),
+				ctrl2: (c22.x, c22.y).into(),
+				to: (e2.x, e2.y).into(),
+			};
+
+			let intersections = path1.cubic_intersections_t(&path2);
+			let intersections: Vec<_> = intersections.into_iter().map(|(s, t)| [s, t]).collect();
+			return intersections;
+		}
+		(PathSegment::Cubic(_, _, _, _), PathSegment::Line(_, _)) => {
+			return cubic_line_intersection(seg0, seg1, eps);
+		}
+		(PathSegment::Line(_, _), PathSegment::Cubic(_, _, _, _)) => {
+			return cubic_line_intersection(seg1, seg0, eps).into_iter().map(|[t1, t0]| [t0, t1]).collect();
+		}
+		_ => (),
+	};
 
 	// https://math.stackexchange.com/questions/20321/how-can-i-tell-when-two-cubic-b%C3%A9zier-curves-intersect
 
@@ -237,6 +247,79 @@ fn calculate_overlap_intersections(seg0: &PathSegment, seg1: &PathSegment, eps:
 	intersections
 }
 
+fn cubic_line_intersection(cubic: &PathSegment, line: &PathSegment, eps: &Epsilons) -> Vec<[f64; 2]> {
+	let (p0, p1, p2, p3) = match cubic {
+		PathSegment::Cubic(p0, p1, p2, p3) => (p0, p1, p2, p3),
+		_ => unreachable!("Expected a cubic Bézier curve"),
+	};
+	let (line_start, line_end) = match line {
+		PathSegment::Line(start, end) => (start, end),
+		_ => unreachable!("Expected a line segment"),
+	};
+
+	// Transform the cubic curve to the line's coordinate system
+	let line_dir = *line_end - *line_start;
+	let line_len = line_dir.length();
+	let line_norm = line_dir / line_len;
+	let perp = DVec2::new(-line_norm.y, line_norm.x);
+
+	let transform = |p: &DVec2| -> DVec2 {
+		let v = *p - *line_start;
+		DVec2::new(v.dot(line_norm), v.dot(perp))
+	};
+
+	let q0 = transform(p0);
+	let q1 = transform(p1);
+	let q2 = transform(p2);
+	let q3 = transform(p3);
+
+	// Now we're looking for roots of the cubic equation in y
+	let a = -q0.y + 3.0 * q1.y - 3.0 * q2.y + q3.y;
+	let b = 3.0 * q0.y - 6.0 * q1.y + 3.0 * q2.y;
+	let c = -3.0 * q0.y + 3.0 * q1.y;
+	let d = q0.y;
+
+	let roots = find_roots_cubic(a, b, c, d);
+
+	let mut intersections = Vec::new();
+
+	let mut process_root = |t: f64| {
+		if (0.0..=1.0).contains(&t) {
+			let x = ((1.0 - t).powi(3) * q0.x + 3.0 * t * (1.0 - t).powi(2) * q1.x + 3.0 * t.powi(2) * (1.0 - t) * q2.x + t.powi(3) * q3.x) / line_len;
+			if (0.0..=1.0).contains(&x) {
+				intersections.push([t, x]);
+			}
+		}
+	};
+
+	match roots {
+		Roots::Three(roots) => {
+			for &t in roots.iter() {
+				process_root(t);
+			}
+		}
+		Roots::Two(roots) => {
+			for &t in roots.iter() {
+				process_root(t);
+			}
+		}
+		Roots::One(roots) => {
+			for &t in roots.iter() {
+				process_root(t);
+			}
+		}
+		_ => {}
+	}
+
+	// Sort intersections by the cubic's t parameter
+	intersections.sort_by(|a, b| a[0].partial_cmp(&b[0]).unwrap());
+
+	// Remove duplicates within epsilon
+	intersections.dedup_by(|a, b| (a[0] - b[0]).abs() < eps.param);
+
+	intersections
+}
+
 fn find_point_on_segment(seg: &PathSegment, point: DVec2, eps: &Epsilons) -> Option<f64> {
 	let start = 0.;
 	let end = 1.;