Move bezier-rs into libraries folder and deploy its interactive docs

libraries/bezier-rs/src/consts.rs

@@ -0,0 +1,28 @@
+// Implementation constants
+
+/// Constant used to determine if `f64`s are equivalent.
+pub const MAX_ABSOLUTE_DIFFERENCE: f64 = 1e-3;
+/// A stricter constant used to determine if `f64`s are equivalent.
+pub const STRICT_MAX_ABSOLUTE_DIFFERENCE: f64 = 1e-6;
+/// Number of distances used in search algorithm for `project`.
+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.;
+
+// Method argument defaults
+
+/// Default `t` value used for the `curve_through_points` functions.
+pub const DEFAULT_T_VALUE: f64 = 0.5;
+/// Default LUT step size in `compute_lookup_table` function.
+pub const DEFAULT_LUT_STEP_SIZE: usize = 10;
+/// Default number of subdivisions used in `length` calculation.
+pub const DEFAULT_LENGTH_SUBDIVISIONS: usize = 1000;
+/// Default step size for `reduce` function.
+pub const DEFAULT_REDUCE_STEP_SIZE: f64 = 0.01;
+
+// SVG constants
+pub const SVG_ARG_CUBIC: &str = "C";
+pub const SVG_ARG_LINEAR: &str = "L";
+pub const SVG_ARG_MOVE: &str = "M";
+pub const SVG_ARG_QUADRATIC: &str = "Q";
+pub const SVG_ARG_CLOSED: &str = "Z";

libraries/bezier-rs/src/structs.rs

@@ -0,0 +1,95 @@
+use glam::DVec2;
+use std::fmt::{Debug, Formatter, Result};
+
+/// Struct to represent optional parameters that can be passed to the `project` function.
+#[derive(Copy, Clone)]
+pub struct ProjectionOptions {
+	/// Size of the lookup table for the initial passthrough. The default value is `20`.
+	pub lut_size: usize,
+	/// Difference used between floating point numbers to be considered as equal. The default value is `0.0001`
+	pub convergence_epsilon: f64,
+	/// Controls the number of iterations needed to consider that minimum distance to have converged. The default value is `3`.
+	pub convergence_limit: usize,
+	/// Controls the maximum total number of iterations to be used. The default value is `10`.
+	pub iteration_limit: usize,
+}
+
+impl Default for ProjectionOptions {
+	fn default() -> Self {
+		ProjectionOptions {
+			lut_size: 20,
+			convergence_epsilon: 1e-4,
+			convergence_limit: 3,
+			iteration_limit: 10,
+		}
+	}
+}
+
+/// Struct used to represent the different strategies for generating arc approximations.
+#[derive(Copy, Clone)]
+pub enum ArcStrategy {
+	/// Start with the greedy strategy of maximizing arc approximations and automatically switch to the divide-and-conquer when the greedy approximations no longer fall within the error bound.
+	Automatic,
+	/// Use the greedy strategy to maximize approximated arcs, despite potentially erroneous arcs.
+	FavorLargerArcs,
+	/// Use the divide-and-conquer strategy that prioritizes correctness over maximal arcs.
+	FavorCorrectness,
+}
+
+/// Struct to represent optional parameters that can be passed to the `arcs` function.
+#[derive(Copy, Clone)]
+pub struct ArcsOptions {
+	/// Determines how the approximated arcs are computed.
+	/// When maximizing the arcs, the algorithm may return incorrect arcs when the curve contains any small loops or segments that look like a very thin "U".
+	/// The enum options behave as follows:
+	/// - `Automatic`: Maximize arcs until an erroneous approximation is found. Compute the arcs of the rest of the curve by first splitting on extremas to ensure no more erroneous cases are encountered.
+	/// - `FavorLargerArcs`: Maximize arcs using the original algorithm from the [Approximating a Bezier curve with circular arcs](https://pomax.github.io/bezierinfo/#arcapproximation) section of Pomax's bezier curve primer. Erroneous arcs are possible.
+	/// - `FavorCorrectness`: Prioritize correctness by first spliting the curve by its extremas and determine the arc approximation of each segment instead.
+	///
+	/// The default value is `Automatic`.
+	pub strategy: ArcStrategy,
+	/// The error used for approximating the arc's fit. The default is `0.5`.
+	pub error: f64,
+	/// The maximum number of segment iterations used as attempts for arc approximations. The default is `100`.
+	pub max_iterations: usize,
+}
+
+impl Default for ArcsOptions {
+	fn default() -> Self {
+		ArcsOptions {
+			strategy: ArcStrategy::Automatic,
+			error: 0.5,
+			max_iterations: 100,
+		}
+	}
+}
+
+/// Struct to represent the circular arc approximation used in the `arcs` bezier function.
+#[derive(Copy, Clone, PartialEq)]
+pub struct CircleArc {
+	/// The center point of the circle.
+	pub center: DVec2,
+	/// The radius of the circle.
+	pub radius: f64,
+	/// The start angle of the circle sector in rad.
+	pub start_angle: f64,
+	/// The end angle of the circle sector in rad.
+	pub end_angle: f64,
+}
+
+impl Debug for CircleArc {
+	fn fmt(&self, f: &mut Formatter<'_>) -> Result {
+		write!(f, "Center: {}, radius: {}, start to end angles: {} to {}", self.center, self.radius, self.start_angle, self.end_angle)
+	}
+}
+
+impl Default for CircleArc {
+	fn default() -> Self {
+		CircleArc {
+			center: DVec2::ZERO,
+			radius: 0.,
+			start_angle: 0.,
+			end_angle: 0.,
+		}
+	}
+}

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

