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 <keavon@keavon.com>

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

@@ -105,7 +105,20 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 
 	/// Returns an iterator of the [Bezier]s along the `Subpath`.
 	pub fn iter(&self) -> SubpathIter<ManipulatorGroupId> {
-		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<ManipulatorGroupId> {
+		SubpathIter {
+			subpath: self,
+			index: 0,
+			is_always_closed: true,
+		}
 	}
 
 	/// Returns a slice of the [ManipulatorGroup]s in the `Subpath`.

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

@@ -30,6 +30,91 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		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<f64>, minimum_separation: Option<f64>) -> 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<f64>, minimum_separation: Option<f64>) -> Option<DVec2> {
+		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 {

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

@@ -29,6 +29,7 @@ unsafe impl<ManipulatorGroupId: crate::Identifier> dyn_any::StaticType for Subpa
 pub struct SubpathIter<'a, ManipulatorGroupId: crate::Identifier> {
 	index: usize,
 	subpath: &'a Subpath<ManipulatorGroupId>,
+	is_always_closed: bool,
 }
 
 impl<ManipulatorGroupId: crate::Identifier> Index<usize> for Subpath<ManipulatorGroupId> {
@@ -55,8 +56,9 @@ impl<ManipulatorGroupId: crate::Identifier> 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,
 			};

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

@@ -116,7 +116,7 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		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<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		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<f64>, minimum_separation: Option<f64>) -> 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.