refactor centroid node

node-graph/gcore/src/math/mod.rs

@@ -1,4 +1,5 @@
 pub mod bbox;
 pub mod math_ext;
+pub mod polynomial;
 pub mod quad;
 pub mod rect;

node-graph/gcore/src/math/polynomial.rs

@@ -0,0 +1,293 @@
+use std::fmt::{self, Display, Formatter};
+use std::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
+
+use kurbo::PathSeg;
+
+/// 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<const N: usize> {
+	coefficients: [f64; N],
+}
+
+impl<const N: usize> Polynomial<N> {
+	/// 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<N> {
+		Polynomial { coefficients }
+	}
+
+	/// Create a polynomial where all its coefficients are zero.
+	pub fn zero() -> Polynomial<N> {
+		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<const M: usize>(&self) -> Option<Polynomial<M>> {
+		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<N> {
+		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<Polynomial<N>> {
+		let mut ans = *self;
+		ans.antiderivative_mut()?;
+		Some(ans)
+	}
+}
+
+impl<const N: usize> Default for Polynomial<N> {
+	fn default() -> Self {
+		Self::zero()
+	}
+}
+
+impl<const N: usize> Display for Polynomial<N> {
+	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<const N: usize> AddAssign<&Polynomial<N>> for Polynomial<N> {
+	fn add_assign(&mut self, rhs: &Polynomial<N>) {
+		self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a += b);
+	}
+}
+
+impl<const N: usize> Add for &Polynomial<N> {
+	type Output = Polynomial<N>;
+
+	fn add(self, other: &Polynomial<N>) -> Polynomial<N> {
+		let mut output = *self;
+		output += other;
+		output
+	}
+}
+
+impl<const N: usize> Neg for &Polynomial<N> {
+	type Output = Polynomial<N>;
+
+	fn neg(self) -> Polynomial<N> {
+		let mut output = *self;
+		output.coefficients.iter_mut().for_each(|x| *x = -*x);
+		output
+	}
+}
+
+impl<const N: usize> Neg for Polynomial<N> {
+	type Output = Polynomial<N>;
+
+	fn neg(mut self) -> Polynomial<N> {
+		self.coefficients.iter_mut().for_each(|x| *x = -*x);
+		self
+	}
+}
+
+impl<const N: usize> SubAssign<&Polynomial<N>> for Polynomial<N> {
+	fn sub_assign(&mut self, rhs: &Polynomial<N>) {
+		self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a -= b);
+	}
+}
+
+impl<const N: usize> Sub for &Polynomial<N> {
+	type Output = Polynomial<N>;
+
+	fn sub(self, other: &Polynomial<N>) -> Polynomial<N> {
+		let mut output = *self;
+		output -= other;
+		output
+	}
+}
+
+impl<const N: usize> MulAssign<&Polynomial<N>> for Polynomial<N> {
+	fn mul_assign(&mut self, rhs: &Polynomial<N>) {
+		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<const N: usize> Mul for &Polynomial<N> {
+	type Output = Polynomial<N>;
+
+	fn mul(self, other: &Polynomial<N>) -> Polynomial<N> {
+		let mut output = *self;
+		output *= other;
+		output
+	}
+}
+
+/// 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 pathseg_to_parametric_polynomial(segment: PathSeg) -> (Polynomial<4>, Polynomial<4>) {
+	match segment {
+		PathSeg::Line(line) => {
+			let term1 = line.p0 - line.p1;
+			(Polynomial::new([line.p0.x, term1.x, 0., 0.]), Polynomial::new([line.p0.y, term1.y, 0., 0.]))
+		}
+		PathSeg::Quad(quad_bez) => {
+			let term1 = 2. * (quad_bez.p1 - quad_bez.p0);
+			let term2 = quad_bez.p0 - 2. * quad_bez.p1.to_vec2() + quad_bez.p2.to_vec2();
+
+			(Polynomial::new([quad_bez.p0.x, term1.x, term2.x, 0.]), Polynomial::new([quad_bez.p0.y, term1.y, term2.y, 0.]))
+		}
+		PathSeg::Cubic(cubic_bez) => {
+			let term1 = 3. * (cubic_bez.p1 - cubic_bez.p0);
+			let term2 = 3. * (cubic_bez.p2 - cubic_bez.p1) - term1;
+			let term3 = cubic_bez.p3 - cubic_bez.p0 - term2 - term1;
+
+			(
+				Polynomial::new([cubic_bez.p0.x, term1.x, term2.x, term3.x]),
+				Polynomial::new([cubic_bez.p0.y, term1.y, term2.y, term3.y]),
+			)
+		}
+	}
+}
+
+#[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");
+	}
+}

node-graph/gcore/src/vector/algorithms/centroid.rs