@@ -0,0 +1,88 @@
+use super::*;
+use crate::consts::*;
+
+/// Functionality relating to core `Subpath` operations, such as constructors and `iter`.
+impl Subpath {
+	/// Create a new `Subpath` using a list of [ManipulatorGroup]s.
+	/// A `Subpath` with less than 2 [ManipulatorGroup]s may not be closed.
+	pub fn new(manipulator_groups: Vec<ManipulatorGroup>, closed: bool) -> Self {
+		assert!(!closed || manipulator_groups.len() > 1, "A closed Subpath must contain more than 1 ManipulatorGroup.");
+		Self { manipulator_groups, closed }
+	}
+
+	/// Create a `Subpath` consisting of 2 manipulator groups from a `Bezier`.
+	pub fn from_bezier(bezier: Bezier) -> Self {
+		Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: bezier.start(),
+					in_handle: None,
+					out_handle: bezier.handle_start(),
+				},
+				ManipulatorGroup {
+					anchor: bezier.end(),
+					in_handle: bezier.handle_end(),
+					out_handle: None,
+				},
+			],
+			false,
+		)
+	}
+
+	/// Returns true if the `Subpath` contains no [ManipulatorGroup].
+	pub fn is_empty(&self) -> bool {
+		self.manipulator_groups.is_empty()
+	}
+
+	/// Returns the number of [ManipulatorGroup]s contained within the `Subpath`.
+	pub fn len(&self) -> usize {
+		self.manipulator_groups.len()
+	}
+
+	/// Returns an iterator of the [Bezier]s along the `Subpath`.
+	pub fn iter(&self) -> SubpathIter {
+		SubpathIter { sub_path: self, index: 0 }
+	}
+
+	/// Returns an SVG representation of the `Subpath`.
+	pub fn to_svg(&self, options: ToSVGOptions) -> String {
+		if self.is_empty() {
+			return String::new();
+		}
+
+		let curve_start_argument = format!("{SVG_ARG_MOVE}{} {}", self[0].anchor.x, self[0].anchor.y);
+		let mut curve_arguments: Vec<String> = self.iter().map(|bezier| bezier.svg_curve_argument()).collect();
+		if self.closed {
+			curve_arguments.push(String::from(SVG_ARG_CLOSED));
+		}
+
+		let anchor_arguments = options.formatted_anchor_arguments();
+		let anchor_circles = self
+			.manipulator_groups
+			.iter()
+			.map(|point| format!(r#"<circle cx="{}" cy="{}" {}/>"#, point.anchor.x, point.anchor.y, anchor_arguments))
+			.collect::<Vec<String>>();
+
+		let handle_point_arguments = options.formatted_handle_point_arguments();
+		let handle_circles: Vec<String> = self
+			.manipulator_groups
+			.iter()
+			.flat_map(|group| [group.in_handle, group.out_handle])
+			.flatten()
+			.map(|handle| format!(r#"<circle cx="{}" cy="{}" {}/>"#, handle.x, handle.y, handle_point_arguments))
+			.collect();
+
+		let handle_pieces: Vec<String> = self.iter().filter_map(|bezier| bezier.svg_handle_line_argument()).collect();
+
+		format!(
+			r#"<path d="{} {}" {}/><path d="{}" {}/>{}{}"#,
+			curve_start_argument,
+			curve_arguments.join(" "),
+			options.formatted_curve_arguments(),
+			handle_pieces.join(" "),
+			options.formatted_handle_line_arguments(),
+			handle_circles.join(""),
+			anchor_circles.join(""),
+		)
+	}
+}

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

@@ -0,0 +1,95 @@
+use super::*;
+
+/// Functionality relating to looking up properties of the `Subpath` or points along the `Subpath`.
+impl Subpath {
+	/// Return the sum of the approximation of the length of each `Bezier` curve along the `Subpath`.
+	/// - `num_subdivisions` - Number of subdivisions used to approximate the curve. The default value is `1000`.
+	pub fn length(&self, num_subdivisions: Option<usize>) -> f64 {
+		self.iter().fold(0., |accumulator, bezier| accumulator + bezier.length(num_subdivisions))
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use crate::Bezier;
+	use glam::DVec2;
+
+	#[test]
+	fn length_quadratic() {
+		let start = DVec2::new(20., 30.);
+		let middle = DVec2::new(80., 90.);
+		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 bezier1 = Bezier::from_quadratic_dvec2(start, handle1, middle);
+		let bezier2 = Bezier::from_quadratic_dvec2(middle, handle2, end);
+		let bezier3 = Bezier::from_quadratic_dvec2(end, handle3, start);
+
+		let mut subpath = Subpath::new(
+			vec![
+				ManipulatorGroup {
+					anchor: start,
+					in_handle: None,
+					out_handle: Some(handle1),
+				},
+				ManipulatorGroup {
+					anchor: middle,
+					in_handle: None,
+					out_handle: Some(handle2),
+				},
+				ManipulatorGroup {
+					anchor: end,
+					in_handle: None,
+					out_handle: Some(handle3),
+				},
+			],
+			false,
+		);
+		assert_eq!(subpath.length(None), bezier1.length(None) + bezier2.length(None));
+
+		subpath.closed = true;
+		assert_eq!(subpath.length(None), bezier1.length(None) + bezier2.length(None) + bezier3.length(None));
+	}
+
+	#[test]
+	fn length_mixed() {
+		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,
+		);
+		assert_eq!(subpath.length(None), linear_bezier.length(None) + quadratic_bezier.length(None));
+
+		subpath.closed = true;
+		assert_eq!(subpath.length(None), linear_bezier.length(None) + quadratic_bezier.length(None) + cubic_bezier.length(None));
+	}
+}

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

@@ -0,0 +1,68 @@
+mod core;
+mod lookup;
+mod structs;
+pub use structs::*;
+
+use crate::Bezier;
+
+use std::ops::{Index, IndexMut};
+
+/// Structure used to represent a path composed of [Bezier] curves.
+pub struct Subpath {
+	manipulator_groups: Vec<ManipulatorGroup>,
+	closed: bool,
+}
+
+/// Iteration structure for iterating across each curve of a `Subpath`, using an intermediate `Bezier` representation.
+pub struct SubpathIter<'a> {
+	index: usize,
+	sub_path: &'a Subpath,
+}
+
+impl Index<usize> for Subpath {
+	type Output = ManipulatorGroup;
+
+	fn index(&self, index: usize) -> &Self::Output {
+		assert!(index < self.len(), "Index out of bounds in trait Index of SubPath.");
+		&self.manipulator_groups[index]
+	}
+}
+
+impl IndexMut<usize> for Subpath {
+	fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+		assert!(index < self.len(), "Index out of bounds in trait IndexMut of SubPath.");
+		&mut self.manipulator_groups[index]
+	}
+}
+
+impl Iterator for SubpathIter<'_> {
+	type Item = Bezier;
+
+	// Returns the Bezier representation of each `Subpath` segment, defined between a pair of adjacent manipulator points.
+	fn next(&mut self) -> Option<Self::Item> {
+		let len = self.sub_path.len() - 1
+			+ match self.sub_path.closed {
+				true => 1,
+				false => 0,
+			};
+		if self.index >= len {
+			return None;
+		}
+		let start_index = self.index;
+		let end_index = (self.index + 1) % self.sub_path.len();
+		self.index += 1;
+
+		let start = self.sub_path[start_index].anchor;
+		let end = self.sub_path[end_index].anchor;
+		let out_handle = self.sub_path[start_index].out_handle;
+		let in_handle = self.sub_path[end_index].in_handle;
+
+		if let (Some(handle1), Some(handle2)) = (out_handle, in_handle) {
+			Some(Bezier::from_cubic_dvec2(start, handle1, handle2, end))
+		} else if let Some(handle) = out_handle.or(in_handle) {
+			Some(Bezier::from_quadratic_dvec2(start, handle, end))
+		} else {
+			Some(Bezier::from_linear_dvec2(start, end))
+		}
+	}
+}

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

