From 3c2fff446545a7ad8e0da3cab879aa9ea86f1664 Mon Sep 17 00:00:00 2001 From: Hannah Li Date: Mon, 4 Jul 2022 19:29:25 -0400 Subject: [PATCH] Implement function to find intersections between a Bezier and a linear line (#708) * Implement line intersection for quadratics, begin work for cubic * Implement line intersection for cubic beziers, add tests for cubic root finding * Rename function and update comments * Minor refactor and adjust comments * Address PR comments --- bezier-rs/docs/interactive-docs/src/App.vue | 39 +++++ .../docs/interactive-docs/wasm/src/lib.rs | 9 ++ bezier-rs/lib/src/lib.rs | 135 +++++++++++++++- bezier-rs/lib/src/utils.rs | 151 +++++++++++++++++- 4 files changed, 322 insertions(+), 12 deletions(-) diff --git a/bezier-rs/docs/interactive-docs/src/App.vue b/bezier-rs/docs/interactive-docs/src/App.vue index 4a9de38b..ac010d3c 100644 --- a/bezier-rs/docs/interactive-docs/src/App.vue +++ b/bezier-rs/docs/interactive-docs/src/App.vue @@ -250,6 +250,45 @@ export default defineComponent({ }); }, }, + { + name: "Rotate", + callback: (canvas: HTMLCanvasElement, bezier: WasmBezierInstance, options: Record): void => { + const context = getContextFromCanvas(canvas); + const rotatedBezier = bezier + .rotate((options.angle * Math.PI) / 180) + .get_points() + .map((p) => JSON.parse(p)); + drawBezier(context, rotatedBezier, null, { curveStrokeColor: COLORS.NON_INTERACTIVE.STROKE_1, radius: 3.5 }); + }, + template: markRaw(SliderExample), + templateOptions: { + sliders: [ + { + variable: "angle", + min: -90, + max: 90, + step: 5, + default: 15, + }, + ], + }, + }, + { + name: "Line Intersection", + callback: (canvas: HTMLCanvasElement, bezier: WasmBezierInstance): void => { + const context = getContextFromCanvas(canvas); + const line = [ + { x: 150, y: 150 }, + { x: 30, y: 30 }, + ]; + const mappedLine = line.map((p) => [p.x, p.y]); + drawLine(context, line[0], line[1], COLORS.NON_INTERACTIVE.STROKE_1); + const intersections: Point[] = bezier.line_intersection(mappedLine).map((p) => JSON.parse(p)); + intersections.forEach((p: Point) => { + drawPoint(context, p, 3, COLORS.NON_INTERACTIVE.STROKE_2); + }); + }, + }, ], }; }, diff --git a/bezier-rs/docs/interactive-docs/wasm/src/lib.rs b/bezier-rs/docs/interactive-docs/wasm/src/lib.rs index 54961d6a..f9a46b54 100644 --- a/bezier-rs/docs/interactive-docs/wasm/src/lib.rs +++ b/bezier-rs/docs/interactive-docs/wasm/src/lib.rs @@ -107,4 +107,13 @@ impl WasmBezier { let local_extrema = self.0.local_extrema(); JsValue::from_serde(&serde_json::to_string(&local_extrema).unwrap()).unwrap() } + + pub fn rotate(&self, angle: f64) -> WasmBezier { + WasmBezier(self.0.rotate(angle)) + } + + pub fn line_intersection(&self, js_points: &JsValue) -> Vec { + let line: [DVec2; 2] = js_points.into_serde().unwrap(); + self.0.line_intersection(line).iter().map(|&p| vec_to_point(&p)).collect::>() + } } diff --git a/bezier-rs/lib/src/lib.rs b/bezier-rs/lib/src/lib.rs index dbc6e0f9..b0d8103a 100644 --- a/bezier-rs/lib/src/lib.rs +++ b/bezier-rs/lib/src/lib.rs @@ -1,4 +1,6 @@ -use glam::DVec2; +//! Bezier-rs: A Bezier Math Library for Rust + +use glam::{DMat2, DVec2}; mod utils; @@ -201,9 +203,7 @@ impl Bezier { /// Calculate the point on the curve based on the `t`-value provided. /// Basis code based off of pseudocode found here: . - pub fn compute(&self, t: f64) -> DVec2 { - assert!((0.0..=1.0).contains(&t)); - + fn unrestricted_compute(&self, t: f64) -> DVec2 { let t_squared = t * t; let one_minus_t = 1.0 - t; let squared_one_minus_t = one_minus_t * one_minus_t; @@ -218,6 +218,13 @@ impl Bezier { } } + /// Calculate the point on the curve based on the `t`-value provided. + /// Expects `t` to be within the inclusive range `[0, 1]`. + pub fn compute(&self, t: f64) -> DVec2 { + assert!((0.0..=1.0).contains(&t)); + self.unrestricted_compute(t) + } + /// Return a selection of equidistant points on the bezier curve. /// If no value is provided for `steps`, then the function will default `steps` to be 10. pub fn compute_lookup_table(&self, steps: Option) -> Vec { @@ -453,15 +460,92 @@ impl Bezier { .try_into() .unwrap() } + + /// Returns a Bezier curve that results from applying the tranformation function to each point in the Bezier. + pub fn apply_transformation(&self, transformation_function: &dyn Fn(DVec2) -> DVec2) -> Bezier { + let transformed_start = transformation_function(self.start); + let transformed_end = transformation_function(self.end); + match self.handles { + BezierHandles::Quadratic { handle } => { + let transformed_handle = transformation_function(handle); + Bezier::from_quadratic_dvec2(transformed_start, transformed_handle, transformed_end) + } + BezierHandles::Cubic { handle_start, handle_end } => { + let transformed_handle_start = transformation_function(handle_start); + let transformed_handle_end = transformation_function(handle_end); + Bezier::from_cubic_dvec2(transformed_start, transformed_handle_start, transformed_handle_end, transformed_end) + } + } + } + + /// Returns a Bezier curve that results from rotating the curve around the origin by the given angle (in radians). + pub fn rotate(&self, angle: f64) -> Bezier { + let rotation_matrix = DMat2::from_angle(angle); + self.apply_transformation(&|point| rotation_matrix.mul_vec2(point)) + } + + /// Returns a Bezier curve that results from translating the curve by the given `DVec2`. + pub fn translate(&self, translation: DVec2) -> Bezier { + self.apply_transformation(&|point| point + translation) + } + + /// Returns a list of points where the provided line segment intersects with the Bezier curve. + /// - `line` - A line segment expected to be received in the format of `[start_point, end_point]`. + pub fn line_intersection(&self, line: [DVec2; 2]) -> Vec { + // Rotate the bezier and the line by the angle that the line makes with the x axis + let slope = line[1] - line[0]; + let angle = slope.angle_between(DVec2::new(1., 0.)); + let rotation_matrix = DMat2::from_angle(angle); + let rotated_bezier = self.apply_transformation(&|point| rotation_matrix.mul_vec2(point)); + let rotated_line = [rotation_matrix.mul_vec2(line[0]), rotation_matrix.mul_vec2(line[1])]; + + // Translate the bezier such that the line becomes aligned on top of the x-axis + let vertical_distance = rotated_line[0].y; + let translated_bezier = rotated_bezier.translate(DVec2::new(0., -vertical_distance)); + + // Compute the roots of the resulting bezier curve + let list_intersection_t = match translated_bezier.handles { + BezierHandles::Quadratic { handle } => { + let a = translated_bezier.start.y - 2. * handle.y + translated_bezier.end.y; + let b = 2. * (handle.y - translated_bezier.start.y); + let c = translated_bezier.start.y; + + let discriminant = b * b - 4. * a * c; + let two_times_a = 2. * a; + + utils::solve_quadratic(discriminant, two_times_a, b, c) + } + BezierHandles::Cubic { handle_start, handle_end } => { + let start_y = translated_bezier.start.y; + let a = -start_y + 3. * handle_start.y - 3. * handle_end.y + translated_bezier.end.y; + let b = 3. * start_y - 6. * handle_start.y + 3. * handle_end.y; + let c = -3. * start_y + 3. * handle_start.y; + let d = start_y; + + utils::solve_cubic(a, b, c, d) + } + }; + let min = line[0].min(line[1]); + let max = line[0].max(line[1]); + let max_abs_diff = 1e-4; + + list_intersection_t + .iter() + .filter(|&&t| utils::f64_approximately_in_range(t, 0., 1., max_abs_diff)) + .map(|&t| self.unrestricted_compute(t)) + .filter(|&point| utils::dvec2_approximately_in_range(point, min, max, max_abs_diff).all()) + .collect::>() + } } #[cfg(test)] mod tests { + use crate::utils; use crate::Bezier; use glam::DVec2; fn compare_points(p1: DVec2, p2: DVec2) -> bool { - p1.abs_diff_eq(p2, 0.001) + utils::dvec2_compare(p1, p2, 1e-3).all() } #[test] @@ -505,4 +589,45 @@ mod tests { let bezier2 = Bezier::from_quadratic_coordinates(0., 0., 0., 100., 100., 100.); assert!(bezier2.project(DVec2::new(100., 0.), 20, 0.0001, 3, 10) == DVec2::new(0., 0.)); } + + #[test] + fn line_intersection_quadratic() { + let p1 = DVec2::new(30., 50.); + let p2 = DVec2::new(140., 30.); + let p3 = DVec2::new(160., 170.); + + // Intersection at edge of curve + let bezier1 = Bezier::from_quadratic_dvec2(p1, p2, p3); + let line1 = [DVec2::new(20., 50.), DVec2::new(40., 50.)]; + let intersections1 = bezier1.line_intersection(line1); + assert!(intersections1.len() == 1); + assert!(compare_points(intersections1[0], p1)); + + // Intersection in the middle of curve + let line2 = [DVec2::new(150., 150.), DVec2::new(30., 30.)]; + let intersections2 = bezier1.line_intersection(line2); + assert!(compare_points(intersections2[0], DVec2::new(47.77355, 47.77354))); + } + + #[test] + fn line_intersection_cubic() { + let p1 = DVec2::new(30., 30.); + let p2 = DVec2::new(60., 140.); + let p3 = DVec2::new(150., 30.); + let p4 = DVec2::new(160., 160.); + + let bezier = Bezier::from_cubic_dvec2(p1, p2, p3, p4); + // Intersection at edge of curve, Discriminant > 0 + let line1 = [DVec2::new(20., 30.), DVec2::new(40., 30.)]; + let intersections1 = bezier.line_intersection(line1); + assert!(intersections1.len() == 1); + assert!(compare_points(intersections1[0], p1)); + + // Intersection at edge and in middle of curve, Discriminant < 0 + let line2 = [DVec2::new(150., 150.), DVec2::new(30., 30.)]; + let intersections2 = bezier.line_intersection(line2); + assert!(intersections2.len() == 2); + assert!(compare_points(intersections2[0], p1)); + assert!(compare_points(intersections2[1], DVec2::new(85.84, 85.84))); + } } diff --git a/bezier-rs/lib/src/utils.rs b/bezier-rs/lib/src/utils.rs index 3bb7587b..616edc7c 100644 --- a/bezier-rs/lib/src/utils.rs +++ b/bezier-rs/lib/src/utils.rs @@ -1,4 +1,5 @@ -use glam::DVec2; +use glam::{BVec2, DVec2}; +use std::f64::consts::PI; /// Helper to perform the computation of a and c, where b is the provided point on the curve. /// Given the correct power of `t` and `(1-t)`, the computation is the same for quadratic and cubic cases. @@ -11,8 +12,8 @@ fn compute_abc_through_points(start_point: DVec2, point_on_curve: DVec2, end_poi [a, point_on_curve, c] } -/// Compute a, b, and c for a quadratic curve that fits the start, end and point on curve at `t`. -/// The definition for the a, b, c points are defined in [the projection identity section](https://pomax.github.io/bezierinfo/#abc) of Pomax's bezier curve primer. +/// Compute `a`, `b`, and `c` for a quadratic curve that fits the start, end and point on curve at `t`. +/// The definition for the `a`, `b`, `c` points are defined in [the projection identity section](https://pomax.github.io/bezierinfo/#abc) of Pomax's bezier curve primer. pub fn compute_abc_for_quadratic_through_points(start_point: DVec2, point_on_curve: DVec2, end_point: DVec2, t: f64) -> [DVec2; 3] { let t_squared = t * t; let one_minus_t = 1. - t; @@ -30,15 +31,16 @@ pub fn compute_abc_for_cubic_through_points(start_point: DVec2, point_on_curve: compute_abc_through_points(start_point, point_on_curve, end_point, t_cubed, cubed_one_minus_t) } +/// Return the index and the value of the closest point in the LUT compared to the provided point. pub fn get_closest_point_in_lut(lut: &[DVec2], point: DVec2) -> (i32, f64) { lut.iter() .enumerate() - .map(|(i, p)| (i as i32, point.distance(*p))) + .map(|(i, p)| (i as i32, point.distance_squared(*p))) .min_by(|x, y| (&(x.1)).partial_cmp(&(y.1)).unwrap()) .unwrap() } -/// Find the roots of the linear equation `ax + b` +/// Find the roots of the linear equation `ax + b`. pub fn solve_linear(a: f64, b: f64) -> Vec { let mut roots = Vec::new(); if a != 0. { @@ -47,8 +49,8 @@ pub fn solve_linear(a: f64, b: f64) -> Vec { roots } -/// Find the roots of the linear equation `ax^2 + bx + c` -/// Precompute the `discriminant` (`b^2 - 4ac`) and `two_times_a` arguments prior to calling this function for efficiency purposes +/// Find the roots of the linear equation `ax^2 + bx + c`. +/// Precompute the `discriminant` (`b^2 - 4ac`) and `two_times_a` arguments prior to calling this function for efficiency purposes. pub fn solve_quadratic(discriminant: f64, two_times_a: f64, b: f64, c: f64) -> Vec { let mut roots = Vec::new(); if two_times_a != 0. { @@ -64,3 +66,138 @@ pub fn solve_quadratic(discriminant: f64, two_times_a: f64, b: f64, c: f64) -> V } roots } + +/// Compute the cube root of a number. +fn cube_root(f: f64) -> f64 { + if f < 0. { + -(-f).powf(1. / 3.) + } else { + f.powf(1. / 3.) + } +} + +/// Solve a cubic of the form `x^3 + px + q`, derivation from: . +pub fn solve_reformatted_cubic(discriminant: f64, a: f64, p: f64, q: f64) -> Vec { + let mut roots = Vec::new(); + if p == 0. { + roots.push(cube_root(-q)); + } else if q == 0. { + if p < 0. { + roots.push((-p).powf(1. / 2.)); + } + } else if discriminant == 0. { + let q_divided_by_2 = q / 2.; + let a_divided_by_3 = a / 3.; + // all roots are real, and 2 are repeated + 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); + } else if discriminant > 0. { + // one real and two imaginary roots + let q_divided_by_2 = q / 2.; + let square_root_discriminant = discriminant.powf(1. / 2.); + roots.push(cube_root(-q_divided_by_2 + square_root_discriminant) - cube_root(q_divided_by_2 + square_root_discriminant) - a / 3.); + } else { + // three real roots + let p_divided_by_3 = p / 3.; + let a_divided_by_3 = a / 3.; + let cube_root_r = (-p_divided_by_3).powf(1. / 2.); + let phi = (-q / (2. * cube_root_r.powi(3))).acos(); + + let two_times_cube_root_r = 2. * cube_root_r; + // three real roots + roots.push(two_times_cube_root_r * (phi / 3.).cos() - a_divided_by_3); + roots.push(two_times_cube_root_r * ((phi + 2. * PI) / 3.).cos() - a_divided_by_3); + roots.push(two_times_cube_root_r * ((phi + 4. * PI) / 3.).cos() - a_divided_by_3); + } + roots +} + +/// Solve a cubic of the form `ax^3 + bx^2 + ct + d`. +pub fn solve_cubic(a: f64, b: f64, c: f64, d: f64) -> Vec { + if a.abs() <= 1e-5 { + if b.abs() <= 1e-5 { + // if both a and b are approximately 0, treat as a linear problem + solve_linear(c, d) + } else { + // if a is approximately 0, treat as a quadratic problem + let discriminant = c * c - 4. * b * d; + solve_quadratic(discriminant, 2. * b, c, d) + } + } else { + let new_a = b / a; + let new_b = c / a; + let new_c = d / a; + + // Refactor cubic to be of the form: a(t^3 + pt + q), derivation from: https://trans4mind.com/personal_development/mathematics/polynomials/cubicAlgebra.htm + let p = (3. * new_b - new_a * new_a) / 3.; + let q = (2. * new_a.powi(3) - 9. * new_a * new_b + 27. * new_c) / 27.; + let discriminant = (p / 3.).powi(3) + (q / 2.).powi(2); + solve_reformatted_cubic(discriminant, new_a, p, q) + } +} + +/// Compare two `f64` numbers with a provided max absolute value difference. +pub fn f64_compare(f1: f64, f2: f64, max_abs_diff: f64) -> bool { + (f1 - f2).abs() < max_abs_diff +} + +/// Determine if an `f64` number is within a given range by using a max absolute value difference comparison. +pub fn f64_approximately_in_range(value: f64, min: f64, max: f64, max_abs_diff: f64) -> bool { + (min..=max).contains(&value) || f64_compare(value, min, max_abs_diff) || f64_compare(value, max, max_abs_diff) +} + +/// Compare the two values in a `DVec2` independently with a provided max absolute value difference. +pub fn dvec2_compare(dv1: DVec2, dv2: DVec2, max_abs_diff: f64) -> BVec2 { + BVec2::new((dv1.x - dv2.x).abs() < max_abs_diff, (dv1.y - dv2.y).abs() < max_abs_diff) +} + +/// Determine if the values in a `DVec2` are within a given range independently by using a max absolute value difference comparison. +pub fn dvec2_approximately_in_range(point: DVec2, min: DVec2, max: DVec2, max_abs_diff: f64) -> BVec2 { + (point.cmpge(min) & point.cmple(max)) | dvec2_compare(point, min, max_abs_diff) | dvec2_compare(point, max, max_abs_diff) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_solve_cubic() { + // discriminant == 0 + let roots1 = solve_cubic(1., 0., 0., 0.); + assert!(roots1.len() == 1); + assert!(roots1[0] == 0.); + + let roots2 = solve_cubic(1., 3., 0., -4.); + assert!(roots2.len() == 2); + assert!(roots2[0] == 1.); + assert!(roots2[1] == -2.); + + // p == 0 + let roots3 = solve_cubic(1., 0., 0., -1.); + assert!(roots3.len() == 1); + assert!(roots3[0] == 1.); + + // discriminant > 0 + let roots4 = solve_cubic(1., 3., 0., 2.); + assert!(roots4.len() == 1); + assert!(f64_compare(roots4[0], -3.196, 1e-3)); + + // discriminant < 0 + let roots5 = solve_cubic(1., 3., 0., -1.); + assert!(roots5.len() == 3); + assert!(f64_compare(roots5[0], 0.532, 1e-3)); + assert!(f64_compare(roots5[1], -2.879, 1e-3)); + assert!(f64_compare(roots5[2], -0.653, 1e-3)); + + // quadratic + let roots6 = solve_cubic(0., 3., 0., -3.); + assert!(roots6.len() == 2); + assert!(roots6[0] == 1.); + assert!(roots6[1] == -1.); + + // linear + let roots7 = solve_cubic(0., 0., 1., -1.); + assert!(roots7.len() == 1); + assert!(roots7[0] == 1.); + } +}