michael/pathfinder
michael
/
pathfinder
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Replace the implicitization algorithm for curve intersection with a 

simpler one based on binary search and the quadratic formula.

This does a better job avoiding floating point error and improves the
rendering of the tiger.
This commit is contained in:
Patrick Walton 2017-09-25 15:58:00 -07:00
parent 2a906569e5
commit 93de8383f9
6 changed files with 92 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
partitioner/src/geometry.rs
View File

@ -73,32 +73,6 @@ pub fn solve_line_t_for_x(x: f32, a: &Point2D<f32>, b: &Point2D<f32>) -> 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<f32>,
p1: &Point2D<f32>,
p2: &Point2D<f32>)
-> 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<f32>,
p1: &Point2D<f32>,
p2: &Point2D<f32>)
-> 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<f32> {
let t = (p0 - p1) / (p0 - 2.0 * p1 + p2);
if t > f32::approx_epsilon() && t < 1.0 - f32::approx_epsilon() {

26
partitioner/src/partitioner.rs
View File

@ -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)
}
}
}

27
path-utils/src/curve.rs
View File

@ -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<f32>; 2],
pub control_point: Point2D<f32>,
@ -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<T>(&self, other: &T) -> Option<Point2D<f32>> where T: Side {
pub fn intersect<T>(&self, other: &T) -> Option<Point2D<f32>> where T: Intersect {
<Curve as Intersect>::intersect(self, other)
}
}

105
path-utils/src/intersection.rs
View File

@ -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>) -> f32;
}
pub(crate) trait Intersect {
fn sample(&self, t: f32) -> Point2D<f32>;
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<T>(&self, other: &T) -> Option<Point2D<f32>> 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<T>(&self, other: &T) -> Option<Point2D<f32>> 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>) -> 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>) -> f32 {
fn l(factor: f32, point: &Point2D<f32>, point_i: &Point2D<f32>, point_j: &Point2D<f32>)
-> 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<f32> {
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<f32> {
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)
}
}

16
path-utils/src/lib.rs
View File

@ -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
}
}

5
path-utils/src/line.rs
View File

@ -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<f32>; 2],
}
@ -42,7 +43,7 @@ impl Line {
}
#[inline]
pub fn intersect<T>(&self, other: &T) -> Option<Point2D<f32>> where T: Side {
pub fn intersect<T>(&self, other: &T) -> Option<Point2D<f32>> where T: Intersect {
<Line as Intersect>::intersect(self, other)
}
}