@@ -0,0 +1,83 @@
+use glam::DVec2;
+
+/// Structure used to represent a single anchor with up to two optional associated handles along a `Subpath`
+pub struct ManipulatorGroup {
+	pub anchor: DVec2,
+	pub in_handle: Option<DVec2>,
+	pub out_handle: Option<DVec2>,
+}
+
+/// Structure to represent optional parameters that can be passed to the `into_svg` function.
+pub struct ToSVGOptions {
+	/// Color of the line segments along the `Subpath`. Defaulted to `black`.
+	pub curve_stroke_color: String,
+	/// Width of the line segments along the `Subpath`. Defaulted to `2.`.
+	pub curve_stroke_width: f64,
+	/// Stroke color outlining circles marking anchors on the `Subpath`. Defaulted to `black`.
+	pub anchor_stroke_color: String,
+	/// Stroke width outlining circles marking anchors on the `Subpath`. Defaulted to `2.`.
+	pub anchor_stroke_width: f64,
+	/// Radius of the circles marking anchors on the `Subpath`. Defaulted to `4.`.
+	pub anchor_radius: f64,
+	/// Fill color of the circles marking anchors on the `Subpath`. Defaulted to `white`.
+	pub anchor_fill: String,
+	/// Color of the line segments connecting anchors to handle points. Defaulted to `gray`.
+	pub handle_line_stroke_color: String,
+	/// Width of the line segments connecting anchors to handle points. Defaulted to `1.`.
+	pub handle_line_stroke_width: f64,
+	/// Stroke color outlining circles marking the handles of `Subpath`. Defaulted to `gray`.
+	pub handle_point_stroke_color: String,
+	/// Stroke color outlining circles marking the handles of `Subpath`. Defaulted to `1.5`.
+	pub handle_point_stroke_width: f64,
+	/// Radius of the circles marking the handles of `Subpath`. Defaulted to `3.`.
+	pub handle_point_radius: f64,
+	/// Fill color of the circles marking the handles of `Subpath`. Defaulted to `white`.
+	pub handle_point_fill: String,
+}
+
+impl ToSVGOptions {
+	/// Combine and format curve styling options for an SVG path.
+	pub(crate) fn formatted_curve_arguments(&self) -> String {
+		format!(r#"stroke="{}" stroke-width="{}" fill="none""#, self.curve_stroke_color, self.curve_stroke_width)
+	}
+
+	/// Combine and format anchor styling options an SVG circle.
+	pub(crate) fn formatted_anchor_arguments(&self) -> String {
+		format!(
+			r#"r="{}", stroke="{}" stroke-width="{}" fill="{}""#,
+			self.anchor_radius, self.anchor_stroke_color, self.anchor_stroke_width, self.anchor_fill
+		)
+	}
+
+	/// Combine and format handle point styling options for an SVG circle.
+	pub(crate) fn formatted_handle_point_arguments(&self) -> String {
+		format!(
+			r#"r="{}", stroke="{}" stroke-width="{}" fill="{}""#,
+			self.handle_point_radius, self.handle_point_stroke_color, self.handle_point_stroke_width, self.handle_point_fill
+		)
+	}
+
+	/// Combine and format handle line styling options an SVG path.
+	pub(crate) fn formatted_handle_line_arguments(&self) -> String {
+		format!(r#"stroke="{}" stroke-width="{}" fill="none""#, self.handle_line_stroke_color, self.handle_line_stroke_width)
+	}
+}
+
+impl Default for ToSVGOptions {
+	fn default() -> Self {
+		ToSVGOptions {
+			curve_stroke_color: String::from("black"),
+			curve_stroke_width: 2.,
+			anchor_stroke_color: String::from("black"),
+			anchor_stroke_width: 2.,
+			anchor_radius: 4.,
+			anchor_fill: String::from("white"),
+			handle_line_stroke_color: String::from("gray"),
+			handle_line_stroke_width: 1.,
+			handle_point_stroke_color: String::from("gray"),
+			handle_point_stroke_width: 1.5,
+			handle_point_radius: 3.,
+			handle_point_fill: String::from("white"),
+		}
+	}
+}

libraries/bezier-rs/src/utils.rs

@@ -0,0 +1,333 @@
+use crate::consts::{MAX_ABSOLUTE_DIFFERENCE, STRICT_MAX_ABSOLUTE_DIFFERENCE};
+
+use glam::{BVec2, DMat2, 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.
+/// Relevant derivation and the definitions of a, b, and c can be found in [the projection identity section](https://pomax.github.io/bezierinfo/#abc) of Pomax's bezier curve primer.
+fn compute_abc_through_points(start_point: DVec2, point_on_curve: DVec2, end_point: DVec2, t_to_nth_power: f64, nth_power_of_one_minus_t: f64) -> [DVec2; 3] {
+	let point_c_ratio = nth_power_of_one_minus_t / (t_to_nth_power + nth_power_of_one_minus_t);
+	let c = point_c_ratio * start_point + (1. - point_c_ratio) * end_point;
+	let ab_bc_ratio = (t_to_nth_power + nth_power_of_one_minus_t - 1.).abs() / (t_to_nth_power + nth_power_of_one_minus_t);
+	let a = point_on_curve + (point_on_curve - c) / ab_bc_ratio;
+	[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.
+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;
+	let squared_one_minus_t = one_minus_t * one_minus_t;
+	compute_abc_through_points(start_point, point_on_curve, end_point, t_squared, squared_one_minus_t)
+}
+
+/// Compute `a`, `b`, and `c` for a cubic 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_cubic_through_points(start_point: DVec2, point_on_curve: DVec2, end_point: DVec2, t: f64) -> [DVec2; 3] {
+	let t_cubed = t * t * t;
+	let one_minus_t = 1. - t;
+	let cubed_one_minus_t = one_minus_t * one_minus_t * one_minus_t;
+
+	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) -> (usize, f64) {
+	lut.iter()
+		.enumerate()
+		.map(|(i, p)| (i, point.distance_squared(*p)))
+		.min_by(|x, y| (&(x.1)).partial_cmp(&(y.1)).unwrap())
+		.unwrap()
+}
+
+// TODO: Use an `Option` return type instead of a `Vec`
+/// Find the roots of the linear equation `ax + b`.
+pub fn solve_linear(a: f64, b: f64) -> Vec<f64> {
+	let mut roots = Vec::new();
+	// There exist roots when `a` is not 0
+	if a.abs() > MAX_ABSOLUTE_DIFFERENCE {
+		roots.push(-b / a);
+	}
+	roots
+}
+
+// TODO: Use an `impl Iterator` return type instead of a `Vec`
+/// 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<f64> {
+	let mut roots = Vec::new();
+	if two_times_a != 0. {
+		if discriminant > 0. {
+			let root_discriminant = discriminant.sqrt();
+			roots.push((-b + root_discriminant) / (two_times_a));
+			roots.push((-b - root_discriminant) / (two_times_a));
+		} else if discriminant == 0. {
+			roots.push(-b / (two_times_a));
+		}
+	} else {
+		roots = solve_linear(b, c);
+	}
+	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.)
+	}
+}
+
+// TODO: Use an `impl Iterator` return type instead of a `Vec`
+/// 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 {
+		// When discriminant is 0 (check for approximation because of floating point errors), all roots are real, and 2 are repeated
+		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);
+	} else if discriminant > 0. {
+		// When discriminant > 0, there is 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 {
+		// Otherwise, discriminant < 0 and there are 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;
+		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
+}
+
+// TODO: Use an `impl Iterator` return type instead of a `Vec`
+/// 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<f64> {
+	if a.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
+		if b.abs() <= STRICT_MAX_ABSOLUTE_DIFFERENCE {
+			// 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)
+	}
+}
+
+/// Determine if two rectangles have any overlap. The rectangles are represented by a pair of coordinates that designate the top left and bottom right corners (in a graphical coordinate system).
+pub fn do_rectangles_overlap(rectangle1: [DVec2; 2], rectangle2: [DVec2; 2]) -> bool {
+	let [bottom_left1, top_right1] = rectangle1;
+	let [bottom_left2, top_right2] = rectangle2;
+
+	top_right1.x >= bottom_left2.x && top_right2.x >= bottom_left1.x && top_right2.y >= bottom_left1.y && top_right1.y >= bottom_left2.y
+}
+
+/// Returns the intersection of two lines. The lines are given by a point on the line and its slope (represented by a vector).
+pub fn line_intersection(point1: DVec2, point1_slope_vector: DVec2, point2: DVec2, point2_slope_vector: DVec2) -> DVec2 {
+	assert!(point1_slope_vector.normalize() != point2_slope_vector.normalize());
+
+	// Find the intersection when the first line is vertical
+	if f64_compare(point1_slope_vector.x, 0., MAX_ABSOLUTE_DIFFERENCE) {
+		let m2 = point2_slope_vector.y / point2_slope_vector.x;
+		let b2 = point2.y - m2 * point2.x;
+		DVec2::new(point1.x, point1.x * m2 + b2)
+	}
+	// Find the intersection when the second line is vertical
+	else if f64_compare(point2_slope_vector.x, 0., MAX_ABSOLUTE_DIFFERENCE) {
+		let m1 = point1_slope_vector.y / point1_slope_vector.x;
+		let b1 = point1.y - m1 * point1.x;
+		DVec2::new(point2.x, point2.x * m1 + b1)
+	}
+	// Find the intersection where neither line is vertical
+	else {
+		let m1 = point1_slope_vector.y / point1_slope_vector.x;
+		let b1 = point1.y - m1 * point1.x;
+		let m2 = point2_slope_vector.y / point2_slope_vector.x;
+		let b2 = point2.y - m2 * point2.x;
+		let intersection_x = (b2 - b1) / (m1 - m2);
+		DVec2::new(intersection_x, intersection_x * m1 + b1)
+	}
+}
+
+/// Check if 3 points are collinear.
+pub fn are_points_collinear(p1: DVec2, p2: DVec2, p3: DVec2) -> bool {
+	let matrix = DMat2::from_cols(p1 - p2, p2 - p3);
+	f64_compare(matrix.determinant() / 2., 0., MAX_ABSOLUTE_DIFFERENCE)
+}
+
+/// Compute the center of the circle that passes through all three provided points. The provided points cannot be collinear.
+pub fn compute_circle_center_from_points(p1: DVec2, p2: DVec2, p3: DVec2) -> Option<DVec2> {
+	if are_points_collinear(p1, p2, p3) {
+		return None;
+	}
+
+	let midpoint_a = p1.lerp(p2, 0.5);
+	let midpoint_b = p2.lerp(p3, 0.5);
+	let midpoint_c = p3.lerp(p1, 0.5);
+
+	let tangent_a = (p1 - p2).perp();
+	let tangent_b = (p2 - p3).perp();
+	let tangent_c = (p3 - p1).perp();
+
+	let intersect_a_b = line_intersection(midpoint_a, tangent_a, midpoint_b, tangent_b);
+	let intersect_b_c = line_intersection(midpoint_b, tangent_b, midpoint_c, tangent_c);
+	let intersect_c_a = line_intersection(midpoint_c, tangent_c, midpoint_a, tangent_a);
+
+	Some((intersect_a_b + intersect_b_c + intersect_c_a) / 3.)
+}
+
+/// 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::*;
+	use crate::consts::MAX_ABSOLUTE_DIFFERENCE;
+
+	/// Compare vectors of `f64`s with a provided max absolute value difference.
+	fn f64_compare_vector(vec1: Vec<f64>, vec2: Vec<f64>, max_abs_diff: f64) -> bool {
+		vec1.len() == vec2.len() && vec1.into_iter().zip(vec2.into_iter()).all(|(a, b)| f64_compare(a, b, max_abs_diff))
+	}
+
+	#[test]
+	fn test_solve_linear() {
+		// Line that is on the x-axis
+		assert!(solve_linear(0., 0.).is_empty());
+		// Line that is parallel to but not on the x-axis
+		assert!(solve_linear(0., 1.).is_empty());
+		// Line with a non-zero slope
+		assert!(solve_linear(2., -8.) == vec![4.]);
+	}
+
+	#[test]
+	fn test_solve_cubic() {
+		// discriminant == 0
+		let roots1 = solve_cubic(1., 0., 0., 0.);
+		assert!(roots1 == vec![0.]);
+
+		let roots2 = solve_cubic(1., 3., 0., -4.);
+		assert!(roots2 == vec![1., -2.]);
+
+		// p == 0
+		let roots3 = solve_cubic(1., 0., 0., -1.);
+		assert!(roots3 == vec![1.]);
+
+		// discriminant > 0
+		let roots4 = solve_cubic(1., 3., 0., 2.);
+		assert!(f64_compare_vector(roots4, vec![-3.196], MAX_ABSOLUTE_DIFFERENCE));
+
+		// discriminant < 0
+		let roots5 = solve_cubic(1., 3., 0., -1.);
+		assert!(f64_compare_vector(roots5, vec![0.532, -2.879, -0.653], MAX_ABSOLUTE_DIFFERENCE));
+
+		// quadratic
+		let roots6 = solve_cubic(0., 3., 0., -3.);
+		assert!(roots6 == vec![1., -1.]);
+
+		// linear
+		let roots7 = solve_cubic(0., 0., 1., -1.);
+		assert!(roots7 == vec![1.]);
+	}
+
+	#[test]
+	fn test_do_rectangles_overlap() {
+		// Rectangles overlap
+		assert!(do_rectangles_overlap([DVec2::new(0., 0.), DVec2::new(20., 20.)], [DVec2::new(10., 10.), DVec2::new(30., 20.)]));
+		// Rectangles share a side
+		assert!(do_rectangles_overlap([DVec2::new(0., 0.), DVec2::new(10., 10.)], [DVec2::new(10., 10.), DVec2::new(30., 30.)]));
+		// Rectangle inside the other
+		assert!(do_rectangles_overlap([DVec2::new(0., 0.), DVec2::new(10., 10.)], [DVec2::new(2., 2.), DVec2::new(6., 4.)]));
+		// No overlap, rectangles are beside each other
+		assert!(!do_rectangles_overlap([DVec2::new(0., 0.), DVec2::new(10., 10.)], [DVec2::new(20., 0.), DVec2::new(30., 10.)]));
+		// No overlap, rectangles are above and below each other
+		assert!(!do_rectangles_overlap([DVec2::new(0., 0.), DVec2::new(10., 10.)], [DVec2::new(0., 20.), DVec2::new(20., 30.)]));
+	}
+
+	#[test]
+	fn test_find_intersection() {
+		// y = 2x + 10
+		// y = 5x + 4
+		// intersect at (2, 14)
+
+		let start1 = DVec2::new(0., 10.);
+		let end1 = DVec2::new(0., 4.);
+		let start_direction1 = DVec2::new(1., 2.);
+		let end_direction1 = DVec2::new(1., 5.);
+		assert!(line_intersection(start1, start_direction1, end1, end_direction1) == DVec2::new(2., 14.));
+
+		// y = x
+		// y = -x + 8
+		// intersect at (4, 4)
+
+		let start2 = DVec2::new(0., 0.);
+		let end2 = DVec2::new(8., 0.);
+		let start_direction2 = DVec2::new(1., 1.);
+		let end_direction2 = DVec2::new(1., -1.);
+		assert!(line_intersection(start2, start_direction2, end2, end_direction2) == DVec2::new(4., 4.));
+	}
+
+	#[test]
+	fn test_are_points_collinear() {
+		assert!(are_points_collinear(DVec2::new(2., 4.), DVec2::new(6., 8.), DVec2::new(4., 6.)));
+		assert!(!are_points_collinear(DVec2::new(1., 4.), DVec2::new(6., 8.), DVec2::new(4., 6.)));
+	}
+
+	#[test]
+	fn test_compute_circle_center_from_points() {
+		// 3/4 of unit circle
+		let center1 = compute_circle_center_from_points(DVec2::new(0., 1.), DVec2::new(-1., 0.), DVec2::new(1., 0.));
+		assert_eq!(center1.unwrap(), DVec2::new(0., 0.));
+		// 1/4 of unit circle
+		let center2 = compute_circle_center_from_points(DVec2::new(-1., 0.), DVec2::new(0., 1.), DVec2::new(1., 0.));
+		assert_eq!(center2.unwrap(), DVec2::new(0., 0.));
+	}
+}