diff --git a/node-graph/gcore/src/vector/algorithms/bezpath_algorithms.rs b/node-graph/gcore/src/vector/algorithms/bezpath_algorithms.rs index c61140e68..6c5b12746 100644 --- a/node-graph/gcore/src/vector/algorithms/bezpath_algorithms.rs +++ b/node-graph/gcore/src/vector/algorithms/bezpath_algorithms.rs @@ -1,10 +1,14 @@ +use super::poisson_disk::poisson_disk_sample; +use crate::vector::PointId; +use bezier_rs::Subpath; +use glam::{DAffine2, DVec2}; +use kurbo::{BezPath, ParamCurve, ParamCurveDeriv, PathSeg, Point, Shape}; + /// Accuracy to find the position on [kurbo::Bezpath]. const POSITION_ACCURACY: f64 = 1e-5; /// Accuracy to find the length of the [kurbo::PathSeg]. pub const PERIMETER_ACCURACY: f64 = 1e-5; -use kurbo::{BezPath, ParamCurve, ParamCurveDeriv, PathSeg, Point, Shape}; - pub fn position_on_bezpath(bezpath: &BezPath, t: f64, euclidian: bool, segments_length: Option<&[f64]>) -> Point { let (segment_index, t) = t_value_to_parametric(bezpath, t, euclidian, segments_length); bezpath.get_seg(segment_index + 1).unwrap().eval(t) @@ -184,3 +188,55 @@ fn bezpath_t_value_to_parametric(bezpath: &kurbo::BezPath, t: BezPathTValue, pre } } } + +/// 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(this: &Subpath, separation_disk_diameter: f64, rng: impl FnMut() -> f64, subpaths: &[(Subpath, [DVec2; 2])], subpath_index: usize) -> Vec { + let Some(bounding_box) = this.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 = this.clone(); + shape.set_closed(true); + shape.apply_transform(DAffine2::from_translation((-offset_x, -offset_y).into())); + + let point_in_shape_checker = |point: DVec2| { + // Check against all paths the point is contained in to compute the correct winding number + let mut number = 0; + for (i, (shape, bb)) in subpaths.iter().enumerate() { + let point = point + bounding_box[0]; + if bb[0].x > point.x || bb[0].y > point.y || bb[1].x < point.x || bb[1].y < point.y { + continue; + } + let winding = shape.winding_order(point); + + if i == subpath_index && winding == 0 { + return false; + } + number += winding; + } + number != 0 + }; + + let square_edges_intersect_shape_checker = |corner1: DVec2, size: f64| { + let corner2 = corner1 + DVec2::splat(size); + this.rectangle_intersections_exist(corner1, corner2) + }; + + let mut points = 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 +} diff --git a/node-graph/gcore/src/vector/algorithms/mod.rs b/node-graph/gcore/src/vector/algorithms/mod.rs index 54fffd0d4..2c8b5833f 100644 --- a/node-graph/gcore/src/vector/algorithms/mod.rs +++ b/node-graph/gcore/src/vector/algorithms/mod.rs @@ -2,3 +2,4 @@ pub mod bezpath_algorithms; mod instance; mod merge_by_distance; pub mod offset_subpath; +mod poisson_disk; diff --git a/node-graph/gcore/src/vector/algorithms/poisson_disk.rs b/node-graph/gcore/src/vector/algorithms/poisson_disk.rs new file mode 100644 index 000000000..cd2d8d45f --- /dev/null +++ b/node-graph/gcore/src/vector/algorithms/poisson_disk.rs @@ -0,0 +1,369 @@ +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" +/// +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 { + 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) && point_in_shape_checker(sub_square + subdivided_size); + // if !square_edges_intersect_shape_checker(sub_square, subdivided_size) { assert_eq!(point_in_shape_checker(sub_square), point_in_shape_checker(sub_square + subdivided_size)); } + // Sometimes this fails so it is necessary to also check the bottom right corner. + + 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: A, b: B) -> impl Iterator +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, + /// 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_checker(corner + square_size) && !square_edges_intersect_shape; + // if !square_edges_intersect_shape { assert_eq!(point_in_shape_checker(corner), point_in_shape_checker(corner + square_size)); } + // Sometimes this fails so it is necessary to also check the bottom right corner. + 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) { + 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 { + // 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 { + // 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, +} + +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 { + 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 { + self.cells.iter().flat_map(|cell| cell.list_cell()).collect() + } +}