From ca01e7eb6e2218783a1ed022eba28a50d72cfef4 Mon Sep 17 00:00:00 2001 From: Patrick Walton Date: Mon, 10 Jul 2017 12:16:58 -0700 Subject: [PATCH] =?UTF-8?q?Subdivide=20cubic=20B=C3=A9ziers=20in=20the=20l?= =?UTF-8?q?egalizer=20if=20not=20monotonically=20increasing=20in=20X.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- partitionfinder/src/geometry.rs | 82 ++++++++++++++++++++++++++++---- partitionfinder/src/legalizer.rs | 55 ++++++++++++++++++++- 2 files changed, 127 insertions(+), 10 deletions(-) diff --git a/partitionfinder/src/geometry.rs b/partitionfinder/src/geometry.rs index ce34fa25..b30c572f 100644 --- a/partitionfinder/src/geometry.rs +++ b/partitionfinder/src/geometry.rs @@ -2,9 +2,35 @@ use euclid::approxeq::ApproxEq; use euclid::{Point2D, Vector2D}; +use std::cmp::Ordering; const NEWTON_RAPHSON_ITERATIONS: u8 = 8; +pub(crate) trait ApproxOrdered { + fn approx_ordered(&self) -> bool; +} + +impl ApproxOrdered for [f32] { + fn approx_ordered(&self) -> bool { + let mut last_ordering = Ordering::Equal; + for window in self.windows(2) { + let (last_value, this_value) = (window[0], window[1]); + let this_ordering = if last_value - this_value < -f32::approx_epsilon() { + Ordering::Less + } else if last_value - this_value > f32::approx_epsilon() { + Ordering::Greater + } else { + Ordering::Equal + }; + if last_ordering != Ordering::Equal && this_ordering != last_ordering { + return false + } + last_ordering = this_ordering + } + true + } +} + // https://stackoverflow.com/a/565282 pub fn line_line_crossing_point(a_p0: &Point2D, a_p1: &Point2D, @@ -67,24 +93,35 @@ fn sample_cubic_bezier(t: f32, p0p1p2.lerp(p1p2p3, t) } -fn sample_cubic_bezier_deriv(t: f32, - p0: &Point2D, - p1: &Point2D, - p2: &Point2D, - p3: &Point2D) - -> Vector2D { +pub fn sample_cubic_bezier_deriv(t: f32, + p0: &Point2D, + p1: &Point2D, + p2: &Point2D, + p3: &Point2D) + -> Vector2D { // https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Cubic_B.C3.A9zier_curves // FIXME(pcwalton): Can this be made faster? let tt = 1.0 - t; return (*p1 - *p0) * 3.0 * tt * tt + (*p2 - *p1) * 6.0 * tt * t + (*p3 - *p2) * 3.0 * t * t } +pub fn sample_cubic_bezier_deriv_deriv(t: f32, + p0: &Point2D, + p1: &Point2D, + p2: &Point2D, + p3: &Point2D) + -> Vector2D { + // https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Cubic_B.C3.A9zier_curves + // FIXME(pcwalton): Can this be made faster? + (*p2 - *p1 * 2.0 + p0.to_vector()).lerp(*p3 - *p2 * 2.0 + p1.to_vector(), t) * 6.0 +} + pub fn solve_line_y_for_x(x: f32, a: &Point2D, b: &Point2D) -> f32 { a.lerp(*b, (x - a.x) / (b.x - a.x)).y } -fn newton_raphson(f: F, dfdx: DFDX, mut x_guess: f32) -> f32 - where F: Fn(f32) -> f32, DFDX: Fn(f32) -> f32 { +pub(crate) fn newton_raphson(f: F, dfdx: DFDX, mut x_guess: f32) -> f32 + where F: Fn(f32) -> f32, DFDX: Fn(f32) -> f32 { for _ in 0..NEWTON_RAPHSON_ITERATIONS { let y = f(x_guess); if y.approx_eq(&0.0) { @@ -115,3 +152,32 @@ pub fn solve_cubic_bezier_y_for_x(x: f32, -> f32 { sample_cubic_bezier(solve_cubic_bezier_t_for_x(x, p0, p1, p2, p3), p0, p1, p2, p3).y } + +#[derive(Clone, Copy, Debug)] +pub struct SubdividedCubicBezier { + pub ap0: Point2D, + pub ap1: Point2D, + pub ap2: Point2D, + pub ap3bp0: Point2D, + pub bp1: Point2D, + pub bp2: Point2D, + pub bp3: Point2D, +} + +impl SubdividedCubicBezier { + pub fn new(t: f32, p0: &Point2D, p1: &Point2D, p2: &Point2D, p3: &Point2D) + -> SubdividedCubicBezier { + let (ap1, p1p2, bp2) = (p0.lerp(*p1, t), p1.lerp(*p2, t), p2.lerp(*p3, t)); + let (ap2, bp1) = (ap1.lerp(p1p2, t), (p1p2.lerp(bp2, t))); + let ap3bp0 = ap2.lerp(bp1, t); + SubdividedCubicBezier { + ap0: *p0, + ap1: ap1, + ap2: ap2, + ap3bp0: ap3bp0, + bp1: bp1, + bp2: bp2, + bp3: *p3, + } + } +} diff --git a/partitionfinder/src/legalizer.rs b/partitionfinder/src/legalizer.rs index ce55e479..6a761c62 100644 --- a/partitionfinder/src/legalizer.rs +++ b/partitionfinder/src/legalizer.rs @@ -1,9 +1,12 @@ // partitionfinder/legalizer.rs use euclid::Point2D; +use geometry::{self, ApproxOrdered, SubdividedCubicBezier}; use std::u32; use {ControlPoints, Endpoint, Subpath}; +const MAX_SUBDIVISIONS: u8 = 16; + pub struct Legalizer { endpoints: Vec, control_points: Vec, @@ -73,8 +76,56 @@ impl Legalizer { point1: &Point2D, point2: &Point2D, endpoint: &Point2D) { - // TODO(pcwalton): Make sure curve points are monotonically increasing in X. de Casteljau - // subdivide if not. + self.bezier_curve_to_subdividing_if_necessary(point1, point2, endpoint, 0) + } + + fn bezier_curve_to_subdividing_if_necessary(&mut self, + point1: &Point2D, + point2: &Point2D, + endpoint: &Point2D, + iteration: u8) { + let last_endpoint_index = self.subpaths + .last() + .expect("`bezier_curve_to()` called with no current_subpath") + .last_endpoint_index; + let point0 = self.endpoints[last_endpoint_index as usize].position; + if iteration >= MAX_SUBDIVISIONS || + [point0.x, point1.x, point2.x, endpoint.x].approx_ordered() { + return self.monotone_bezier_curve_to(point1, point2, endpoint) + } + + let t = geometry::newton_raphson(|t| { + geometry::sample_cubic_bezier_deriv(t, + &point0, + point1, + point2, + endpoint).x + }, + |t| { + geometry::sample_cubic_bezier_deriv_deriv(t, + &point0, + point1, + point2, + endpoint).x + }, + 0.5); + + let subdivision = SubdividedCubicBezier::new(t, &point0, point1, point2, endpoint); + + self.bezier_curve_to_subdividing_if_necessary(&subdivision.ap1, + &subdivision.ap2, + &subdivision.ap3bp0, + iteration + 1); + self.bezier_curve_to_subdividing_if_necessary(&subdivision.bp1, + &subdivision.bp2, + &subdivision.bp3, + iteration + 1); + } + + fn monotone_bezier_curve_to(&mut self, + point1: &Point2D, + point2: &Point2D, + endpoint: &Point2D) { self.subpaths .last_mut() .expect("`bezier_curve_to()` called with no current subpath")