Bezier-rs: Add parametric evaluate and line intersect to subpath (#852)

* add slider to subpath component + change evaluate to take an enum

Co-authored-by: Rob Nadal <RobNadal@users.noreply.github.com>

wip - add intersect to subpath, TODO fix bug

Co-authored-by: Rob Nadal <RobNadal@users.noreply.github.com>

add unit tests to subpath intersections

stress, testing

Co-authored-by: Hannah Li <hannahli2010@gmail.com>

* add parametric eval impl to subpath

* add line intersection to subpath

* Uncomment and #[ignore] disabled tests

* Reorder a few imports

* change subpath:eval slider to radio button

* fixed bug with solve_cubic, fixed unit tests, improved intersection accuracy

* fix failing test

Co-authored-by: Hannah Li <hannahli2010@gmail.com>
Co-authored-by: Keavon Chambers <keavon@keavon.com>

libraries/bezier-rs/src/bezier/solvers.rs

@@ -222,11 +222,36 @@ impl Bezier {
 		}
 	}
 
+	// TODO: Use an `impl Iterator` return type instead of a `Vec`
+	/// Returns a list of filtered `t` values that correspond to intersection points between the current bezier curve and the provided one
+	/// such that the difference between adjacent `t` values in sorted order is greater than some minimum seperation value. If the difference
+	/// between 2 adjacent `t` values is lesss than the minimum difference, the filtering takes the larger `t` value and discards the smaller `t` value.
+	/// The returned `t` values are with respect to the current bezier, not the provided parameter.
+	/// If the provided curve is linear, then zero intersection points will be returned along colinear segments.
+	/// - `error` - For intersections where the provided bezier is non-linear, `error` defines the threshold for bounding boxes to be considered an intersection point.
+	/// - `minimum_seperation` - The minimum difference between adjacent `t` values in sorted order
+	pub fn intersections(&self, other: &Bezier, error: Option<f64>, minimum_seperation: Option<f64>) -> Vec<f64> {
+		// TODO: Consider using the `intersections_between_vectors_of_curves` helper function here
+		// Otherwise, use bounding box to determine intersections
+		let mut intersection_t_values = self.unfiltered_intersections(other, error);
+		intersection_t_values.sort_by(|a, b| a.partial_cmp(b).unwrap());
+
+		// println!("<<<<< intersection_t_values :: {:?}", intersection_t_values);
+
+		intersection_t_values.iter().fold(Vec::new(), |mut accumulator, t| {
+			if !accumulator.is_empty() && (accumulator.last().unwrap() - t).abs() < minimum_seperation.unwrap_or(MIN_SEPERATION_VALUE) {
+				accumulator.pop();
+			}
+			accumulator.push(*t);
+			accumulator
+		})
+	}
+
 	// TODO: Use an `impl Iterator` return type instead of a `Vec`
 	/// Returns a list of `t` values that correspond to intersection points between the current bezier curve and the provided one. The returned `t` values are with respect to the current bezier, not the provided parameter.
 	/// If the provided curve is linear, then zero intersection points will be returned along colinear segments.
 	/// - `error` - For intersections where the provided bezier is non-linear, `error` defines the threshold for bounding boxes to be considered an intersection point.
