Attribute-based vector format refactor (#1624)

* Initial vector format structure

* Click targets

* Code review pass

* Remove subpaths from vector data

* Morph node & vector node tests

* Insignificant change

---------

Co-authored-by: Keavon Chambers <keavon@keavon.com>
This commit is contained in:
0HyperCube 2024-03-09 18:27:30 +00:00 committed by GitHub
parent c8ea9e05a6
commit 218e9675fd
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
30 changed files with 991 additions and 436 deletions

View file

@ -1,4 +1,4 @@
use crate::utils::{f64_compare, TValue, TValueType};
use crate::utils::{TValue, TValueType};
use super::*;
@ -21,44 +21,57 @@ impl Bezier {
return 1.;
}
let mut low = 0.;
let mut mid = 0.5;
let mut high = 1.;
match self.handles {
BezierHandles::Linear => euclidean_t,
BezierHandles::Quadratic { handle } => {
// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
fn recurse(a0: DVec2, a1: DVec2, a2: DVec2, level: u8, desired_len: f64) -> (f64, f64) {
let lower = a0.distance(a2);
let upper = a0.distance(a1) + a1.distance(a2);
if level >= 8 {
let approx_len = (lower + upper) / 2.;
return (approx_len, desired_len / approx_len);
}
// The euclidean t-value input generally correlates with the parametric t-value result.
// So we can assume a low t-value has a short length from the start of the curve, and a high t-value has a short length from the end of the curve.
// We'll use a strategy where we measure from either end of the curve depending on which side is closer than thus more likely to be proximate to the sought parametric t-value.
// This allows us to use fewer segments to approximate the curve, which usually won't go much beyond half the curve.
let result_likely_closer_to_start = euclidean_t < 0.5;
// If the curve is near either end, we need even fewer segments to approximate the curve with reasonable accuracy.
// A point that's likely near the center is the worst case where we need to use up to half the predefined number of max subdivisions.
let subdivisions_proportional_to_likely_length = ((euclidean_t - 0.5).abs() * DEFAULT_LENGTH_SUBDIVISIONS as f64).round().max(1.) as usize;
let b1 = 0.5 * (a0 + a1);
let c1 = 0.5 * (a1 + a2);
let b2 = 0.5 * (b1 + c1);
let (first_len, t) = recurse(a0, b1, b2, level + 1, desired_len);
if first_len > desired_len {
return (first_len, t * 0.5);
}
let (second_len, t) = recurse(b2, c1, a2, level + 1, desired_len - first_len);
(first_len + second_len, t * 0.5 + 0.5)
}
recurse(self.start, handle, self.end, 0, total_length * euclidean_t).1
}
BezierHandles::Cubic { handle_start, handle_end } => {
// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
fn recurse(a0: DVec2, a1: DVec2, a2: DVec2, a3: DVec2, level: u8, desired_len: f64) -> (f64, f64) {
let lower = a0.distance(a3);
let upper = a0.distance(a1) + a1.distance(a2) + a2.distance(a3);
if level >= 8 {
let approx_len = (lower + upper) / 2.;
return (approx_len, desired_len / approx_len);
}
// Binary search for the parametric t-value that corresponds to the euclidean distance ratio by trimming the curve between the start and the tested parametric t-value during each iteration of the search.
while low < high {
mid = (low + high) / 2.;
// We can search from the curve start to the sought point, or from the sought point to the curve end, depending on which side is likely closer to the result.
let current_length = if result_likely_closer_to_start {
let trimmed = self.trim(TValue::Parametric(0.), TValue::Parametric(mid));
trimmed.length(Some(subdivisions_proportional_to_likely_length))
} else {
let trimmed = self.trim(TValue::Parametric(mid), TValue::Parametric(1.));
let trimmed_length = trimmed.length(Some(subdivisions_proportional_to_likely_length));
total_length - trimmed_length
};
let current_euclidean_t = current_length / total_length;
if f64_compare(current_euclidean_t, euclidean_t, error) {
break;
} else if current_euclidean_t < euclidean_t {
low = mid;
} else {
high = mid;
let b1 = 0.5 * (a0 + a1);
let t0 = 0.5 * (a1 + a2);
let c1 = 0.5 * (a2 + a3);
let b2 = 0.5 * (b1 + t0);
let c2 = 0.5 * (t0 + c1);
let b3 = 0.5 * (b2 + c2);
let (first_len, t) = recurse(a0, b1, b2, b3, level + 1, desired_len);
if first_len > desired_len {
return (first_len, t * 0.5);
}
let (second_len, t) = recurse(b3, c2, c1, a3, level + 1, desired_len - first_len);
(first_len + second_len, t * 0.5 + 0.5)
}
recurse(self.start, handle_start, handle_end, self.end, 0, total_length * euclidean_t).1
}
}
mid
.clamp(0., 1.)
}
/// Convert a [TValue] to a parametric `t`-value.
@ -109,133 +122,86 @@ impl Bezier {
/// Return a selection of equidistant points on the bezier curve.
/// If no value is provided for `steps`, then the function will default `steps` to be 10.
/// <iframe frameBorder="0" width="100%" height="350px" src="https://graphite.rs/libraries/bezier-rs#bezier/lookup-table/solo" title="Lookup-Table Demo"></iframe>
pub fn compute_lookup_table(&self, steps: Option<usize>, tvalue_type: Option<TValueType>) -> Vec<DVec2> {
pub fn compute_lookup_table(&self, steps: Option<usize>, tvalue_type: Option<TValueType>) -> impl Iterator<Item = DVec2> + '_ {
let steps = steps.unwrap_or(DEFAULT_LUT_STEP_SIZE);
let tvalue_type = tvalue_type.unwrap_or(TValueType::Parametric);
(0..=steps)
.map(|t| {
let tvalue = match tvalue_type {
TValueType::Parametric => TValue::Parametric(t as f64 / steps as f64),
TValueType::Euclidean => TValue::Euclidean(t as f64 / steps as f64),
};
self.evaluate(tvalue)
})
.collect()
(0..=steps).map(move |t| {
let tvalue = match tvalue_type {
TValueType::Parametric => TValue::Parametric(t as f64 / steps as f64),
TValueType::Euclidean => TValue::Euclidean(t as f64 / steps as f64),
};
self.evaluate(tvalue)
})
}
/// Return an approximation of the length of the bezier curve.
/// - `num_subdivisions` - Number of subdivisions used to approximate the curve. The default value is 1000.
/// - `tolerance` - Tolerance used to approximate the curve.
/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/length/solo" title="Length Demo"></iframe>
pub fn length(&self, num_subdivisions: Option<usize>) -> f64 {
pub fn length(&self, tolerance: Option<f64>) -> f64 {
match self.handles {
BezierHandles::Linear => (self.start - self.end).length(),
_ => {
// Code example from <https://gamedev.stackexchange.com/questions/5373/moving-ships-between-two-planets-along-a-bezier-missing-some-equations-for-acce/5427#5427>.
BezierHandles::Quadratic { handle } => {
// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
fn recurse(a0: DVec2, a1: DVec2, a2: DVec2, tolerance: f64, level: u8) -> f64 {
let lower = a0.distance(a2);
let upper = a0.distance(a1) + a1.distance(a2);
if upper - lower <= 2. * tolerance || level >= 8 {
return (lower + upper) / 2.;
}
// We will use an approximate approach where we split the curve into many subdivisions
// and calculate the euclidean distance between the two endpoints of the subdivision
let lookup_table = self.compute_lookup_table(Some(num_subdivisions.unwrap_or(DEFAULT_LENGTH_SUBDIVISIONS)), Some(TValueType::Parametric));
let approx_curve_length: f64 = lookup_table.windows(2).map(|points| (points[1] - points[0]).length()).sum();
let b1 = 0.5 * (a0 + a1);
let c1 = 0.5 * (a1 + a2);
let b2 = 0.5 * (b1 + c1);
recurse(a0, b1, b2, 0.5 * tolerance, level + 1) + recurse(b2, c1, a2, 0.5 * tolerance, level + 1)
}
recurse(self.start, handle, self.end, tolerance.unwrap_or_default(), 0)
}
BezierHandles::Cubic { handle_start, handle_end } => {
// Use Casteljau subdivision, noting that the length is more than the straight line distance from start to end but less than the straight line distance through the handles
fn recurse(a0: DVec2, a1: DVec2, a2: DVec2, a3: DVec2, tolerance: f64, level: u8) -> f64 {
let lower = a0.distance(a3);
let upper = a0.distance(a1) + a1.distance(a2) + a2.distance(a3);
if upper - lower <= 2. * tolerance || level >= 8 {
return (lower + upper) / 2.;
}
approx_curve_length
let b1 = 0.5 * (a0 + a1);
let t0 = 0.5 * (a1 + a2);
let c1 = 0.5 * (a2 + a3);
let b2 = 0.5 * (b1 + t0);
let c2 = 0.5 * (t0 + c1);
let b3 = 0.5 * (b2 + c2);
recurse(a0, b1, b2, b3, 0.5 * tolerance, level + 1) + recurse(b3, c2, c1, a3, 0.5 * tolerance, level + 1)
}
recurse(self.start, handle_start, handle_end, self.end, tolerance.unwrap_or_default(), 0)
}
}
}
/// Returns the parametric `t`-value that corresponds to the closest point on the curve to the provided point.
/// Uses a searching algorithm akin to binary search that can be customized using the optional [ProjectionOptions] struct.
/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/project/solo" title="Project Demo"></iframe>
pub fn project(&self, point: DVec2, options: Option<ProjectionOptions>) -> f64 {
let options = options.unwrap_or_default();
let ProjectionOptions {
lut_size,
convergence_epsilon,
convergence_limit,
iteration_limit,
} = options;
pub fn project(&self, point: DVec2) -> f64 {
let sbasis = crate::symmetrical_basis::to_symmetrical_basis_pair(*self);
let derivative = sbasis.derivative();
let dd = (sbasis - point).dot(&derivative);
let roots = dd.roots();
// TODO: Consider optimizations from precomputing useful values, or using the GPU
// First find the closest point from the results of a lookup table
let lut = self.compute_lookup_table(Some(lut_size), Some(TValueType::Parametric));
let (minimum_position, minimum_distance) = utils::get_closest_point_in_lut(&lut, point);
let mut closest = 0.;
let mut min_dist_squared = self.evaluate(TValue::Parametric(0.)).distance_squared(point);
// Get the t values to the left and right of the closest result in the lookup table
let lut_size_f64 = lut_size as f64;
let minimum_position_f64 = minimum_position as f64;
let mut left_t = (minimum_position_f64 - 1.).max(0.) / lut_size_f64;
let mut right_t = (minimum_position_f64 + 1.).min(lut_size_f64) / lut_size_f64;
// Perform a finer search by finding closest t from 5 points between [left_t, right_t] inclusive
// Choose new left_t and right_t for a smaller range around the closest t and repeat the process
let mut final_t = left_t;
let mut distance;
// Increment minimum_distance to ensure that the distance < minimum_distance comparison will be true for at least one iteration
let mut new_minimum_distance = minimum_distance + 1.;
// Maintain the previous distance to identify convergence
let mut previous_distance;
// Counter to limit the number of iterations
let mut iteration_count = 0;
// Counter to identify how many iterations have had a similar result. Used for convergence test
let mut convergence_count = 0;
// Store calculated distances to minimize unnecessary recomputations
let mut distances: [f64; NUM_DISTANCES] = [
point.distance(lut[(minimum_position as i64 - 1).max(0) as usize]),
0.,
0.,
0.,
point.distance(lut[lut_size.min(minimum_position + 1)]),
];
while left_t <= right_t && convergence_count < convergence_limit && iteration_count < iteration_limit {
previous_distance = new_minimum_distance;
let step = (right_t - left_t) / (NUM_DISTANCES as f64 - 1.);
let mut iterator_t = left_t;
let mut target_index = 0;
// Iterate through first 4 points and will handle the right most point later
for (step_index, table_distance) in distances.iter_mut().enumerate().take(4) {
// Use previously computed distance for the left most point, and compute new values for the others
if step_index == 0 {
distance = *table_distance;
} else {
distance = point.distance(self.evaluate(TValue::Parametric(iterator_t)));
*table_distance = distance;
}
if distance < new_minimum_distance {
new_minimum_distance = distance;
target_index = step_index;
final_t = iterator_t
}
iterator_t += step;
}
// Check right most edge separately since step may not perfectly add up to it (floating point errors)
if distances[NUM_DISTANCES - 1] < new_minimum_distance {
new_minimum_distance = distances[NUM_DISTANCES - 1];
final_t = right_t;
}
// Update left_t and right_t to be the t values (final_t +/- step), while handling the edges (i.e. if final_t is 0, left_t will be 0 instead of -step)
// Ensure that the t values never exceed the [0, 1] range
left_t = (final_t - step).max(0.);
right_t = (final_t + step).min(1.);
// Re-use the corresponding computed distances (target_index is the index corresponding to final_t)
// Since target_index is a u_size, can't subtract one if it is zero
distances[0] = distances[if target_index == 0 { 0 } else { target_index - 1 }];
distances[NUM_DISTANCES - 1] = distances[(target_index + 1).min(NUM_DISTANCES - 1)];
iteration_count += 1;
// update count for consecutive iterations of similar minimum distances
if previous_distance - new_minimum_distance < convergence_epsilon {
convergence_count += 1;
} else {
convergence_count = 0;
for time in roots {
let distance = self.evaluate(TValue::Parametric(time)).distance_squared(point);
if distance < min_dist_squared {
closest = time;
min_dist_squared = distance;
}
}
final_t
if self.evaluate(TValue::Parametric(1.)).distance_squared(point) < min_dist_squared {
closest = 1.;
}
closest
}
}
@ -259,11 +225,11 @@ mod tests {
#[test]
fn test_compute_lookup_table() {
let bezier1 = Bezier::from_quadratic_coordinates(10., 10., 30., 30., 50., 10.);
let lookup_table1 = bezier1.compute_lookup_table(Some(2), Some(TValueType::Parametric));
let lookup_table1 = bezier1.compute_lookup_table(Some(2), Some(TValueType::Parametric)).collect::<Vec<_>>();
assert_eq!(lookup_table1, vec![bezier1.start(), bezier1.evaluate(TValue::Parametric(0.5)), bezier1.end()]);
let bezier2 = Bezier::from_cubic_coordinates(10., 10., 30., 30., 70., 70., 90., 10.);
let lookup_table2 = bezier2.compute_lookup_table(Some(4), Some(TValueType::Parametric));
let lookup_table2 = bezier2.compute_lookup_table(Some(4), Some(TValueType::Parametric)).collect::<Vec<_>>();
assert_eq!(
lookup_table2,
vec![
@ -296,10 +262,10 @@ mod tests {
#[test]
fn test_project() {
let bezier1 = Bezier::from_cubic_coordinates(4., 4., 23., 45., 10., 30., 56., 90.);
assert_eq!(bezier1.project(DVec2::ZERO, None), 0.);
assert_eq!(bezier1.project(DVec2::new(100., 100.), None), 1.);
assert_eq!(bezier1.project(DVec2::ZERO), 0.);
assert_eq!(bezier1.project(DVec2::new(100., 100.)), 1.);
let bezier2 = Bezier::from_quadratic_coordinates(0., 0., 0., 100., 100., 100.);
assert_eq!(bezier2.project(DVec2::new(100., 0.), None), 0.);
assert_eq!(bezier2.project(DVec2::new(100., 0.)), 0.);
}
}

View file

@ -57,20 +57,12 @@ impl Bezier {
/// Get the coordinates of the bezier segment's first handle point. This represents the only handle in a quadratic segment.
pub fn handle_start(&self) -> Option<DVec2> {
match self.handles {
BezierHandles::Linear => None,
BezierHandles::Quadratic { handle } => Some(handle),
BezierHandles::Cubic { handle_start, .. } => Some(handle_start),
}
self.handles.start()
}
/// Get the coordinates of the second handle point. This will return `None` for a quadratic segment.
pub fn handle_end(&self) -> Option<DVec2> {
match self.handles {
BezierHandles::Linear { .. } => None,
BezierHandles::Quadratic { .. } => None,
BezierHandles::Cubic { handle_end, .. } => Some(handle_end),
}
self.handles.end()
}
/// Get an iterator over the coordinates of all points in a vector.

View file

@ -14,7 +14,7 @@ use glam::DVec2;
use std::fmt::{Debug, Formatter, Result};
/// Representation of the handle point(s) in a bezier segment.
#[derive(Copy, Clone, PartialEq)]
#[derive(Copy, Clone, PartialEq, Debug)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum BezierHandles {
Linear,
@ -31,10 +31,55 @@ pub enum BezierHandles {
handle_end: DVec2,
},
}
impl std::hash::Hash for BezierHandles {
fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
std::mem::discriminant(self).hash(state);
match self {
BezierHandles::Linear => {}
BezierHandles::Quadratic { handle } => handle.to_array().map(|v| v.to_bits()).hash(state),
BezierHandles::Cubic { handle_start, handle_end } => [handle_start, handle_end].map(|handle| handle.to_array().map(|v| v.to_bits())).hash(state),
}
}
}
impl BezierHandles {
pub fn is_cubic(&self) -> bool {
matches!(self, Self::Cubic { .. })
}
/// Get the coordinates of the bezier segment's first handle point. This represents the only handle in a quadratic segment.
pub fn start(&self) -> Option<DVec2> {
match *self {
BezierHandles::Cubic { handle_start, .. } | BezierHandles::Quadratic { handle: handle_start } => Some(handle_start),
_ => None,
}
}
/// Get the coordinates of the second handle point. This will return `None` for a quadratic segment.
pub fn end(&self) -> Option<DVec2> {
match *self {
BezierHandles::Cubic { handle_end, .. } => Some(handle_end),
_ => None,
}
}
/// Returns a Bezier curve that results from applying the transformation function to each handle point in the Bezier.
#[must_use]
pub fn apply_transformation(&self, transformation_function: impl Fn(DVec2) -> DVec2) -> Self {
match *self {
BezierHandles::Linear => Self::Linear,
BezierHandles::Quadratic { handle } => {
let handle = transformation_function(handle);
Self::Quadratic { handle }
}
BezierHandles::Cubic { handle_start, handle_end } => {
let handle_start = transformation_function(handle_start);
let handle_end = transformation_function(handle_end);
Self::Cubic { handle_start, handle_end }
}
}
}
}
#[cfg(feature = "dyn-any")]

View file

@ -1,30 +1,6 @@
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 {
Self {
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 {

View file

@ -105,19 +105,10 @@ impl Bezier {
/// Returns a Bezier curve that results from applying the transformation function to each point in the Bezier.
pub fn apply_transformation(&self, transformation_function: impl Fn(DVec2) -> DVec2) -> Bezier {
let transformed_start = transformation_function(self.start);
let transformed_end = transformation_function(self.end);
match self.handles {
BezierHandles::Linear => Bezier::from_linear_dvec2(transformed_start, transformed_end),
BezierHandles::Quadratic { handle } => {
let transformed_handle = transformation_function(handle);
Bezier::from_quadratic_dvec2(transformed_start, transformed_handle, transformed_end)
}
BezierHandles::Cubic { handle_start, handle_end } => {
let transformed_handle_start = transformation_function(handle_start);
let transformed_handle_end = transformation_function(handle_end);
Bezier::from_cubic_dvec2(transformed_start, transformed_handle_start, transformed_handle_end, transformed_end)
}
Self {
start: transformation_function(self.start),
end: transformation_function(self.end),
handles: self.handles.apply_transformation(transformation_function),
}
}
@ -315,12 +306,12 @@ impl Bezier {
BezierHandles::Linear => Bezier::from_linear_dvec2(transformed_start, transformed_end),
BezierHandles::Quadratic { handle: _ } => unreachable!(),
BezierHandles::Cubic { handle_start, handle_end } => {
let handle_start_closest_t = intermediate.project(handle_start, None);
let handle_start_closest_t = intermediate.project(handle_start);
let handle_start_scale_distance = (1. - handle_start_closest_t) * start_distance + handle_start_closest_t * end_distance;
let transformed_handle_start =
utils::scale_point_from_direction_vector(handle_start, intermediate.normal(TValue::Parametric(handle_start_closest_t)), false, handle_start_scale_distance);
let handle_end_closest_t = intermediate.project(handle_start, None);
let handle_end_closest_t = intermediate.project(handle_start);
let handle_end_scale_distance = (1. - handle_end_closest_t) * start_distance + handle_end_closest_t * end_distance;
let transformed_handle_end = utils::scale_point_from_direction_vector(handle_end, intermediate.normal(TValue::Parametric(handle_end_closest_t)), false, handle_end_scale_distance);
Bezier::from_cubic_dvec2(transformed_start, transformed_handle_start, transformed_handle_end, transformed_end)
@ -810,7 +801,7 @@ mod tests {
.iter()
.map(|t| {
let offset_point = offset_segment.evaluate(TValue::Parametric(*t));
let closest_point_t = bezier.project(offset_point, None);
let closest_point_t = bezier.project(offset_point);
let closest_point = bezier.evaluate(TValue::Parametric(closest_point_t));
let actual_distance = offset_point.distance(closest_point);

View file

@ -4,8 +4,6 @@
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.;
/// Minimum allowable separation between adjacent `t` values when calculating curve intersections
@ -19,8 +17,6 @@ pub const DEFAULT_EUCLIDEAN_ERROR_BOUND: f64 = 0.001;
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;

View file

@ -8,6 +8,7 @@ use std::fmt::Write;
impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
/// Create a new `Subpath` using a list of [ManipulatorGroup]s.
/// A `Subpath` with less than 2 [ManipulatorGroup]s may not be closed.
#[track_caller]
pub fn new(manipulator_groups: Vec<ManipulatorGroup<ManipulatorGroupId>>, closed: bool) -> Self {
assert!(!closed || manipulator_groups.len() > 1, "A closed Subpath must contain more than 1 ManipulatorGroup.");
Self { manipulator_groups, closed }
@ -276,61 +277,71 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
// Number of points = number of points to find handles for
let len_points = points.len();
// matrix coefficients a, b and c (see https://mathworld.wolfram.com/CubicSpline.html)
// because the 'a' coefficients are all 1 they need not be stored
// this algorithm does a variation of the above algorithm.
// Instead of using the traditional cubic: a + bt + ct^2 + dt^3, we use the bezier cubic.
let mut b = vec![DVec2::new(4., 4.); len_points];
b[0] = DVec2::new(2., 2.);
b[len_points - 1] = DVec2::new(2., 2.);
let mut c = vec![DVec2::new(1., 1.); len_points];
// 'd' is the the second point in a cubic bezier, which is what we solve for
let mut d = vec![DVec2::ZERO; len_points];
d[0] = DVec2::new(2. * points[1].x + points[0].x, 2. * points[1].y + points[0].y);
d[len_points - 1] = DVec2::new(3. * points[len_points - 1].x, 3. * points[len_points - 1].y);
for idx in 1..(len_points - 1) {
d[idx] = DVec2::new(4. * points[idx].x + 2. * points[idx + 1].x, 4. * points[idx].y + 2. * points[idx + 1].y);
}
// Solve with Thomas algorithm (see https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)
// do row operations to eliminate `a` coefficients
c[0] /= -b[0];
d[0] /= -b[0];
#[allow(clippy::assign_op_pattern)]
for i in 1..len_points {
b[i] += c[i - 1];
// for some reason the below line makes the borrow checker mad
//d[i] += d[i-1]
d[i] = d[i] + d[i - 1];
c[i] /= -b[i];
d[i] /= -b[i];
}
// at this point b[i] == -a[i + 1], a[i] == 0,
// do row operations to eliminate 'c' coefficients and solve
d[len_points - 1] *= -1.;
#[allow(clippy::assign_op_pattern)]
for i in (0..len_points - 1).rev() {
d[i] = d[i] - (c[i] * d[i + 1]);
d[i] *= -1.; //d[i] /= b[i]
}
let out_handles = solve_spline_first_handle(&points);
let mut subpath = Subpath::new(Vec::new(), false);
// given the second point in the n'th cubic bezier, the third point is given by 2 * points[n+1] - b[n+1].
// to find 'handle1_pos' for the n'th point we need the n-1 cubic bezier
subpath.manipulator_groups.push(ManipulatorGroup::new(points[0], None, Some(d[0])));
subpath.manipulator_groups.push(ManipulatorGroup::new(points[0], None, Some(out_handles[0])));
for i in 1..len_points - 1 {
subpath.manipulator_groups.push(ManipulatorGroup::new(points[i], Some(2. * points[i] - d[i]), Some(d[i])));
subpath
.manipulator_groups
.push(ManipulatorGroup::new(points[i], Some(2. * points[i] - out_handles[i]), Some(out_handles[i])));
}
subpath
.manipulator_groups
.push(ManipulatorGroup::new(points[len_points - 1], Some(2. * points[len_points - 1] - d[len_points - 1]), None));
.push(ManipulatorGroup::new(points[len_points - 1], Some(2. * points[len_points - 1] - out_handles[len_points - 1]), None));
subpath
}
}
pub fn solve_spline_first_handle(points: &[DVec2]) -> Vec<DVec2> {
let len_points = points.len();
// matrix coefficients a, b and c (see https://mathworld.wolfram.com/CubicSpline.html)
// because the 'a' coefficients are all 1 they need not be stored
// this algorithm does a variation of the above algorithm.
// Instead of using the traditional cubic: a + bt + ct^2 + dt^3, we use the bezier cubic.
let mut b = vec![DVec2::new(4., 4.); len_points];
b[0] = DVec2::new(2., 2.);
b[len_points - 1] = DVec2::new(2., 2.);
let mut c = vec![DVec2::new(1., 1.); len_points];
// 'd' is the the second point in a cubic bezier, which is what we solve for
let mut d = vec![DVec2::ZERO; len_points];
d[0] = DVec2::new(2. * points[1].x + points[0].x, 2. * points[1].y + points[0].y);
d[len_points - 1] = DVec2::new(3. * points[len_points - 1].x, 3. * points[len_points - 1].y);
for idx in 1..(len_points - 1) {
d[idx] = DVec2::new(4. * points[idx].x + 2. * points[idx + 1].x, 4. * points[idx].y + 2. * points[idx + 1].y);
}
// Solve with Thomas algorithm (see https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)
// do row operations to eliminate `a` coefficients
c[0] /= -b[0];
d[0] /= -b[0];
#[allow(clippy::assign_op_pattern)]
for i in 1..len_points {
b[i] += c[i - 1];
// for some reason the below line makes the borrow checker mad
//d[i] += d[i-1]
d[i] = d[i] + d[i - 1];
c[i] /= -b[i];
d[i] /= -b[i];
}
// at this point b[i] == -a[i + 1], a[i] == 0,
// do row operations to eliminate 'c' coefficients and solve
d[len_points - 1] *= -1.;
#[allow(clippy::assign_op_pattern)]
for i in (0..len_points - 1).rev() {
d[i] = d[i] - (c[i] * d[i + 1]);
d[i] *= -1.; //d[i] /= b[i]
}
d
}

View file

@ -1,7 +1,6 @@
use super::*;
use crate::consts::{DEFAULT_EUCLIDEAN_ERROR_BOUND, DEFAULT_LUT_STEP_SIZE};
use crate::utils::{SubpathTValue, TValue, TValueType};
use crate::ProjectionOptions;
use glam::DVec2;
/// Functionality relating to looking up properties of the `Subpath` or points along the `Subpath`.
@ -25,10 +24,10 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
}
/// 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`.
/// - `tolerance` - Tolerance used to approximate the curve.
/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#subpath/length/solo" title="Length Demo"></iframe>
pub fn length(&self, num_subdivisions: Option<usize>) -> f64 {
self.iter().map(|bezier| bezier.length(num_subdivisions)).sum()
pub fn length(&self, tolerance: Option<f64>) -> f64 {
self.iter().map(|bezier| bezier.length(tolerance)).sum()
}
/// Converts from a subpath (composed of multiple segments) to a point along a certain segment represented.
@ -98,9 +97,8 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
}
/// Returns the segment index and `t` value that corresponds to the closest point on the curve to the provided point.
/// Uses a searching algorithm akin to binary search that can be customized using the [ProjectionOptions] structure.
/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#subpath/project/solo" title="Project Demo"></iframe>
pub fn project(&self, point: DVec2, options: Option<ProjectionOptions>) -> Option<(usize, f64)> {
pub fn project(&self, point: DVec2) -> Option<(usize, f64)> {
if self.is_empty() {
return None;
}
@ -109,7 +107,7 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
let (index, (_, project_t)) = self
.iter()
.map(|bezier| {
let project_t = bezier.project(point, options);
let project_t = bezier.project(point);
(bezier.evaluate(TValue::Parametric(project_t)).distance(point), project_t)
})
.enumerate()

View file

@ -4,6 +4,7 @@ mod manipulators;
mod solvers;
mod structs;
mod transform;
pub use core::*;
pub use structs::*;
use crate::Bezier;

View file

@ -296,8 +296,8 @@ impl<ManipulatorGroupId: crate::Identifier> Subpath<ManipulatorGroupId> {
let start_tangent = second_bezier.non_normalized_tangent(0.);
// Compute an average unit vector, weighing the segments by a rough estimation of their relative size.
let segment1_len = first_bezier.length(Some(5));
let segment2_len = second_bezier.length(Some(5));
let segment1_len = first_bezier.length(None);
let segment2_len = second_bezier.length(None);
let average_unit_tangent = (end_tangent.normalize() * segment1_len + start_tangent.normalize() * segment2_len) / (segment1_len + segment2_len);
// Adjust start and end handles to fit the average tangent

View file

@ -90,11 +90,6 @@ pub fn compute_abc_for_cubic_through_points(start_point: DVec2, point_on_curve:
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).total_cmp(&(y.1))).unwrap()
}
/// Find the roots of the linear equation `ax + b`.
pub fn solve_linear(a: f64, b: f64) -> [Option<f64>; 3] {
// There exist roots when `a` is not 0

View file

@ -259,10 +259,13 @@ impl_type!(
);
#[cfg(feature = "std")]
use std::sync::*;
use std::{
collections::{HashMap, HashSet},
sync::*,
};
#[cfg(feature = "std")]
impl_type!(Once, Mutex<T>, RwLock<T>);
impl_type!(Once, Mutex<T>, RwLock<T>, HashSet<T>, HashMap<K, V>);
#[cfg(feature = "rc")]
use std::rc::Rc;