Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions sparse_strips/vello_common/src/flatten.rs
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ pub use crate::flatten_simd::FlattenCtx;

/// The flattening tolerance.
const TOL: f64 = 0.25;
pub(crate) const TOL_2: f64 = TOL * TOL;

/// A point.
#[derive(Clone, Copy, Debug, PartialEq)]
Expand Down
104 changes: 75 additions & 29 deletions sparse_strips/vello_common/src/flatten_simd.rs
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,15 @@
//! well as some code that was copied from kurbo, which is needed to reimplement the
//! full `flatten` method.

use crate::flatten::TOL_2;
#[cfg(not(feature = "std"))]
use crate::kurbo::common::FloatFuncs as _;
use crate::kurbo::{CubicBez, ParamCurve, PathEl, Point, QuadBez};
use alloc::vec;
use alloc::vec::Vec;
use bytemuck::{Pod, Zeroable};
use fearless_simd::*;

#[cfg(not(feature = "std"))]
use crate::kurbo::common::FloatFuncs as _;

// Unlike kurbo, which takes a closure with a callback for outputting the lines, we use a trait
// instead. The reason is that this way the callback can be inlined, which is not possible with
// a closure and turned out to have a noticeable overhead.
Expand Down Expand Up @@ -49,38 +49,84 @@ pub(crate) fn flatten<S: Simd>(
}
PathEl::QuadTo(p1, p2) => {
if let Some(p0) = last_pt {
let q = QuadBez::new(p0, p1, p2);
let params = q.estimate_subdiv(sqrt_tol);
let n = ((0.5 * params.val / sqrt_tol).ceil() as usize).max(1);
let step = 1.0 / (n as f64);
for i in 1..n {
let u = (i as f64) * step;
let t = q.determine_subdiv_t(&params, u);
let p = q.eval(t);
callback.callback(PathEl::LineTo(p));
// An upper bound on the shortest distance of any point on the quadratic Bezier
// curve to the line segment [p0, p2] is 1/2 of the maximum of the
// endpoint-to-control-point distances.
//
// The derivation is similar to that for the cubic Bezier (see below). In
// short:
//
// q(t) = B0(t) p0 + B1(t) p1 + B2(t) p2
// dist(q(t), [p0, p1]) <= B1(t) dist(p1, [p0, p1])
// = 2 (1-t)t dist(p1, [p0, p1]).
//
// The maximum occurs at t=1/2, hence
// max(dist(q(t), [p0, p1] <= 1/2 dist(p1, [p0, p1])).
//
// A cheap upper bound for dist(p1, [p0, p1]) is max(dist(p1, p0), dist(p1, p2)).
//
// The following takes the square to elide the square root of the Euclidean
// distance.
if f64::max((p1 - p0).hypot2(), (p1 - p2).hypot2()) <= 4. * TOL_2 {
callback.callback(PathEl::LineTo(p2));
} else {
let q = QuadBez::new(p0, p1, p2);
let params = q.estimate_subdiv(sqrt_tol);
let n = ((0.5 * params.val / sqrt_tol).ceil() as usize).max(1);
let step = 1.0 / (n as f64);
for i in 1..n {
let u = (i as f64) * step;
let t = q.determine_subdiv_t(&params, u);
let p = q.eval(t);
callback.callback(PathEl::LineTo(p));
}
callback.callback(PathEl::LineTo(p2));
}
callback.callback(PathEl::LineTo(p2));
}
last_pt = Some(p2);
}
PathEl::CurveTo(p1, p2, p3) => {
if let Some(p0) = last_pt {
let c = CubicBez::new(p0, p1, p2, p3);
let max = simd.vectorize(
#[inline(always)]
|| {
flatten_cubic_simd(
simd,
c,
flatten_ctx,
tolerance as f32,
&mut flattened_cubics,
)
},
);

for p in &flattened_cubics[1..max] {
callback.callback(PathEl::LineTo(Point::new(p.x as f64, p.y as f64)));
// An upper bound on the shortest distance of any point on the cubic Bezier
// curve to the line segment [p0, p3] is 3/4 of the maximum of the
// endpoint-to-control-point distances.
//
// With Bernstein weights Bi(t), we have
// c(t) = B0(t) p0 + B1(t) p1 + B2(t) p2 + B3(t) p3
// with t from 0 to 1 (inclusive).
//
// Through convexivity of the Euclidean distance function and the line segment,
// we have
// dist(c(t), [p0, p3]) <= B1(t) dist(p1, [p0, p3]) + B2(t) dist(p2, [p0, p3])
// <= B1(t) ||p1-p0|| + B2(t) ||p2-p3||
// <= (B1(t) + B2(t)) max(||p1-p0||, ||p2-p3|||)
// = 3 ((1-t)t^2 + (1-t)^2t) max(||p1-p0||, ||p2-p3||).
//
// The inner polynomial has its maximum of 1/4 at t=1/2, hence
// max(dist(c(t), [p0, p3])) <= 3/4 max(||p1-p0||, ||p2-p3||).
//
// The following takes the square to elide the square root of the Euclidean
// distance.
if f64::max((p0 - p1).hypot2(), (p3 - p2).hypot2()) <= 16. / 9. * TOL_2 {
callback.callback(PathEl::LineTo(p3));
} else {
let c = CubicBez::new(p0, p1, p2, p3);
let max = simd.vectorize(
#[inline(always)]
|| {
flatten_cubic_simd(
simd,
c,
flatten_ctx,
tolerance as f32,
&mut flattened_cubics,
)
},
);

for p in &flattened_cubics[1..max] {
callback.callback(PathEl::LineTo(Point::new(p.x as f64, p.y as f64)));
}
}
}
last_pt = Some(p3);
Expand Down
Loading