-	pub fn intersections(&self, other: &Bezier, error: Option<f64>) -> Vec<f64> {
+	fn unfiltered_intersections(&self, other: &Bezier, error: Option<f64>) -> Vec<f64> {
 		let error = error.unwrap_or(0.5);
 		if other.handles == BezierHandles::Linear {
 			// Rotate the bezier and the line by the angle that the line makes with the x axis
@@ -295,7 +320,7 @@ impl Bezier {
 		let segment_pairs = subcurves1.iter().flat_map(move |(curve1, curve1_t_pair)| {
 			subcurves2
 				.iter()
-				.filter_map(move |(curve2, curve2_t_pair)| utils::do_rectangles_overlap(curve1.bounding_box(), curve2.bounding_box()).then(|| (curve1, curve1_t_pair, curve2, curve2_t_pair)))
+				.filter_map(move |(curve2, curve2_t_pair)| utils::do_rectangles_overlap(curve1.bounding_box(), curve2.bounding_box()).then_some((curve1, curve1_t_pair, curve2, curve2_t_pair)))
 		});
 		segment_pairs
 			.flat_map(|(curve1, curve1_t_pair, curve2, curve2_t_pair)| curve1.intersections_between_subcurves(curve1_t_pair.clone(), curve2, curve2_t_pair.clone(), error))
@@ -563,13 +588,13 @@ mod tests {
 		// Intersection at edge of curve
 		let bezier = Bezier::from_linear_dvec2(p1, p2);
 		let line1 = Bezier::from_linear_coordinates(20., 60., 70., 60.);
-		let intersections1 = bezier.intersections(&line1, None);
+		let intersections1 = bezier.intersections(&line1, None, None);
 		assert!(intersections1.len() == 1);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections1[0])), DVec2::new(30., 60.)));
 
 		// Intersection in the middle of curve
 		let line2 = Bezier::from_linear_coordinates(150., 150., 30., 30.);
-		let intersections2 = bezier.intersections(&line2, None);
+		let intersections2 = bezier.intersections(&line2, None, None);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections2[0])), DVec2::new(96., 96.)));
 	}
 
@@ -582,13 +607,13 @@ mod tests {
 		// Intersection at edge of curve
 		let bezier = Bezier::from_quadratic_dvec2(p1, p2, p3);
 		let line1 = Bezier::from_linear_coordinates(20., 50., 40., 50.);
-		let intersections1 = bezier.intersections(&line1, None);
+		let intersections1 = bezier.intersections(&line1, None, None);
 		assert!(intersections1.len() == 1);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections1[0])), p1));
 
 		// Intersection in the middle of curve
 		let line2 = Bezier::from_linear_coordinates(150., 150., 30., 30.);
-		let intersections2 = bezier.intersections(&line2, None);
+		let intersections2 = bezier.intersections(&line2, None, None);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections2[0])), DVec2::new(47.77355, 47.77354)));
 	}
 
@@ -602,30 +627,63 @@ mod tests {
 		let bezier = Bezier::from_cubic_dvec2(p1, p2, p3, p4);
 		// Intersection at edge of curve, Discriminant > 0
 		let line1 = Bezier::from_linear_coordinates(20., 30., 40., 30.);
-		let intersections1 = bezier.intersections(&line1, None);
+		let intersections1 = bezier.intersections(&line1, None, None);
 		assert!(intersections1.len() == 1);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections1[0])), p1));
 
 		// Intersection at edge and in middle of curve, Discriminant < 0
 		let line2 = Bezier::from_linear_coordinates(150., 150., 30., 30.);
-		let intersections2 = bezier.intersections(&line2, None);
+		let intersections2 = bezier.intersections(&line2, None, None);
 		assert!(intersections2.len() == 2);
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections2[0])), p1));
 		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections2[1])), DVec2::new(85.84, 85.84)));
 	}
 