@@ -0,0 +1,152 @@
+use glam::DVec2;
+use kurbo::{BezPath, CubicBez, ParamCurve, PathEl, PathSeg, QuadBez, Vec2};
+
+use crate::math::polynomial::pathseg_to_parametric_polynomial;
+
+use super::{contants::MAX_ABSOLUTE_DIFFERENCE, intersection::bezpath_all_self_intersections};
+
+/// Return the approximation of the length centroid, together with the length, of the `BezPath`.
+///
+/// The length centroid is the center of mass for the arc length of the solid shape's perimeter.
+/// An infinitely thin wire forming the subpath's closed shape would balance at this point.
+///
+/// It will return `None` if Bezpath has no segments.
+/// - `accuracy` is used to approximate the curve.
+/// - `always_closed` to consider the BezPath as closed always.
+pub fn bezpath_length_centroid_and_length(mut bezpath: BezPath, accuracy: Option<f64>, always_closed: bool) -> Option<(Vec2, f64)> {
+	// TODO: Take the Bezpath as a reference instead of value to avoid allocation. Presently we do it so we can close the path.
+
+	if !bezpath.elements().last().is_some_and(|element| *element == PathEl::ClosePath) && always_closed {
+		bezpath.close_path();
+	}
+
+	bezpath
+		.segments()
+		.map(|segment| pathseg_length_centroid_and_length(segment, accuracy))
+		.map(|(centroid, length)| (centroid * length, length))
+		.reduce(|(centroid_part1, length1), (centroid_part2, length2)| (centroid_part1 + centroid_part2, length1 + length2))
+		.map(|(centroid_part, length)| (centroid_part / length, length))
+}
+
+/// Return an approximation of the length centroid, together with the length, of the bezier curve.
+///
+/// The length centroid is the center of mass for the arc length of the Bezier segment.
+/// An infinitely thin wire forming the Bezier segment's shape would balance at this point.
+///
+/// - `accuracy` is used to approximate the curve.
+fn pathseg_length_centroid_and_length(segment: PathSeg, accuracy: Option<f64>) -> (Vec2, f64) {
+	match segment {
+		PathSeg::Line(line) => ((line.start().to_vec2() + line.end().to_vec2()) / 2., (line.start().to_vec2() - line.end().to_vec2()).length()),
+		PathSeg::Quad(quad_bez) => {
+			let QuadBez { p0, p1, p2 } = quad_bez;
+			// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
+			fn recurse(a0: Vec2, a1: Vec2, a2: Vec2, accuracy: f64, level: u8) -> (f64, Vec2) {
+				let lower = (a2 - a1).length();
+				let upper = (a1 - a0).length() + (a2 - a1).length();
+				if upper - lower <= 2. * accuracy || level >= 8 {
+					let length = (lower + upper) / 2.;
+					return (length, length * (a0 + a1 + a2) / 3.);
+				}
+
+				let b1 = 0.5 * (a0 + a1);
+				let c1 = 0.5 * (a1 + a2);
+				let b2 = 0.5 * (b1 + c1);
+
+				let (length1, centroid_part1) = recurse(a0, b1, b2, 0.5 * accuracy, level + 1);
+				let (length2, centroid_part2) = recurse(b2, c1, a2, 0.5 * accuracy, level + 1);
+				(length1 + length2, centroid_part1 + centroid_part2)
+			}
+
+			let (length, centroid_parts) = recurse(p0.to_vec2(), p1.to_vec2(), p2.to_vec2(), accuracy.unwrap_or_default(), 0);
+			(centroid_parts / length, length)
+		}
+		PathSeg::Cubic(cubic_bez) => {
+			let CubicBez { p0, p1, p2, p3 } = cubic_bez;
+
+			// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
+			fn recurse(a0: Vec2, a1: Vec2, a2: Vec2, a3: Vec2, accuracy: f64, level: u8) -> (f64, Vec2) {
+				let lower = (a3 - a0).length();
+				let upper = (a1 - a0).length() + (a2 - a1).length() + (a3 - a2).length();
+				if upper - lower <= 2. * accuracy || level >= 8 {
+					let length = (lower + upper) / 2.;
+					return (length, length * (a0 + a1 + a2 + a3) / 4.);
+				}
+
+				let b1 = 0.5 * (a0 + a1);
+				let t0 = 0.5 * (a1 + a2);
+				let c1 = 0.5 * (a2 + a3);
+				let b2 = 0.5 * (b1 + t0);
+				let c2 = 0.5 * (t0 + c1);
+				let b3 = 0.5 * (b2 + c2);
+
+				let (length1, centroid_part1) = recurse(a0, b1, b2, b3, 0.5 * accuracy, level + 1);
+				let (length2, centroid_part2) = recurse(b3, c2, c1, a3, 0.5 * accuracy, level + 1);
+				(length1 + length2, centroid_part1 + centroid_part2)
+			}
+
+			let (length, centroid_parts) = recurse(p0.to_vec2(), p1.to_vec2(), p2.to_vec2(), p3.to_vec2(), accuracy.unwrap_or_default(), 0);
+			(centroid_parts / length, length)
+		}
+	}
+}
+
+/// Return the area centroid, together with the area, of the `BezPath` always considering it as a closed subpath. The area will always be a positive value.
+///
+/// The area centroid is the center of mass for the area of a solid shape's interior.
+/// An infinitely flat material forming the subpath's closed shape would balance at this point.
+///
+/// It will return `None` if no manipulator is present. If the area is less than `error`, it will return `Some((DVec2::NAN, 0.))`.
+///
+/// Because the calculation of area and centroid 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 bezpath_area_centroid_and_area(mut bezpath: BezPath, error: Option<f64>, minimum_separation: Option<f64>) -> Option<(DVec2, f64)> {
+	let all_intersections = bezpath_all_self_intersections(bezpath.clone(), error, minimum_separation);
+	let mut current_sign: f64 = 1.;
+
+	bezpath.close_path();
+
+	let (x_sum, y_sum, area) = bezpath
+		.segments()
+		.enumerate()
+		.map(|(index, segment)| {
+			let (f_x, f_y) = pathseg_to_parametric_polynomial(segment);
+			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))?;
+
+	if area.abs() < error.unwrap_or(MAX_ABSOLUTE_DIFFERENCE) {
+		return Some((DVec2::NAN, 0.));
+	}
+
+	Some((DVec2::new(x_sum / area, y_sum / area), area.abs()))
+}

