Add Poisson-disk sampling node and Bezier-rs 0.4 release (#1586)

* Add Poisson-disk sampling node and Bezier-rs 0.4 release

* Additional optimizations

* More performance optimizations with help from 0Hypercube

* Add comments

Cargo.lock

@@ -580,7 +580,7 @@ checksum = "9d297deb1925b89f2ccc13d7635fa0714f12c87adce1c75356b39ca9b7178567"
 
 [[package]]
 name = "bezier-rs"
-version = "0.3.0"
+version = "0.4.0"
 dependencies = [
  "dyn-any",
  "glam",
@@ -2329,6 +2329,7 @@ dependencies = [
  "node-macro",
  "num-derive",
  "num-traits",
+ "rand 0.8.5",
  "rand_chacha 0.3.1",
  "rustybuzz 0.8.0",
  "serde",

editor/src/messages/portfolio/document/node_graph/node_graph_message_handler/document_node_types.rs

@@ -2706,6 +2706,18 @@ fn static_nodes() -> Vec<DocumentNodeDefinition> {
 			properties: node_properties::sample_points_properties,
 			..Default::default()
 		},
+		DocumentNodeDefinition {
+			name: "Poisson-Disk Points",
+			category: "Vector",
+			implementation: NodeImplementation::proto("graphene_core::vector::PoissonDiskPoints<_>"),
+			inputs: vec![
+				DocumentInputType::value("Vector Data", TaggedValue::VectorData(graphene_core::vector::VectorData::empty()), true),
+				DocumentInputType::value("Separation Disk Diameter", TaggedValue::F32(10.), false),
+			],
+			outputs: vec![DocumentOutputType::new("Vector", FrontendGraphDataType::Subpath)],
+			properties: node_properties::poisson_disk_points_properties,
+			..Default::default()
+		},
 		DocumentNodeDefinition {
 			name: "Splines from Points",
 			category: "Vector",
@@ -2723,7 +2735,7 @@ fn static_nodes() -> Vec<DocumentNodeDefinition> {
 				DocumentInputType::value("Source", TaggedValue::VectorData(graphene_core::vector::VectorData::empty()), true),
 				DocumentInputType::value("Target", TaggedValue::VectorData(graphene_core::vector::VectorData::empty()), true),
 				DocumentInputType::value("Start Index", TaggedValue::U32(0), false),
-				DocumentInputType::value("Time", TaggedValue::F64(0.5), false),
+				DocumentInputType::value("Time", TaggedValue::F32(0.5), false),
 			],
 			outputs: vec![DocumentOutputType::new("Vector", FrontendGraphDataType::Subpath)],
 			manual_composition: Some(concrete!(Footprint)),

editor/src/messages/portfolio/document/node_graph/node_graph_message_handler/node_properties.rs

@@ -251,7 +251,7 @@ fn vec_dvec2_input(document_node: &DocumentNode, node_id: NodeId, index: usize,
 			.filter(|x| !x.is_empty())
 			.map(|x| x.parse::<f64>().ok())
 			.collect::<Option<Vec<_>>>()
-			.map(|numbers| numbers.chunks_exact(2).map(|vals| DVec2::new(vals[0], vals[1])).collect())
+			.map(|numbers| numbers.chunks_exact(2).map(|values| DVec2::new(values[0], values[1])).collect())
 			.map(TaggedValue::VecDVec2)
 	};
 
@@ -2055,6 +2055,19 @@ pub fn sample_points_properties(document_node: &DocumentNode, node_id: NodeId, _
 	]
 }
 
+pub fn poisson_disk_points_properties(document_node: &DocumentNode, node_id: NodeId, _context: &mut NodePropertiesContext) -> Vec<LayoutGroup> {
+	let spacing = number_widget(
+		document_node,
+		node_id,
+		1,
+		"Separation Disk Diameter",
+		NumberInput::default().min(0.01).mode_range().range_min(Some(1.)).range_max(Some(100.)),
+		true,
+	);
+
+	vec![LayoutGroup::Row { widgets: spacing }]
+}
+
 pub fn morph_properties(document_node: &DocumentNode, node_id: NodeId, _context: &mut NodePropertiesContext) -> Vec<LayoutGroup> {
 	let start_index = number_widget(document_node, node_id, 2, "Start Index", NumberInput::default().min(0.), true);
 	let time = number_widget(document_node, node_id, 3, "Time", NumberInput::default().min(0.).max(1.).mode_range(), true);

editor/src/messages/tool/tool_messages/select_tool.rs

@@ -8,7 +8,6 @@ use crate::messages::portfolio::document::overlays::utility_types::OverlayContex
 use crate::messages::portfolio::document::utility_types::document_metadata::LayerNodeIdentifier;
 use crate::messages::portfolio::document::utility_types::misc::{AlignAggregate, AlignAxis, FlipAxis};
 use crate::messages::portfolio::document::utility_types::transformation::Selected;
-use crate::messages::tool;
 use crate::messages::tool::common_functionality::graph_modification_utils::is_layer_fed_by_node_of_name;
 use crate::messages::tool::common_functionality::pivot::Pivot;
 use crate::messages::tool::common_functionality::snapping::{self, SnapCandidatePoint, SnapConstraint, SnapData, SnapManager, SnappedPoint};

libraries/bezier-rs/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "bezier-rs"
-version = "0.3.0"
+version = "0.4.0"
 rust-version = "1.66.0"
 edition = "2021"
 authors = ["Graphite Authors <contact@graphite.rs>"]

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

@@ -127,10 +127,10 @@ impl Bezier {
 	pub fn write_curve_argument(&self, svg: &mut String) -> std::fmt::Result {
 		match self.handles {
 			BezierHandles::Linear => svg.push_str(SVG_ARG_LINEAR),
-			BezierHandles::Quadratic { handle } => write!(svg, "{SVG_ARG_QUADRATIC}{},{}", handle.x, handle.y)?,
-			BezierHandles::Cubic { handle_start, handle_end } => write!(svg, "{SVG_ARG_CUBIC}{},{} {},{}", handle_start.x, handle_start.y, handle_end.x, handle_end.y)?,
+			BezierHandles::Quadratic { handle } => write!(svg, "{SVG_ARG_QUADRATIC}{:.6},{:.6}", handle.x, handle.y)?,
+			BezierHandles::Cubic { handle_start, handle_end } => write!(svg, "{SVG_ARG_CUBIC}{:.6},{:.6} {:.6},{:.6}", handle_start.x, handle_start.y, handle_end.x, handle_end.y)?,
 		}
-		write!(svg, " {},{}", self.end.x, self.end.y)
+		write!(svg, " {:.6},{:.6}", self.end.x, self.end.y)
 	}
 
 	/// Return the string argument used to create the lines connecting handles to endpoints in an SVG `path`
@@ -138,17 +138,17 @@ impl Bezier {
 		match self.handles {
 			BezierHandles::Linear => None,
 			BezierHandles::Quadratic { handle } => {
-				let handle_line = format!("{SVG_ARG_LINEAR}{} {}", handle.x, handle.y);
+				let handle_line = format!("{SVG_ARG_LINEAR}{:.6} {:.6}", handle.x, handle.y);
 				Some(format!(
-					"{SVG_ARG_MOVE}{} {} {handle_line} {SVG_ARG_MOVE}{} {} {handle_line}",
+					"{SVG_ARG_MOVE}{:.6} {:.6} {handle_line} {SVG_ARG_MOVE}{:.6} {:.6} {handle_line}",
 					self.start.x, self.start.y, self.end.x, self.end.y
 				))
 			}
 			BezierHandles::Cubic { handle_start, handle_end } => {
-				let handle_start_line = format!("{SVG_ARG_LINEAR}{} {}", handle_start.x, handle_start.y);
+				let handle_start_line = format!("{SVG_ARG_LINEAR}{:.6} {:.6}", handle_start.x, handle_start.y);
 				let handle_end_line = format!("{SVG_ARG_LINEAR}{} {}", handle_end.x, handle_end.y);
 				Some(format!(
-					"{SVG_ARG_MOVE}{} {} {handle_start_line} {SVG_ARG_MOVE}{} {} {handle_end_line}",
+					"{SVG_ARG_MOVE}{:.6} {:.6} {handle_start_line} {SVG_ARG_MOVE}{:.6} {:.6} {handle_end_line}",
 					self.start.x, self.start.y, self.end.x, self.end.y
 				))
 			}

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

@@ -181,6 +181,15 @@ impl Bezier {
 		[endpoints_min, endpoints_max]
 	}
 
+	/// Return the min and max corners that represent the bounding box enclosing this Bezier's two anchor points and any handles.
+	pub fn bounding_box_of_anchors_and_handles(&self) -> [DVec2; 2] {
+		match self.handles {
+			BezierHandles::Linear => [self.start.min(self.end), self.start.max(self.end)],
+			BezierHandles::Quadratic { handle } => [self.start.min(self.end).min(handle), self.start.max(self.end).max(handle)],
+			BezierHandles::Cubic { handle_start, handle_end } => [self.start.min(self.end).min(handle_start).min(handle_end), self.start.max(self.end).max(handle_start).max(handle_end)],
+		}
+	}
+
 	/// Returns `true` if the bounding box of the bezier is contained entirely within a rectangle defined by its minimum and maximum corners.
 	pub fn is_contained_within(&self, min_corner: DVec2, max_corner: DVec2) -> bool {
 		let [bounding_box_min, bounding_box_max] = self.bounding_box();
@@ -261,7 +270,7 @@ impl Bezier {
 	/// The list of `t`-values returned are filtered such that they fall within the range `[0, 1]`.
 	/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/inflections/solo" title="Inflections Demo"></iframe>
 	pub fn inflections(&self) -> Vec<f64> {
-		self.unrestricted_inflections().into_iter().filter(|&t| t > 0. && t < 1.).collect::<Vec<f64>>()
+		self.unrestricted_inflections().filter(|&t| t > 0. && t < 1.).collect::<Vec<f64>>()
 	}
 
 	/// Implementation of the algorithm to find curve intersections by iterating on bounding boxes.
@@ -343,22 +352,21 @@ impl Bezier {
 			let line_directional_vector = other.end - other.start;
 			let angle = line_directional_vector.angle_between(DVec2::new(0., 1.));
 			let rotation_matrix = DMat2::from_angle(angle);
-			let rotated_bezier = self.apply_transformation(|point| rotation_matrix.mul_vec2(point));
-			let rotated_line = [rotation_matrix.mul_vec2(other.start), rotation_matrix.mul_vec2(other.end)];
+			let rotated_bezier = self.apply_transformation(|point| rotation_matrix * point);
 
 			// Translate the bezier such that the line becomes aligned on top of the x-axis
-			let vertical_distance = rotated_line[0].x;
+			let vertical_distance = (rotation_matrix * other.start).x;
 			let translated_bezier = rotated_bezier.translate(DVec2::new(-vertical_distance, 0.));
 
 			// Compute the roots of the resulting bezier curve
 			let list_intersection_t = translated_bezier.find_tvalues_for_x(0.);
 
-			let min = other.start.min(other.end);
-			let max = other.start.max(other.end);
+			// Calculate line's bounding box
+			let [min_corner, max_corner] = other.bounding_box_of_anchors_and_handles();
 
 			return list_intersection_t
 				// Accept the t value if it is approximately in [0, 1] and if the corresponding coordinates are within the range of the linear line
-				.filter(|&t| utils::dvec2_approximately_in_range(self.unrestricted_parametric_evaluate(t), min, max, MAX_ABSOLUTE_DIFFERENCE).all())
+				.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::<Vec<f64>>();
@@ -369,6 +377,59 @@ impl Bezier {
 		self.intersections_between_subcurves(0. ..1., other, 0. ..1., error).iter().map(|t_values| t_values[0]).collect()
 	}
 
+	/// 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.)
+	/// If this needs to be called frequently with a line of the same rotation angle, consider instead using [`line_test_crossings_prerotated`] and moving this function's setup code into your own logic before the repeated call.
+	pub fn line_test_crossings(&self, point_on_line: DVec2, direction_vector: DVec2) -> impl Iterator<Item = f64> + '_ {
+		// Rotate the bezier and the line by the angle that the line makes with the x axis
+		let angle = direction_vector.angle_between(DVec2::new(0., 1.));
+		let rotation_matrix = DMat2::from_angle(angle);
+		let rotated_bezier = self.apply_transformation(|point| rotation_matrix * point);
+
+		self.line_test_crossings_prerotated(point_on_line, rotation_matrix, rotated_bezier)
+	}
+
+	/// Returns a list of `t` values that correspond to points on this Bezier segment where they intersect with the given infinite line.
+	/// This version of the function is for better performance when calling it frequently without needing to change the rotation between each call.
+	/// If that isn't important, use [`line_test_crossings`] which wraps this and provides an easier interface by taking a line rotation vector.
+	/// Instead, this version requires a rotation matrix for the line's rotation and a version of this Bezier segment that has had its rotation already applied.
+	pub fn line_test_crossings_prerotated(&self, point_on_line: DVec2, rotation_matrix: DMat2, rotated_bezier: Self) -> impl Iterator<Item = f64> + '_ {
+		// Translate the bezier such that the line becomes aligned on top of the x-axis
+		let vertical_distance = (rotation_matrix.x_axis.x * point_on_line.x) + (rotation_matrix.y_axis.x * point_on_line.y);
+		let translated_bezier = rotated_bezier.translate(DVec2::new(-vertical_distance, 0.));
+
+		// Compute the roots of the resulting bezier curve
+		translated_bezier.find_tvalues_for_x(0.)
+	}
+
+	/// Returns a list of `t` values that correspond to points on this Bezier segment where they intersect with the given ray. (`ray_direction` does not need to be normalized.)
+	/// If this needs to be called frequently with a ray of the same rotation angle, consider instead using [`ray_test_crossings_prerotated`] and moving this function's setup code into your own logic before the repeated call.
+	pub fn ray_test_crossings(&self, ray_start: DVec2, ray_direction: DVec2) -> impl Iterator<Item = f64> + '_ {
+		// Rotate the bezier and the line by the angle that the line makes with the x axis
+		let angle = ray_direction.angle_between(DVec2::new(0., 1.));
+		let rotation_matrix = DMat2::from_angle(angle);
+		let rotated_bezier = self.apply_transformation(|point| rotation_matrix * point);
+
+		self.ray_test_crossings_prerotated(ray_start, rotation_matrix, rotated_bezier)
+	}
+
+	/// Returns a list of `t` values that correspond to points on this Bezier segment where they intersect with the given infinite ray.
+	/// This version of the function is for better performance when calling it frequently without needing to change the rotation between each call.
+	/// If that isn't important, use [`ray_test_crossings`] which wraps this and provides an easier interface by taking a ray direction vector.
+	/// Instead, this version requires a rotation matrix for the ray's rotation and a version of this Bezier segment that has had its rotation already applied.
+	pub fn ray_test_crossings_prerotated(&self, ray_start: DVec2, rotation_matrix: DMat2, rotated_bezier: Self) -> impl Iterator<Item = f64> + '_ {
+		// Intersection t-values include those beyond the [0-1] range where the segment's ends extend through the X-axis
+		let intersection_t_values_on_rotated_bezier = self.line_test_crossings_prerotated(ray_start, rotation_matrix, rotated_bezier);
+
+		intersection_t_values_on_rotated_bezier
+			// Accept the t value if it is approximately in [0, 1] and if the corresponding coordinates are within the range of the linear line
+			.filter(move |&t| {
+				let point = self.unrestricted_parametric_evaluate(t);
+				// Ensure the returned value is within the correct range
+				let in_bounds = point.cmpge(ray_start) | utils::dvec2_compare(point, ray_start, MAX_ABSOLUTE_DIFFERENCE);
+				in_bounds.x && in_bounds.y
+			})
+	}
+
 	/// Helper function to compute intersections between lists of subcurves.
 	/// This function uses the algorithm implemented in `intersections_between_subcurves`.
 	fn intersections_between_vectors_of_curves(subcurves1: &[(Bezier, Range<f64>)], subcurves2: &[(Bezier, Range<f64>)], error: f64) -> Vec<[f64; 2]> {

libraries/bezier-rs/src/lib.rs

@@ -4,6 +4,7 @@ pub(crate) mod compare;
 
 mod bezier;
 mod consts;
+mod poisson_disk;
 mod subpath;
 mod symmetrical_basis;
 mod utils;

libraries/bezier-rs/src/poisson_disk.rs

@@ -0,0 +1,364 @@
+use core::f64;
+use glam::DVec2;
+
+const DEEPEST_SUBDIVISION_LEVEL_BEFORE_DISCARDING: usize = 8;
+
+/// Fast (O(n) with respect to time and memory) algorithm for generating a maximal set of points using Poisson-disk sampling.
+/// Based on the paper:
+/// "Poisson Disk Point Sets by Hierarchical Dart Throwing"
+/// <https://scholarsarchive.byu.edu/facpub/237/>
+pub fn poisson_disk_sample(
+	width: f64,
+	height: f64,
+	diameter: f64,
+	point_in_shape_checker: impl Fn(DVec2) -> bool,
+	square_edges_intersect_shape_checker: impl Fn(DVec2, f64) -> bool,
+	rng: impl FnMut() -> f64,
+) -> Vec<DVec2> {
+	let mut rng = rng;
+	let diameter_squared = diameter.powi(2);
+
+	// Initialize a place to store the generated points within a spatial acceleration structure
+	let mut points_grid = AccelerationGrid::new(width, height, diameter);
+
+	// Pick a grid size for the base-level domain that's as large as possible, while also:
+	// - Dividing into an integer number of cells across the dartboard domain, to avoid wastefully throwing darts beyond the width and height of the dartboard domain
+	// - Being fully covered by the radius around a dart thrown anywhere in its area, where the worst-case is a corner which has a distance of sqrt(2) to the opposite corner
+	let greater_dimension = width.max(height);
+	let base_level_grid_size = greater_dimension / (greater_dimension * std::f64::consts::SQRT_2 / (diameter / 2.)).ceil();
+
+	// Initialize the problem by including all base-level squares in the active list since they're all part of the yet-to-be-targetted dartboard domain
+	let base_level = ActiveListLevel::new_filled(base_level_grid_size, width, height, &point_in_shape_checker, &square_edges_intersect_shape_checker);
+	// In the future, if necessary, this could be turned into a fixed-length array with worst-case length `f64::MANTISSA_DIGITS`
+	let mut active_list_levels = vec![base_level];
+
+	// Loop until all active squares have been processed, meaning all of the dartboard domain has been checked
+	while active_list_levels.iter().any(|active_list| active_list.not_empty()) {
+		// Randomly pick a square in the dartboard domain, with probability proportional to its area
+		let (active_square_level, active_square_index_in_level) = target_active_square(&active_list_levels, &mut rng);
+
+		// The level contains the list of all active squares at this target square's subdivision depth
+		let level = &mut active_list_levels[active_square_level];
+
+		// Take the targetted active square out of the list and get its size
+		let active_square = level.take_square(active_square_index_in_level);
+		let active_square_size = level.square_size();
+
+		// Skip this target square if it's within range of any current points, since more nearby points could have been added after this square was included in the active list
+		if !square_not_covered_by_poisson_points(active_square.top_left_corner(), active_square_size / 2., diameter_squared, &points_grid) {
+			continue;
+		}
+
+		// Throw a dart by picking a random point within this target square
+		let point = {
+			let active_top_left_corner = active_square.top_left_corner();
+			let x = active_top_left_corner.x + rng() * active_square_size;
+			let y = active_top_left_corner.y + rng() * active_square_size;
+			(x, y).into()
+		};
+
+		// If the dart hit a valid spot, save that point (we're now permanently done with this target square's region)
+		if point_not_covered_by_poisson_points(point, diameter_squared, &points_grid) {
+			// Silently reject the point if it lies outside the shape
+			if active_square.fully_in_shape() || point_in_shape_checker(point) {
+				points_grid.insert(point);
+			}
+		}
+		// Otherwise, subdivide this target square and add valid sub-squares back to the active list for later targetting
+		else {
+			// Discard any targetable domain smaller than this limited number of subdivision levels since it's too small to matter
+			let next_level_deeper_level = active_square_level + 1;
+			if next_level_deeper_level > DEEPEST_SUBDIVISION_LEVEL_BEFORE_DISCARDING {
+				continue;
+			}
+
+			// If necessary for the following step, add another layer of depth to store squares at the next subdivision level
+			if active_list_levels.len() <= next_level_deeper_level {
+				active_list_levels.push(ActiveListLevel::new(active_square_size / 2.))
+			}
+
+			// Get the list of active squares at the level of depth beneath this target square's level
+			let next_level_deeper = &mut active_list_levels[next_level_deeper_level];
+
+			// Subdivide this target square into four sub-squares; running out of numerical precision will make this terminate at very small scales
+			let subdivided_size = active_square_size / 2.;
+			let active_top_left_corner = active_square.top_left_corner();
+			let subdivided = [
+				active_top_left_corner + DVec2::new(0., 0.),
+				active_top_left_corner + DVec2::new(subdivided_size, 0.),
+				active_top_left_corner + DVec2::new(0., subdivided_size),
+				active_top_left_corner + DVec2::new(subdivided_size, subdivided_size),
+			];
+
+			// Add the sub-squares which aren't within the radius of a nearby point to the sub-level's active list
+			let half_subdivided_size = subdivided_size / 2.;
+			let new_sub_squares = subdivided.into_iter().filter_map(|sub_square| {
+				// Any sub-squares within the radius of a nearby point are filtered out
+				if !square_not_covered_by_poisson_points(sub_square, half_subdivided_size, diameter_squared, &points_grid) {
+					return None;
+				}
+
+				// Fully inside the shape
+				if active_square.fully_in_shape() {
+					Some(ActiveSquare::new(sub_square, true))
+				}
+				// Intersecting the shape's border
+				else {
+					// The sub-square is fully inside the shape if its top-left corner is inside and its edges don't intersect the shape border
+					let sub_square_fully_inside_shape = !square_edges_intersect_shape_checker(sub_square, subdivided_size) && point_in_shape_checker(sub_square);
+					Some(ActiveSquare::new(sub_square, sub_square_fully_inside_shape))
+				}
+			});
+			next_level_deeper.add_squares(new_sub_squares);
+		}
+	}
+
+	points_grid.final_points()
+}
+
+/// Randomly pick a square in the dartboard domain, with probability proportional to its area.
+/// Returns a tuple with the subdivision level depth and the square index at that depth.
+fn target_active_square(active_list_levels: &[ActiveListLevel], rng: &mut impl FnMut() -> f64) -> (usize, usize) {
+	let active_squares_total_area: f64 = active_list_levels.iter().map(|active_list| active_list.total_area()).sum();
+	let mut index_into_area = rng() * active_squares_total_area;
+
+	for (level, active_list_level) in active_list_levels.iter().enumerate() {
+		let subtracted = index_into_area - active_list_level.total_area();
+		if subtracted > 0. {
+			index_into_area = subtracted;
+			continue;
+		}
+
+		let active_square_index_in_level = (index_into_area / active_list_levels[level].square_area()).floor() as usize;
+		return (level, active_square_index_in_level);
+	}
+
+	panic!("index_into_area couldn't be be mapped to a square in any level of the active lists");
+}
+
+fn point_not_covered_by_poisson_points(point: DVec2, diameter_squared: f64, points_grid: &AccelerationGrid) -> bool {
+	points_grid.nearby_points(point).all(|nearby_point| {
+		let x_separation = nearby_point.x - point.x;
+		let y_separation = nearby_point.y - point.y;
+
+		x_separation.powi(2) + y_separation.powi(2) > diameter_squared
+	})
+}
+
+fn square_not_covered_by_poisson_points(point: DVec2, half_square_size: f64, diameter_squared: f64, points_grid: &AccelerationGrid) -> bool {
+	let square_center_x = point.x + half_square_size;
+	let square_center_y = point.y + half_square_size;
+
+	points_grid.nearby_points(point).all(|nearby_point| {
+		let x_distance = (square_center_x - nearby_point.x).abs() + half_square_size;
+		let y_distance = (square_center_y - nearby_point.y).abs() + half_square_size;
+
+		x_distance.powi(2) + y_distance.powi(2) > diameter_squared
+	})
+}
+
+#[inline(always)]
+fn cartesian_product<A, B>(a: A, b: B) -> impl Iterator<Item = (A::Item, B::Item)>
+where
+	A: Iterator + Clone,
+	B: Iterator + Clone,
+	A::Item: Clone,
+	B::Item: Clone,
+{
+	a.flat_map(move |i| (b.clone().map(move |j| (i.clone(), j))))
+}
+
+/// A square (represented by its top left corner position and width/height of `square_size`) that is currently a candidate for targetting by the dart throwing process.
+/// The positive sign bit encodes if the square is contained entirely within the masking shape, or negative if it's outside or intersects the shape path.
+pub struct ActiveSquare(DVec2);
+
+impl ActiveSquare {
+	pub fn new(top_left_corner: DVec2, fully_in_shape: bool) -> Self {
+		Self(if fully_in_shape { top_left_corner } else { -top_left_corner })
+	}
+
+	pub fn top_left_corner(&self) -> DVec2 {
+		self.0.abs()
+	}
+
+	pub fn fully_in_shape(&self) -> bool {
+		self.0.x.is_sign_positive()
+	}
+}
+
+pub struct ActiveListLevel {
+	/// List of all subdivided squares of the same size that are currently candidates for targetting by the dart throwing process
+	active_squares: Vec<ActiveSquare>,
+	/// Width and height of the squares in this level of subdivision
+	square_size: f64,
+	/// Current sum of the area in all active squares in this subdivision level
+	total_area: f64,
+}
+
+impl ActiveListLevel {
+	#[inline(always)]
+	pub fn new(square_size: f64) -> Self {
+		Self {
+			active_squares: Vec::new(),
+			square_size,
+			total_area: 0.,
+		}
+	}
+
+	#[inline(always)]
+	pub fn new_filled(square_size: f64, width: f64, height: f64, point_in_shape_checker: impl Fn(DVec2) -> bool, square_edges_intersect_shape_checker: impl Fn(DVec2, f64) -> bool) -> Self {
+		// These should divide evenly but rounding is to protect against small numerical imprecision errors
+		let x_squares = (width / square_size).round() as usize;
+		let y_squares = (height / square_size).round() as usize;
+
+		// Populate each square with its top-left corner coordinate
+		let active_squares: Vec<_> = cartesian_product(0..x_squares, 0..y_squares)
+			.filter_map(|(x, y)| {
+				let corner = (x as f64 * square_size, y as f64 * square_size).into();
+
+				let point_in_shape = point_in_shape_checker(corner);
+				let square_edges_intersect_shape = square_edges_intersect_shape_checker(corner, square_size);
+				let square_not_outside_shape = point_in_shape || square_edges_intersect_shape;
+				let square_in_shape = point_in_shape && !square_edges_intersect_shape;
+
+				square_not_outside_shape.then_some(ActiveSquare::new(corner, square_in_shape))
+			})
+			.collect();
+
+		// Sum every square's area to get the total
+		let total_area = square_size.powi(2) * active_squares.len() as f64;
+
+		Self {
+			active_squares,
+			square_size,
+			total_area,
+		}
+	}
+
+	#[must_use]
+	#[inline(always)]
+	pub fn take_square(&mut self, active_square_index: usize) -> ActiveSquare {
+		let targetted_square = self.active_squares.swap_remove(active_square_index);
+		self.total_area = self.square_size.powi(2) * self.active_squares.len() as f64;
+		targetted_square
+	}
+
+	#[inline(always)]
+	pub fn add_squares(&mut self, new_squares: impl Iterator<Item = ActiveSquare>) {
+		for new_square in new_squares {
+			self.active_squares.push(new_square);
+		}
+		self.total_area = self.square_size.powi(2) * self.active_squares.len() as f64;
+	}
+
+	#[inline(always)]
+	pub fn square_size(&self) -> f64 {
+		self.square_size
+	}
+
+	#[inline(always)]
+	pub fn square_area(&self) -> f64 {
+		self.square_size.powi(2)
+	}
+
+	#[inline(always)]
+	pub fn total_area(&self) -> f64 {
+		self.total_area
+	}
+
+	#[inline(always)]
+	pub fn not_empty(&self) -> bool {
+		!self.active_squares.is_empty()
+	}
+}
+
+#[derive(Clone, Default)]
+pub struct PointsList {
+	// The worst-case number of points in a 3x3 grid is 16 (one at each intersection of the four gridlines per axis)
+	storage_slots: [DVec2; 16],
+	length: usize,
+}
+
+impl PointsList {
+	#[inline(always)]
+	pub fn push(&mut self, point: DVec2) {
+		self.storage_slots[self.length] = point;
+		self.length += 1;
+	}
+
+	#[inline(always)]
+	pub fn list_cell_and_neighbors(&self) -> impl Iterator<Item = DVec2> {
+		// The negative bit is used to store whether a point belongs to a neighboring cell
+		self.storage_slots.into_iter().take(self.length).map(|point| (point.x.abs(), point.y.abs()).into())
+	}
+
+	#[inline(always)]
+	pub fn list_cell(&self) -> impl Iterator<Item = DVec2> {
+		// The negative bit is used to store whether a point belongs to a neighboring cell
+		self.storage_slots
+			.into_iter()
+			.take(self.length)
+			.filter(|point| point.x.is_sign_positive() && point.y.is_sign_positive())
+	}
+}
+
+pub struct AccelerationGrid {
+	size: f64,
+	dimension_x: usize,
+	dimension_y: usize,
+	cells: Vec<PointsList>,
+}
+
+impl AccelerationGrid {
+	#[inline(always)]
+	pub fn new(width: f64, height: f64, size: f64) -> Self {
+		let dimension_x = (width / size).ceil() as usize + 1;
+		let dimension_y = (height / size).ceil() as usize + 1;
+
+		Self {
+			size,
+			dimension_x,
+			dimension_y,
+			cells: vec![PointsList::default(); dimension_x * dimension_y],
+		}
+	}
+
+	#[inline(always)]
+	pub fn insert(&mut self, point: DVec2) {
+		let x = (point.x / self.size).floor() as usize;
+		let y = (point.y / self.size).floor() as usize;
+
+		// Insert this point at this cell and the surrounding cells in a 3x3 patch
+		for (x_offset, y_offset) in cartesian_product((-1)..=1, (-1)..=1) {
+			// Avoid going negative
+			let (x, y) = (x as isize + x_offset, y as isize + y_offset);
+			if x < 0 || y < 0 {
+				continue;
+			}
+			// Avoid going beyond the width or height
+			let (x, y) = (x as usize, y as usize);
+			if x > self.dimension_x - 1 || y > self.dimension_y - 1 {
+				continue;
+			}
+
+			// Get the cell corresponding to the (x, y) index
+			let cell = &mut self.cells[y * self.dimension_x + x];
+
+			// Store the given point in this grid cell, and use the negative bit to indicate if this belongs to a neighboring cell
+			cell.push(if x_offset == 0 && y_offset == 0 { point } else { -point });
+		}
+	}
+
+	#[inline(always)]
+	pub fn nearby_points(&self, point: DVec2) -> impl Iterator<Item = DVec2> {
+		let x = (point.x / self.size).floor() as usize;
+		let y = (point.y / self.size).floor() as usize;
+
+		self.cells[y * self.dimension_x + x].list_cell_and_neighbors()
+	}
+
+	#[inline(always)]
+	pub fn final_points(&self) -> Vec<DVec2> {
+		self.cells.iter().flat_map(|cell| cell.list_cell()).collect()
+	}
+}

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

@@ -150,7 +150,7 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 			return Ok(());
 		}
 		let start = transform.transform_point2(self[0].anchor);
-		write!(svg, "{SVG_ARG_MOVE}{},{}", start.x, start.y)?;
+		write!(svg, "{SVG_ARG_MOVE}{:.6},{:.6}", start.x, start.y)?;
 		for bezier in self.iter() {
 			bezier.apply_transformation(|pos| transform.transform_point2(pos)).write_curve_argument(svg)?;
 			svg.push(' ');

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

@@ -3,7 +3,7 @@ use crate::consts::MAX_ABSOLUTE_DIFFERENCE;
 use crate::utils::{compute_circular_subpath_details, line_intersection, SubpathTValue};
 use crate::TValue;
 
-use glam::{DMat2, DVec2};
+use glam::{DAffine2, DMat2, DVec2};
 use std::f64::consts::PI;
 
 impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
@@ -23,6 +23,10 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 	/// - `error`: an optional f64 value to provide an error bound
 	/// - `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.
+	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-linear/solo" title="Intersection Demo"></iframe>
+	///
+	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-quadratic/solo" title="Intersection Demo"></iframe>
+	///
 	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-cubic/solo" title="Intersection Demo"></iframe>
 	pub fn intersections(&self, other: &Bezier, error: Option<f64>, minimum_separation: Option<f64>) -> Vec<(usize, f64)> {
 		self.iter()
@@ -35,20 +39,73 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 	/// This function expects the following:
 	/// - other: a [Bezier] curve to check intersections against
 	/// - error: an optional f64 value to provide an error bound
-	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-cubic/solo" title="Intersection Demo"></iframe>
 	pub fn subpath_intersections(&self, other: &Subpath<ManipulatorGroupId>, error: Option<f64>, minimum_separation: Option<f64>) -> Vec<(usize, f64)> {
 		let mut intersection_t_values: Vec<(usize, f64)> = other.iter().flat_map(|bezier| self.intersections(&bezier, error, minimum_separation)).collect();
 		intersection_t_values.sort_by(|a, b| a.partial_cmp(b).unwrap());
 		intersection_t_values
 	}
 
+	/// Returns how many times a given ray intersects with this subpath. (`ray_direction` does not need to be normalized.)
+	/// If this needs to be called frequently with a ray of the same rotation angle, consider instead using [`ray_test_crossings_count_prerotated`].
+	pub fn ray_test_crossings_count(&self, ray_start: DVec2, ray_direction: DVec2) -> usize {
+		self.iter().map(|bezier| bezier.ray_test_crossings(ray_start, ray_direction).count()).sum()
+	}
+
+	/// Returns how many times a given ray intersects with this subpath. (`ray_direction` does not need to be normalized.)
+	/// This version of the function is for better performance when calling it frequently without needing to change the rotation between each call.
+	/// If that isn't important, use [`ray_test_crossings_count`] which provides an easier interface by taking a ray direction vector.
+	/// Instead, this version requires a rotation matrix for the ray's rotation and a prerotated version of this subpath that has had its rotation applied.
+	pub fn ray_test_crossings_count_prerotated(&self, ray_start: DVec2, rotation_matrix: DMat2, rotated_subpath: &Self) -> usize {
+		self.iter()
+			.zip(rotated_subpath.iter())
+			.map(|(bezier, rotated_bezier)| bezier.ray_test_crossings_prerotated(ray_start, rotation_matrix, rotated_bezier).count())
+			.sum()
+	}
+
+	/// Returns true if the given point is inside this subpath. Open paths are NOT automatically closed so you'll need to call `set_closed(true)` before calling this.
+	/// Self-intersecting subpaths use the `evenodd` fill rule for checking in/outside-ness: <https://developer.mozilla.org/en-US/docs/Web/SVG/Attribute/fill-rule>.
+	/// If this needs to be called frequently, consider instead using [`point_inside_prerotated`] and moving this function's setup code into your own logic before the repeated call.
+	pub fn point_inside(&self, point: DVec2) -> bool {
+		// The directions use prime numbers to reduce the likelihood of running across two anchor points simultaneously
+		const SIN_13DEG: f64 = 0.22495105434;
+		const COS_13DEG: f64 = 0.97437006478;
+		const DIRECTION1: DVec2 = DVec2::new(SIN_13DEG, COS_13DEG);
+		const DIRECTION2: DVec2 = DVec2::new(-COS_13DEG, -SIN_13DEG);
+
+		// We (inefficiently) check for odd crossings in two directions and make sure they agree to reduce how often anchor points cause a double-increment
+		let test1 = self.ray_test_crossings_count(point, DIRECTION1) % 2 == 1;
+		let test2 = self.ray_test_crossings_count(point, DIRECTION2) % 2 == 1;
+
+		test1 && test2
+	}
+
+	/// Returns true if the given point is inside this subpath. Open paths are NOT automatically closed so you'll need to call `set_closed(true)` before calling this.
+	/// Self-intersecting subpaths use the `evenodd` fill rule for checking in/outside-ness: <https://developer.mozilla.org/en-US/docs/Web/SVG/Attribute/fill-rule>.
+	/// This version of the function is for better performance when calling it frequently because it lets the caller precompute the rotations once instead of every call.
+	/// If that isn't important, use [`point_inside`] which provides an easier interface.
+	/// Instead, this version requires a pair of rotation matrices for the ray's rotation and a pair of prerotated versions of this subpath.
+	/// They should face in different directions that are unlikely to align in the real world. Consider using the following rotations:
+	/// ```rs
+	/// const SIN_13DEG: f64 = 0.22495105434;
+	/// const COS_13DEG: f64 = 0.97437006478;
+	/// const DIRECTION1: DVec2 = DVec2::new(SIN_13DEG, COS_13DEG);
+	/// const DIRECTION2: DVec2 = DVec2::new(-COS_13DEG, -SIN_13DEG);
+	/// ```
+	pub fn point_inside_prerotated(&self, point: DVec2, rotation_matrix1: DMat2, rotation_matrix2: DMat2, rotated_subpath1: &Self, rotated_subpath2: &Self) -> bool {
+		// We (inefficiently) check for odd crossings in two directions and make sure they agree to reduce how often anchor points cause a double-increment
+		let test1 = self.ray_test_crossings_count_prerotated(point, rotation_matrix1, rotated_subpath1) % 2 == 1;
+		let test2 = self.ray_test_crossings_count_prerotated(point, rotation_matrix2, rotated_subpath2) % 2 == 1;
+
+		test1 && test2
+	}
+
 	/// Returns a list of `t` values that correspond to the self intersection points of the subpath. 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.
 	/// - `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.
-	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/self-intersect/solo" title="Self-Intersection Demo"></iframe>
+	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-self/solo" title="Self-Intersection Demo"></iframe>
 	pub fn 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);
@@ -68,6 +125,79 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		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.
+	/// Expects the following:
+	/// - `corner1`: any corner of the axis-aligned box to intersect with
+	/// - `corner2`: the corner opposite to `corner1`
+	/// - `error`: an optional f64 value to provide an error bound
+	/// - `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.
+	/// <iframe frameBorder="0" width="100%" height="375px" src="https://graphite.rs/libraries/bezier-rs#subpath/intersect-rectangle/solo" title="Intersection Demo"></iframe>
+	pub fn rectangle_intersections(&self, corner1: DVec2, corner2: DVec2, error: Option<f64>, minimum_separation: Option<f64>) -> Vec<(usize, f64)> {
+		[
+			Bezier::from_linear_coordinates(corner1.x, corner1.y, corner2.x, corner1.y),
+			Bezier::from_linear_coordinates(corner2.x, corner1.y, corner2.x, corner2.y),
+			Bezier::from_linear_coordinates(corner2.x, corner2.y, corner1.x, corner2.y),
+			Bezier::from_linear_coordinates(corner1.x, corner2.y, corner1.x, corner1.y),
+		]
+		.iter()
+		.flat_map(|bezier| self.intersections(bezier, error, minimum_separation))
+		.collect()
+	}
+
+	/// Checks if any intersections exist between this subpath and the four edges of the rectangle defined by the top-left `corner1` and bottom-right `corner2`.
+	/// This is faster than calling [`rectangle_intersections`]`.len()` because it short-circuits as soon as an intersection is found.
+	pub fn rectangle_intersections_exist(&self, corner1: DVec2, corner2: DVec2) -> bool {
+		let rotate_by_90deg = |point| DMat2::from_angle(std::f64::consts::FRAC_PI_2) * point;
+
+		for bezier in self.iter() {
+			// Check that the two bounding boxes don't intersect, since we can avoid doing intersection's cubic root finding in that case
+			let [bezier_corner1, bezier_corner2] = bezier.bounding_box_of_anchors_and_handles();
+			if !(((corner1.x <= bezier_corner1.x) && (bezier_corner1.x <= corner2.x) || (corner1.x <= bezier_corner2.x) && (bezier_corner2.x <= corner2.x))
+				&& corner1.y <= bezier_corner2.y
+				&& corner2.y >= bezier_corner1.y
+				|| ((corner1.y <= bezier_corner1.y) && (bezier_corner1.y <= corner2.y) || (corner1.y <= bezier_corner2.y) && (bezier_corner2.y <= corner2.y))
+					&& corner1.x <= bezier_corner2.x
+					&& corner2.x >= bezier_corner1.x)
+			{
+				continue;
+			}
+
+			// Original rotation axis
+			if bezier.line_test_crossings_prerotated(corner1, DMat2::IDENTITY, bezier).any(|intersection_point| {
+				let (_, y) = bezier.unrestricted_parametric_evaluate(intersection_point).into();
+				y >= corner1.y && y <= corner2.y
+			}) {
+				return true;
+			}
+			if bezier.line_test_crossings_prerotated(corner2, DMat2::IDENTITY, bezier).any(|intersection_point| {
+				let (_, y) = bezier.unrestricted_parametric_evaluate(intersection_point).into();
+				y >= corner1.y && y <= corner2.y
+			}) {
+				return true;
+			}
+
+			// Perpendicular to original rotation axis
+			let rotated_bezier = bezier.apply_transformation(rotate_by_90deg);
+			if bezier.line_test_crossings_prerotated(corner1, DMat2::IDENTITY, rotated_bezier).any(|intersection_point| {
+				let (x, _) = bezier.unrestricted_parametric_evaluate(intersection_point).into();
+				x >= corner1.x && x <= corner2.x
+			}) {
+				return true;
+			}
+			if bezier.line_test_crossings_prerotated(corner2, DMat2::IDENTITY, rotated_bezier).any(|intersection_point| {
+				let (x, _) = bezier.unrestricted_parametric_evaluate(intersection_point).into();
+				x >= corner1.x && x <= corner2.x
+			}) {
+				return true;
+			}
+		}
+
+		false
+	}
+
 	/// Returns a normalized unit vector representing the tangent on the subpath based on the parametric `t`-value provided.
 	/// <iframe frameBorder="0" width="100%" height="350px" src="https://graphite.rs/libraries/bezier-rs#subpath/tangent/solo" title="Tangent Demo"></iframe>
 	pub fn tangent(&self, t: SubpathTValue) -> DVec2 {
@@ -137,6 +267,55 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		self.iter().map(|bezier| bezier.winding(target_point)).sum::<i32>() != 0
 	}
 
+	/// Randomly places points across the filled surface of this subpath (which is assumed to be closed).
+	/// The `separation_disk_diameter` determines the minimum distance between all points from one another.
+	/// Conceptually, this works by "throwing a dart" at the subpath's bounding box and keeping the dart only if:
+	/// - It's inside the shape
+	/// - It's not closer than `separation_disk_diameter` to any other point from a previous accepted dart throw
+	/// This repeats until accepted darts fill all possible areas between one another.
+	///
+	/// While the conceptual process described above asymptotically slows down and is never guaranteed to produce a maximal set in finite time,
+	/// this is implemented with an algorithm that produces a maximal set in O(n) time. The slowest part is actually checking if points are inside the subpath shape.
+	pub fn poisson_disk_points(&self, separation_disk_diameter: f64, rng: impl FnMut() -> f64) -> Vec<DVec2> {
+		let Some(bounding_box) = self.bounding_box() else { return Vec::new() };
+		let (offset_x, offset_y) = bounding_box[0].into();
+		let (width, height) = (bounding_box[1] - bounding_box[0]).into();
+
+		// TODO: Optimize the following code and make it more robust
+
+		let mut shape = self.clone();
+		shape.set_closed(true);
+		shape.apply_transform(DAffine2::from_translation((-offset_x, -offset_y).into()));
+
+		const SIN_13DEG: f64 = 0.22495105434;
+		const COS_13DEG: f64 = 0.97437006478;
+		let rotated_subpath = |ray_direction: DVec2| {
+			// Rotate the bezier and the line by the angle that the line makes with the x axis
+			let angle = ray_direction.angle_between(DVec2::new(0., 1.));
+			let rotation_matrix = DMat2::from_angle(angle);
+
+			let mut prerotated = shape.clone();
+			prerotated.apply_transform(DAffine2::from_angle(angle));
+			(rotation_matrix, prerotated)
+		};
+		// The directions use prime numbers to reduce the likelihood of running across two anchor points simultaneously
+		let (matrix1, prerotated1) = rotated_subpath(DVec2::new(SIN_13DEG, COS_13DEG));
+		let (matrix2, prerotated2) = rotated_subpath(DVec2::new(-COS_13DEG, -SIN_13DEG));
+		let point_in_shape_checker = |point: DVec2| shape.point_inside_prerotated(point, matrix1, matrix2, &prerotated1, &prerotated2);
+
+		let square_edges_intersect_shape_checker = |corner1: DVec2, size: f64| {
+			let corner2 = corner1 + DVec2::splat(size);
+			self.rectangle_intersections_exist(corner1, corner2)
+		};
+
+		let mut points = crate::poisson_disk::poisson_disk_sample(width, height, separation_disk_diameter, point_in_shape_checker, square_edges_intersect_shape_checker, rng);
+		for point in &mut points {
+			point.x += offset_x;
+			point.y += offset_y;
+		}
+		points
+	}
+
 	/// Returns the manipulator point that is needed for a miter join if it is possible.
 	/// - `miter_limit`: Defines a limit for the ratio between the miter length and the stroke width.
 	/// Alternatively, this can be interpreted as limiting the angle that the miter can form.
@@ -222,7 +401,7 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
 		compute_circular_subpath_details(left, arc_point, right, center, None)
 	}
 
-	/// Returns the two manipulator groups that create a sqaure cap between the end of `self` and the beginning of `other`.
+	/// Returns the two manipulator groups that create a square cap between the end of `self` and the beginning of `other`.
 	pub(crate) fn square_cap(&self, other: &Subpath<ManipulatorGroupId>) -> [ManipulatorGroup<ManipulatorGroupId>; 2] {
 		let left = self.manipulator_groups[self.len() - 1].anchor;
 		let right = other.manipulator_groups[0].anchor;

libraries/bezier-rs/src/utils.rs

@@ -125,9 +125,9 @@ pub fn solve_quadratic(discriminant: f64, two_times_a: f64, b: f64, c: f64) -> [
 /// Compute the cube root of a number.
 fn cube_root(f: f64) -> f64 {
 	if f < 0. {
-		-(-f).powf(1. / 3.)
+		-(-f).cbrt()
 	} else {
-		f.powf(1. / 3.)
+		f.cbrt()
 	}
 }
 
@@ -149,14 +149,14 @@ pub fn solve_reformatted_cubic(discriminant: f64, a: f64, p: f64, q: f64) -> [Op
 	} 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.);
+		let square_root_discriminant = discriminant.sqrt();
 
 		roots[0] = Some(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 cube_root_r = (-p_divided_by_3).sqrt();
 		let phi = (-q / (2. * cube_root_r.powi(3))).acos();
 
 		let two_times_cube_root_r = 2. * cube_root_r;
@@ -271,8 +271,8 @@ pub fn dvec2_compare(a: DVec2, b: DVec2, max_abs_diff: f64) -> BVec2 {
 }
 
 /// 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)
+pub fn dvec2_approximately_in_range(point: DVec2, min_corner: DVec2, max_corner: DVec2, max_abs_diff: f64) -> BVec2 {
+	(point.cmpge(min_corner) & point.cmple(max_corner)) | dvec2_compare(point, min_corner, max_abs_diff) | dvec2_compare(point, max_corner, max_abs_diff)
 }
 
 /// Calculate a new position for a point given its original position, a unit vector in the desired direction, and a distance to move it by.

node-graph/gcore/Cargo.toml

@@ -52,16 +52,13 @@ image = { workspace = true, optional = true, default-features = false, features
 	"png",
 ] }
 specta = { workspace = true, optional = true }
-
 rustybuzz = { workspace = true, optional = true }
-
 num-derive = { workspace = true }
 num-traits = { workspace = true, default-features = false, features = ["i128"] }
-
-
 wasm-bindgen = { workspace = true, optional = true }
 js-sys = { workspace = true, optional = true }
 usvg = { workspace = true }
+rand = { workspace = true, default-features = false, features = ["std_rng"] }
 
 [dependencies.web-sys]
 workspace = true

node-graph/gcore/src/raster/color.rs

@@ -483,7 +483,7 @@ impl Color {
 		if luminance <= 0.008856 {
 			(luminance * 903.3) / 100.
 		} else {
-			(luminance.powf(1. / 3.) * 116. - 16.) / 100.
+			(luminance.cbrt() * 116. - 16.) / 100.
 		}
 	}
 

node-graph/gcore/src/vector/vector_nodes.rs

@@ -7,6 +7,7 @@ use core::future::Future;
 
 use bezier_rs::{Subpath, SubpathTValue, TValue};
 use glam::{DAffine2, DVec2};
+use rand::{Rng, SeedableRng};
 
 #[derive(Debug, Clone, Copy)]
 pub struct SetFillNode<FillType, SolidColor, GradientType, Start, End, Transform, Positions> {
@@ -276,6 +277,30 @@ async fn sample_points<FV: Future<Output = VectorData>, FL: Future<Output = Vec<
 	vector_data
 }
 
+#[derive(Debug, Clone, Copy)]
+pub struct PoissonDiskPoints<SeparationDiskDiameter> {
+	separation_disk_diameter: SeparationDiskDiameter,
+}
+
+#[node_macro::node_fn(PoissonDiskPoints)]
+fn poisson_disk_points(mut vector_data: VectorData, separation_disk_diameter: f32) -> VectorData {
+	let mut rng = rand::rngs::StdRng::seed_from_u64(0);
+	for subpath in &mut vector_data.subpaths.iter_mut() {
+		if subpath.manipulator_groups().len() < 3 {
+			continue;
+		}
+
+		subpath.apply_transform(vector_data.transform);
+
+		let points = subpath.poisson_disk_points(separation_disk_diameter as f64, || rng.gen::<f64>()).into_iter().map(|point| point.into());
+		*subpath = Subpath::from_anchors(points, false);
+
+		subpath.apply_transform(vector_data.transform.inverse());
+	}
+
+	vector_data
+}
+
 #[derive(Debug, Clone, Copy)]
 pub struct LengthsOfSegmentsOfSubpaths;
 
@@ -323,8 +348,10 @@ async fn morph<SourceFuture: Future<Output = VectorData>, TargetFuture: Future<O
 	source: impl Node<Footprint, Output = SourceFuture>,
 	target: impl Node<Footprint, Output = TargetFuture>,
 	start_index: u32,
-	time: f64,
+	time: f32,
 ) -> VectorData {
+	let time = time as f64;
+
 	let mut source = self.source.eval(footprint).await;
 	let mut target = self.target.eval(footprint).await;
 

node-graph/gstd/Cargo.toml

@@ -28,19 +28,19 @@ wayland = []
 
 [dependencies]
 fastnoise-lite = { workspace = true }
-rand = { workspace = true, features = [
+rand = { workspace = true, default-features = false, features = [
 	"alloc",
 	"small_rng",
-], default-features = false }
+] }
 rand_chacha = { workspace = true }
 autoquant = { git = "https://github.com/truedoctor/autoquant", optional = true, features = [
 	"fitting",
 ] }
-graphene-core = { path = "../gcore", features = [
+graphene-core = { path = "../gcore", default-features = false, features = [
 	"std",
 	"serde",
 	"alloc",
-], default-features = false }
+] }
 dyn-any = { path = "../../libraries/dyn-any", features = ["derive"] }
 graph-craft = { path = "../graph-craft", features = ["serde"] }
 vulkan-executor = { path = "../vulkan-executor", optional = true }

node-graph/interpreted-executor/src/node_registry.rs

@@ -735,9 +735,10 @@ fn node_registry() -> HashMap<ProtoNodeIdentifier, HashMap<NodeIOTypes, NodeCons
 		async_node!(graphene_core::vector::CopyToPoints<_, _>, input: Footprint, output: VectorData, fn_params: [Footprint => VectorData, Footprint => VectorData]),
 		async_node!(graphene_core::vector::CopyToPoints<_, _>, input: Footprint, output: GraphicGroup, fn_params: [Footprint => VectorData, Footprint => GraphicGroup]),
 		async_node!(graphene_core::vector::SamplePoints<_, _, _, _, _, _>, input: Footprint, output: VectorData, fn_params: [Footprint => VectorData, () => f32, () => f32, () => f32, () => bool, Footprint => Vec<Vec<f64>>]),
+		register_node!(graphene_core::vector::PoissonDiskPoints<_>, input: VectorData, params: [f32]),
 		register_node!(graphene_core::vector::LengthsOfSegmentsOfSubpaths, input: VectorData, params: []),
 		register_node!(graphene_core::vector::SplinesFromPointsNode, input: VectorData, params: []),
-		async_node!(graphene_core::vector::MorphNode<_, _, _, _>, input: Footprint, output: VectorData, fn_params: [Footprint => VectorData, Footprint => VectorData, () => u32, () => f64]),
+		async_node!(graphene_core::vector::MorphNode<_, _, _, _>, input: Footprint, output: VectorData, fn_params: [Footprint => VectorData, Footprint => VectorData, () => u32, () => f32]),
 		register_node!(graphene_core::vector::generator_nodes::CircleGenerator<_>, input: (), params: [f32]),
 		register_node!(graphene_core::vector::generator_nodes::EllipseGenerator<_, _>, input: (), params: [f32, f32]),
 		register_node!(graphene_core::vector::generator_nodes::RectangleGenerator<_, _>, input: (), params: [f32, f32]),

website/other/bezier-rs-demos/src/features/bezier-features.ts

@@ -409,7 +409,7 @@ const bezierFeatures = {
 		})(),
 	},
 	"intersect-linear": {
-		name: "Intersect (Line Segment)",
+		name: "Intersect (Linear Segment)",
 		callback: (bezier: WasmBezierInstance): string => {
 			const line = [
 				[45, 30],
@@ -419,7 +419,7 @@ const bezierFeatures = {
 		},
 	},
 	"intersect-quadratic": {
-		name: "Intersect (Quadratic)",
+		name: "Intersect (Quadratic Segment)",
 		callback: (bezier: WasmBezierInstance, options: Record<string, number>): string => {
 			const quadratic = [
 				[45, 80],
@@ -435,7 +435,7 @@ const bezierFeatures = {
 		},
 	},
 	"intersect-cubic": {
-		name: "Intersect (Cubic)",
+		name: "Intersect (Cubic Segment)",
 		callback: (bezier: WasmBezierInstance, options: Record<string, number>): string => {
 			const cubic = [
 				[65, 20],

website/other/bezier-rs-demos/src/features/subpath-features.ts

@@ -1,4 +1,4 @@
-import { capOptions, joinOptions, tSliderOptions, subpathTValueVariantOptions, intersectionErrorOptions, minimumSeparationOptions } from "@/utils/options";
+import { capOptions, joinOptions, tSliderOptions, subpathTValueVariantOptions, intersectionErrorOptions, minimumSeparationOptions, separationDiskDiameter } from "@/utils/options";
 import type { SubpathCallback, SubpathInputOption, WasmSubpathInstance } from "@/utils/types";
 import { SUBPATH_T_VALUE_VARIANTS } from "@/utils/types";
 
@@ -59,12 +59,17 @@ const subpathFeatures = {
 		name: "Bounding Box",
 		callback: (subpath: WasmSubpathInstance): string => subpath.bounding_box(),
 	},
+	"poisson-disk-points": {
+		name: "Poisson-Disk Points",
+		callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.poisson_disk_points(options.separation_disk_diameter),
+		inputOptions: [separationDiskDiameter],
+	},
 	inflections: {
 		name: "Inflections",
 		callback: (subpath: WasmSubpathInstance): string => subpath.inflections(),
 	},
 	"intersect-linear": {
-		name: "Intersect (Line Segment)",
+		name: "Intersect (Linear Segment)",
 		callback: (subpath: WasmSubpathInstance, options: Record<string, number>): string =>
 			subpath.intersect_line_segment(
 				[
@@ -105,11 +110,24 @@ const subpathFeatures = {
 			),
 		inputOptions: [intersectionErrorOptions, minimumSeparationOptions],
 	},
-	"self-intersect": {
-		name: "Self Intersect",
+	"intersect-self": {
+		name: "Intersect (Self)",
 		callback: (subpath: WasmSubpathInstance, options: Record<string, number>): string => subpath.self_intersections(options.error, options.minimum_separation),
 		inputOptions: [intersectionErrorOptions, minimumSeparationOptions],
 	},
+	"intersect-rectangle": {
+		name: "Intersect (Rectangle)",
+		callback: (subpath: WasmSubpathInstance, options: Record<string, number>): string =>
+			subpath.intersect_rectangle(
+				[
+					[75, 50],
+					[175, 150],
+				],
+				options.error,
+				options.minimum_separation,
+			),
+		inputOptions: [intersectionErrorOptions, minimumSeparationOptions],
+	},
 	curvature: {
 		name: "Curvature",
 		callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.curvature(options.t, SUBPATH_T_VALUE_VARIANTS[options.TVariant]),

website/other/bezier-rs-demos/src/utils/options.ts

@@ -1,11 +1,11 @@
 import { BEZIER_T_VALUE_VARIANTS, CAP_VARIANTS, JOIN_VARIANTS, SUBPATH_T_VALUE_VARIANTS } from "@/utils/types";
 
 export const tSliderOptions = {
+	variable: "t",
 	min: 0,
 	max: 1,
 	step: 0.01,
 	default: 0.5,
-	variable: "t",
 };
 
 export const errorOptions = {
@@ -32,6 +32,14 @@ export const intersectionErrorOptions = {
 	default: 0.02,
 };
 
+export const separationDiskDiameter = {
+	variable: "separation_disk_diameter",
+	min: 2.5,
+	max: 25,
+	step: 0.1,
+	default: 5,
+};
+
 export const bezierTValueVariantOptions = {
 	variable: "TVariant",
 	default: 0,

website/other/bezier-rs-demos/wasm/src/subpath.rs

@@ -4,6 +4,7 @@ use crate::utils::{parse_cap, parse_join};
 use bezier_rs::{Bezier, ManipulatorGroup, Subpath, SubpathTValue, TValueType};
 
 use glam::DVec2;
+use js_sys::Math;
 use std::fmt::Write;
 use wasm_bindgen::prelude::*;
 
@@ -171,7 +172,7 @@ impl WasmSubpath {
 			None => wrap_svg_tag(subpath_svg),
 			Some(bounding_box) => {
 				let content = format!(
-					"{subpath_svg}<rect x={} y ={} width=\"{}\" height=\"{}\" style=\"fill:{NONE};stroke:{RED};stroke-width:1\" />",
+					"{subpath_svg}<rect x={} y={} width=\"{}\" height=\"{}\" style=\"fill:{NONE};stroke:{RED};stroke-width:1\" />",
 					bounding_box[0].x,
 					bounding_box[0].y,
 					bounding_box[1].x - bounding_box[0].x,
@@ -182,6 +183,21 @@ impl WasmSubpath {
 		}
 	}
 
+	pub fn poisson_disk_points(&self, separation_disk_diameter: f64) -> String {
+		let r = separation_disk_diameter / 2.;
+
+		let subpath_svg = self.to_default_svg();
+		let points = self.0.poisson_disk_points(separation_disk_diameter, Math::random);
+
+		let points_style = format!("<style class=\"poisson\">style.poisson ~ circle {{ fill: {RED}; opacity: 0.25; }}</style>");
+		let content = points
+			.iter()
+			.map(|point| format!("<circle cx=\"{}\" cy=\"{}\" r=\"{r}\" />", point.x, point.y))
+			.collect::<Vec<_>>()
+			.join("");
+		wrap_svg_tag(format!("{subpath_svg}{points_style}{content}"))
+	}
+
 	pub fn inflections(&self) -> String {
 		let inflections: Vec<f64> = self.0.inflections();
 
@@ -342,6 +358,37 @@ impl WasmSubpath {
 		wrap_svg_tag(format!("{subpath_svg}{self_intersections_svg}"))
 	}
 
+	pub fn intersect_rectangle(&self, js_points: JsValue, error: f64, minimum_separation: f64) -> String {
+		let points: [DVec2; 2] = serde_wasm_bindgen::from_value(js_points).unwrap();
+
+		let subpath_svg = self.to_default_svg();
+
+		let mut rectangle_svg = String::new();
+		[
+			Bezier::from_linear_coordinates(points[0].x, points[0].y, points[1].x, points[0].y),
+			Bezier::from_linear_coordinates(points[1].x, points[0].y, points[1].x, points[1].y),
+			Bezier::from_linear_coordinates(points[1].x, points[1].y, points[0].x, points[1].y),
+			Bezier::from_linear_coordinates(points[0].x, points[1].y, points[0].x, points[0].y),
+		]
+		.iter()
+		.for_each(|line| line.to_svg(&mut rectangle_svg, CURVE_ATTRIBUTES.to_string().replace(BLACK, RED), String::new(), String::new(), String::new()));
+
+		let intersections_svg = self
+			.0
+			.rectangle_intersections(points[0], points[1], Some(error), Some(minimum_separation))
+			.iter()
+			.map(|(segment_index, intersection_t)| {
+				let point = self.0.evaluate(SubpathTValue::Parametric {
+					segment_index: *segment_index,
+					t: *intersection_t,
+				});
+				draw_circle(point, 4., RED, 1.5, WHITE)
+			})
+			.fold(String::new(), |acc, item| format!("{acc}{item}"));
+
+		wrap_svg_tag(format!("{subpath_svg}{rectangle_svg}{intersections_svg}"))
+	}
+
 	pub fn curvature(&self, t: f64, t_variant: String) -> String {
 		let subpath = self.to_default_svg();
 		let t = parse_t_variant(&t_variant, t);