+	#[test]
+	fn test_intersect_curve_cubic_anchor_handle_overlap() {
+		// M31 94 C40 40 107 107 106 106
+
+		let p1 = DVec2::new(31., 94.);
+		let p2 = DVec2::new(40., 40.);
+		let p3 = DVec2::new(107., 107.);
+		let p4 = DVec2::new(106., 106.);
+		let bezier = Bezier::from_cubic_dvec2(p1, p2, p3, p4);
+
+		let line = Bezier::from_linear_coordinates(150., 150., 20., 20.);
+		let intersections = bezier.intersections(&line, None, None);
+
+		assert_eq!(intersections.len(), 1);
+		assert!(compare_points(bezier.evaluate(ComputeType::Parametric(intersections[0])), p4));
+	}
+
+	#[test]
+	fn test_intersect_curve_cubic_edge_case() {
+		// M34 107 C40 40 120 120 102 29
+
+		let p1 = DVec2::new(34., 107.);
+		let p2 = DVec2::new(40., 40.);
+		let p3 = DVec2::new(120., 120.);
+		let p4 = DVec2::new(102., 29.);
+		let bezier = Bezier::from_cubic_dvec2(p1, p2, p3, p4);
+
+		let line = Bezier::from_linear_coordinates(150., 150., 20., 20.);
+		let intersections = bezier.intersections(&line, None, None);
+
+		assert_eq!(intersections.len(), 1);
+	}
+
 	#[test]
 	fn test_intersect_curve() {
 		let bezier1 = Bezier::from_cubic_coordinates(30., 30., 60., 140., 150., 30., 160., 160.);
 		let bezier2 = Bezier::from_quadratic_coordinates(175., 140., 20., 20., 120., 20.);
 
-		let intersections = bezier1.intersections(&bezier2, None);
-		let intersections2 = bezier2.intersections(&bezier1, None);
-		assert!(compare_vec_of_points(
-			intersections.iter().map(|&t| bezier1.evaluate(ComputeType::Parametric(t))).collect(),
-			intersections2.iter().map(|&t| bezier2.evaluate(ComputeType::Parametric(t))).collect(),
-			2.
-		));
+		let intersections1 = bezier1.intersections(&bezier2, None, None);
+		let intersections2 = bezier2.intersections(&bezier1, None, None);
+
+		let intersections1_points: Vec<DVec2> = intersections1.iter().map(|&t| bezier1.evaluate(ComputeType::Parametric(t))).collect();
+		let intersections2_points: Vec<DVec2> = intersections2.iter().map(|&t| bezier2.evaluate(ComputeType::Parametric(t))).rev().collect();
+
+		assert!(compare_vec_of_points(intersections1_points, intersections2_points, 2.));
 	}
 
 	#[test]

libraries/bezier-rs/src/consts.rs

@@ -8,6 +8,8 @@ pub const STRICT_MAX_ABSOLUTE_DIFFERENCE: f64 = 1e-6;
 pub const NUM_DISTANCES: usize = 5;
 /// Maximum allowed angle that the normal of the `start` or `end` point can make with the normal of the corresponding handle for a curve to be considered scalable/simple.
 pub const SCALABLE_CURVE_MAX_ENDPOINT_NORMAL_ANGLE: f64 = std::f64::consts::PI / 3.;
+/// Minimum allowable separation between adjacent `t` values when calculating curve intersections
+pub const MIN_SEPERATION_VALUE: f64 = 5. * 1e-3;
 
 // Method argument defaults
 

libraries/bezier-rs/src/subpath/core.rs

@@ -1,5 +1,6 @@
 use super::*;
 use crate::consts::*;
+
 use std::fmt::Write;
 
 /// Functionality relating to core `Subpath` operations, such as constructors and `iter`.
@@ -40,6 +41,15 @@ impl Subpath {
 		self.manipulator_groups.len()
 	}
 