node-graph/gcore/src/vector/algorithms/contants.rs

@@ -1,6 +1,4 @@
 /// Minimum allowable separation between adjacent `t` values when calculating curve intersections
 pub const MIN_SEPARATION_VALUE: f64 = 5. * 1e-3;
-
 /// Constant used to determine if `f64`s are equivalent.
-#[cfg(test)]
 pub const MAX_ABSOLUTE_DIFFERENCE: f64 = 1e-3;

node-graph/gcore/src/vector/algorithms/intersection.rs

@@ -1,4 +1,4 @@
-use super::contants::MIN_SEPARATION_VALUE;
+use super::contants::{MAX_ABSOLUTE_DIFFERENCE, MIN_SEPARATION_VALUE};
 use kurbo::{BezPath, DEFAULT_ACCURACY, ParamCurve, PathSeg, Shape};
 
 /// Calculates the intersection points the bezpath has with a given segment and returns a list of `(usize, f64)` tuples,
@@ -125,6 +125,62 @@ pub fn filtered_all_segment_intersections(segment1: PathSeg, segment2: PathSeg,
 	})
 }
 
+// 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 pathseg_self_intersections(segment: PathSeg, accuracy: Option<f64>, minimum_separation: Option<f64>) -> Vec<(f64, f64)> {
+	let (first_half, second_half) = segment.subdivide();
+	let mut intersection_t_values = segment_intersections(first_half, second_half, accuracy);
+
+	intersection_t_values.sort_by(|a, b| (a.0 + a.1).partial_cmp(&(b.0 + b.1)).unwrap());
+
+	intersection_t_values.iter().filter(|(t1, t2)| !(*t1 == 0.5 && *t2 == 0.)).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 `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 bezpath_all_self_intersections(mut bezpath: BezPath, error: Option<f64>, minimum_separation: Option<f64>) -> Vec<(usize, f64)> {
+	// TODO: Take the Bezpath as a reference instead of value to avoid allocation. Presently we do it so we can close the path.
+	bezpath.close_path();
+
+	let mut intersections_vec = Vec::new();
+	let err = error.unwrap_or(MAX_ABSOLUTE_DIFFERENCE);
+	let num_curves = bezpath.segments().count();
+	// TODO: optimization opportunity - this for-loop currently compares all intersections with all curve-segments in the subpath collection
+	bezpath.segments().enumerate().for_each(|(i, other)| {
+		intersections_vec.extend(pathseg_self_intersections(other, error, minimum_separation).iter().flat_map(|value| [(i, value.0), (i, value.1)]));
+		bezpath.segments().enumerate().skip(i + 1).for_each(|(j, curve)| {
+			intersections_vec.extend(
+				filtered_all_segment_intersections(curve, 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
+}
+
 #[cfg(test)]
 mod tests {
 	use super::{bezpath_and_segment_intersections, filtered_segment_intersections};

node-graph/gcore/src/vector/algorithms/mod.rs

@@ -1,4 +1,5 @@
 pub mod bezpath_algorithms;
+pub mod centroid;
 mod contants;
 pub mod instance;
 pub mod intersection;