From beb88d280c59bb1017855367657a049d0ecf787b Mon Sep 17 00:00:00 2001 From: Elbert Ronnie <103196773+elbertronnie@users.noreply.github.com> Date: Tue, 30 Apr 2024 07:49:12 +0530 Subject: [PATCH] Bezier-rs: Add calculations for area and centroid of subpaths (#1729) * add code to calculate area and centroid of subpath * change library to poly_it * add a demo of area and centroid * modify algorithm to consider negetive area as positive * add code for manipulating polynomials in bezier-rs * remove `poly_it` dependency and use custom Polynomial * formatting floats to skip last zero * add debug mechanism * collect both intersection points instead of one * fix test and cargo fmt * apply minimum separation filtering in self_intersection and use better endpoint filtering algorithm * remove debug mechanism and cargo fmt * consider the subpath as closed for intersection calculation * add documentation for polynomial.rs * impl display for Polynomial * make area always positive * add missing docs * fix test and cargo fmt * Naming/formatting code review --------- Co-authored-by: Keavon Chambers --- libraries/bezier-rs/src/bezier/solvers.rs | 97 ++++++- libraries/bezier-rs/src/lib.rs | 1 + libraries/bezier-rs/src/polynomial.rs | 264 ++++++++++++++++++ libraries/bezier-rs/src/subpath/core.rs | 15 +- libraries/bezier-rs/src/subpath/lookup.rs | 151 ++++++++++ libraries/bezier-rs/src/subpath/mod.rs | 4 +- libraries/bezier-rs/src/subpath/solvers.rs | 32 ++- .../src/features/bezier-features.ts | 4 +- .../src/features/subpath-features.ts | 10 + .../other/bezier-rs-demos/wasm/src/bezier.rs | 4 +- .../other/bezier-rs-demos/wasm/src/subpath.rs | 10 + 11 files changed, 577 insertions(+), 15 deletions(-) create mode 100644 libraries/bezier-rs/src/polynomial.rs diff --git a/libraries/bezier-rs/src/bezier/solvers.rs b/libraries/bezier-rs/src/bezier/solvers.rs index 6330d935..849a7fab 100644 --- a/libraries/bezier-rs/src/bezier/solvers.rs +++ b/libraries/bezier-rs/src/bezier/solvers.rs @@ -1,4 +1,5 @@ use super::*; +use crate::polynomial::Polynomial; use crate::utils::{solve_cubic, solve_quadratic, TValue}; use crate::{to_symmetrical_basis_pair, SymmetricalBasis}; @@ -40,6 +41,31 @@ impl Bezier { de_casteljau_points } + /// Returns two [`Polynomial`]s representing the parametric equations for x and y coordinates of the bezier curve respectively. + /// The domain of both the equations are from t=0.0 representing the start and t=1.0 representing the end of the bezier curve. + pub fn parametric_polynomial(&self) -> (Polynomial<4>, Polynomial<4>) { + match self.handles { + BezierHandles::Linear => { + let term1 = self.end - self.start; + + (Polynomial::new([self.start.x, term1.x, 0., 0.]), Polynomial::new([self.start.y, term1.y, 0., 0.])) + } + BezierHandles::Quadratic { handle } => { + let term1 = 2. * (handle - self.start); + let term2 = self.start - 2. * handle + self.end; + + (Polynomial::new([self.start.x, term1.x, term2.x, 0.]), Polynomial::new([self.start.y, term1.y, term2.y, 0.])) + } + BezierHandles::Cubic { handle_start, handle_end } => { + let term1 = 3. * (handle_start - self.start); + let term2 = 3. * (handle_end - handle_start) - term1; + let term3 = self.end - self.start - term2 - term1; + + (Polynomial::new([self.start.x, term1.x, term2.x, term3.x]), Polynomial::new([self.start.y, term1.y, term2.y, term3.y])) + } + } + } + /// Returns a [Bezier] representing the derivative of the original curve. /// - This function returns `None` for a linear segment. /// @@ -337,10 +363,36 @@ impl Bezier { let mut intersection_t_values = self.unfiltered_intersections(other, error); intersection_t_values.sort_by(|a, b| a.partial_cmp(b).unwrap()); - intersection_t_values.iter().fold(Vec::new(), |mut accumulator, t| { + intersection_t_values.iter().map(|x| x[0]).fold(Vec::new(), |mut accumulator, t| { if !accumulator.is_empty() && (accumulator.last().unwrap() - t).abs() < minimum_separation.unwrap_or(MIN_SEPARATION_VALUE) { accumulator.pop(); } + accumulator.push(t); + accumulator + }) + } + + // TODO: Use an `impl Iterator` return type instead of a `Vec` + /// Returns a list of pairs of filtered parametric `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 separation value. If the difference + /// between 2 adjacent `t` values is less than the minimum difference, the filtering takes the larger `t` value and discards the smaller `t` value. + /// The first value in pair is with respect to the current bezier and the second value in pair is with respect to 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_separation` - The minimum difference between adjacent `t` values in sorted order + pub fn all_intersections(&self, other: &Bezier, error: Option, minimum_separation: Option) -> Vec<[f64; 2]> { + // 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[0] + a[1]).partial_cmp(&(b[0] + b[1])).unwrap()); + + intersection_t_values.iter().fold(Vec::new(), |mut accumulator, t| { + if !accumulator.is_empty() + && (accumulator.last().unwrap()[0] - t[0]).abs() < minimum_separation.unwrap_or(MIN_SEPARATION_VALUE) + && (accumulator.last().unwrap()[1] - t[1]).abs() < minimum_separation.unwrap_or(MIN_SEPARATION_VALUE) + { + accumulator.pop(); + } accumulator.push(*t); accumulator }) @@ -350,7 +402,7 @@ impl Bezier { /// 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. - fn unfiltered_intersections(&self, other: &Bezier, error: Option) -> Vec { + pub fn unfiltered_intersections(&self, other: &Bezier, error: Option) -> Vec<[f64; 2]> { 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 @@ -363,6 +415,9 @@ impl Bezier { let vertical_distance = (rotation_matrix * other.start).x; let translated_bezier = rotated_bezier.translate(DVec2::new(-vertical_distance, 0.)); + let y_start = (rotation_matrix * other.start).y; + let y_end = (rotation_matrix * other.end).y; + // Compute the roots of the resulting bezier curve let list_intersection_t = translated_bezier.find_tvalues_for_x(0.); @@ -374,12 +429,17 @@ impl Bezier { .filter(|&t| utils::dvec2_approximately_in_range(self.unrestricted_parametric_evaluate(t), min_corner, max_corner, MAX_ABSOLUTE_DIFFERENCE).all()) // Ensure the returned value is within the correct range .map(|t| t.clamp(0., 1.)) - .collect::>(); + .map(|t| { + let y = translated_bezier.evaluate(TValue::Parametric(t)).y; + let other_t = (y-y_start)/(y_end-y_start); + [t, other_t] + }) + .collect::>(); } // TODO: Consider using the `intersections_between_vectors_of_curves` helper function here // Otherwise, use bounding box to determine intersections - self.intersections_between_subcurves(0. ..1., other, 0. ..1., error).iter().map(|t_values| t_values[0]).collect() + self.intersections_between_subcurves(0. ..1., other, 0. ..1., error).to_vec() } /// Returns a list of `t` values that correspond to points on this Bezier segment where they intersect with the given line. (`direction_vector` does not need to be normalized.) @@ -452,7 +512,7 @@ impl Bezier { /// Returns a list of parametric `t` values that correspond to the self intersection points of the current bezier curve. For each intersection point, the returned `t` value is the smaller of the two that correspond to the point. /// - `error` - For intersections with non-linear beziers, `error` defines the threshold for bounding boxes to be considered an intersection point. /// - pub fn self_intersections(&self, error: Option) -> Vec<[f64; 2]> { + fn unfiltered_self_intersections(&self, error: Option) -> Vec<[f64; 2]> { if self.handles == BezierHandles::Linear || matches!(self.handles, BezierHandles::Quadratic { .. }) { return vec![]; } @@ -482,6 +542,27 @@ impl Bezier { .collect() } + // TODO: Use an `impl Iterator` return type instead of a `Vec` + /// Returns a list of parametric `t` values that correspond to the self intersection points of the current bezier curve. For each intersection point, the returned `t` value is the smaller of the two that correspond to the point. + /// If the difference between 2 adjacent `t` values is less than the minimum difference, the filtering takes the larger `t` value and discards the smaller `t` value. + /// - `error` - For intersections with non-linear beziers, `error` defines the threshold for bounding boxes to be considered an intersection point. + /// - `minimum_separation` - The minimum difference between adjacent `t` values in sorted order + pub fn self_intersections(&self, error: Option, minimum_separation: Option) -> Vec<[f64; 2]> { + let mut intersection_t_values = self.unfiltered_self_intersections(error); + intersection_t_values.sort_by(|a, b| (a[0] + a[1]).partial_cmp(&(b[0] + b[1])).unwrap()); + + intersection_t_values.iter().fold(Vec::new(), |mut accumulator, t| { + if !accumulator.is_empty() + && (accumulator.last().unwrap()[0] - t[0]).abs() < minimum_separation.unwrap_or(MIN_SEPARATION_VALUE) + && (accumulator.last().unwrap()[1] - t[1]).abs() < minimum_separation.unwrap_or(MIN_SEPARATION_VALUE) + { + accumulator.pop(); + } + accumulator.push(*t); + accumulator + }) + } + /// Returns a list of parametric `t` values that correspond to the intersection points between the curve and a rectangle defined by opposite corners. /// pub fn rectangle_intersections(&self, corner1: DVec2, corner2: DVec2) -> Vec { @@ -1062,13 +1143,13 @@ mod tests { #[test] fn test_intersect_with_self() { let bezier = Bezier::from_cubic_coordinates(160., 180., 170., 10., 30., 90., 180., 140.); - let intersections = bezier.self_intersections(Some(0.5)); + let intersections = bezier.self_intersections(Some(0.5), None); assert!(compare_vec_of_points( intersections.iter().map(|&t| bezier.evaluate(TValue::Parametric(t[0]))).collect(), intersections.iter().map(|&t| bezier.evaluate(TValue::Parametric(t[1]))).collect(), 2. )); - assert!(Bezier::from_linear_coordinates(160., 180., 170., 10.).self_intersections(None).is_empty()); - assert!(Bezier::from_quadratic_coordinates(160., 180., 170., 10., 30., 90.).self_intersections(None).is_empty()); + assert!(Bezier::from_linear_coordinates(160., 180., 170., 10.).self_intersections(None, None).is_empty()); + assert!(Bezier::from_quadratic_coordinates(160., 180., 170., 10., 30., 90.).self_intersections(None, None).is_empty()); } } diff --git a/libraries/bezier-rs/src/lib.rs b/libraries/bezier-rs/src/lib.rs index 05cc9738..1630e28f 100644 --- a/libraries/bezier-rs/src/lib.rs +++ b/libraries/bezier-rs/src/lib.rs @@ -5,6 +5,7 @@ pub(crate) mod compare; mod bezier; mod consts; mod poisson_disk; +mod polynomial; mod subpath; mod symmetrical_basis; mod utils; diff --git a/libraries/bezier-rs/src/polynomial.rs b/libraries/bezier-rs/src/polynomial.rs new file mode 100644 index 00000000..4e8fdb20 --- /dev/null +++ b/libraries/bezier-rs/src/polynomial.rs @@ -0,0 +1,264 @@ +use std::fmt::{self, Display, Formatter}; +use std::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; + +/// A struct that represents a polynomial with a maximum degree of `N-1`. +/// +/// It provides basic mathematical operations for polynomials like addition, multiplication, differentiation, integration, etc. +#[derive(Copy, Clone, Debug, PartialEq)] +pub struct Polynomial { + coefficients: [f64; N], +} + +impl Polynomial { + /// Create a new polynomial from the coefficients given in the array. + /// + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. + pub fn new(coefficients: [f64; N]) -> Polynomial { + Polynomial { coefficients } + } + + /// Create a polynomial where all its coefficients are zero. + pub fn zero() -> Polynomial { + Polynomial { coefficients: [0.; N] } + } + + /// Return an immutable reference to the coefficients. + /// + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. + pub fn coefficients(&self) -> &[f64; N] { + &self.coefficients + } + + /// Return a mutable reference to the coefficients. + /// + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. + pub fn coefficients_mut(&mut self) -> &mut [f64; N] { + &mut self.coefficients + } + + /// Evaluate the polynomial at `value`. + pub fn eval(&self, value: f64) -> f64 { + self.coefficients.iter().rev().copied().reduce(|acc, x| acc * value + x).unwrap() + } + + /// Return the same polynomial but with a different maximum degree of `M-1`.\ + /// + /// Returns `None` if the polynomial cannot fit in the specified size. + pub fn as_size(&self) -> Option> { + let mut coefficients = [0.; M]; + + if M >= N { + coefficients[..N].copy_from_slice(&self.coefficients); + } else if self.coefficients.iter().rev().take(N - M).all(|&x| x == 0.) { + coefficients.copy_from_slice(&self.coefficients[..M]) + } else { + return None; + } + + Some(Polynomial { coefficients }) + } + + /// Computes the derivative in place. + pub fn derivative_mut(&mut self) { + self.coefficients.iter_mut().enumerate().for_each(|(index, x)| *x *= index as f64); + self.coefficients.rotate_left(1); + } + + /// Computes the antiderivative at `C = 0` in place. + /// + /// Returns `None` if the polynomial is not big enough to accommodate the extra degree. + pub fn antiderivative_mut(&mut self) -> Option<()> { + if self.coefficients[N - 1] != 0. { + return None; + } + self.coefficients.rotate_right(1); + self.coefficients.iter_mut().enumerate().skip(1).for_each(|(index, x)| *x /= index as f64); + Some(()) + } + + /// Computes the polynomial's derivative. + pub fn derivative(&self) -> Polynomial { + let mut ans = *self; + ans.derivative_mut(); + ans + } + + /// Computes the antiderivative at `C = 0`. + /// + /// Returns `None` if the polynomial is not big enough to accommodate the extra degree. + pub fn antiderivative(&self) -> Option> { + let mut ans = *self; + ans.antiderivative_mut()?; + Some(ans) + } +} + +impl Default for Polynomial { + fn default() -> Self { + Self::zero() + } +} + +impl Display for Polynomial { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + let mut first = true; + for (index, coefficient) in self.coefficients.iter().enumerate().rev().filter(|(_, &coefficient)| coefficient != 0.) { + if first { + first = false; + } else { + f.write_str(" + ")? + } + + coefficient.fmt(f)?; + if index == 0 { + continue; + } + f.write_str("x")?; + if index == 1 { + continue; + } + f.write_str("^")?; + index.fmt(f)?; + } + + Ok(()) + } +} + +impl AddAssign<&Polynomial> for Polynomial { + fn add_assign(&mut self, rhs: &Polynomial) { + self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a += b); + } +} + +impl Add for &Polynomial { + type Output = Polynomial; + + fn add(self, other: &Polynomial) -> Polynomial { + let mut output = *self; + output += other; + output + } +} + +impl Neg for &Polynomial { + type Output = Polynomial; + + fn neg(self) -> Polynomial { + let mut output = *self; + output.coefficients.iter_mut().for_each(|x| *x = -*x); + output + } +} + +impl Neg for Polynomial { + type Output = Polynomial; + + fn neg(mut self) -> Polynomial { + self.coefficients.iter_mut().for_each(|x| *x = -*x); + self + } +} + +impl SubAssign<&Polynomial> for Polynomial { + fn sub_assign(&mut self, rhs: &Polynomial) { + self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a -= b); + } +} + +impl Sub for &Polynomial { + type Output = Polynomial; + + fn sub(self, other: &Polynomial) -> Polynomial { + let mut output = *self; + output -= other; + output + } +} + +impl MulAssign<&Polynomial> for Polynomial { + fn mul_assign(&mut self, rhs: &Polynomial) { + for i in (0..N).rev() { + self.coefficients[i] = self.coefficients[i] * rhs.coefficients[0]; + for j in 0..i { + self.coefficients[i] += self.coefficients[j] * rhs.coefficients[i - j]; + } + } + } +} + +impl Mul for &Polynomial { + type Output = Polynomial; + + fn mul(self, other: &Polynomial) -> Polynomial { + let mut output = *self; + output *= other; + output + } +} + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn evaluation() { + let p = Polynomial::new([1., 2., 3.]); + + assert_eq!(p.eval(1.), 6.); + assert_eq!(p.eval(2.), 17.); + } + + #[test] + fn size_change() { + let p1 = Polynomial::new([1., 2., 3.]); + let p2 = Polynomial::new([1., 2., 3., 0.]); + + assert_eq!(p1.as_size(), Some(p2)); + assert_eq!(p2.as_size(), Some(p1)); + + assert_eq!(p2.as_size::<2>(), None); + } + + #[test] + fn addition_and_subtaction() { + let p1 = Polynomial::new([1., 2., 3.]); + let p2 = Polynomial::new([4., 5., 6.]); + + let addition = Polynomial::new([5., 7., 9.]); + let subtraction = Polynomial::new([-3., -3., -3.]); + + assert_eq!(&p1 + &p2, addition); + assert_eq!(&p1 - &p2, subtraction); + } + + #[test] + fn multiplication() { + let p1 = Polynomial::new([1., 2., 3.]).as_size().unwrap(); + let p2 = Polynomial::new([4., 5., 6.]).as_size().unwrap(); + + let multiplication = Polynomial::new([4., 13., 28., 27., 18.]); + + assert_eq!(&p1 * &p2, multiplication); + } + + #[test] + fn derivative_and_antiderivative() { + let mut p = Polynomial::new([1., 2., 3.]); + let p_deriv = Polynomial::new([2., 6., 0.]); + + assert_eq!(p.derivative(), p_deriv); + + p.coefficients_mut()[0] = 0.; + assert_eq!(p_deriv.antiderivative().unwrap(), p); + + assert_eq!(p.antiderivative(), None); + } + + #[test] + fn display() { + let p = Polynomial::new([1., 2., 0., 3.]); + + assert_eq!(format!("{:.2}", p), "3.00x^3 + 2.00x + 1.00"); + } +} diff --git a/libraries/bezier-rs/src/subpath/core.rs b/libraries/bezier-rs/src/subpath/core.rs index bd59d5a9..ce17cda4 100644 --- a/libraries/bezier-rs/src/subpath/core.rs +++ b/libraries/bezier-rs/src/subpath/core.rs @@ -105,7 +105,20 @@ impl Subpath { /// Returns an iterator of the [Bezier]s along the `Subpath`. pub fn iter(&self) -> SubpathIter { - SubpathIter { subpath: self, index: 0 } + SubpathIter { + subpath: self, + index: 0, + is_always_closed: false, + } + } + + /// Returns an iterator of the [Bezier]s along the `Subpath` always considering it as a closed subpath. + pub fn iter_closed(&self) -> SubpathIter { + SubpathIter { + subpath: self, + index: 0, + is_always_closed: true, + } } /// Returns a slice of the [ManipulatorGroup]s in the `Subpath`. diff --git a/libraries/bezier-rs/src/subpath/lookup.rs b/libraries/bezier-rs/src/subpath/lookup.rs index 24e050ec..9c09b53f 100644 --- a/libraries/bezier-rs/src/subpath/lookup.rs +++ b/libraries/bezier-rs/src/subpath/lookup.rs @@ -30,6 +30,91 @@ impl Subpath { self.iter().map(|bezier| bezier.length(tolerance)).sum() } + /// Return the area enclosed by the `Subpath` always considering it as a closed subpath. It will always give a positive value. + /// + /// Because the calculation of area for self-intersecting path requires finding the intersections, the following parameters are used: + /// - `error` - For intersections with non-linear beziers, `error` defines the threshold for bounding boxes to be considered an intersection point. + /// - `minimum_separation`: the minimum difference two adjacent `t`-values must have when comparing adjacent `t`-values in sorted order. + /// If the comparison condition is not satisfied, the function takes the larger `t`-value of the two + /// + /// **NOTE**: if an intersection were to occur within an `error` distance away from an anchor point, the algorithm will filter that intersection out. + pub fn area(&self, error: Option, minimum_separation: Option) -> f64 { + let all_intersections = self.all_self_intersections(error, minimum_separation); + let mut current_sign: f64 = 1.; + + let area: f64 = self + .iter_closed() + .enumerate() + .map(|(index, bezier)| { + let (f_x, f_y) = bezier.parametric_polynomial(); + let (f_x, mut f_y) = (f_x.as_size::<7>().unwrap(), f_y.as_size::<7>().unwrap()); + f_y.derivative_mut(); + f_y *= &f_x; + f_y.antiderivative_mut(); + + let mut curve_sum = -current_sign * f_y.eval(0.); + for (_, t) in all_intersections.iter().filter(|(i, _)| *i == index) { + curve_sum += 2. * current_sign * f_y.eval(*t); + current_sign *= -1.; + } + curve_sum += current_sign * f_y.eval(1.); + curve_sum + }) + .sum(); + + area.abs() + } + + /// Return the centroid of the `Subpath` always considering it as a closed subpath. + /// It will return `None` if no manipulator is present. + /// + /// Because the calculation of area for self-intersecting path requires finding the intersections, the following parameters are used: + /// - `error` - For intersections with non-linear beziers, `error` defines the threshold for bounding boxes to be considered an intersection point. + /// - `minimum_separation`: the minimum difference two adjacent `t`-values must have when comparing adjacent `t`-values in sorted order. + /// If the comparison condition is not satisfied, the function takes the larger `t`-value of the two + /// + /// **NOTE**: if an intersection were to occur within an `error` distance away from an anchor point, the algorithm will filter that intersection out. + pub fn centroid(&self, error: Option, minimum_separation: Option) -> Option { + let all_intersections = self.all_self_intersections(error, minimum_separation); + let mut current_sign: f64 = 1.; + + let (x_sum, y_sum, area) = self + .iter_closed() + .enumerate() + .map(|(index, bezier)| { + let (f_x, f_y) = bezier.parametric_polynomial(); + let (f_x, f_y) = (f_x.as_size::<10>().unwrap(), f_y.as_size::<10>().unwrap()); + let f_y_prime = f_y.derivative(); + let f_x_prime = f_x.derivative(); + let f_xy = &f_x * &f_y; + + let mut x_part = &f_xy * &f_x_prime; + let mut y_part = &f_xy * &f_y_prime; + let mut area_part = &f_x * &f_y_prime; + x_part.antiderivative_mut(); + y_part.antiderivative_mut(); + area_part.antiderivative_mut(); + + let mut curve_sum_x = -current_sign * x_part.eval(0.); + let mut curve_sum_y = -current_sign * y_part.eval(0.); + let mut curve_sum_area = -current_sign * area_part.eval(0.); + for (_, t) in all_intersections.iter().filter(|(i, _)| *i == index) { + curve_sum_x += 2. * current_sign * x_part.eval(*t); + curve_sum_y += 2. * current_sign * y_part.eval(*t); + curve_sum_area += 2. * current_sign * area_part.eval(*t); + current_sign *= -1.; + } + curve_sum_x += current_sign * x_part.eval(1.); + curve_sum_y += current_sign * y_part.eval(1.); + curve_sum_area += current_sign * area_part.eval(1.); + + (-curve_sum_x, curve_sum_y, curve_sum_area) + }) + .reduce(|(x1, y1, area1), (x2, y2, area2)| (x1 + x2, y1 + y2, area1 + area2))?; + + Some(DVec2::new(x_sum / area, y_sum / area)) + } + /// Converts from a subpath (composed of multiple segments) to a point along a certain segment represented. /// The returned tuple represents the segment index and the `t` value along that segment. /// Both the input global `t` value and the output `t` value are in euclidean space, meaning there is a constant rate of change along the arc length. @@ -208,6 +293,72 @@ mod tests { assert_eq!(subpath.length(None), linear_bezier.length(None) + quadratic_bezier.length(None) + cubic_bezier.length(None)); } + #[test] + fn area() { + let start = DVec2::new(0., 0.); + let end = DVec2::new(1., 1.); + let handle = DVec2::new(0., 1.); + + let mut subpath = Subpath::new( + vec![ + ManipulatorGroup { + anchor: start, + in_handle: None, + out_handle: Some(handle), + id: EmptyId, + }, + ManipulatorGroup { + anchor: end, + in_handle: None, + out_handle: None, + id: EmptyId, + }, + ], + false, + ); + + let expected_area = 1. / 3.; + let epsilon = 0.00001; + + assert!((subpath.area(Some(0.001), Some(0.001)) - expected_area).abs() < epsilon); + + subpath.closed = true; + assert!((subpath.area(Some(0.001), Some(0.001)) - expected_area).abs() < epsilon); + } + + #[test] + fn centroid() { + let start = DVec2::new(0., 0.); + let end = DVec2::new(1., 1.); + let handle = DVec2::new(0., 1.); + + let mut subpath = Subpath::new( + vec![ + ManipulatorGroup { + anchor: start, + in_handle: None, + out_handle: Some(handle), + id: EmptyId, + }, + ManipulatorGroup { + anchor: end, + in_handle: None, + out_handle: None, + id: EmptyId, + }, + ], + false, + ); + + let expected_centroid = DVec2::new(0.4, 0.6); + let epsilon = 0.00001; + + assert!(subpath.centroid(Some(0.001), Some(0.001)).unwrap().abs_diff_eq(expected_centroid, epsilon)); + + subpath.closed = true; + assert!(subpath.centroid(Some(0.001), Some(0.001)).unwrap().abs_diff_eq(expected_centroid, epsilon)); + } + #[test] fn t_value_to_parametric_global_parametric_open_subpath() { let mock_manipulator_group = ManipulatorGroup { diff --git a/libraries/bezier-rs/src/subpath/mod.rs b/libraries/bezier-rs/src/subpath/mod.rs index 78e5bca8..854022bb 100644 --- a/libraries/bezier-rs/src/subpath/mod.rs +++ b/libraries/bezier-rs/src/subpath/mod.rs @@ -29,6 +29,7 @@ unsafe impl dyn_any::StaticType for Subpa pub struct SubpathIter<'a, ManipulatorGroupId: crate::Identifier> { index: usize, subpath: &'a Subpath, + is_always_closed: bool, } impl Index for Subpath { @@ -55,8 +56,9 @@ impl Iterator for SubpathIter<'_, Manipul if self.subpath.is_empty() { return None; } + let closed = if self.is_always_closed { true } else { self.subpath.closed }; let len = self.subpath.len() - 1 - + match self.subpath.closed { + + match closed { true => 1, false => 0, }; diff --git a/libraries/bezier-rs/src/subpath/solvers.rs b/libraries/bezier-rs/src/subpath/solvers.rs index d2a26aaf..a88f1f0b 100644 --- a/libraries/bezier-rs/src/subpath/solvers.rs +++ b/libraries/bezier-rs/src/subpath/solvers.rs @@ -116,7 +116,7 @@ impl Subpath { let err = error.unwrap_or(MAX_ABSOLUTE_DIFFERENCE); // TODO: optimization opportunity - this for-loop currently compares all intersections with all curve-segments in the subpath collection self.iter().enumerate().for_each(|(i, other)| { - intersections_vec.extend(other.self_intersections(error).iter().map(|value| (i, value[0]))); + intersections_vec.extend(other.self_intersections(error, minimum_separation).iter().map(|value| (i, value[0]))); self.iter().enumerate().skip(i + 1).for_each(|(j, curve)| { intersections_vec.extend( curve @@ -130,6 +130,36 @@ impl Subpath { intersections_vec } + /// Returns a list of `t` values that correspond to all the self intersection points of the subpath always considering it as a closed subpath. The index and `t` value of both will be returned that corresponds to a point. + /// The points will be sorted based on their index and `t` repsectively. + /// - `error` - For intersections with non-linear beziers, `error` defines the threshold for bounding boxes to be considered an intersection point. + /// - `minimum_separation`: the minimum difference two adjacent `t`-values must have when comparing adjacent `t`-values in sorted order. + /// If the comparison condition is not satisfied, the function takes the larger `t`-value of the two + /// + /// **NOTE**: if an intersection were to occur within an `error` distance away from an anchor point, the algorithm will filter that intersection out. + pub fn all_self_intersections(&self, error: Option, minimum_separation: Option) -> Vec<(usize, f64)> { + let mut intersections_vec = Vec::new(); + let err = error.unwrap_or(MAX_ABSOLUTE_DIFFERENCE); + let num_curves = self.len(); + // TODO: optimization opportunity - this for-loop currently compares all intersections with all curve-segments in the subpath collection + self.iter_closed().enumerate().for_each(|(i, other)| { + intersections_vec.extend(other.self_intersections(error, minimum_separation).iter().flat_map(|value| [(i, value[0]), (i, value[1])])); + self.iter_closed().enumerate().skip(i + 1).for_each(|(j, curve)| { + intersections_vec.extend( + curve + .all_intersections(&other, error, minimum_separation) + .iter() + .filter(|&value| (j != i + 1 || value[0] > err || (1. - value[1]) > err) && (j != num_curves - 1 || i != 0 || value[1] > err || (1. - value[0]) > err)) + .flat_map(|value| [(j, value[0]), (i, value[1])]), + ); + }); + }); + + intersections_vec.sort_by(|a, b| a.partial_cmp(b).unwrap()); + + intersections_vec + } + /// Calculates the intersection points the subpath has with a given rectangle and returns a list of `(usize, f64)` tuples, /// where the `usize` represents the index of the curve in the subpath, and the `f64` represents the `t`-value local to /// that curve where the intersection occurred. diff --git a/website/other/bezier-rs-demos/src/features/bezier-features.ts b/website/other/bezier-rs-demos/src/features/bezier-features.ts index d824cc98..72b63e11 100644 --- a/website/other/bezier-rs-demos/src/features/bezier-features.ts +++ b/website/other/bezier-rs-demos/src/features/bezier-features.ts @@ -453,7 +453,7 @@ const bezierFeatures = { }, "intersect-self": { name: "Intersect (Self)", - callback: (bezier: WasmBezierInstance, options: Record): string => bezier.intersect_self(options.error), + callback: (bezier: WasmBezierInstance, options: Record): string => bezier.intersect_self(options.error, options.minimum_separation), demoOptions: { Linear: { disabled: true, @@ -462,7 +462,7 @@ const bezierFeatures = { disabled: true, }, Cubic: { - inputOptions: [errorOptions], + inputOptions: [errorOptions, minimumSeparationOptions], customPoints: [ [160, 180], [170, 10], diff --git a/website/other/bezier-rs-demos/src/features/subpath-features.ts b/website/other/bezier-rs-demos/src/features/subpath-features.ts index f12e0b39..004c31a7 100644 --- a/website/other/bezier-rs-demos/src/features/subpath-features.ts +++ b/website/other/bezier-rs-demos/src/features/subpath-features.ts @@ -16,6 +16,16 @@ const subpathFeatures = { name: "Length", callback: (subpath: WasmSubpathInstance): string => subpath.length(), }, + area: { + name: "Area", + callback: (subpath: WasmSubpathInstance, options: Record, _: undefined): string => subpath.area(options.error, options.minimum_separation), + inputOptions: [intersectionErrorOptions, minimumSeparationOptions], + }, + centroid: { + name: "Centroid", + callback: (subpath: WasmSubpathInstance, options: Record, _: undefined): string => subpath.centroid(options.error, options.minimum_separation), + inputOptions: [intersectionErrorOptions, minimumSeparationOptions], + }, evaluate: { name: "Evaluate", callback: (subpath: WasmSubpathInstance, options: Record, _: undefined): string => subpath.evaluate(options.t, SUBPATH_T_VALUE_VARIANTS[options.TVariant]), diff --git a/website/other/bezier-rs-demos/wasm/src/bezier.rs b/website/other/bezier-rs-demos/wasm/src/bezier.rs index 6d0e28e6..094c00f8 100644 --- a/website/other/bezier-rs-demos/wasm/src/bezier.rs +++ b/website/other/bezier-rs-demos/wasm/src/bezier.rs @@ -498,11 +498,11 @@ impl WasmBezier { } /// The wrapped return type is `Vec<[f64; 2]>`. - pub fn intersect_self(&self, error: f64) -> String { + pub fn intersect_self(&self, error: f64, minimum_separation: f64) -> String { let bezier_curve_svg = self.get_bezier_path(); let intersect_self_svg = self .0 - .self_intersections(Some(error)) + .self_intersections(Some(error), Some(minimum_separation)) .iter() .map(|intersection_t| { let point = &self.0.evaluate(TValue::Parametric(intersection_t[0])); diff --git a/website/other/bezier-rs-demos/wasm/src/subpath.rs b/website/other/bezier-rs-demos/wasm/src/subpath.rs index b115e685..9b5685af 100644 --- a/website/other/bezier-rs-demos/wasm/src/subpath.rs +++ b/website/other/bezier-rs-demos/wasm/src/subpath.rs @@ -92,6 +92,16 @@ impl WasmSubpath { wrap_svg_tag(format!("{}{}", self.to_default_svg(), length_text)) } + pub fn area(&self, error: f64, minimum_separation: f64) -> String { + let area_text = draw_text(format!("Area: {}", self.0.area(Some(error), Some(minimum_separation))), 5., 193., BLACK); + wrap_svg_tag(format!("{}{}", self.to_default_svg(), area_text)) + } + + pub fn centroid(&self, error: f64, minimum_separation: f64) -> String { + let point_text = draw_circle(self.0.centroid(Some(error), Some(minimum_separation)).unwrap(), 4., RED, 1.5, WHITE); + wrap_svg_tag(format!("{}{}", self.to_default_svg(), point_text)) + } + pub fn evaluate(&self, t: f64, t_variant: String) -> String { let t = parse_t_variant(&t_variant, t); let point = self.0.evaluate(t);