+	/// Returns the number of segments contained within the `Subpath`.
+	pub fn len_segments(&self) -> usize {
+		let mut number_of_curves = self.len();
+		if !self.closed {
+			number_of_curves -= 1
+		}
+		number_of_curves
+	}
+
 	/// Returns an iterator of the [Bezier]s along the `Subpath`.
 	pub fn iter(&self) -> SubpathIter {
 		SubpathIter { sub_path: self, index: 0 }

libraries/bezier-rs/src/subpath/mod.rs

@@ -1,5 +1,6 @@
 mod core;
 mod lookup;
+mod solvers;
 mod structs;
 pub use structs::*;
 

libraries/bezier-rs/src/subpath/solvers.rs

@@ -0,0 +1,405 @@
+use super::*;
+use crate::{consts::MIN_SEPERATION_VALUE, ComputeType};
+
+use glam::DVec2;
+
+impl Subpath {
+	/// Calculate the point on the subpath based on the parametric `t`-value provided.
+	/// Expects `t` to be within the inclusive range `[0, 1]`.
+	pub fn evaluate(&self, t: ComputeType) -> DVec2 {
+		match t {
+			ComputeType::Parametric(t) => {
+				assert!((0.0..=1.).contains(&t));
+
+				let number_of_curves = self.len_segments() as f64;
+				let scaled_t = t * number_of_curves;
+
+				let target_curve_index = scaled_t.floor() as i32;
+				let target_curve_t = scaled_t % 1.;
+
+				if let Some(curve) = self.iter().nth(target_curve_index as usize) {
+					curve.evaluate(ComputeType::Parametric(target_curve_t))
+				} else {
+					self.iter().last().unwrap().evaluate(ComputeType::Parametric(1.))
+				}
+			}
+			// TODO: change this implementation to Euclidean compute
+			ComputeType::Euclidean(_t) => self.iter().next().unwrap().evaluate(ComputeType::Parametric(0.)),
+			ComputeType::EuclideanWithinError { t: _, epsilon: _ } => todo!(),
+		}
+	}
+
+	/// Calculates the intersection points the subpath has with a given line and returns a list of parameteric `t`-values.
+	/// This function expects the following:
+	/// - other: a [Bezier] curve to check intersections against
+	/// - error: an optional f64 value to provide an error bound
+	pub fn intersections(&self, other: &Bezier, error: Option<f64>, minimum_seperation: Option<f64>) -> Vec<f64> {
+		// TODO: account for either euclidean or parametric type
+		let number_of_curves = self.len_segments() as f64;
+		let intersection_t_values: Vec<f64> = self
+			.iter()
+			.enumerate()
+			.flat_map(|(index, bezier)| {
+				bezier
+					.intersections(other, error, minimum_seperation)
+					.into_iter()
+					.map(|t| ((index as f64) + t) / number_of_curves)
+					.collect::<Vec<f64>>()
+			})
+			.collect();
+
+		intersection_t_values.iter().fold(Vec::new(), |mut accumulator, t| {
+			if !accumulator.is_empty() && (accumulator.last().unwrap() - t).abs() < minimum_seperation.unwrap_or(MIN_SEPERATION_VALUE) {
+				accumulator.pop();
+			}
+			accumulator.push(*t);
+			accumulator
+		});
+
+		intersection_t_values
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use crate::Bezier;
+	use glam::DVec2;
+
+	use crate::consts::MAX_ABSOLUTE_DIFFERENCE;
+	use crate::utils;
+
+	fn normalize_t(n: i64, t: f64) -> f64 {
+		t * (n as f64) % 1.
+	}
+
+	#[test]
+	fn evaluate_one_subpath_curve() {
+		let start = DVec2::new(20., 30.);
+		let end = DVec2::new(60., 45.);
+		let handle = DVec2::new(75., 85.);
+
+		let bezier = Bezier::from_quadratic_dvec2(start, handle, end);
+		let subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: start,
+					in_handle: None,
+					out_handle: Some(handle),
+				},
+				ManipulatorGroup {
+					anchor: end,
+					in_handle: None,
+					out_handle: Some(handle),
+				},
+			],
+			false,
+		);
+
+		let t0 = 0.;
+		assert_eq!(subpath.evaluate(ComputeType::Parametric(t0)), bezier.evaluate(ComputeType::Parametric(t0)));
+
+		let t1 = 0.25;
+		assert_eq!(subpath.evaluate(ComputeType::Parametric(t1)), bezier.evaluate(ComputeType::Parametric(t1)));
+
+		let t2 = 0.50;
+		assert_eq!(subpath.evaluate(ComputeType::Parametric(t2)), bezier.evaluate(ComputeType::Parametric(t2)));
+
+		let t3 = 1.;
+		assert_eq!(subpath.evaluate(ComputeType::Parametric(t3)), bezier.evaluate(ComputeType::Parametric(t3)));
+	}
+
+	#[test]
+	fn evaluate_multiple_subpath_curves() {
+		let start = DVec2::new(20., 30.);
+		let middle = DVec2::new(70., 70.);
+		let end = DVec2::new(60., 45.);
+		let handle1 = DVec2::new(75., 85.);
+		let handle2 = DVec2::new(40., 30.);
+		let handle3 = DVec2::new(10., 10.);
+
+		let linear_bezier = Bezier::from_linear_dvec2(start, middle);
+		let quadratic_bezier = Bezier::from_quadratic_dvec2(middle, handle1, end);
+		let cubic_bezier = Bezier::from_cubic_dvec2(end, handle2, handle3, start);
+
+		let mut subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: start,
+					in_handle: Some(handle3),
+					out_handle: None,
+				},
+				ManipulatorGroup {
+					anchor: middle,
+					in_handle: None,
+					out_handle: Some(handle1),
+				},
+				ManipulatorGroup {
+					anchor: end,
+					in_handle: None,
+					out_handle: Some(handle2),
+				},
+			],
+			false,
+		);
+
+		// Test open subpath
+
+		let mut n = (subpath.len() as i64) - 1;
+
+		let t0 = 0.;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t0)),
+			linear_bezier.evaluate(ComputeType::Parametric(normalize_t(n, t0))),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		let t1 = 0.25;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t1)),
+			linear_bezier.evaluate(ComputeType::Parametric(normalize_t(n, t1))),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		let t2 = 0.50;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t2)),
+			quadratic_bezier.evaluate(ComputeType::Parametric(normalize_t(n, t2))),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		let t3 = 0.75;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t3)),
+			quadratic_bezier.evaluate(ComputeType::Parametric(normalize_t(n, t3))),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		let t4 = 1.0;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t4)),
+			quadratic_bezier.evaluate(ComputeType::Parametric(1.)),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		// Test closed subpath
+
+		subpath.closed = true;
+		n = subpath.len() as i64;
+
+		let t5 = 2. / 3.;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t5)),
+			cubic_bezier.evaluate(ComputeType::Parametric(normalize_t(n, t5))),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		let t6 = 1.;
+		assert!(utils::dvec2_compare(
+			subpath.evaluate(ComputeType::Parametric(t6)),
+			cubic_bezier.evaluate(ComputeType::Parametric(1.)),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+	}
+
+	#[test]
+	fn intersection_linear_multiple_subpath_curves_test_one() {
+		// M 35 125 C 40 40 120 120 43 43 Q 175 90 145 150 Q 70 185 35 125 Z
+
+		let cubic_start = DVec2::new(35., 125.);
+		let cubic_handle_1 = DVec2::new(40., 40.);
+		let cubic_handle_2 = DVec2::new(120., 120.);
+		let cubic_end = DVec2::new(43., 43.);
+
+		let quadratic_1_handle = DVec2::new(175., 90.);
+		let quadratic_end = DVec2::new(145., 150.);
+
+		let quadratic_2_handle = DVec2::new(70., 185.);
+
+		let cubic_bezier = Bezier::from_cubic_dvec2(cubic_start, cubic_handle_1, cubic_handle_2, cubic_end);
+		let quadratic_bezier_1 = Bezier::from_quadratic_dvec2(cubic_end, quadratic_1_handle, quadratic_end);
+
+		let subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: cubic_start,
+					in_handle: None,
+					out_handle: Some(cubic_handle_1),
+				},
+				ManipulatorGroup {
+					anchor: cubic_end,
+					in_handle: Some(cubic_handle_2),
+					out_handle: None,
+				},
+				ManipulatorGroup {
+					anchor: quadratic_end,
+					in_handle: Some(quadratic_1_handle),
+					out_handle: Some(quadratic_2_handle),
+				},
+			],
+			true,
+		);
+
+		let line = Bezier::from_linear_coordinates(150., 150., 20., 20.);
+
+		let cubic_intersections = cubic_bezier.intersections(&line, None, None);
+		let quadratic_1_intersections = quadratic_bezier_1.intersections(&line, None, None);
+		let subpath_intersections = subpath.intersections(&line, None, None);
+
+		assert!(utils::dvec2_compare(
+			cubic_bezier.evaluate(ComputeType::Parametric(cubic_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[0])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		assert!(utils::dvec2_compare(
+			quadratic_bezier_1.evaluate(ComputeType::Parametric(quadratic_1_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[1])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		assert!(utils::dvec2_compare(
+			quadratic_bezier_1.evaluate(ComputeType::Parametric(quadratic_1_intersections[1])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[2])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+	}
+
+	#[test]
+	fn intersection_linear_multiple_subpath_curves_test_two() {
+		// M34 107 C40 40 120 120 102 29 Q175 90 129 171 Q70 185 34 107 Z
+		// M150 150 L 20 20
+
+		let cubic_start = DVec2::new(34., 107.);
+		let cubic_handle_1 = DVec2::new(40., 40.);
+		let cubic_handle_2 = DVec2::new(120., 120.);
+		let cubic_end = DVec2::new(102., 29.);
+
+		let quadratic_1_handle = DVec2::new(175., 90.);
+		let quadratic_end = DVec2::new(129., 171.);
+
+		let quadratic_2_handle = DVec2::new(70., 185.);
+
+		let cubic_bezier = Bezier::from_cubic_dvec2(cubic_start, cubic_handle_1, cubic_handle_2, cubic_end);
+		let quadratic_bezier_1 = Bezier::from_quadratic_dvec2(cubic_end, quadratic_1_handle, quadratic_end);
+
+		let subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: cubic_start,
+					in_handle: None,
+					out_handle: Some(cubic_handle_1),
+				},
+				ManipulatorGroup {
+					anchor: cubic_end,
+					in_handle: Some(cubic_handle_2),
+					out_handle: None,
+				},
+				ManipulatorGroup {
+					anchor: quadratic_end,
+					in_handle: Some(quadratic_1_handle),
+					out_handle: Some(quadratic_2_handle),
+				},
+			],
+			true,
+		);
+
+		let line = Bezier::from_linear_coordinates(150., 150., 20., 20.);
+
+		let cubic_intersections = cubic_bezier.intersections(&line, None, None);
+		let quadratic_1_intersections = quadratic_bezier_1.intersections(&line, None, None);
+		let subpath_intersections = subpath.intersections(&line, None, None);
+
+		assert!(utils::dvec2_compare(
+			cubic_bezier.evaluate(ComputeType::Parametric(cubic_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[0])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		assert!(utils::dvec2_compare(
+			quadratic_bezier_1.evaluate(ComputeType::Parametric(quadratic_1_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[1])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+	}
+
+	#[test]
+	fn intersection_linear_multiple_subpath_curves_test_three() {
+		// M35 125 C40 40 120 120 44 44 Q175 90 145 150 Q70 185 35 125 Z
+
+		let cubic_start = DVec2::new(35., 125.);
+		let cubic_handle_1 = DVec2::new(40., 40.);
+		let cubic_handle_2 = DVec2::new(120., 120.);
+		let cubic_end = DVec2::new(44., 44.);
+
+		let quadratic_1_handle = DVec2::new(175., 90.);
+		let quadratic_end = DVec2::new(145., 150.);
+
+		let quadratic_2_handle = DVec2::new(70., 185.);
+
+		let cubic_bezier = Bezier::from_cubic_dvec2(cubic_start, cubic_handle_1, cubic_handle_2, cubic_end);
+		let quadratic_bezier_1 = Bezier::from_quadratic_dvec2(cubic_end, quadratic_1_handle, quadratic_end);
+
+		let subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: cubic_start,
+					in_handle: None,
+					out_handle: Some(cubic_handle_1),
+				},
+				ManipulatorGroup {
+					anchor: cubic_end,
+					in_handle: Some(cubic_handle_2),
+					out_handle: None,
+				},
+				ManipulatorGroup {
+					anchor: quadratic_end,
+					in_handle: Some(quadratic_1_handle),
+					out_handle: Some(quadratic_2_handle),
+				},
+			],
+			true,
+		);
+
+		let line = Bezier::from_linear_coordinates(150., 150., 20., 20.);
+
+		let cubic_intersections = cubic_bezier.intersections(&line, None, None);
+		let quadratic_1_intersections = quadratic_bezier_1.intersections(&line, None, None);
+		let subpath_intersections = subpath.intersections(&line, None, None);
+
+		assert!(utils::dvec2_compare(
+			cubic_bezier.evaluate(ComputeType::Parametric(cubic_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[0])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		assert!(utils::dvec2_compare(
+			quadratic_bezier_1.evaluate(ComputeType::Parametric(quadratic_1_intersections[0])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[1])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+
+		assert!(utils::dvec2_compare(
+			quadratic_bezier_1.evaluate(ComputeType::Parametric(quadratic_1_intersections[1])),
+			subpath.evaluate(ComputeType::Parametric(subpath_intersections[2])),
+			MAX_ABSOLUTE_DIFFERENCE
+		)
+		.all());
+	}
+
+	// TODO: add more intersection tests
+}

libraries/bezier-rs/src/utils.rs

@@ -1,4 +1,4 @@
-use crate::consts::{MAX_ABSOLUTE_DIFFERENCE, STRICT_MAX_ABSOLUTE_DIFFERENCE};
+use crate::consts::{MAX_ABSOLUTE_DIFFERENCE, MIN_SEPERATION_VALUE, STRICT_MAX_ABSOLUTE_DIFFERENCE};
 
 use glam::{BVec2, DMat2, DVec2};
 use std::f64::consts::PI;
@@ -90,21 +90,17 @@ fn cube_root(f: f64) -> f64 {
 /// Solve a cubic of the form `x^3 + px + q`, derivation from: <https://trans4mind.com/personal_development/mathematics/polynomials/cubicAlgebra.htm>.
 pub fn solve_reformatted_cubic(discriminant: f64, a: f64, p: f64, q: f64) -> Vec<f64> {
 	let mut roots = Vec::new();
-	if p.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
-		// Handle when p is approximately 0
-		roots.push(cube_root(-q));
-	} else if q.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
-		// Handle when q is approximately 0
-		if p < 0. {
-			roots.push((-p).powf(1. / 2.));
-		}
-	} else if discriminant.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
+	if discriminant.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
 		// When discriminant is 0 (check for approximation because of floating point errors), all roots are real, and 2 are repeated
+		// filter out repeated roots (ie. roots whose distance is less than some epsilon)
 		let q_divided_by_2 = q / 2.;
 		let a_divided_by_3 = a / 3.;
-
-		roots.push(2. * cube_root(-q_divided_by_2) - a_divided_by_3);
-		roots.push(cube_root(q_divided_by_2) - a_divided_by_3);
+		let root_1 = 2. * cube_root(-q_divided_by_2) - a_divided_by_3;
+		let root_2 = cube_root(q_divided_by_2) - a_divided_by_3;
+		if (root_1 - root_2).abs() > MIN_SEPERATION_VALUE {
+			roots.push(root_1);
+		}
+		roots.push(root_2);
 	} else if discriminant > 0. {
 		// When discriminant > 0, there is one real and two imaginary roots
 		let q_divided_by_2 = q / 2.;
@@ -139,6 +135,7 @@ pub fn solve_cubic(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
 			solve_quadratic(discriminant, 2. * b, c, d)
 		}
 	} else {
+		// convert at^3 + bt^2 + ct + d ==> t^3 + a't^2 + b't + c'
 		let new_a = b / a;
 		let new_b = c / a;
 		let new_c = d / a;