diff --git a/partitioner/src/geometry.rs b/partitioner/src/geometry.rs index 1884a13c..004b79ae 100644 --- a/partitioner/src/geometry.rs +++ b/partitioner/src/geometry.rs @@ -73,32 +73,6 @@ pub fn solve_line_t_for_x(x: f32, a: &Point2D, b: &Point2D) -> f32 { } } -// Use the Citardauq Formula to avoid precision problems. -// -// https://math.stackexchange.com/a/311397 -pub fn solve_quadratic_bezier_t_for_x(x: f32, - p0: &Point2D, - p1: &Point2D, - p2: &Point2D) - -> f32 { - let (p0x, p1x, p2x, x) = (p0.x as f64, p1.x as f64, p2.x as f64, x as f64); - - let a = p0x - 2.0 * p1x + p2x; - let b = -2.0 * p0x + 2.0 * p1x; - let c = p0x - x; - - let t = 2.0 * c / (-b - (b * b - 4.0 * a * c).sqrt()); - t.max(0.0).min(1.0) as f32 -} - -pub fn solve_quadratic_bezier_y_for_x(x: f32, - p0: &Point2D, - p1: &Point2D, - p2: &Point2D) - -> f32 { - Curve::new(p0, p1, p2).sample(solve_quadratic_bezier_t_for_x(x, p0, p1, p2)).y -} - fn quadratic_bezier_axis_inflection_point(p0: f32, p1: f32, p2: f32) -> Option { let t = (p0 - p1) / (p0 - 2.0 * p1 + p2); if t > f32::approx_epsilon() && t < 1.0 - f32::approx_epsilon() { diff --git a/partitioner/src/partitioner.rs b/partitioner/src/partitioner.rs index 5848b5f2..b64f9284 100644 --- a/partitioner/src/partitioner.rs +++ b/partitioner/src/partitioner.rs @@ -303,6 +303,7 @@ impl<'a> Partitioner<'a> { } fn sort_active_edge_list_and_emit_self_intersections(&mut self, endpoint_index: u32) { + let x = self.endpoints[endpoint_index as usize].position.x; loop { let mut swapped = false; for lower_active_edge_index in 1..(self.active_edges.len() as u32) { @@ -310,7 +311,7 @@ impl<'a> Partitioner<'a> { if self.active_edges_are_ordered(upper_active_edge_index, lower_active_edge_index, - endpoint_index) { + x) { continue } @@ -435,7 +436,7 @@ impl<'a> Partitioner<'a> { fn active_edges_are_ordered(&mut self, prev_active_edge_index: u32, next_active_edge_index: u32, - reference_endpoint_index: u32) + x: f32) -> bool { let prev_active_edge = &self.active_edges[prev_active_edge_index as usize]; let next_active_edge = &self.active_edges[next_active_edge_index as usize]; @@ -447,11 +448,10 @@ impl<'a> Partitioner<'a> { // TODO(pcwalton): See if we can speed this up. It's trickier than it seems, due to path // self intersection! - let reference_endpoint = &self.endpoints[reference_endpoint_index as usize]; - let prev_active_edge_y = self.solve_active_edge_y_for_x(reference_endpoint.position.x, - prev_active_edge); - let next_active_edge_y = self.solve_active_edge_y_for_x(reference_endpoint.position.x, - next_active_edge); + let prev_active_edge_t = self.solve_active_edge_t_for_x(x, prev_active_edge); + let next_active_edge_t = self.solve_active_edge_t_for_x(x, next_active_edge); + let prev_active_edge_y = self.sample_active_edge(prev_active_edge_t, prev_active_edge).y; + let next_active_edge_y = self.sample_active_edge(next_active_edge_t, next_active_edge).y; prev_active_edge_y <= next_active_edge_y } @@ -846,10 +846,7 @@ impl<'a> Partitioner<'a> { -> (SubdividedActiveEdge, SubdividedActiveEdge) { let curve = subdivision.to_curve(&self.b_vertex_positions) .expect("subdivide_active_edge_again_at_x(): not a curve!"); - let t = geometry::solve_quadratic_bezier_t_for_x(x, - &curve.endpoints[0], - &curve.control_point, - &curve.endpoints[1]); + let t = curve.solve_t_for_x(x); self.subdivide_active_edge_again_at_t(subdivision, t, bottom) } @@ -922,10 +919,9 @@ impl<'a> Partitioner<'a> { } control_point_vertex_index => { let control_point = &self.b_vertex_positions[control_point_vertex_index as usize]; - geometry::solve_quadratic_bezier_t_for_x(x, - left_vertex_position, - control_point, - right_endpoint_position) + Curve::new(left_vertex_position, + control_point, + right_endpoint_position).solve_t_for_x(x) } } } diff --git a/path-utils/src/curve.rs b/path-utils/src/curve.rs index 4bdfc10c..23c0a472 100644 --- a/path-utils/src/curve.rs +++ b/path-utils/src/curve.rs @@ -15,9 +15,10 @@ use euclid::Point2D; use std::f32; use PathSegment; -use intersection::{Intersect, Side}; +use intersection::Intersect; use line::Line; +#[derive(Clone, Copy, PartialEq, Debug)] pub struct Curve { pub endpoints: [Point2D; 2], pub control_point: Point2D, @@ -62,6 +63,28 @@ impl Curve { (inflection_point_x, inflection_point_y) } + /// Uses the Citardauq Formula to avoid precision problems. + /// + /// https://math.stackexchange.com/a/311397 + pub fn solve_t_for_x(&self, x: f32) -> f32 { + let p0x = self.endpoints[0].x as f64; + let p1x = self.control_point.x as f64; + let p2x = self.endpoints[1].x as f64; + let x = x as f64; + + let a = p0x - 2.0 * p1x + p2x; + let b = -2.0 * p0x + 2.0 * p1x; + let c = p0x - x; + + let t = 2.0 * c / (-b - (b * b - 4.0 * a * c).sqrt()); + t.max(0.0).min(1.0) as f32 + } + + #[inline] + pub fn solve_y_for_x(&self, x: f32) -> f32 { + self.sample(self.solve_t_for_x(x)).y + } + #[inline] pub fn baseline(&self) -> Line { Line::new(&self.endpoints[0], &self.endpoints[1]) @@ -81,7 +104,7 @@ impl Curve { } #[inline] - pub fn intersect(&self, other: &T) -> Option> where T: Side { + pub fn intersect(&self, other: &T) -> Option> where T: Intersect { ::intersect(self, other) } } diff --git a/path-utils/src/intersection.rs b/path-utils/src/intersection.rs index 6e61463d..5a6ba1c5 100644 --- a/path-utils/src/intersection.rs +++ b/path-utils/src/intersection.rs @@ -15,88 +15,73 @@ use euclid::Point2D; use curve::Curve; use line::Line; -use {det2x2, det3x3, lerp}; +use {lerp, sign}; -pub trait Side { - fn side(&self, point: &Point2D) -> f32; -} - -pub(crate) trait Intersect { - fn sample(&self, t: f32) -> Point2D; +pub trait Intersect { + fn min_x(&self) -> f32; + fn max_x(&self) -> f32; + fn solve_y_for_x(&self, t: f32) -> f32; /// Requires that any curves be monotonic. (See the `monotonic` module for that.) /// /// This should work for line segments, but it is inefficient. /// - /// See T.W. Sederberg, "Computer Aided Geometric Design Course Notes" § 17.8. - fn intersect(&self, other: &T) -> Option> where T: Side { - let (mut t_min, mut t_max) = (0.0, 1.0); - let mut iteration = 0; - while t_max - t_min > f32::approx_epsilon() { - let (p_min, p_max) = (self.sample(t_min), self.sample(t_max)); + /// This algorithm used to be smarter (based on implicitization) but floating point error + /// forced the adoption of this simpler, but slower, technique. Improvements are welcome. + fn intersect(&self, other: &T) -> Option> where T: Intersect { + let mut min_x = f32::max(self.min_x(), other.min_x()); + let mut max_x = f32::min(self.max_x(), other.max_x()); - let side_min = other.side(&p_min).signum(); - let side_max = other.side(&p_max).signum(); - if iteration == 0 && side_min == side_max { - return None - } + while max_x - min_x > f32::approx_epsilon() { + let mid_x = lerp(min_x, max_x, 0.5); + let min_sign = sign(self.solve_y_for_x(min_x) - other.solve_y_for_x(min_x)); + let mid_sign = sign(self.solve_y_for_x(mid_x) - other.solve_y_for_x(mid_x)); + let max_sign = sign(self.solve_y_for_x(max_x) - other.solve_y_for_x(max_x)); - let t_mid = lerp(t_min, t_max, 0.5); - let p_mid = self.sample(t_mid); - let side_mid = other.side(&p_mid).signum(); - - if side_mid == side_min { - t_min = t_mid - } else if side_mid == side_max { - t_max = t_mid + if min_sign == mid_sign && mid_sign != max_sign { + min_x = mid_x + } else if min_sign != mid_sign && mid_sign == max_sign { + max_x = min_x } else { break } - - iteration += 1 } - Some(self.sample(lerp(t_min, t_max, 0.5))) - } -} - -impl Side for Line { - #[inline] - fn side(&self, point: &Point2D) -> f32 { - Line::side(self, point) - } -} - -impl Side for Curve { - /// See T.W. Sederberg, "Computer Aided Geometric Design Course Notes" § 17.6.1. - fn side(&self, point: &Point2D) -> f32 { - fn l(factor: f32, point: &Point2D, point_i: &Point2D, point_j: &Point2D) - -> f32 { - factor * det3x3(&[ - point.x, point.y, 1.0, - point_i.x, point_i.y, 1.0, - point_j.x, point_j.y, 1.0, - ]) - } - - let l20 = l(1.0 * 1.0, point, &self.endpoints[1], &self.endpoints[0]); - det2x2(&[ - l(2.0 * 1.0, point, &self.endpoints[1], &self.control_point), l20, - l20, l(2.0 * 1.0, point, &self.control_point, &self.endpoints[0]), - ]) + let mid_x = lerp(min_x, max_x, 0.5); + Some(Point2D::new(mid_x, self.solve_y_for_x(mid_x))) } } impl Intersect for Line { #[inline] - fn sample(&self, t: f32) -> Point2D { - Line::sample(self, t) + fn min_x(&self) -> f32 { + f32::min(self.endpoints[0].x, self.endpoints[1].x) + } + + #[inline] + fn max_x(&self) -> f32 { + f32::max(self.endpoints[0].x, self.endpoints[1].x) + } + + #[inline] + fn solve_y_for_x(&self, x: f32) -> f32 { + self.sample((x - self.endpoints[0].x) / (self.endpoints[1].x - self.endpoints[0].x)).y } } impl Intersect for Curve { #[inline] - fn sample(&self, t: f32) -> Point2D { - Curve::sample(self, t) + fn min_x(&self) -> f32 { + f32::min(self.endpoints[0].x, self.endpoints[1].x) + } + + #[inline] + fn max_x(&self) -> f32 { + f32::max(self.endpoints[0].x, self.endpoints[1].x) + } + + #[inline] + fn solve_y_for_x(&self, x: f32) -> f32 { + Curve::solve_y_for_x(self, x) } } diff --git a/path-utils/src/lib.rs b/path-utils/src/lib.rs index 8c1eac00..e2f75d9b 100644 --- a/path-utils/src/lib.rs +++ b/path-utils/src/lib.rs @@ -208,12 +208,12 @@ pub(crate) fn lerp(a: f32, b: f32, t: f32) -> f32 { } #[inline] -pub(crate) fn det2x2(m: &[f32; 4]) -> f32 { - m[0] * m[3] - m[1] * m[2] -} - -pub(crate) fn det3x3(m: &[f32; 9]) -> f32 { - m[0] * det2x2(&[m[4], m[5], m[7], m[8]]) - - m[1] * det2x2(&[m[3], m[5], m[6], m[8]]) + - m[2] * det2x2(&[m[3], m[4], m[6], m[7]]) +pub(crate) fn sign(x: f32) -> f32 { + if x < 0.0 { + -1.0 + } else if x > 0.0 { + 1.0 + } else { + x + } } diff --git a/path-utils/src/line.rs b/path-utils/src/line.rs index a8213f83..70dbe066 100644 --- a/path-utils/src/line.rs +++ b/path-utils/src/line.rs @@ -12,8 +12,9 @@ use euclid::{Point2D, Vector2D}; -use intersection::{Intersect, Side}; +use intersection::Intersect; +#[derive(Clone, Copy, PartialEq, Debug)] pub struct Line { pub endpoints: [Point2D; 2], } @@ -42,7 +43,7 @@ impl Line { } #[inline] - pub fn intersect(&self, other: &T) -> Option> where T: Side { + pub fn intersect(&self, other: &T) -> Option> where T: Intersect { ::intersect(self, other) } }