From 6e0ec0f5fd57b00c56202ec8d0d42981d55e555f Mon Sep 17 00:00:00 2001 From: Simon Cozens Date: Thu, 16 Sep 2021 09:49:33 +0100 Subject: [PATCH 01/20] =?UTF-8?q?Borrow=20Lyon=E2=80=99s=20implementation?= =?UTF-8?q?=20of=20curve-curve=20intersection?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/common.rs | 10 + src/cubicbez.rs | 134 ++++++++- src/curve_intersections.rs | 603 +++++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + src/param_curve.rs | 18 ++ 5 files changed, 763 insertions(+), 3 deletions(-) create mode 100644 src/curve_intersections.rs diff --git a/src/common.rs b/src/common.rs index 3788a5fd..6a78a2a1 100644 --- a/src/common.rs +++ b/src/common.rs @@ -43,6 +43,16 @@ impl FloatExt for f32 { } } +/// Order two things into minimum and maximum +#[inline] +pub fn min_max(a: T, b: T) -> (T, T) { + if a < b { + (a, b) + } else { + (b, a) + } +} + /// Find real roots of cubic equation. /// /// The implementation is not (yet) fully robust, but it does handle the case diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 104996eb..4ddc27fa 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -6,11 +6,12 @@ use crate::MAX_EXTREMA; use crate::{Line, QuadSpline, Vec2}; use arrayvec::ArrayVec; -use crate::common::solve_quadratic; use crate::common::GAUSS_LEGENDRE_COEFFS_9; +use crate::common::{min_max, solve_cubic, solve_quadratic}; use crate::{ - Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature, - ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, QuadBez, Rect, Shape, + Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, + ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, + QuadBez, Rect, Shape, }; const MAX_SPLINE_SPLIT: usize = 100; @@ -305,6 +306,13 @@ impl CubicBez { pub fn is_nan(&self) -> bool { self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan() || self.p3.is_nan() } + + /// Is this cubic Bezier curve a line? + #[inline] + fn is_linear(&self, accuracy: f64) -> bool { + self.baseline().nearest(self.p1, accuracy).distance_sq <= accuracy + && self.baseline().nearest(self.p2, accuracy).distance_sq <= accuracy + } } /// An iterator for cubic beziers. @@ -527,6 +535,101 @@ impl ParamCurveExtrema for CubicBez { } } +// Note that the line is unbounded here! +fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { + let vec2 = l.p1 - l.p0; + let a = -vec2.y; + let b = vec2.x; + let c = -(a * l.start().x + b * l.start().y); + a * p.x + b * p.y + c +} + +impl ParamCurveBezierClipping for CubicBez { + fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b, c, d) = self.parameters(); + solve_cubic(d.x - x, c.x, b.x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b, c, d) = self.parameters(); + solve_cubic(d.y - y, c.y, b.y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self) -> (f64, f64) { + let baseline = self.baseline(); + let (d1, d2) = min_max( + signed_distance_from_ray_to_point(&baseline, self.p1), + signed_distance_from_ray_to_point(&baseline, self.p2), + ); + let factor = if (d1 * d2) > 0.0 { + 3.0 / 4.0 + } else { + 4.0 / 9.0 + }; + + let d_min = factor * f64::min(d1, 0.0); + let d_max = factor * f64::max(d2, 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.p1); + let d2 = signed_distance_from_ray_to_point(l, self.p2); + let d3 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0 / 3.0, d1); + let p2 = Point::new(2.0 / 3.0, d2); + let p3 = Point::new(1.0, d3); + // Compute the vertical signed distance of p1 and p2 from [p0, p3]. + let dist1 = d1 - (2.0 * d0 + d3) / 3.0; + let dist2 = d2 - (d0 + 2.0 * d3) / 3.0; + + // Compute the hull assuming p1 is on top - we'll switch later if needed. + let mut hull = if dist1 * dist2 < 0.0 { + // p1 and p2 lie on opposite sides of [p0, p3], so the hull is a quadrilateral: + (vec![p0, p1, p3], vec![p0, p2, p3]) + } else { + // p1 and p2 lie on the same side of [p0, p3]. The hull can be a triangle or a + // quadrilateral, and [p0, p3] is part of the hull. The hull is a triangle if the vertical + // distance of one of the middle points p1, p2 is <= half the vertical distance of the + // other middle point. + let dist1 = dist1.abs(); + let dist2 = dist2.abs(); + if dist1 >= 2.0 * dist2 { + (vec![p0, p1, p3], vec![p0, p3]) + } else if dist2 >= 2.0 * dist1 { + (vec![p0, p2, p3], vec![p0, p3]) + } else { + (vec![p0, p1, p2, p3], vec![p0, p3]) + } + }; + + // Flip the hull if needed: + if dist1 < 0.0 || (dist1 == 0.0 && dist2 < 0.0) { + hull = (hull.1, hull.0); + } + + hull + } +} + impl Mul for Affine { type Output = CubicBez; @@ -594,6 +697,7 @@ mod tests { cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez, }; + use arrayvec::{Array, ArrayVec}; #[test] fn cubicbez_deriv() { @@ -873,4 +977,28 @@ mod tests { converted[0].points()[2].distance(Point::new(88639.0 / 90.0, 52584.0 / 90.0)) < 0.0001 ); } + + use crate::param_curve::ParamCurveBezierClipping; + #[test] + fn solve_t_for_xy() { + fn verify>(mut roots: ArrayVec, expected: &[f64]) { + assert_eq!(expected.len(), roots.len()); + let epsilon = 1e-6; + roots.sort_by(|a, b| a.partial_cmp(b).unwrap()); + + for i in 0..expected.len() { + assert!((roots[i] - expected[i]).abs() < epsilon); + } + } + + let curve = CubicBez::new((0.0, 0.0), (0.0, 8.0), (10.0, 8.0), (10.0, 0.0)); + verify(curve.solve_t_for_x(5.0), &[0.5]); + verify(curve.solve_t_for_y(6.0), &[0.5]); + + { + let curve = CubicBez::new((0.0, 10.0), (0.0, 10.0), (10.0, 10.0), (10.0, 10.0)); + + verify(curve.solve_t_for_y(10.0), &[]); + } + } } diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs new file mode 100644 index 00000000..6dd9c6b1 --- /dev/null +++ b/src/curve_intersections.rs @@ -0,0 +1,603 @@ +use crate::ParamCurveBezierClipping; +use crate::{ParamCurve, ParamCurveExtrema}; +use crate::{Point, Rect}; +use arrayvec::ArrayVec; +use std::ops::Range; + +pub fn curve_curve_intersections( + curve1: &T, + curve2: &U, +) -> ArrayVec<[(f64, f64); 9]> { + let mut av = ArrayVec::new(); + add_curve_intersections( + curve1, + curve2, + &(0.0..1.0), + &(0.0..1.0), + &mut av, + false, + 0, + 0, + curve1, + curve2, + ); + av +} + +// This function implements the main bézier clipping algorithm by recursively subdividing curve1 and +// curve2 in to smaller and smaller portions of the original curves with the property that one of +// the curves intersects the fat line of the other curve at each stage. +// +// curve1 and curve2 at each stage are sub-bézier curves of the original curves; flip tells us +// whether curve1 at a given stage is a subcurve of the original curve1 or the original curve2; +// similarly for curve2. domain1 and domain2 shrink (or stay the same) at each stage and describe +// which subdomain of an original curve the current curve1 and curve2 correspond to. (The domains of +// curve1 and curve2 are 0..1 at every stage.) +#[allow(clippy::too_many_arguments)] +fn add_curve_intersections( + curve1: &T, + curve2: &U, + domain1: &Range, + domain2: &Range, + intersections: &mut ArrayVec<[(f64, f64); 9]>, + flip: bool, + mut recursion_count: u32, + mut call_count: u32, + orig_curve1: &T, + orig_curve2: &U, +) -> u32 { + call_count += 1; + recursion_count += 1; + if call_count >= 4096 || recursion_count >= 60 { + return call_count; + } + + let epsilon = 1e-9; + + if (domain2.start - domain2.end).abs() < f64::EPSILON || curve2.arclen(epsilon) == 0.0 { + add_point_curve_intersection( + curve2, + /* point is curve1 */ false, + curve1, + domain2, + domain1, + intersections, + flip, + ); + return call_count; + } else if curve2.start() == curve2.end() { + // There's no curve2 baseline to fat-line against (and we'll (debug) crash if we try with + // the current implementation), so split curve2 and try again. + let new_2_curves = orig_curve2.subsegment(domain2.clone()).subdivide(); + let domain2_mid = (domain2.start + domain2.end) * 0.5; + call_count = add_curve_intersections( + curve1, + &new_2_curves.0, + domain1, + &(domain2.start..domain2_mid), + intersections, + flip, + recursion_count, + call_count, + orig_curve1, + orig_curve2, + ); + call_count = add_curve_intersections( + curve1, + &new_2_curves.1, + domain1, + &(domain2_mid..domain2.end), + intersections, + flip, + recursion_count, + call_count, + orig_curve1, + orig_curve2, + ); + return call_count; + } + + // (Don't call this before checking for point curves: points are inexact and can lead to false + // negatives here.) + if !rectangles_overlap(&curve1.bounding_box(), &curve2.bounding_box()) { + return call_count; + } + + let (t_min_clip, t_max_clip) = match restrict_curve_to_fat_line(curve1, curve2) { + Some((min, max)) => (min, max), + None => return call_count, + }; + + // t_min_clip and t_max_clip are (0, 1)-based, so project them back to get the new restricted + // range: + let new_domain1 = + &(domain_value_at_t(&domain1, t_min_clip)..domain_value_at_t(&domain1, t_max_clip)); + + if (domain2.end - domain2.start).max(new_domain1.end - new_domain1.start) < epsilon { + let t1 = (new_domain1.start + new_domain1.end) * 0.5; + let t2 = (domain2.start + domain2.end) * 0.5; + // There's an unfortunate tendency for curve2 endpoints that end near (but not all + // that near) to the interior of curve1 to register as intersections, so try to avoid + // that. (We could be discarding a legitimate intersection here.) + let end_eps = 1e-3; + if (t2 < end_eps || t2 > 1.0 - end_eps) + && (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot() > 5.0 + { + return call_count; + } + add_intersection(t1, orig_curve1, t2, orig_curve2, flip, intersections); + return call_count; + } + + // Reduce curve1 to the part that might intersect curve2. + let curve1 = &orig_curve1.subsegment(new_domain1.clone()); + + // (Note: it's possible for new_domain1 to have become a point, even if + // t_min_clip < t_max_clip. It's also possible for curve1 to not be a point even if new_domain1 + // is a point (but then curve1 will be very small).) + if (new_domain1.start - new_domain1.end).abs() < f64::EPSILON || curve1.arclen(epsilon) == 0.0 { + add_point_curve_intersection( + curve1, + /* point is curve1 */ true, + curve2, + new_domain1, + domain2, + intersections, + flip, + ); + return call_count; + } + + // If the new range is still 80% or more of the old range, subdivide and try again. + if t_max_clip - t_min_clip > 0.8 { + // Subdivide the curve which has converged the least. + if new_domain1.end - new_domain1.start > domain2.end - domain2.start { + let new_1_curves = curve1.subdivide(); + let new_domain1_mid = (new_domain1.start + new_domain1.end) * 0.5; + call_count = add_curve_intersections( + curve2, + &new_1_curves.0, + domain2, + &(new_domain1.start..new_domain1_mid), + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + ); + call_count = add_curve_intersections( + curve2, + &new_1_curves.1, + domain2, + &(new_domain1_mid..new_domain1.end), + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + ); + } else { + let new_2_curves = orig_curve2.subsegment(domain2.clone()).subdivide(); + let domain2_mid = (domain2.start + domain2.end) * 0.5; + call_count = add_curve_intersections( + &new_2_curves.0, + curve1, + &(domain2.start..domain2_mid), + new_domain1, + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + ); + call_count = add_curve_intersections( + &new_2_curves.1, + curve1, + &(domain2_mid..domain2.end), + new_domain1, + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + ); + } + } else { + // Iterate. + if domain2.end - domain2.start >= epsilon { + call_count = add_curve_intersections( + curve2, + curve1, + domain2, + new_domain1, + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + ); + } else { + // The interval on curve2 is already tight enough, so just continue iterating on curve1. + call_count = add_curve_intersections( + curve1, + curve2, + new_domain1, + domain2, + intersections, + flip, + recursion_count, + call_count, + orig_curve1, + orig_curve2, + ); + } + } + + call_count +} + +fn add_point_curve_intersection( + pt_curve: &T, + pt_curve_is_curve1: bool, + curve: &U, + pt_domain: &Range, + curve_domain: &Range, + intersections: &mut ArrayVec<[(f64, f64); 9]>, + flip: bool, +) { + let pt = pt_curve.start(); + // We assume pt is curve1 when we add intersections below. + let flip = if pt_curve_is_curve1 { flip } else { !flip }; + + // Generally speaking |curve| will be quite small at this point, so see if we can get away with + // just sampling here. + + let epsilon = epsilon_for_point(pt); + let pt_t = (pt_domain.start + pt_domain.end) * 0.5; + + let curve_t = { + let mut t_for_min = 0.0; + let mut min_dist_sq = epsilon; + let tenths = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; + for &t in tenths.iter() { + let d = (pt - curve.eval(t)).hypot2(); + if d < min_dist_sq { + t_for_min = t; + min_dist_sq = d; + } + } + + if (min_dist_sq - epsilon).abs() < f64::EPSILON { + -1.0 + } else { + domain_value_at_t(curve_domain, t_for_min) + } + }; + + if (curve_t - -1.0).abs() > f64::EPSILON { + add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections); + return; + } + + // If sampling didn't work, try a different approach. + let results = point_curve_intersections(pt, curve, epsilon); + for t in results { + let curve_t = domain_value_at_t(curve_domain, t); + add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections); + } +} + +fn point_curve_intersections( + pt: Point, + curve: &T, + epsilon: f64, +) -> ArrayVec<[f64; 9]> { + let mut result = ArrayVec::new(); + + // (If both endpoints are epsilon close, we only return 0.0.) + if (pt - curve.start()).hypot2() < epsilon { + result.push(0.0); + return result; + } + if (pt - curve.end()).hypot2() < epsilon { + result.push(1.0); + return result; + } + + let curve_x_t_params = curve.solve_t_for_x(pt.x); + let curve_y_t_params = curve.solve_t_for_y(pt.y); + // We want to coalesce parameters representing the same intersection from the x and y + // directions, but the parameter calculations aren't very accurate, so give a little more + // leeway there (TODO: this isn't perfect, as you might expect - the dupes that pass here are + // currently being detected in add_intersection). + let param_eps = 10.0 * epsilon; + for params in [curve_x_t_params, curve_y_t_params].iter() { + for t in params { + let t = *t; + if (pt - curve.eval(t)).hypot2() > epsilon { + continue; + } + let mut already_found_t = false; + for u in &result { + if f64::abs(t - *u) < param_eps { + already_found_t = true; + break; + } + } + if !already_found_t { + result.push(t); + } + } + } + + if !result.is_empty() { + return result; + } + + // The remaining case is if pt is within epsilon of an interior point of curve, but not within + // the x-range or y-range of the curve (which we already checked) - for example if curve is a + // horizontal line that extends beyond its endpoints, and pt is just outside an end of the line; + // or if the curve has a cusp in one of the corners of its convex hull and pt is + // diagonally just outside the hull. This is a rare case (could we even ignore it?). + #[inline] + fn maybe_add( + t: f64, + pt: Point, + curve: &T, + epsilon: f64, + result: &mut ArrayVec<[f64; 9]>, + ) -> bool { + if (curve.eval(t) - pt).hypot2() < epsilon { + result.push(t); + return true; + } + false + } + + for ex in curve.extrema() { + maybe_add(ex, pt, curve, epsilon, &mut result); + } + + result +} + +// If we're comparing distances between samples of curves, our epsilon should depend on how big the +// points we're comparing are. This function returns an epsilon appropriate for the size of pt. +fn epsilon_for_point(pt: Point) -> f64 { + let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); + if max <= 9.0 { + 0.001 + } else if max <= 99.0 { + 0.01 + } else if max <= 999.0 { + 0.1 + } else if max <= 9_999.0 { + 0.25 + } else if max <= 999_999.0 { + 0.5 + } else { + 1.0 + } +} + +fn add_intersection( + t1: f64, + orig_curve1: &T, + t2: f64, + orig_curve2: &U, + flip: bool, + intersections: &mut ArrayVec<[(f64, f64); 9]>, +) { + let (t1, t2) = if flip { (t2, t1) } else { (t1, t2) }; + // (This should probably depend in some way on how large our input coefficients are.) + let epsilon = 1e-3; + // Discard endpoint/endpoint intersections. + let t1_is_an_endpoint = t1 < epsilon || t1 > 1.0 - epsilon; + let t2_is_an_endpoint = t2 < epsilon || t2 > 1.0 - epsilon; + if t1_is_an_endpoint && t2_is_an_endpoint { + return; + } + + // We can get repeated intersections when we split a curve at an intersection point, or when + // two curves intersect at a point where the curves are very close together, or when the fat + // line process breaks down. + for i in 0..intersections.len() { + let (old_t1, old_t2) = intersections[i]; + // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" + // duplicate intersection, so always keep the one that minimizes sample distance. + if (t1 - old_t1).abs() < epsilon && (t2 - old_t2).abs() < epsilon { + let cur_dist = (orig_curve1.eval(old_t1) - orig_curve2.eval(old_t2)).hypot2(); + let new_dist = (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot2(); + if new_dist < cur_dist { + intersections[i] = (t1, t2); + } + return; + } + } + + if intersections.len() < 9 { + intersections.push((t1, t2)); + } +} + +// Returns an interval (t_min, t_max) with the property that for parameter values outside that +// interval, curve1 is guaranteed to not intersect curve2; uses the fat line of curve2 as its basis +// for the guarantee. (See the Sederberg document for what's going on here.) +fn restrict_curve_to_fat_line< + T: ParamCurveBezierClipping + ParamCurve + ParamCurveExtrema, + U: ParamCurveBezierClipping + ParamCurve + ParamCurveExtrema, +>( + curve1: &T, + curve2: &U, +) -> Option<(f64, f64)> { + // TODO: Consider clipping against the perpendicular fat line as well (recommended by + // Sederberg). + // TODO: The current algorithm doesn't handle the (rare) case where curve1 and curve2 are + // overlapping lines. + + let baseline2 = curve2.baseline(); + let (mut top, mut bottom) = curve1.convex_hull_from_line(&baseline2); + let (d_min, d_max) = curve2.fat_line_min_max(); + + clip_convex_hull_to_fat_line(&mut top, &mut bottom, d_min, d_max) +} + +// Returns the min and max values at which the convex hull enters the fat line min/max offset lines. +fn clip_convex_hull_to_fat_line( + hull_top: &mut Vec, + hull_bottom: &mut Vec, + d_min: f64, + d_max: f64, +) -> Option<(f64, f64)> { + // Walk from the left corner of the convex hull until we enter the fat line limits: + let t_clip_min = walk_convex_hull_start_to_fat_line(&hull_top, &hull_bottom, d_min, d_max)?; + + // Now walk from the right corner of the convex hull until we enter the fat line limits - to + // walk right to left we just reverse the order of the hull vertices, so that hull_top and + // hull_bottom start at the right corner now: + hull_top.reverse(); + hull_bottom.reverse(); + let t_clip_max = walk_convex_hull_start_to_fat_line(&hull_top, &hull_bottom, d_min, d_max)?; + + Some((t_clip_min, t_clip_max)) +} + +// Walk the edges of the convex hull until you hit a fat line offset value, starting from the +// (first vertex in hull_top_vertices == first vertex in hull_bottom_vertices). +fn walk_convex_hull_start_to_fat_line( + hull_top_vertices: &[Point], + hull_bottom_vertices: &[Point], + d_min: f64, + d_max: f64, +) -> Option { + let start_corner = hull_top_vertices[0]; + + if start_corner.y < d_min { + walk_convex_hull_edges_to_fat_line(hull_top_vertices, true, d_min) + } else if start_corner.y > d_max { + walk_convex_hull_edges_to_fat_line(hull_bottom_vertices, false, d_max) + } else { + Some(start_corner.x) + } +} + +// Do the actual walking, starting from the first vertex of hull_vertices. +fn walk_convex_hull_edges_to_fat_line( + hull_vertices: &[Point], + vertices_are_for_top: bool, + threshold: f64, +) -> Option { + for i in 0..hull_vertices.len() - 1 { + let p = hull_vertices[i]; + let q = hull_vertices[i + 1]; + if (vertices_are_for_top && q.y >= threshold) || (!vertices_are_for_top && q.y <= threshold) + { + return if (q.y - threshold).abs() < f64::EPSILON { + Some(q.x) + } else { + Some(p.x + (threshold - p.y) * (q.x - p.x) / (q.y - p.y)) + }; + } + } + // All points of the hull are outside the threshold: + None +} + +#[inline] +// Return the point of domain corresponding to the point t, 0 <= t <= 1. +fn domain_value_at_t(domain: &Range, t: f64) -> f64 { + domain.start + (domain.end - domain.start) * t +} + +#[inline] +// Rect.intersects doesn't count edge/corner intersections, this version does. +fn rectangles_overlap(r1: &Rect, r2: &Rect) -> bool { + r1.origin().x <= r2.origin().x + r2.size().width + && r2.origin().x <= r1.origin().x + r1.size().width + && r1.origin().y <= r2.origin().y + r2.size().height + && r2.origin().y <= r1.origin().y + r1.size().height +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::CubicBez; + + fn do_test( + curve1: &T, + curve2: &U, + count: usize, + ) { + let arr1 = curve_curve_intersections(curve1, curve2); + assert_eq!(arr1.len(), count); + let arr2 = curve_curve_intersections(curve2, curve1); + assert_eq!(arr2.len(), count); + } + + #[test] + fn test_cubic_cubic_intersections() { + do_test( + &CubicBez::new((0.0, 0.0), (0.0, 1.0), (0.0, 1.0), (1.0, 1.0)), + &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), + 1, + ); + do_test( + &CubicBez::new((48.0, 84.0), (104.0, 176.0), (190.0, 37.0), (121.0, 75.0)), + &CubicBez::new((68.0, 145.0), (74.0, 6.0), (143.0, 197.0), (138.0, 55.0)), + 4, + ); + do_test( + &CubicBez::new((0.0, 0.0), (0.5, 1.0), (0.5, 1.0), (1.0, 0.0)), + &CubicBez::new((0.0, 1.0), (0.5, 0.0), (0.5, 0.0), (1.0, 1.0)), + 2, + ); + do_test( + &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), + &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), + 9, + ); + + // (A previous version of the code was returning two practically identical + // intersection points here.) + do_test( + &CubicBez::new( + (718133.1363092018, 673674.987999388), + (-53014.13135835016, 286988.87959900266), + (-900630.1880107201, -7527.6889376943), + (417822.48349384824, -149039.14932848653), + ), + &CubicBez::new( + (924715.3309247112, 719414.5221912428), + (965365.9679664494, -563421.3040676294), + (273552.85484064696, 643090.0890117711), + (-113963.134524995, 732017.9466050486), + ), + 1, + ); + + // On these curves the algorithm runs to a state at which the new clipped domain1 becomes a + // point even though t_min_clip < t_max_clip (because domain1 was small enough to begin with + // relative to the small distance between t_min_clip and t_max_clip), and the new curve1 is not + // a point (it's split off the old curve1 using t_min_clip < t_max_clip). + do_test( + &CubicBez::new( + (423394.5967598548, -91342.7434613118), + (333212.450870987, 225564.45711810607), + (668108.668469816, -626100.8367380127), + (-481885.0610437216, 893767.5320803947), + ), + &CubicBez::new( + (-484505.2601961801, -222621.44229855016), + (22432.829984141514, -944727.7102144773), + (-433294.66549074976, -168018.60431004688), + (567688.5977972192, 13975.09633399453), + ), + 3, + ); + } +} diff --git a/src/lib.rs b/src/lib.rs index d164e31f..a8dd7ba2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -84,6 +84,7 @@ mod bezpath; mod circle; pub mod common; mod cubicbez; +mod curve_intersections; mod ellipse; mod insets; mod line; diff --git a/src/param_curve.rs b/src/param_curve.rs index 7bed03e7..bf06d6cd 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -1,5 +1,6 @@ //! A trait for curves parametrized by a scalar. +use crate::Line; use std::ops::Range; use arrayvec::ArrayVec; @@ -213,3 +214,20 @@ pub trait ParamCurveExtrema: ParamCurve { bbox } } + +/// A parameterized curve that can be used in the Bezier clipping algorithm +pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveExtrema + ParamCurveArclen { + /// Returns a line from the curve's start point to its end point + fn baseline(&self) -> Line { + Line::new(self.start(), self.end()) + } + + /// Find the time `t` at which the curve has the given x value + fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]>; + /// Find the time `t` at which the curve has the given x value + fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]>; + /// Returns the upper and lower convex hull + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); + /// Returns the minimum and maximum distances of the "fat line" enclosing this curve + fn fat_line_min_max(&self) -> (f64, f64); +} From f52accd5afd8df483d3e363cb4457773ba56c547 Mon Sep 17 00:00:00 2001 From: Simon Cozens Date: Thu, 16 Sep 2021 10:15:26 +0100 Subject: [PATCH 02/20] Make clippy less sad --- src/curve_intersections.rs | 7 ++++--- src/lib.rs | 1 + 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 6dd9c6b1..10c65e61 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -4,6 +4,7 @@ use crate::{Point, Rect}; use arrayvec::ArrayVec; use std::ops::Range; +/// Compute the intersections between two Bézier curves pub fn curve_curve_intersections( curve1: &T, curve2: &U, @@ -111,7 +112,7 @@ fn add_curve_intersections Option<(f64, f64)> { // Walk from the left corner of the convex hull until we enter the fat line limits: - let t_clip_min = walk_convex_hull_start_to_fat_line(&hull_top, &hull_bottom, d_min, d_max)?; + let t_clip_min = walk_convex_hull_start_to_fat_line(hull_top, hull_bottom, d_min, d_max)?; // Now walk from the right corner of the convex hull until we enter the fat line limits - to // walk right to left we just reverse the order of the hull vertices, so that hull_top and // hull_bottom start at the right corner now: hull_top.reverse(); hull_bottom.reverse(); - let t_clip_max = walk_convex_hull_start_to_fat_line(&hull_top, &hull_bottom, d_min, d_max)?; + let t_clip_max = walk_convex_hull_start_to_fat_line(hull_top, hull_bottom, d_min, d_max)?; Some((t_clip_min, t_clip_max)) } diff --git a/src/lib.rs b/src/lib.rs index a8dd7ba2..e2111e1f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -107,6 +107,7 @@ pub use crate::arc::*; pub use crate::bezpath::*; pub use crate::circle::*; pub use crate::cubicbez::*; +pub use crate::curve_intersections::curve_curve_intersections; pub use crate::ellipse::*; pub use crate::insets::*; pub use crate::line::*; From 4aeaeeba28b5716311d6fe970feeefb6be80102f Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Tue, 1 Feb 2022 08:10:23 -0700 Subject: [PATCH 03/20] Squashed in curve-curve-intersections-point-along: * Added Epilog Laser/myself to the authors file * Added tests for real.rs and fixed bug * Renamed function for clarity, made public * Reworded comment for clarity * Fixed testing for curves that have multiple points at a t-value * Made testing multiple t values along curve easier * Fixed last commit * Added crate for comparing real values --- AUTHORS | 1 + src/curve_intersections.rs | 112 +++++++++++++++++++++------------- src/lib.rs | 3 +- src/real.rs | 120 +++++++++++++++++++++++++++++++++++++ 4 files changed, 193 insertions(+), 43 deletions(-) create mode 100644 src/real.rs diff --git a/AUTHORS b/AUTHORS index b4fa084d..08eca596 100644 --- a/AUTHORS +++ b/AUTHORS @@ -5,3 +5,4 @@ # of contributors, see the revision history in source control. Raph Levien Nicolas Silva +Epilog Laser (Tyler Scott) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 10c65e61..966ed64b 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,6 +1,7 @@ use crate::ParamCurveBezierClipping; use crate::{ParamCurve, ParamCurveExtrema}; use crate::{Point, Rect}; +use crate::real::*; use arrayvec::ArrayVec; use std::ops::Range; @@ -55,7 +56,7 @@ fn add_curve_intersections= epsilon { + if !real_is_equal(domain2.end, domain2.start) { call_count = add_curve_intersections( curve2, curve1, @@ -251,7 +252,7 @@ fn add_point_curve_intersection, flip: bool, ) { - let pt = pt_curve.start(); + let pt = pt_curve.eval(0.5); // We assume pt is curve1 when we add intersections below. let flip = if pt_curve_is_curve1 { flip } else { !flip }; @@ -273,39 +274,38 @@ fn add_point_curve_intersection f64::EPSILON { + if !real_is_equal(curve_t, -1.0) { add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections); return; } // If sampling didn't work, try a different approach. - let results = point_curve_intersections(pt, curve, epsilon); + let results = t_along_curve_for_point(pt, curve); for t in results { let curve_t = domain_value_at_t(curve_domain, t); add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections); } } -fn point_curve_intersections( +pub fn t_along_curve_for_point( pt: Point, curve: &T, - epsilon: f64, ) -> ArrayVec<[f64; 9]> { let mut result = ArrayVec::new(); - // (If both endpoints are epsilon close, we only return 0.0.) - if (pt - curve.start()).hypot2() < epsilon { + // If both endpoints are approximately close, we only return 0.0. + if point_is_equal(pt, curve.start()) { result.push(0.0); return result; } - if (pt - curve.end()).hypot2() < epsilon { + if point_is_equal(pt, curve.end()) { result.push(1.0); return result; } @@ -316,11 +316,11 @@ fn point_curve_intersections( // directions, but the parameter calculations aren't very accurate, so give a little more // leeway there (TODO: this isn't perfect, as you might expect - the dupes that pass here are // currently being detected in add_intersection). - let param_eps = 10.0 * epsilon; + let param_eps = 10.0 * epsilon_for_point(pt); for params in [curve_x_t_params, curve_y_t_params].iter() { for t in params { let t = *t; - if (pt - curve.eval(t)).hypot2() > epsilon { + if !point_is_equal(pt, curve.eval(t)) { continue; } let mut already_found_t = false; @@ -340,20 +340,21 @@ fn point_curve_intersections( return result; } - // The remaining case is if pt is within epsilon of an interior point of curve, but not within - // the x-range or y-range of the curve (which we already checked) - for example if curve is a - // horizontal line that extends beyond its endpoints, and pt is just outside an end of the line; - // or if the curve has a cusp in one of the corners of its convex hull and pt is - // diagonally just outside the hull. This is a rare case (could we even ignore it?). + // The remaining case is if pt is approximately equal to an interior point + // of curve, but not within the x-range or y-range of the curve (which we + // already checked) due to floating point errors - for example if curve is + // a horizontal line that extends beyond its endpoints, and pt is on the + // extrema, but just barely outside the x-y limits; or if the curve has a + // cusp in one of the corners of its convex hull and pt is diagonally just + // outside the hull. #[inline] fn maybe_add( t: f64, pt: Point, curve: &T, - epsilon: f64, result: &mut ArrayVec<[f64; 9]>, ) -> bool { - if (curve.eval(t) - pt).hypot2() < epsilon { + if point_is_equal(curve.eval(t), pt) { result.push(t); return true; } @@ -361,31 +362,12 @@ fn point_curve_intersections( } for ex in curve.extrema() { - maybe_add(ex, pt, curve, epsilon, &mut result); + maybe_add(ex, pt, curve, &mut result); } result } -// If we're comparing distances between samples of curves, our epsilon should depend on how big the -// points we're comparing are. This function returns an epsilon appropriate for the size of pt. -fn epsilon_for_point(pt: Point) -> f64 { - let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); - if max <= 9.0 { - 0.001 - } else if max <= 99.0 { - 0.01 - } else if max <= 999.0 { - 0.1 - } else if max <= 9_999.0 { - 0.25 - } else if max <= 999_999.0 { - 0.5 - } else { - 1.0 - } -} - fn add_intersection( t1: f64, orig_curve1: &T, @@ -535,13 +517,48 @@ mod tests { count: usize, ) { let arr1 = curve_curve_intersections(curve1, curve2); - assert_eq!(arr1.len(), count); let arr2 = curve_curve_intersections(curve2, curve1); + assert_eq!(arr1.len(), count); assert_eq!(arr2.len(), count); } + fn test_t_along_curve( + curve: &T, + t_test: f64, + ) { + let pt_test = curve.eval(t_test); + let t_values = t_along_curve_for_point(pt_test, curve); + assert!(!t_values.is_empty()); + + let mut found_t = false; + for t in &t_values { + let pt = curve.eval(*t); + + if real_is_equal(*t, t_test) + { + found_t = true; + } + + println!("t for pt on curve: {} {}, expected {} {}", t, pt, t_test, pt_test); + assert!(point_is_equal(pt, pt_test)) + } + assert!(found_t) + } + + fn test_t( + curve: &T, + ) { + let tenths = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; + for &t in tenths.iter() { + test_t_along_curve(curve, t); + } + } + #[test] fn test_cubic_cubic_intersections() { + test_t(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0))); + test_t(&CubicBez::new((0.0, 0.0), (-1.0, 0.0), (1.0, 0.0), (1.0, 0.0))); + do_test( &CubicBez::new((0.0, 0.0), (0.0, 1.0), (0.0, 1.0), (1.0, 1.0)), &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), @@ -562,6 +579,17 @@ mod tests { &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, ); + // THESE TESTS FAIL, WORKING ON FIXES + // do_test( + // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + // 2, + // ); + // do_test( + // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + // &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), + // 2, + // ); // (A previous version of the code was returning two practically identical // intersection points here.) diff --git a/src/lib.rs b/src/lib.rs index e2111e1f..d857507e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -71,7 +71,7 @@ #![forbid(unsafe_code)] #![deny(missing_docs, clippy::trivially_copy_pass_by_ref)] -#![warn(broken_intra_doc_links)] +#![warn(rustdoc::broken_intra_doc_links)] #![allow( clippy::unreadable_literal, clippy::many_single_char_names, @@ -101,6 +101,7 @@ mod size; mod svg; mod translate_scale; mod vec2; +mod real; pub use crate::affine::*; pub use crate::arc::*; diff --git a/src/real.rs b/src/real.rs new file mode 100644 index 00000000..9b933150 --- /dev/null +++ b/src/real.rs @@ -0,0 +1,120 @@ +//! A utility for comparing real values + +use crate::{Point}; + +/// If we're comparing numbers, our epsilon should depend on how big the number +/// is. This function returns an epsilon appropriate for the size the largest +/// number. +pub fn epsilon_for_value(num: f64) -> f64 { + // These values have been obtained experimentally and can be changed if + // necessary. + (num.abs() * 0.0000000001).max(f64::EPSILON) +} + +/// If we're comparing distances between samples of curves, our epsilon should +/// depend on how big the points we're comparing are. This function returns an +/// epsilon appropriate for the size of pt. +pub fn epsilon_for_point(pt: Point) -> f64 { + let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); + epsilon_for_value(max) +} + +/// Compare if two real values are approximately equal +pub fn real_is_equal(num1: f64, num2: f64) -> bool +{ + if num1 == f64::INFINITY || num2 == f64::INFINITY + { + return num1 == num2; // Do exact comparison + } + + // Do approximate comparison + (num1 - num2).abs() <= epsilon_for_value(num1.max(num2)) +} + +/// Compare if a real value is approximately zero +pub fn real_is_zero(num1: f64) -> bool +{ + real_is_equal(num1, 0.) +} + +/// Compare if num1 < num2 and not approximately equal +pub fn real_lt(num1: f64, num2: f64) -> bool +{ + (num1 < num2) && !real_is_equal(num1, num2) +} + +/// Compare if num1 < num2 or approximately equal +pub fn real_lte(num1: f64, num2: f64) -> bool +{ + (num1 < num2) || real_is_equal(num1, num2) +} + +/// Compare if num1 > num2 and not approximately equal +pub fn real_gt(num1: f64, num2: f64) -> bool +{ + (num1 > num2) && !real_is_equal(num1, num2) +} + +/// Compare if num1 > num2 or approximately equal +pub fn real_gte(num1: f64, num2: f64) -> bool +{ + (num1 > num2) || real_is_equal(num1, num2) +} + +/// Compare if two points are approximately equal +pub fn point_is_equal(pt1: Point, pt2: Point) -> bool +{ + real_is_equal(pt1.x, pt2.x) && real_is_equal(pt1.y, pt2.y) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_real_comparisons() { + // Equals + assert!(real_is_equal(1., 1.)); + assert!(real_is_equal(1000000., 1000000.)); + assert!(real_is_equal(f64::INFINITY, f64::INFINITY)); + assert!(!real_is_equal(f64::INFINITY, f64::NEG_INFINITY)); + assert!(real_is_equal(f64::EPSILON, f64::EPSILON)); + assert!(!real_is_equal(f64::EPSILON, -f64::EPSILON)); + assert!(real_is_equal(f64::EPSILON, 0.)); + assert!(real_is_equal(0., f64::EPSILON)); + + // Less than + assert!(real_lte(1., 1.)); + assert!(real_lte(1000000., 1000000.)); + assert!(real_lte(f64::INFINITY, f64::INFINITY)); + assert!(!real_lt(1., 1.)); + assert!(!real_lt(1000000., 1000000.)); + assert!(!real_lt(f64::INFINITY, f64::INFINITY)); + assert!(real_lt(0., f64::INFINITY)); + assert!(real_lt(f64::NEG_INFINITY, 0.)); + assert!(real_lt(-1., 1.)); + assert!(real_lt(f64::NEG_INFINITY, f64::INFINITY)); + + // Greater than + assert!(real_gte(1., 1.)); + assert!(real_gte(1000000., 1000000.)); + assert!(real_gte(f64::INFINITY, f64::INFINITY)); + assert!(!real_gt(1., 1.)); + assert!(!real_gt(1000000., 1000000.)); + assert!(!real_gt(f64::INFINITY, f64::INFINITY)); + assert!(real_gt(f64::INFINITY, 0.)); + assert!(real_gt(0., f64::NEG_INFINITY)); + assert!(real_gt(1., -1.)); + assert!(real_gt(f64::INFINITY, f64::NEG_INFINITY)); + } + + #[test] + fn test_point_comparisons() { + assert!(point_is_equal(Point::new(0., 0.), Point::new(0., 0.))); + assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.))); + assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY))); + assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1.))); + assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY))); + assert!(point_is_equal(Point::new(0., 0.), Point::new(f64::EPSILON, f64::EPSILON))); + } +} \ No newline at end of file From d5d0be89dd6257a3cf464d19594d77c94cf40dc6 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Tue, 1 Feb 2022 08:16:49 -0700 Subject: [PATCH 04/20] Squashed in curve-curve-intersection-cubic-quad-line: * Added tests * Enabled intersections between lines and other curve types * Added tests for intersections between quadratic bezier curves * Moved ParamCurveBezierClipping into new file --- src/common.rs | 21 +++- src/cubicbez.rs | 103 +------------------ src/curve_intersections.rs | 32 +++++- src/lib.rs | 2 + src/line.rs | 8 ++ src/param_curve_clipping.rs | 200 ++++++++++++++++++++++++++++++++++++ src/quadbez.rs | 19 +++- 7 files changed, 279 insertions(+), 106 deletions(-) create mode 100644 src/param_curve_clipping.rs diff --git a/src/common.rs b/src/common.rs index 6a78a2a1..a7151382 100644 --- a/src/common.rs +++ b/src/common.rs @@ -118,12 +118,8 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec<[f64; 2]> { let sc1 = c1 * c2.recip(); if !sc0.is_finite() || !sc1.is_finite() { // c2 is zero or very small, treat as linear eqn - let root = -c0 / c1; - if root.is_finite() { + for root in solve_linear(c0, c1) { result.push(root); - } else if c0 == 0.0 && c1 == 0.0 { - // Degenerate case - result.push(0.0); } return result; } @@ -158,6 +154,21 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec<[f64; 2]> { result } +/// Find real roots of linear equation. +/// +/// Return values of x for which c0 + c1 x = 0. +pub fn solve_linear(c0: f64, c1: f64) -> ArrayVec<[f64; 1]> { + let mut result = ArrayVec::new(); + let root = -c0 / c1; + if root.is_finite() { + result.push(root); + } else if c0 == 0.0 && c1 == 0.0 { + // Degenerate case + result.push(0.0); + } + result +} + /// Solve an arbitrary function for a zero-crossing. /// /// This uses the [ITP method], as described in the paper diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 4ddc27fa..50540de1 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -7,7 +7,7 @@ use crate::{Line, QuadSpline, Vec2}; use arrayvec::ArrayVec; use crate::common::GAUSS_LEGENDRE_COEFFS_9; -use crate::common::{min_max, solve_cubic, solve_quadratic}; +use crate::common::{solve_quadratic}; use crate::{ Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, @@ -235,7 +235,9 @@ impl CubicBez { }) } - fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) { + /// Get the parameters such that the curve can be represented by the following formula: + /// B(t) = a*t^3 + b*t^2 + c*t + d + pub fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) { let c = (self.p1 - self.p0) * 3.0; let b = (self.p2 - self.p1) * 3.0 - c; let d = self.p0.to_vec2(); @@ -309,7 +311,7 @@ impl CubicBez { /// Is this cubic Bezier curve a line? #[inline] - fn is_linear(&self, accuracy: f64) -> bool { + pub fn is_linear(&self, accuracy: f64) -> bool { self.baseline().nearest(self.p1, accuracy).distance_sq <= accuracy && self.baseline().nearest(self.p2, accuracy).distance_sq <= accuracy } @@ -535,101 +537,6 @@ impl ParamCurveExtrema for CubicBez { } } -// Note that the line is unbounded here! -fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { - let vec2 = l.p1 - l.p0; - let a = -vec2.y; - let b = vec2.x; - let c = -(a * l.start().x + b * l.start().y); - a * p.x + b * p.y + c -} - -impl ParamCurveBezierClipping for CubicBez { - fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { - if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { - return ArrayVec::new(); - } - let (a, b, c, d) = self.parameters(); - solve_cubic(d.x - x, c.x, b.x, a.x) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { - if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { - return ArrayVec::new(); - } - - let (a, b, c, d) = self.parameters(); - solve_cubic(d.y - y, c.y, b.y, a.y) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - - fn fat_line_min_max(&self) -> (f64, f64) { - let baseline = self.baseline(); - let (d1, d2) = min_max( - signed_distance_from_ray_to_point(&baseline, self.p1), - signed_distance_from_ray_to_point(&baseline, self.p2), - ); - let factor = if (d1 * d2) > 0.0 { - 3.0 / 4.0 - } else { - 4.0 / 9.0 - }; - - let d_min = factor * f64::min(d1, 0.0); - let d_max = factor * f64::max(d2, 0.0); - - (d_min, d_max) - } - - fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { - let d0 = signed_distance_from_ray_to_point(l, self.start()); - let d1 = signed_distance_from_ray_to_point(l, self.p1); - let d2 = signed_distance_from_ray_to_point(l, self.p2); - let d3 = signed_distance_from_ray_to_point(l, self.end()); - - let p0 = Point::new(0.0, d0); - let p1 = Point::new(1.0 / 3.0, d1); - let p2 = Point::new(2.0 / 3.0, d2); - let p3 = Point::new(1.0, d3); - // Compute the vertical signed distance of p1 and p2 from [p0, p3]. - let dist1 = d1 - (2.0 * d0 + d3) / 3.0; - let dist2 = d2 - (d0 + 2.0 * d3) / 3.0; - - // Compute the hull assuming p1 is on top - we'll switch later if needed. - let mut hull = if dist1 * dist2 < 0.0 { - // p1 and p2 lie on opposite sides of [p0, p3], so the hull is a quadrilateral: - (vec![p0, p1, p3], vec![p0, p2, p3]) - } else { - // p1 and p2 lie on the same side of [p0, p3]. The hull can be a triangle or a - // quadrilateral, and [p0, p3] is part of the hull. The hull is a triangle if the vertical - // distance of one of the middle points p1, p2 is <= half the vertical distance of the - // other middle point. - let dist1 = dist1.abs(); - let dist2 = dist2.abs(); - if dist1 >= 2.0 * dist2 { - (vec![p0, p1, p3], vec![p0, p3]) - } else if dist2 >= 2.0 * dist1 { - (vec![p0, p2, p3], vec![p0, p3]) - } else { - (vec![p0, p1, p2, p3], vec![p0, p3]) - } - }; - - // Flip the hull if needed: - if dist1 < 0.0 || (dist1 == 0.0 && dist2 < 0.0) { - hull = (hull.1, hull.0); - } - - hull - } -} - impl Mul for Affine { type Output = CubicBez; diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 966ed64b..1874ea07 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -509,7 +509,7 @@ fn rectangles_overlap(r1: &Rect, r2: &Rect) -> bool { #[cfg(test)] mod tests { use super::*; - use crate::CubicBez; + use crate::{CubicBez, QuadBez, Line}; fn do_test( curve1: &T, @@ -579,6 +579,36 @@ mod tests { &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, ); + do_test( + &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), + &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), + 3, + ); + do_test( + &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), + &QuadBez::new((0.0, 0.25), (0.5, 0.75), (1.0, 0.25)), + 1, + ); + do_test( + &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), + &QuadBez::new((0.0, 0.25), (0.5, -0.25), (1.0, 0.25)), + 2, + ); + do_test( + &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), + &Line::new((0.0, 0.5), (1.0, 0.25)), + 2, + ); + do_test( + &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), + &Line::new((0.0, 0.5), (1.0, 0.5)), + 1, + ); + do_test( + &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), + &Line::new((0.0, 0.5), (1.0, 0.25)), + 3, + ); // THESE TESTS FAIL, WORKING ON FIXES // do_test( // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), diff --git a/src/lib.rs b/src/lib.rs index d857507e..7de0a0e8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -90,6 +90,7 @@ mod insets; mod line; mod mindist; mod param_curve; +mod param_curve_clipping; mod point; mod quadbez; mod quadspline; @@ -113,6 +114,7 @@ pub use crate::ellipse::*; pub use crate::insets::*; pub use crate::line::*; pub use crate::param_curve::*; +pub use crate::param_curve_clipping::*; pub use crate::point::*; pub use crate::quadbez::*; pub use crate::quadspline::*; diff --git a/src/line.rs b/src/line.rs index 62b97c42..9375c09a 100644 --- a/src/line.rs +++ b/src/line.rs @@ -48,6 +48,14 @@ impl Line { Some(other.p0 + cd * h) } + /// Get the parameters such that the curve can be represented by the following formula: + /// B(t) = a*t + b + pub fn parameters(&self) -> (Vec2, Vec2) { + let b = self.p0.to_vec2(); + let a = self.p1 - self.p0; + (a, b) + } + /// Is this line finite? #[inline] pub fn is_finite(self) -> bool { diff --git a/src/param_curve_clipping.rs b/src/param_curve_clipping.rs new file mode 100644 index 00000000..639fc966 --- /dev/null +++ b/src/param_curve_clipping.rs @@ -0,0 +1,200 @@ +//! A trait for clipping parametrized curves. + +use crate::common::{min_max, solve_cubic, solve_quadratic, solve_linear}; +use crate::{CubicBez, QuadBez, Line, Point, ParamCurve, ParamCurveBezierClipping}; +use arrayvec::ArrayVec; + +// Note that the line is unbounded here! +fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { + let vec2 = l.p1 - l.p0; + let a = -vec2.y; + let b = vec2.x; + let c = -(a * l.start().x + b * l.start().y); + a * p.x + b * p.y + c +} + +impl ParamCurveBezierClipping for Line { + fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + if (self.p0.x - self.p1.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b) = self.parameters(); + solve_linear(b.x - x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + if (self.p0.y - self.p1.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b) = self.parameters(); + solve_linear(b.y - y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self) -> (f64, f64) { + (0., 0.) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0, d1); + + // The hull is simply the line itself + (vec![p0, p1], vec![p0, p1]) + } +} + +impl ParamCurveBezierClipping for QuadBez { + fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.x - self.p2.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b, c) = self.parameters(); + solve_quadratic(c.x - x, b.x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.y - self.p2.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b, c) = self.parameters(); + solve_quadratic(c.y - y, b.y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self) -> (f64, f64) { + let baseline = self.baseline(); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let factor = 1.0 / 2.0; + + let d_min = factor * f64::min(d1, 0.0); + let d_max = factor * f64::max(d1, 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.p1); + let d2 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0 / 2.0, d1); + let p2 = Point::new(1.0, d2); + // Compute the vertical signed distance of p1 and p2 from [p0, p3]. + let dist1 = d1 - (d0 + d2) / 2.0; + + // Compute the hull assuming p1 is on top - we'll switch later if needed. + let mut hull = (vec![p0, p1, p2], vec![p0, p2]); + + // Flip the hull if needed: + if dist1 < 0.0 { + hull = (hull.1, hull.0); + } + + hull + } +} + +impl ParamCurveBezierClipping for CubicBez { + fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b, c, d) = self.parameters(); + solve_cubic(d.x - x, c.x, b.x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b, c, d) = self.parameters(); + solve_cubic(d.y - y, c.y, b.y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self) -> (f64, f64) { + let baseline = self.baseline(); + let (d1, d2) = min_max( + signed_distance_from_ray_to_point(&baseline, self.p1), + signed_distance_from_ray_to_point(&baseline, self.p2), + ); + let factor = if (d1 * d2) > 0.0 { + 3.0 / 4.0 + } else { + 4.0 / 9.0 + }; + + let d_min = factor * f64::min(d1, 0.0); + let d_max = factor * f64::max(d2, 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.p1); + let d2 = signed_distance_from_ray_to_point(l, self.p2); + let d3 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0 / 3.0, d1); + let p2 = Point::new(2.0 / 3.0, d2); + let p3 = Point::new(1.0, d3); + // Compute the vertical signed distance of p1 and p2 from [p0, p3]. + let dist1 = d1 - (2.0 * d0 + d3) / 3.0; + let dist2 = d2 - (d0 + 2.0 * d3) / 3.0; + + // Compute the hull assuming p1 is on top - we'll switch later if needed. + let mut hull = if dist1 * dist2 < 0.0 { + // p1 and p2 lie on opposite sides of [p0, p3], so the hull is a quadrilateral: + (vec![p0, p1, p3], vec![p0, p2, p3]) + } else { + // p1 and p2 lie on the same side of [p0, p3]. The hull can be a triangle or a + // quadrilateral, and [p0, p3] is part of the hull. The hull is a triangle if the vertical + // distance of one of the middle points p1, p2 is <= half the vertical distance of the + // other middle point. + let dist1 = dist1.abs(); + let dist2 = dist2.abs(); + if dist1 >= 2.0 * dist2 { + (vec![p0, p1, p3], vec![p0, p3]) + } else if dist2 >= 2.0 * dist1 { + (vec![p0, p2, p3], vec![p0, p3]) + } else { + (vec![p0, p1, p2, p3], vec![p0, p3]) + } + }; + + // Flip the hull if needed: + if dist1 < 0.0 || (dist1 == 0.0 && dist2 < 0.0) { + hull = (hull.1, hull.0); + } + + hull + } +} diff --git a/src/quadbez.rs b/src/quadbez.rs index 5759136c..080c9d40 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -7,9 +7,9 @@ use arrayvec::ArrayVec; use crate::common::solve_cubic; use crate::MAX_EXTREMA; use crate::{ - Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, + Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, - Rect, Shape, + Rect, Shape, Vec2 }; /// A single quadratic Bézier segment. @@ -92,6 +92,15 @@ impl QuadBez { (u - params.u0) * params.uscale } + /// Get the parameters such that the curve can be represented by the following formula: + /// B(t) = a*t^2 + b*t + c + pub fn parameters(&self) -> (Vec2, Vec2, Vec2) { + let c = self.p0.to_vec2(); + let b = (self.p1 - self.p0) * 2.0; + let a = c - self.p1.to_vec2() * 2.0 + self.p2.to_vec2(); + (a, b, c) + } + /// Is this quadratic Bezier curve finite? #[inline] pub fn is_finite(&self) -> bool { @@ -103,6 +112,12 @@ impl QuadBez { pub fn is_nan(&self) -> bool { self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan() } + + /// Is this quadratic Bezier curve a line? + #[inline] + pub fn is_linear(&self, accuracy: f64) -> bool { + self.baseline().nearest(self.p1, accuracy).distance_sq <= accuracy + } } /// An iterator for quadratic beziers. From 1375912c9ab798c1d659849f9fe66e3d302ad26a Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Mon, 7 Feb 2022 14:08:41 -0700 Subject: [PATCH 05/20] Squashed in curve-curve-intersection-overlap: * Added more flexibility to root finding to account for edge cases * Added function that checks if two curves perfectly overlap for some range * Added flags to curve intersection methods --- Cargo.toml | 1 + src/curve_intersections.rs | 319 ++++++++++++++++++++++++++++++++----- src/lib.rs | 1 + src/param_curve.rs | 2 +- src/real.rs | 18 ++- 5 files changed, 296 insertions(+), 45 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 63a57e0d..375dfeff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,6 +15,7 @@ features = ["mint", "serde"] [dependencies] arrayvec = "0.5.1" +bitflags = "1.3" [dependencies.mint] version = "0.5.1" diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 1874ea07..4d991c2f 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -5,12 +5,27 @@ use crate::real::*; use arrayvec::ArrayVec; use std::ops::Range; +bitflags::bitflags! { + pub struct CurveIntersectionFlags: u32 { + const NONE = 0b00000000; + const KEEP_ENDPOINT_INTERSECTIONS = 0b00000001; + const KEEP_DUPLICATE_INTERSECTIONS = 0b00000010; + } +} + /// Compute the intersections between two Bézier curves pub fn curve_curve_intersections( curve1: &T, curve2: &U, + flags: CurveIntersectionFlags, ) -> ArrayVec<[(f64, f64); 9]> { let mut av = ArrayVec::new(); + + let overlaps = check_for_overlap(curve1, curve2, flags, &mut av); + if overlaps { + return av; + } + add_curve_intersections( curve1, curve2, @@ -22,10 +37,129 @@ pub fn curve_curve_intersections( + curve1: &T, + curve2: &U, + flags: CurveIntersectionFlags, + intersections: &mut ArrayVec<[(f64, f64); 9]>, +) -> bool { + // If the two curves overlap, there are two main cases we need to account for: + // 1) One curve is completely along some portion of the other curve, thus + // both endpoints of the first curve lie along the other curve. + // ex. curve 1: O------------------O + // curve 2: O---O + // 2) The two curves overlap, thus one endpoint from each curve lies + // along the other. + // ex. curve 1: O------------------O + // curve 2: O-----------O + + // Use this function to get the t-value on another curve that corresponds to the t-value on this curve + #[inline] + fn t_value_on_other( + t_this: f64, + curve_this: &T, + curve_other: &U, + ) -> (bool, f64) { + let t_values_other = point_is_on_curve(curve_this.eval(t_this), curve_other); + if !t_values_other.0 { + return (false, 0.0); + } + + assert!(t_values_other.1.len() >= 1); // This has to be true or we really messed up! + if t_values_other.1.len() == 1 { + return (true, t_values_other.1[0]) + } + + // Things get kind of complicated if we have more than one t value + // here. The best solution I can think of is to compare the derivatives + // at each point. + let curve_deriv_at_t_this = curve_this.deriv().eval(t_this); + let curve_deriv_other = curve_other.deriv(); + + // Find the t-value with the closest derivative + let mut deriv = Point::new(f64::MAX, f64::MAX); + let mut index = 0; + for (i, t) in t_values_other.1.iter().enumerate() { + let curve_deriv_at_t_other = curve_deriv_other.eval(*t); + let d = (curve_deriv_at_t_other - deriv).hypot2(); + let d_compare = (curve_deriv_at_t_this - deriv).hypot2(); + if d < d_compare { + deriv = curve_deriv_at_t_other; + index = i; + } + } + (true, t_values_other.1[index]) + } + + // First, test which endpoints, if any lie along the other curve + let start1_on_2 = t_value_on_other(0., curve1, curve2); + let end1_on_2 = t_value_on_other(1., curve1, curve2); + let start2_on_1 = t_value_on_other(0., curve2, curve1); + let end2_on_1 = t_value_on_other(1., curve2, curve1); + + // At least 2 endpoints must lie on the other curve for us to proceed + let count = if start1_on_2.0 { 1 } else { 0 } + + if end1_on_2. 0 { 1 } else { 0 } + + if start2_on_1.0 { 1 } else { 0 } + + if end2_on_1.0 { 1 } else { 0 }; + if count < 2 { + return false; // No overlap + } + + // Check a few intermediate points to make sure that it's not just the + // endpoints that lie on top of the other curve. I'll make the claim that + // two curves can intersect at 4 points, but not be co-linear (or co-curved?). + // This could happen if two bezier with the same start and end point, one + // curved and one looped, intersected where the second curve intersects itself. + // Therefore, it's safe if we check 5 positions in total. Since we've already + // checked the endpoints, we just need to check 3 intermediate positions. + + fn get_projection(t_value: (bool, f64), default_value: f64) -> f64 { + if t_value.0 { t_value.1 } else { default_value } + } + + // Determine the range of curve1 that lies along curve2 + let overlap_along_curve1 = ( + if start1_on_2.0 { 0. } else { (get_projection(start2_on_1, 1.)).min(get_projection(end2_on_1, 1.)) }, + if end1_on_2.0 { 1. } else { (get_projection(start2_on_1, 0.)).max(get_projection(end2_on_1, 0.)) }, + ); + assert!(overlap_along_curve1.0 <= overlap_along_curve1.1); + + // Determine the range of curve2 that lies along curve1 such that + // overlap_along_curve2.0 corresponds to the same position in space as + // overlap_along_curve1.0. + // Note: It's not guaranteed that overlap_along_curve2.0 < overlap_along_curve2.1 + let overlap_along_curve2 = ( + if start1_on_2.0 { start1_on_2.1 } + else { let p1 = get_projection(start2_on_1, 1.); let p2 = get_projection(end2_on_1, 1.); if p1 < p2 { 0. } else { 1. } }, + if end1_on_2.0 { end1_on_2.1 } + else { let p1 = get_projection(start2_on_1, 0.); let p2 = get_projection(end2_on_1, 0.); if p1 < p2 { 1. } else { 0. } }, + ); + + // Check the 3 intermediate locations + let interval1 = (overlap_along_curve1.1 - overlap_along_curve1.0) / 4.; + let interval2 = (overlap_along_curve2.1 - overlap_along_curve2.0) / 4.; + for n in 1..3 { + if !point_is_equal(curve1.eval(overlap_along_curve1.0 + (n as f64) * interval1), curve2.eval(overlap_along_curve2.0 + (n as f64) * interval2)) { + return false; + } + } + + if flags.contains(CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS) { + // Add the endpoints as intersections + intersections.push((overlap_along_curve1.0, overlap_along_curve2.0)); + intersections.push((overlap_along_curve1.1, overlap_along_curve2.1)); + } + + true // The two curves overlap! +} + // This function implements the main bézier clipping algorithm by recursively subdividing curve1 and // curve2 in to smaller and smaller portions of the original curves with the property that one of // the curves intersects the fat line of the other curve at each stage. @@ -47,6 +181,7 @@ fn add_curve_intersections u32 { call_count += 1; recursion_count += 1; @@ -65,6 +200,7 @@ fn add_curve_intersections 1.0 - end_eps) - && (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot() > 5.0 + if !flags.contains(CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS) + && (real_lte(t2, 0.0) || real_gte(t2, 1.0)) + && point_is_equal(orig_curve1.eval(t1), orig_curve2.eval(t2)) { return call_count; } - add_intersection(t1, orig_curve1, t2, orig_curve2, flip, intersections); + add_intersection(t1, orig_curve1, t2, orig_curve2, flip, intersections, flags); return call_count; } @@ -146,6 +284,7 @@ fn add_curve_intersections, intersections: &mut ArrayVec<[(f64, f64); 9]>, flip: bool, + flags: CurveIntersectionFlags, ) { let pt = pt_curve.eval(0.5); // We assume pt is curve1 when we add intersections below. @@ -282,7 +428,7 @@ fn add_point_curve_intersection( pt: Point, curve: &T, -) -> ArrayVec<[f64; 9]> { +) -> ArrayVec<[f64; 2]> { let mut result = ArrayVec::new(); // If both endpoints are approximately close, we only return 0.0. @@ -320,12 +466,15 @@ pub fn t_along_curve_for_point( for params in [curve_x_t_params, curve_y_t_params].iter() { for t in params { let t = *t; - if !point_is_equal(pt, curve.eval(t)) { + + // Note: use point_is_approx below for a little more wiggle room + // since our root finding algorithm is perfect + if !point_is_approx(pt, curve.eval(t), 2.) { continue; } let mut already_found_t = false; for u in &result { - if f64::abs(t - *u) < param_eps { + if real_is_equal(t, *u) { already_found_t = true; break; } @@ -352,7 +501,7 @@ pub fn t_along_curve_for_point( t: f64, pt: Point, curve: &T, - result: &mut ArrayVec<[f64; 9]>, + result: &mut ArrayVec<[f64; 2]>, ) -> bool { if point_is_equal(curve.eval(t), pt) { result.push(t); @@ -368,6 +517,19 @@ pub fn t_along_curve_for_point( result } +fn point_is_on_curve( + pt: Point, + curve: &T, +) -> (bool, ArrayVec<[f64; 2]>) { + let t_values = t_along_curve_for_point(pt, curve); + for t in &t_values { + if point_is_equal(pt, curve.eval(*t)) { + return (true, t_values.clone()); // Point lies along curve + } + } + (false, t_values.clone()) // Point is not along curve +} + fn add_intersection( t1: f64, orig_curve1: &T, @@ -375,31 +537,40 @@ fn add_intersection( orig_curve2: &U, flip: bool, intersections: &mut ArrayVec<[(f64, f64); 9]>, + flags: CurveIntersectionFlags, ) { let (t1, t2) = if flip { (t2, t1) } else { (t1, t2) }; - // (This should probably depend in some way on how large our input coefficients are.) - let epsilon = 1e-3; + // Discard endpoint/endpoint intersections. - let t1_is_an_endpoint = t1 < epsilon || t1 > 1.0 - epsilon; - let t2_is_an_endpoint = t2 < epsilon || t2 > 1.0 - epsilon; - if t1_is_an_endpoint && t2_is_an_endpoint { - return; + if !flags.contains(CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS) + { + let t1_is_an_endpoint = real_lte(t1, 0.0) || real_gte(t1, 1.0); + let t2_is_an_endpoint = real_lte(t2, 0.0) || real_gte(t2, 1.0); + if t1_is_an_endpoint && t2_is_an_endpoint { + return; + } } - // We can get repeated intersections when we split a curve at an intersection point, or when - // two curves intersect at a point where the curves are very close together, or when the fat - // line process breaks down. - for i in 0..intersections.len() { - let (old_t1, old_t2) = intersections[i]; - // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" - // duplicate intersection, so always keep the one that minimizes sample distance. - if (t1 - old_t1).abs() < epsilon && (t2 - old_t2).abs() < epsilon { - let cur_dist = (orig_curve1.eval(old_t1) - orig_curve2.eval(old_t2)).hypot2(); - let new_dist = (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot2(); - if new_dist < cur_dist { - intersections[i] = (t1, t2); + if !flags.contains(CurveIntersectionFlags::KEEP_DUPLICATE_INTERSECTIONS) + { + // (This should probably depend in some way on how large our input coefficients are.) + let epsilon = 1e-3; + + // We can get repeated intersections when we split a curve at an intersection point, or when + // two curves intersect at a point where the curves are very close together, or when the fat + // line process breaks down. + for i in 0..intersections.len() { + let (old_t1, old_t2) = intersections[i]; + // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" + // duplicate intersection, so always keep the one that minimizes sample distance. + if (t1 - old_t1).abs() < epsilon && (t2 - old_t2).abs() < epsilon { + let cur_dist = (orig_curve1.eval(old_t1) - orig_curve2.eval(old_t2)).hypot2(); + let new_dist = (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot2(); + if new_dist < cur_dist { + intersections[i] = (t1, t2); + } + return; } - return; } } @@ -515,9 +686,10 @@ mod tests { curve1: &T, curve2: &U, count: usize, + flags: CurveIntersectionFlags, ) { - let arr1 = curve_curve_intersections(curve1, curve2); - let arr2 = curve_curve_intersections(curve2, curve1); + let arr1 = curve_curve_intersections(curve1, curve2, flags); + let arr2 = curve_curve_intersections(curve2, curve1, flags); assert_eq!(arr1.len(), count); assert_eq!(arr2.len(), count); } @@ -554,6 +726,39 @@ mod tests { } } + fn test_overlapping( + curve: &T, + t1: f64, + t2: f64, + ) { + assert!(t1 < t2); + + do_test( + &curve.subsegment(0.0..t2), + &curve.subsegment(t1..1.0), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + do_test( + &curve.subsegment(0.0..t2), + &curve.subsegment(1.0..t1), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + do_test( + &curve.subsegment(0.0..1.0), + &curve.subsegment(t1..t2), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + do_test( + &curve.subsegment(0.0..1.0), + &curve.subsegment(t2..t1), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + } + #[test] fn test_cubic_cubic_intersections() { test_t(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0))); @@ -563,63 +768,91 @@ mod tests { &CubicBez::new((0.0, 0.0), (0.0, 1.0), (0.0, 1.0), (1.0, 1.0)), &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), 1, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &CubicBez::new((48.0, 84.0), (104.0, 176.0), (190.0, 37.0), (121.0, 75.0)), &CubicBez::new((68.0, 145.0), (74.0, 6.0), (143.0, 197.0), (138.0, 55.0)), 4, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &CubicBez::new((0.0, 0.0), (0.5, 1.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.0, 1.0), (0.5, 0.0), (0.5, 0.0), (1.0, 1.0)), 2, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), 3, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, 0.75), (1.0, 0.25)), 1, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, -0.25), (1.0, 0.25)), 2, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 2, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.5)), 1, + CurveIntersectionFlags::NONE, // Standard algorithm ); do_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 3, + CurveIntersectionFlags::NONE, // Standard algorithm + ); + + // These tests intersect at the endpoints which is only sometimes useful to know about. + do_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + do_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), + 0, + CurveIntersectionFlags::NONE, // Discard endpoint intersections ); - // THESE TESTS FAIL, WORKING ON FIXES - // do_test( - // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), - // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), - // 2, - // ); - // do_test( - // &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), - // &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), - // 2, - // ); + do_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), + 2, + CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + do_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), + 0, + CurveIntersectionFlags::NONE, // Discard endpoint intersections + ); + + // Test curves that lie exactly along one-another + test_overlapping(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), 0.2, 0.8); // (A previous version of the code was returning two practically identical // intersection points here.) @@ -637,6 +870,7 @@ mod tests { (-113963.134524995, 732017.9466050486), ), 1, + CurveIntersectionFlags::NONE, // Standard algorithm ); // On these curves the algorithm runs to a state at which the new clipped domain1 becomes a @@ -657,6 +891,7 @@ mod tests { (567688.5977972192, 13975.09633399453), ), 3, + CurveIntersectionFlags::NONE, // Standard algorithm ); } } diff --git a/src/lib.rs b/src/lib.rs index 7de0a0e8..f131e667 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -126,3 +126,4 @@ pub use crate::size::*; pub use crate::svg::*; pub use crate::translate_scale::*; pub use crate::vec2::*; +pub use crate::real::*; diff --git a/src/param_curve.rs b/src/param_curve.rs index bf06d6cd..0c903813 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -216,7 +216,7 @@ pub trait ParamCurveExtrema: ParamCurve { } /// A parameterized curve that can be used in the Bezier clipping algorithm -pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveExtrema + ParamCurveArclen { +pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen { /// Returns a line from the curve's start point to its end point fn baseline(&self) -> Line { Line::new(self.start(), self.end()) diff --git a/src/real.rs b/src/real.rs index 9b933150..b8d06d44 100644 --- a/src/real.rs +++ b/src/real.rs @@ -21,6 +21,13 @@ pub fn epsilon_for_point(pt: Point) -> f64 { /// Compare if two real values are approximately equal pub fn real_is_equal(num1: f64, num2: f64) -> bool +{ + real_is_approx(num1, num2, 1.) +} + +/// Compare if two real values are approximately equal to a smaller degree than +/// real_is_equal +pub fn real_is_approx(num1: f64, num2: f64, scale: f64) -> bool { if num1 == f64::INFINITY || num2 == f64::INFINITY { @@ -28,7 +35,7 @@ pub fn real_is_equal(num1: f64, num2: f64) -> bool } // Do approximate comparison - (num1 - num2).abs() <= epsilon_for_value(num1.max(num2)) + (num1 - num2).abs() <= (scale * epsilon_for_value(num1.max(num2))) } /// Compare if a real value is approximately zero @@ -64,7 +71,14 @@ pub fn real_gte(num1: f64, num2: f64) -> bool /// Compare if two points are approximately equal pub fn point_is_equal(pt1: Point, pt2: Point) -> bool { - real_is_equal(pt1.x, pt2.x) && real_is_equal(pt1.y, pt2.y) + point_is_approx(pt1, pt2, 1.) +} + +/// Compare if two points are approximately equal to a smaller degree than +/// point_is_equal +pub fn point_is_approx(pt1: Point, pt2: Point, scale: f64) -> bool +{ + real_is_approx(pt1.x, pt2.x, scale) && real_is_approx(pt1.y, pt2.y, scale) } #[cfg(test)] From dc5b41a15a10e622f850c15674d96fa4019379f6 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 24 Feb 2022 09:00:50 -0700 Subject: [PATCH 06/20] Squashed in curve-curve-intersection-perpendicular-fatline: * Now checking perpendicular fat line, fixed bugs --- src/common.rs | 10 ------ src/curve_intersections.rs | 58 ++++++++++++++++++++++++++----- src/param_curve.rs | 2 +- src/param_curve_clipping.rs | 69 ++++++++++++++++++++++++++++--------- src/quadbez.rs | 5 +++ 5 files changed, 108 insertions(+), 36 deletions(-) diff --git a/src/common.rs b/src/common.rs index a7151382..fdd80caa 100644 --- a/src/common.rs +++ b/src/common.rs @@ -43,16 +43,6 @@ impl FloatExt for f32 { } } -/// Order two things into minimum and maximum -#[inline] -pub fn min_max(a: T, b: T) -> (T, T) { - if a < b { - (a, b) - } else { - (b, a) - } -} - /// Find real roots of cubic equation. /// /// The implementation is not (yet) fully robust, but it does handle the case diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 4d991c2f..c49f19b5 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,6 +1,6 @@ use crate::ParamCurveBezierClipping; use crate::{ParamCurve, ParamCurveExtrema}; -use crate::{Point, Rect}; +use crate::{Point, Rect, Affine}; use crate::real::*; use arrayvec::ArrayVec; use std::ops::Range; @@ -185,7 +185,7 @@ fn add_curve_intersections u32 { call_count += 1; recursion_count += 1; - if call_count >= 4096 || recursion_count >= 60 { + if call_count >= 4096 || recursion_count >= 128 { return call_count; } @@ -243,15 +243,28 @@ fn add_curve_intersections (min, max), - None => return call_count, + // Clip to the normal fatline and the perpendicular fatline and limit our search to the + // intersection of the two intervals. + let clip_norm = restrict_curve_to_fat_line(curve1, curve2); + let clip_perp = restrict_curve_to_perpendicular_fat_line(curve1, curve2); + let (t_min_clip, t_max_clip) = match clip_norm { + Some((min_norm, max_norm)) => { + match clip_perp { + Some((min_perp, max_perp)) => (min_norm.max(min_perp), max_norm.min(max_perp)), + None => (min_norm, max_norm), + } + } + None => { + match clip_perp { + Some((min_perp, max_perp)) => (min_perp, max_perp), + None => return call_count, + } + }, }; // t_min_clip and t_max_clip are (0, 1)-based, so project them back to get the new restricted // range: - let new_domain1 = - &(domain_value_at_t(domain1, t_min_clip)..domain_value_at_t(domain1, t_max_clip)); + let new_domain1 = &(domain_value_at_t(domain1, t_min_clip)..domain_value_at_t(domain1, t_max_clip)); if (domain2.end - domain2.start).max(new_domain1.end - new_domain1.start) < epsilon { let t1 = (new_domain1.start + new_domain1.end) * 0.5; @@ -462,7 +475,6 @@ pub fn t_along_curve_for_point( // directions, but the parameter calculations aren't very accurate, so give a little more // leeway there (TODO: this isn't perfect, as you might expect - the dupes that pass here are // currently being detected in add_intersection). - let param_eps = 10.0 * epsilon_for_point(pt); for params in [curve_x_t_params, curve_y_t_params].iter() { for t in params { let t = *t; @@ -596,7 +608,35 @@ fn restrict_curve_to_fat_line< let baseline2 = curve2.baseline(); let (mut top, mut bottom) = curve1.convex_hull_from_line(&baseline2); - let (d_min, d_max) = curve2.fat_line_min_max(); + let (d_min, d_max) = curve2.fat_line_min_max(&baseline2); + + clip_convex_hull_to_fat_line(&mut top, &mut bottom, d_min, d_max) +} + +// Returns an interval (t_min, t_max) with the property that for parameter values outside that +// interval, curve1 is guaranteed to not intersect curve2; uses the perpendicular fatline of +// curve2 as its basis for the guarantee. We use check the perpendicular fatline as well as the +// regular fatline because sometimes it converges to the result faster. +fn restrict_curve_to_perpendicular_fat_line< + T: ParamCurveBezierClipping + ParamCurve + ParamCurveExtrema, + U: ParamCurveBezierClipping + ParamCurve + ParamCurveExtrema, +>( + curve1: &T, + curve2: &U, +) -> Option<(f64, f64)> { + use std::f64::consts::PI; + + // TODO: Consider clipping against the perpendicular fat line as well (recommended by + // Sederberg). + // TODO: The current algorithm doesn't handle the (rare) case where curve1 and curve2 are + // overlapping lines. + + let baseline2 = curve2.baseline(); + let center_baseline2 = baseline2.bounding_box().center().to_vec2(); + let affine_baseline2 = Affine::translate(center_baseline2) * Affine::rotate(PI / 2.) * Affine::translate(-center_baseline2); + let baseline2 = affine_baseline2 * baseline2; + let (mut top, mut bottom) = curve1.convex_hull_from_line(&baseline2); + let (d_min, d_max) = curve2.fat_line_min_max(&baseline2); clip_convex_hull_to_fat_line(&mut top, &mut bottom, d_min, d_max) } diff --git a/src/param_curve.rs b/src/param_curve.rs index 0c903813..ac037647 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -229,5 +229,5 @@ pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveDeriv + ParamCurveExt /// Returns the upper and lower convex hull fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); /// Returns the minimum and maximum distances of the "fat line" enclosing this curve - fn fat_line_min_max(&self) -> (f64, f64); + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64); } diff --git a/src/param_curve_clipping.rs b/src/param_curve_clipping.rs index 639fc966..1d3c18a6 100644 --- a/src/param_curve_clipping.rs +++ b/src/param_curve_clipping.rs @@ -1,6 +1,6 @@ //! A trait for clipping parametrized curves. -use crate::common::{min_max, solve_cubic, solve_quadratic, solve_linear}; +use crate::common::{solve_cubic, solve_quadratic, solve_linear}; use crate::{CubicBez, QuadBez, Line, Point, ParamCurve, ParamCurveBezierClipping}; use arrayvec::ArrayVec; @@ -10,6 +10,19 @@ fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { let a = -vec2.y; let b = vec2.x; let c = -(a * l.start().x + b * l.start().y); + + let len = (a * a + b * b).sqrt(); + let a = a / len; + let b = b / len; + let c = c / len; + + if a.is_infinite() || b.is_infinite() || c.is_infinite() + { + // Can't compute distance from zero-length line, so return distance + // from p0 instead + return (p - l.p0).hypot(); + } + a * p.x + b * p.y + c } @@ -38,8 +51,16 @@ impl ParamCurveBezierClipping for Line { .collect() } - fn fat_line_min_max(&self) -> (f64, f64) { - (0., 0.) + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + + // Calculate min and max distances + let d_min = f64::min(d0.min(d1), 0.0); + let d_max = f64::max(d0.max(d1), 0.0); + + (d_min, d_max) } fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { @@ -79,13 +100,24 @@ impl ParamCurveBezierClipping for QuadBez { .collect() } - fn fat_line_min_max(&self) -> (f64, f64) { - let baseline = self.baseline(); + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); + + // Shrink d1 to a more accurate number, but test the distance compared + // to d0 and d3 in case we're looking at a perpendicular fat line let factor = 1.0 / 2.0; + //let factor = 1.0; + let d1_1 = (d1 - d0) * factor + d0; + let d1_2 = (d1 - d2) * factor + d2; + // let d1_1 = d1; + // let d1_2 = d1; - let d_min = factor * f64::min(d1, 0.0); - let d_max = factor * f64::max(d1, 0.0); + // Calculate min and max distances + let d_min = f64::min(d0.min(d1_1).min(d1_2).min(d2), 0.0); + let d_max = f64::max(d0.max(d1_1).max(d1_2).max(d2), 0.0); (d_min, d_max) } @@ -138,20 +170,25 @@ impl ParamCurveBezierClipping for CubicBez { .collect() } - fn fat_line_min_max(&self) -> (f64, f64) { - let baseline = self.baseline(); - let (d1, d2) = min_max( - signed_distance_from_ray_to_point(&baseline, self.p1), - signed_distance_from_ray_to_point(&baseline, self.p2), - ); - let factor = if (d1 * d2) > 0.0 { + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); + let d3 = signed_distance_from_ray_to_point(&baseline, self.p3); + + // Shrink d1 and d2 to more accurate numbers + let factor: f64 = if ((d1 - d0) * (d2 - d3)) > 0.0 { 3.0 / 4.0 } else { 4.0 / 9.0 }; + let d1 = (d1 - d0) * factor + d0; + let d2 = (d2 - d3) * factor + d3; - let d_min = factor * f64::min(d1, 0.0); - let d_max = factor * f64::max(d2, 0.0); + // Calculate min and max distances + let d_min = f64::min(d0.min(d1).min(d2).min(d3), 0.0); + let d_max = f64::max(d0.max(d1).max(d2).max(d3), 0.0); (d_min, d_max) } diff --git a/src/quadbez.rs b/src/quadbez.rs index 080c9d40..2abaa349 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -279,6 +279,11 @@ impl ParamCurveArclen for QuadBez { let c2 = 2.0 * c.sqrt(); let ba_c2 = b * a2 + c2; + if a2.is_infinite() { + // The arclength is zero + return 0. + } + let v0 = 0.25 * a2 * a2 * b * (2.0 * sabc - c2) + sabc; // TODO: justify and fine-tune this exact constant. if ba_c2 < 1e-13 { From f4cf1911e4b4b28058f9b15aed172059f07266eb Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Wed, 2 Mar 2022 11:32:04 -0700 Subject: [PATCH 07/20] Split KEEP_ENDPOINT_INTERSECTIONS into individual flags for each endpoint --- src/curve_intersections.rs | 172 ++++++++++++++++++++++++++----------- 1 file changed, 120 insertions(+), 52 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index c49f19b5..4fdc9f64 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -7,9 +7,18 @@ use std::ops::Range; bitflags::bitflags! { pub struct CurveIntersectionFlags: u32 { - const NONE = 0b00000000; - const KEEP_ENDPOINT_INTERSECTIONS = 0b00000001; - const KEEP_DUPLICATE_INTERSECTIONS = 0b00000010; + const NONE = 0b00000000; + const KEEP_CURVE1_T0_INTERSECTION = 0b00000001; + const KEEP_CURVE1_T1_INTERSECTION = 0b00000010; + const KEEP_CURVE2_T0_INTERSECTION = 0b00000100; + const KEEP_CURVE2_T1_INTERSECTION = 0b00001000; + const KEEP_CURVE1_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE1_T0_INTERSECTION.bits + | Self::KEEP_CURVE1_T1_INTERSECTION.bits; + const KEEP_CURVE2_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE2_T0_INTERSECTION.bits + | Self::KEEP_CURVE2_T1_INTERSECTION.bits; + const KEEP_ALL_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE1_ENDPOINT_INTERSECTIONS.bits + | Self::KEEP_CURVE2_ENDPOINT_INTERSECTIONS.bits; + const KEEP_DUPLICATE_INTERSECTIONS = 0b00010000; } } @@ -151,10 +160,25 @@ fn check_for_overlap( } } - if flags.contains(CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS) { - // Add the endpoints as intersections - intersections.push((overlap_along_curve1.0, overlap_along_curve2.0)); - intersections.push((overlap_along_curve1.1, overlap_along_curve2.1)); + let overlap_along_curve2_is_flipped = overlap_along_curve2.1 < overlap_along_curve2.0; + let overlap_along_curve2_0_is_start = !overlap_along_curve2_is_flipped && start2_on_1.0; + let overlap_along_curve2_1_is_start = overlap_along_curve2_is_flipped && start2_on_1.0; + let overlap_along_curve2_0_is_end = overlap_along_curve2_is_flipped && end2_on_1.0; + let overlap_along_curve2_1_is_end = !overlap_along_curve2_is_flipped && end2_on_1.0; + if flags.intersects(CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS) { + // Add the endpoints as intersections if necessary + if (!start1_on_2.0 && (!overlap_along_curve2_0_is_start && !overlap_along_curve2_0_is_end)) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && start1_on_2.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && overlap_along_curve2_0_is_start) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && overlap_along_curve2_0_is_end) { + intersections.push((overlap_along_curve1.0, overlap_along_curve2.0)); + } + if (!end1_on_2.0 && (!overlap_along_curve2_1_is_start && !overlap_along_curve2_1_is_end)) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && end1_on_2.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && overlap_along_curve2_1_is_start) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && overlap_along_curve2_1_is_end) { + intersections.push((overlap_along_curve1.1, overlap_along_curve2.1)); + } } true // The two curves overlap! @@ -269,17 +293,14 @@ fn add_curve_intersections( ) { let (t1, t2) = if flip { (t2, t1) } else { (t1, t2) }; - // Discard endpoint/endpoint intersections. - if !flags.contains(CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS) - { - let t1_is_an_endpoint = real_lte(t1, 0.0) || real_gte(t1, 1.0); - let t2_is_an_endpoint = real_lte(t2, 0.0) || real_gte(t2, 1.0); - if t1_is_an_endpoint && t2_is_an_endpoint { - return; - } + // Discard endpoint/endpoint intersections if desired + let keep_intersection = + (!real_lte(t1, 0.0) && !real_gte(t1, 1.0) && !real_lte(t2, 0.0) && !real_gte(t2, 1.0)) + || flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && real_lte(t1, 0.0) + || flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && real_gte(t1, 1.0) + || flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && real_lte(t2, 0.0) + || flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && real_gte(t2, 1.0); + if !keep_intersection { + return; } if !flags.contains(CurveIntersectionFlags::KEEP_DUPLICATE_INTERSECTIONS) @@ -722,7 +744,17 @@ mod tests { use super::*; use crate::{CubicBez, QuadBez, Line}; - fn do_test( + fn do_single_test( + curve1: &T, + curve2: &U, + count: usize, + flags: CurveIntersectionFlags, + ) { + let arr1 = curve_curve_intersections(curve1, curve2, flags); + assert_eq!(arr1.len(), count); + } + + fn do_double_test( curve1: &T, curve2: &U, count: usize, @@ -751,7 +783,6 @@ mod tests { found_t = true; } - println!("t for pt on curve: {} {}, expected {} {}", t, pt, t_test, pt_test); assert!(point_is_equal(pt, pt_test)) } assert!(found_t) @@ -773,29 +804,54 @@ mod tests { ) { assert!(t1 < t2); - do_test( + do_double_test( &curve.subsegment(0.0..t2), &curve.subsegment(t1..1.0), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections ); - do_test( + do_double_test( &curve.subsegment(0.0..t2), &curve.subsegment(1.0..t1), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections ); - do_test( + do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t1..t2), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections ); - do_test( + do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t2..t1), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + ); + + do_single_test( + &curve.subsegment(0.0..1.0), + &curve.subsegment(t1..t2), + 1, + CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection + ); + do_single_test( + &curve.subsegment(t1..t2), + &curve.subsegment(0.0..1.0), + 0, + CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection + ); + do_single_test( + &curve.subsegment(0.0..1.0), + &curve.subsegment(t1..t2), + 1, + CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection + ); + do_single_test( + &curve.subsegment(t1..t2), + &curve.subsegment(0.0..1.0), + 0, + CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection ); } @@ -804,61 +860,61 @@ mod tests { test_t(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0))); test_t(&CubicBez::new((0.0, 0.0), (-1.0, 0.0), (1.0, 0.0), (1.0, 0.0))); - do_test( + do_double_test( &CubicBez::new((0.0, 0.0), (0.0, 1.0), (0.0, 1.0), (1.0, 1.0)), &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), 1, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &CubicBez::new((48.0, 84.0), (104.0, 176.0), (190.0, 37.0), (121.0, 75.0)), &CubicBez::new((68.0, 145.0), (74.0, 6.0), (143.0, 197.0), (138.0, 55.0)), 4, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &CubicBez::new((0.0, 0.0), (0.5, 1.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.0, 1.0), (0.5, 0.0), (0.5, 0.0), (1.0, 1.0)), 2, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), 3, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, 0.75), (1.0, 0.25)), 1, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, -0.25), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.5)), 1, CurveIntersectionFlags::NONE, // Standard algorithm ); - do_test( + do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 3, @@ -866,37 +922,49 @@ mod tests { ); // These tests intersect at the endpoints which is only sometimes useful to know about. - do_test( + do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections ); - do_test( + do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections ); - do_test( + do_double_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), + 1, + CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION, // Include curve1 t=0 intersection + ); + do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 2, - CurveIntersectionFlags::KEEP_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections ); - do_test( + do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections ); + do_double_test( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), + 1, + CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION, // Include curve1 t=1 intersection + ); // Test curves that lie exactly along one-another test_overlapping(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), 0.2, 0.8); // (A previous version of the code was returning two practically identical // intersection points here.) - do_test( + do_double_test( &CubicBez::new( (718133.1363092018, 673674.987999388), (-53014.13135835016, 286988.87959900266), @@ -917,7 +985,7 @@ mod tests { // point even though t_min_clip < t_max_clip (because domain1 was small enough to begin with // relative to the small distance between t_min_clip and t_max_clip), and the new curve1 is not // a point (it's split off the old curve1 using t_min_clip < t_max_clip). - do_test( + do_double_test( &CubicBez::new( (423394.5967598548, -91342.7434613118), (333212.450870987, 225564.45711810607), From 84400efa974204044314a0fa992712331e6378f8 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 3 Mar 2022 10:33:38 -0700 Subject: [PATCH 08/20] Added ability to define how much accuracy is needed --- src/curve_intersections.rs | 211 +++++++++++++++++++++++-------------- src/real.rs | 31 ++++-- 2 files changed, 155 insertions(+), 87 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 4fdc9f64..ab690cfd 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -27,14 +27,21 @@ pub fn curve_curve_intersections ArrayVec<[(f64, f64); 9]> { let mut av = ArrayVec::new(); - let overlaps = check_for_overlap(curve1, curve2, flags, &mut av); + // Make sure we don't use too small of an accuracy + let bounds = curve1.bounding_box().union(curve2.bounding_box()); + let accuracy = accuracy.max(epsilon_for_point(Point::new(bounds.max_x(), bounds.max_y()))); + + // Check if both curves lie directly on each other for some span + let overlaps = check_for_overlap(curve1, curve2, flags, accuracy, &mut av); if overlaps { return av; } + // Find intersections add_curve_intersections( curve1, curve2, @@ -47,6 +54,7 @@ pub fn curve_curve_intersections( curve1: &T, curve2: &U, flags: CurveIntersectionFlags, + accuracy: f64, intersections: &mut ArrayVec<[(f64, f64); 9]>, ) -> bool { // If the two curves overlap, there are two main cases we need to account for: @@ -74,8 +83,9 @@ fn check_for_overlap( t_this: f64, curve_this: &T, curve_other: &U, + accuracy: f64, ) -> (bool, f64) { - let t_values_other = point_is_on_curve(curve_this.eval(t_this), curve_other); + let t_values_other = point_is_on_curve(curve_this.eval(t_this), curve_other, accuracy); if !t_values_other.0 { return (false, 0.0); } @@ -107,10 +117,10 @@ fn check_for_overlap( } // First, test which endpoints, if any lie along the other curve - let start1_on_2 = t_value_on_other(0., curve1, curve2); - let end1_on_2 = t_value_on_other(1., curve1, curve2); - let start2_on_1 = t_value_on_other(0., curve2, curve1); - let end2_on_1 = t_value_on_other(1., curve2, curve1); + let start1_on_2 = t_value_on_other(0., curve1, curve2, accuracy); + let end1_on_2 = t_value_on_other(1., curve1, curve2, accuracy); + let start2_on_1 = t_value_on_other(0., curve2, curve1, accuracy); + let end2_on_1 = t_value_on_other(1., curve2, curve1, accuracy); // At least 2 endpoints must lie on the other curve for us to proceed let count = if start1_on_2.0 { 1 } else { 0 } + @@ -155,7 +165,9 @@ fn check_for_overlap( let interval1 = (overlap_along_curve1.1 - overlap_along_curve1.0) / 4.; let interval2 = (overlap_along_curve2.1 - overlap_along_curve2.0) / 4.; for n in 1..3 { - if !point_is_equal(curve1.eval(overlap_along_curve1.0 + (n as f64) * interval1), curve2.eval(overlap_along_curve2.0 + (n as f64) * interval2)) { + let curve1_pt = curve1.eval(overlap_along_curve1.0 + (n as f64) * interval1); + let curve2_pt = curve2.eval(overlap_along_curve2.0 + (n as f64) * interval2); + if !point_is_equal(curve1_pt, curve2_pt, accuracy) { return false; } } @@ -206,6 +218,7 @@ fn add_curve_intersections u32 { call_count += 1; recursion_count += 1; @@ -213,9 +226,9 @@ fn add_curve_intersections, intersections: &mut ArrayVec<[(f64, f64); 9]>, flip: bool, + orig_curve1: &T, + orig_curve2: &U, flags: CurveIntersectionFlags, + accuracy: f64, ) { - let pt = pt_curve.eval(0.5); // We assume pt is curve1 when we add intersections below. let flip = if pt_curve_is_curve1 { flip } else { !flip }; - // Generally speaking |curve| will be quite small at this point, so see if we can get away with - // just sampling here. - - let epsilon = epsilon_for_point(pt); + // Get the mid-point + let pt = pt_curve.eval(0.5); let pt_t = (pt_domain.start + pt_domain.end) * 0.5; - let curve_t = { - let mut t_for_min = 0.0; - let mut min_dist_sq = epsilon; - let tenths = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; - for &t in tenths.iter() { - let d = (pt - curve.eval(t)).hypot2(); - if d < min_dist_sq { - t_for_min = t; - min_dist_sq = d; - } - } - - if real_is_zero(min_dist_sq) { - -1.0 - } else { - domain_value_at_t(curve_domain, t_for_min) - } - }; - - if !real_is_equal(curve_t, -1.0) { - add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections, flags); + // Generally speaking |curve| will be quite small at this point, so see if we can get away with + // just finding one of the points along the curve (even though there may be multiple). + let results = t_along_curve_for_point(pt, curve, accuracy, false); + let curve_t = results.iter() + .fold((-1., f64::MAX), |(t, d_sq), result| { + let result_pos = curve.eval(*result); + let result_d_sq = (pt - result_pos).hypot2(); + if result_d_sq < d_sq { return (*result, result_d_sq); } + (t, d_sq) + }).0; + + if real_is_equal(curve_t, -1.) { + // No intersection found return; } - // If sampling didn't work, try a different approach. - let results = t_along_curve_for_point(pt, curve); - for t in results { - let curve_t = domain_value_at_t(curve_domain, t); - add_intersection(pt_t, pt_curve, curve_t, curve, flip, intersections, flags); - } + let curve_t = domain_value_at_t(curve_domain, curve_t); // convert to proper domain + add_intersection(pt_t, orig_curve1, curve_t, orig_curve2, flip, intersections, flags); } pub fn t_along_curve_for_point( pt: Point, curve: &T, -) -> ArrayVec<[f64; 2]> { + accuracy: f64, + near_pts_only: bool, +) -> ArrayVec<[f64; 4]> { let mut result = ArrayVec::new(); // If both endpoints are approximately close, we only return 0.0. - if point_is_equal(pt, curve.start()) { + if point_is_equal(pt, curve.start(), accuracy) { result.push(0.0); return result; } - if point_is_equal(pt, curve.end()) { + if point_is_equal(pt, curve.end(), accuracy) { result.push(1.0); return result; } - let curve_x_t_params = curve.solve_t_for_x(pt.x); - let curve_y_t_params = curve.solve_t_for_y(pt.y); // We want to coalesce parameters representing the same intersection from the x and y // directions, but the parameter calculations aren't very accurate, so give a little more // leeway there (TODO: this isn't perfect, as you might expect - the dupes that pass here are // currently being detected in add_intersection). + let curve_x_t_params = curve.solve_t_for_x(pt.x); + let curve_y_t_params = curve.solve_t_for_y(pt.y); for params in [curve_x_t_params, curve_y_t_params].iter() { for t in params { let t = *t; - // Note: use point_is_approx below for a little more wiggle room - // since our root finding algorithm is perfect - if !point_is_approx(pt, curve.eval(t), 2.) { + if near_pts_only && !point_is_equal(pt, curve.eval(t), 1e-3) { continue; } + let mut already_found_t = false; for u in &result { if real_is_equal(t, *u) { @@ -534,9 +552,10 @@ pub fn t_along_curve_for_point( t: f64, pt: Point, curve: &T, - result: &mut ArrayVec<[f64; 2]>, + accuracy: f64, + result: &mut ArrayVec<[f64; 4]>, ) -> bool { - if point_is_equal(curve.eval(t), pt) { + if point_is_equal(curve.eval(t), pt, accuracy) { result.push(t); return true; } @@ -544,7 +563,7 @@ pub fn t_along_curve_for_point( } for ex in curve.extrema() { - maybe_add(ex, pt, curve, &mut result); + maybe_add(ex, pt, curve, accuracy, &mut result); } result @@ -553,10 +572,11 @@ pub fn t_along_curve_for_point( fn point_is_on_curve( pt: Point, curve: &T, -) -> (bool, ArrayVec<[f64; 2]>) { - let t_values = t_along_curve_for_point(pt, curve); + accuracy: f64, +) -> (bool, ArrayVec<[f64; 4]>) { + let t_values = t_along_curve_for_point(pt, curve, accuracy, true); for t in &t_values { - if point_is_equal(pt, curve.eval(*t)) { + if point_is_equal(pt, curve.eval(*t), accuracy) { return (true, t_values.clone()); // Point lies along curve } } @@ -572,6 +592,8 @@ fn add_intersection( intersections: &mut ArrayVec<[(f64, f64); 9]>, flags: CurveIntersectionFlags, ) { + // Note: orig_curve1 and orig_curve2 are actually swapped with flip == true + let (t1_flipped, t2_flipped) = (t1, t2); let (t1, t2) = if flip { (t2, t1) } else { (t1, t2) }; // Discard endpoint/endpoint intersections if desired @@ -587,20 +609,20 @@ fn add_intersection( if !flags.contains(CurveIntersectionFlags::KEEP_DUPLICATE_INTERSECTIONS) { - // (This should probably depend in some way on how large our input coefficients are.) - let epsilon = 1e-3; + let (pt1, pt2) = (orig_curve1.eval(t1_flipped), orig_curve2.eval(t2_flipped)); // We can get repeated intersections when we split a curve at an intersection point, or when // two curves intersect at a point where the curves are very close together, or when the fat // line process breaks down. for i in 0..intersections.len() { let (old_t1, old_t2) = intersections[i]; + let (old_t1_flipped, old_t2_flipped) = if flip { (old_t2, old_t1) } else { (old_t1, old_t2) }; + let (old_pt1, old_pt2) = (orig_curve1.eval(old_t1_flipped), orig_curve2.eval(old_t2_flipped)); + // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" // duplicate intersection, so always keep the one that minimizes sample distance. - if (t1 - old_t1).abs() < epsilon && (t2 - old_t2).abs() < epsilon { - let cur_dist = (orig_curve1.eval(old_t1) - orig_curve2.eval(old_t2)).hypot2(); - let new_dist = (orig_curve1.eval(t1) - orig_curve2.eval(t2)).hypot2(); - if new_dist < cur_dist { + if point_is_equal(pt1, old_pt1, 1e-3) && point_is_equal(pt2, old_pt2, 1e-3) { + if (pt1 - pt2).hypot2() < (old_pt1 - old_pt2).hypot2() { intersections[i] = (t1, t2); } return; @@ -749,8 +771,9 @@ mod tests { curve2: &U, count: usize, flags: CurveIntersectionFlags, + accuracy: f64, ) { - let arr1 = curve_curve_intersections(curve1, curve2, flags); + let arr1 = curve_curve_intersections(curve1, curve2, flags, accuracy); assert_eq!(arr1.len(), count); } @@ -759,9 +782,10 @@ mod tests { curve2: &U, count: usize, flags: CurveIntersectionFlags, + accuracy: f64, ) { - let arr1 = curve_curve_intersections(curve1, curve2, flags); - let arr2 = curve_curve_intersections(curve2, curve1, flags); + let arr1 = curve_curve_intersections(curve1, curve2, flags, accuracy); + let arr2 = curve_curve_intersections(curve2, curve1, flags, accuracy); assert_eq!(arr1.len(), count); assert_eq!(arr2.len(), count); } @@ -769,21 +793,22 @@ mod tests { fn test_t_along_curve( curve: &T, t_test: f64, + accuracy: f64, ) { let pt_test = curve.eval(t_test); - let t_values = t_along_curve_for_point(pt_test, curve); + let t_values = t_along_curve_for_point(pt_test, curve, accuracy, true); assert!(!t_values.is_empty()); let mut found_t = false; for t in &t_values { let pt = curve.eval(*t); - if real_is_equal(*t, t_test) + if real_is_equal(*t, t_test) // ^^^ REPLACE WITH CMP { found_t = true; } - assert!(point_is_equal(pt, pt_test)) + assert!(point_is_equal(pt, pt_test, 1e-3)) } assert!(found_t) } @@ -793,7 +818,7 @@ mod tests { ) { let tenths = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; for &t in tenths.iter() { - test_t_along_curve(curve, t); + test_t_along_curve(curve, t, 0.); } } @@ -801,6 +826,7 @@ mod tests { curve: &T, t1: f64, t2: f64, + accuracy: f64, ) { assert!(t1 < t2); @@ -809,24 +835,28 @@ mod tests { &curve.subsegment(t1..1.0), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + accuracy ); do_double_test( &curve.subsegment(0.0..t2), &curve.subsegment(1.0..t1), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + accuracy ); do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t1..t2), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + accuracy ); do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t2..t1), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + accuracy ); do_single_test( @@ -834,24 +864,28 @@ mod tests { &curve.subsegment(t1..t2), 1, CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection + accuracy ); do_single_test( &curve.subsegment(t1..t2), &curve.subsegment(0.0..1.0), 0, CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection + accuracy ); do_single_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t1..t2), 1, CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection + accuracy ); do_single_test( &curve.subsegment(t1..t2), &curve.subsegment(0.0..1.0), 0, CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection + accuracy ); } @@ -865,60 +899,70 @@ mod tests { &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), 1, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((48.0, 84.0), (104.0, 176.0), (190.0, 37.0), (121.0, 75.0)), &CubicBez::new((68.0, 145.0), (74.0, 6.0), (143.0, 197.0), (138.0, 55.0)), 4, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.5, 1.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.0, 1.0), (0.5, 0.0), (0.5, 0.0), (1.0, 1.0)), 2, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), 3, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, 0.75), (1.0, 0.25)), 1, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, -0.25), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.5)), 1, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 3, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); // These tests intersect at the endpoints which is only sometimes useful to know about. @@ -927,40 +971,51 @@ mod tests { &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 1, CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION, // Include curve1 t=0 intersection + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 1, CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION, // Include curve1 t=1 intersection + 0., // Tightest possible accuracy ); // Test curves that lie exactly along one-another - test_overlapping(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), 0.2, 0.8); + test_overlapping( + &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), + 0.2, + 0.8, + 0. + ); // (A previous version of the code was returning two practically identical // intersection points here.) @@ -979,6 +1034,7 @@ mod tests { ), 1, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); // On these curves the algorithm runs to a state at which the new clipped domain1 becomes a @@ -1000,6 +1056,7 @@ mod tests { ), 3, CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy ); } } diff --git a/src/real.rs b/src/real.rs index b8d06d44..96b6ca08 100644 --- a/src/real.rs +++ b/src/real.rs @@ -2,6 +2,11 @@ use crate::{Point}; +pub fn get_epsilon(num: f64, accuracy: f64) -> f64 { + // Use the provided number as the epsilon unless it's too small + accuracy.abs().max(epsilon_for_value(num)) +} + /// If we're comparing numbers, our epsilon should depend on how big the number /// is. This function returns an epsilon appropriate for the size the largest /// number. @@ -69,16 +74,16 @@ pub fn real_gte(num1: f64, num2: f64) -> bool } /// Compare if two points are approximately equal -pub fn point_is_equal(pt1: Point, pt2: Point) -> bool +pub fn point_is_equal(pt1: Point, pt2: Point, accuracy: f64) -> bool { - point_is_approx(pt1, pt2, 1.) + point_is_approx(pt1, pt2, accuracy, 1.) } /// Compare if two points are approximately equal to a smaller degree than /// point_is_equal -pub fn point_is_approx(pt1: Point, pt2: Point, scale: f64) -> bool +pub fn point_is_approx(pt1: Point, pt2: Point, accuracy: f64, scale: f64) -> bool { - real_is_approx(pt1.x, pt2.x, scale) && real_is_approx(pt1.y, pt2.y, scale) + (pt2 - pt1).hypot() <= get_epsilon((pt1.x).max(pt1.y).max(pt2.x).max(pt2.y), accuracy * scale) } #[cfg(test)] @@ -124,11 +129,17 @@ mod tests { #[test] fn test_point_comparisons() { - assert!(point_is_equal(Point::new(0., 0.), Point::new(0., 0.))); - assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.))); - assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY))); - assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1.))); - assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY))); - assert!(point_is_equal(Point::new(0., 0.), Point::new(f64::EPSILON, f64::EPSILON))); + assert!(point_is_equal(Point::new(0., 0.), Point::new(0., 0.), 0.)); + assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.), 0.)); + assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), 0.)); + assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1.), 0.)); + assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), 0.)); + + let epsilon = 0.1; + assert!(point_is_equal(Point::new(0., 0.), Point::new(0. + epsilon, 0.), epsilon)); + assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000. + epsilon), epsilon)); + assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), epsilon)); + assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1. + epsilon), epsilon)); + assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), epsilon)); } } \ No newline at end of file From 7398c236d09af28cb34d1a9dc52a4818bdbb04b3 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 3 Mar 2022 16:08:33 -0700 Subject: [PATCH 09/20] Removed real comparisons, only point comparisons remailn --- src/curve_intersections.rs | 41 +++++------ src/lib.rs | 2 - src/point.rs | 47 ++++++++++++ src/real.rs | 145 ------------------------------------- 4 files changed, 67 insertions(+), 168 deletions(-) delete mode 100644 src/real.rs diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index ab690cfd..5285400a 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,7 +1,6 @@ use crate::ParamCurveBezierClipping; use crate::{ParamCurve, ParamCurveExtrema}; use crate::{Point, Rect, Affine}; -use crate::real::*; use arrayvec::ArrayVec; use std::ops::Range; @@ -33,7 +32,7 @@ pub fn curve_curve_intersections( for n in 1..3 { let curve1_pt = curve1.eval(overlap_along_curve1.0 + (n as f64) * interval1); let curve2_pt = curve2.eval(overlap_along_curve2.0 + (n as f64) * interval2); - if !point_is_equal(curve1_pt, curve2_pt, accuracy) { + if !Point::is_near(curve1_pt, curve2_pt, accuracy) { return false; } } @@ -228,7 +227,7 @@ fn add_curve_intersections accuracy { call_count = add_curve_intersections( curve2, curve1, @@ -482,7 +481,7 @@ fn add_point_curve_intersection( let mut result = ArrayVec::new(); // If both endpoints are approximately close, we only return 0.0. - if point_is_equal(pt, curve.start(), accuracy) { + if Point::is_near(pt, curve.start(), accuracy) { result.push(0.0); return result; } - if point_is_equal(pt, curve.end(), accuracy) { + if Point::is_near(pt, curve.end(), accuracy) { result.push(1.0); return result; } @@ -519,13 +518,13 @@ pub fn t_along_curve_for_point( for t in params { let t = *t; - if near_pts_only && !point_is_equal(pt, curve.eval(t), 1e-3) { + if near_pts_only && !Point::is_near(pt, curve.eval(t), 1e-3) { continue; } let mut already_found_t = false; for u in &result { - if real_is_equal(t, *u) { + if t.eq(u) { already_found_t = true; break; } @@ -555,7 +554,7 @@ pub fn t_along_curve_for_point( accuracy: f64, result: &mut ArrayVec<[f64; 4]>, ) -> bool { - if point_is_equal(curve.eval(t), pt, accuracy) { + if Point::is_near(curve.eval(t), pt, accuracy) { result.push(t); return true; } @@ -576,7 +575,7 @@ fn point_is_on_curve( ) -> (bool, ArrayVec<[f64; 4]>) { let t_values = t_along_curve_for_point(pt, curve, accuracy, true); for t in &t_values { - if point_is_equal(pt, curve.eval(*t), accuracy) { + if Point::is_near(pt, curve.eval(*t), accuracy) { return (true, t_values.clone()); // Point lies along curve } } @@ -598,11 +597,11 @@ fn add_intersection( // Discard endpoint/endpoint intersections if desired let keep_intersection = - (!real_lte(t1, 0.0) && !real_gte(t1, 1.0) && !real_lte(t2, 0.0) && !real_gte(t2, 1.0)) - || flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && real_lte(t1, 0.0) - || flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && real_gte(t1, 1.0) - || flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && real_lte(t2, 0.0) - || flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && real_gte(t2, 1.0); + (t1 > 0.0 && t1 < 1.0 && t2 > 0.0 && t2 < 1.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && t1 <= 0.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && t1 >= 1.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && t2 <= 0.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && t2 >= 1.0); if !keep_intersection { return; } @@ -621,7 +620,7 @@ fn add_intersection( // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" // duplicate intersection, so always keep the one that minimizes sample distance. - if point_is_equal(pt1, old_pt1, 1e-3) && point_is_equal(pt2, old_pt2, 1e-3) { + if Point::is_near(pt1, old_pt1, 1e-3) && Point::is_near(pt2, old_pt2, 1e-3) { if (pt1 - pt2).hypot2() < (old_pt1 - old_pt2).hypot2() { intersections[i] = (t1, t2); } @@ -803,12 +802,12 @@ mod tests { for t in &t_values { let pt = curve.eval(*t); - if real_is_equal(*t, t_test) // ^^^ REPLACE WITH CMP + if (*t - t_test) <= 1e-9 { found_t = true; } - assert!(point_is_equal(pt, pt_test, 1e-3)) + assert!(Point::is_near(pt, pt_test, 1e-3)) } assert!(found_t) } diff --git a/src/lib.rs b/src/lib.rs index f131e667..cb897da7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -102,7 +102,6 @@ mod size; mod svg; mod translate_scale; mod vec2; -mod real; pub use crate::affine::*; pub use crate::arc::*; @@ -126,4 +125,3 @@ pub use crate::size::*; pub use crate::svg::*; pub use crate::translate_scale::*; pub use crate::vec2::*; -pub use crate::real::*; diff --git a/src/point.rs b/src/point.rs index 7c94016e..0c6bc975 100644 --- a/src/point.rs +++ b/src/point.rs @@ -169,6 +169,37 @@ impl Point { pub fn is_nan(self) -> bool { self.x.is_nan() || self.y.is_nan() } + + /// If we're comparing distances between samples of curves, our epsilon should + /// depend on how big the points we're comparing are. This function returns an + /// epsilon appropriate for the size of pt. + pub fn epsilon(pt: Point, accuracy: f64) -> f64 { + // Get max number in the point + let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); + + // These values have been obtained experimentally and can be changed if + // necessary. + accuracy.max(max.abs() * 0.0000000001).max(f64::EPSILON) + } + + /// Compare if two points are approximately equal + pub fn is_near(pt1: Point, pt2: Point, accuracy: f64) -> bool + { + Point::is_approx(pt1, pt2, accuracy, 1.) + } + + /// Compare if two points are approximately equal to some scaled accuracy + pub fn is_approx(pt1: Point, pt2: Point, accuracy: f64, scale: f64) -> bool + { + if pt1.x.is_infinite() || pt1.y.is_infinite() || pt2.x.is_infinite() || pt2.y.is_infinite() + { + return pt1 == pt2; + } + + let accuracy_scaled = accuracy * scale; + let epsilon = Point::epsilon(pt1, accuracy_scaled).max(Point::epsilon(pt2, accuracy_scaled)); + (pt2 - pt1).hypot() <= epsilon + } } impl From<(f64, f64)> for Point { @@ -309,6 +340,22 @@ mod tests { let p = Point::new(0.12345, 9.87654); assert_eq!(format!("{:.2}", p), "(0.12, 9.88)"); } + + #[test] + fn test_point_comparisons() { + assert!(Point::is_near(Point::new(0., 0.), Point::new(0., 0.), 0.)); + assert!(Point::is_near(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.), 0.)); + assert!(Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), 0.)); + assert!(!Point::is_near(Point::new(0., 0.), Point::new(1., 1.), 0.)); + assert!(!Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), 0.)); + + let epsilon = 0.1; + assert!(Point::is_near(Point::new(0., 0.), Point::new(0. + epsilon, 0.), epsilon)); + assert!(Point::is_near(Point::new(1000000., 1000000.), Point::new(1000000., 1000000. + epsilon), epsilon)); + assert!(Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), epsilon)); + assert!(!Point::is_near(Point::new(0., 0.), Point::new(1., 1. + epsilon), epsilon)); + assert!(!Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), epsilon)); + } } #[cfg(feature = "mint")] diff --git a/src/real.rs b/src/real.rs deleted file mode 100644 index 96b6ca08..00000000 --- a/src/real.rs +++ /dev/null @@ -1,145 +0,0 @@ -//! A utility for comparing real values - -use crate::{Point}; - -pub fn get_epsilon(num: f64, accuracy: f64) -> f64 { - // Use the provided number as the epsilon unless it's too small - accuracy.abs().max(epsilon_for_value(num)) -} - -/// If we're comparing numbers, our epsilon should depend on how big the number -/// is. This function returns an epsilon appropriate for the size the largest -/// number. -pub fn epsilon_for_value(num: f64) -> f64 { - // These values have been obtained experimentally and can be changed if - // necessary. - (num.abs() * 0.0000000001).max(f64::EPSILON) -} - -/// If we're comparing distances between samples of curves, our epsilon should -/// depend on how big the points we're comparing are. This function returns an -/// epsilon appropriate for the size of pt. -pub fn epsilon_for_point(pt: Point) -> f64 { - let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); - epsilon_for_value(max) -} - -/// Compare if two real values are approximately equal -pub fn real_is_equal(num1: f64, num2: f64) -> bool -{ - real_is_approx(num1, num2, 1.) -} - -/// Compare if two real values are approximately equal to a smaller degree than -/// real_is_equal -pub fn real_is_approx(num1: f64, num2: f64, scale: f64) -> bool -{ - if num1 == f64::INFINITY || num2 == f64::INFINITY - { - return num1 == num2; // Do exact comparison - } - - // Do approximate comparison - (num1 - num2).abs() <= (scale * epsilon_for_value(num1.max(num2))) -} - -/// Compare if a real value is approximately zero -pub fn real_is_zero(num1: f64) -> bool -{ - real_is_equal(num1, 0.) -} - -/// Compare if num1 < num2 and not approximately equal -pub fn real_lt(num1: f64, num2: f64) -> bool -{ - (num1 < num2) && !real_is_equal(num1, num2) -} - -/// Compare if num1 < num2 or approximately equal -pub fn real_lte(num1: f64, num2: f64) -> bool -{ - (num1 < num2) || real_is_equal(num1, num2) -} - -/// Compare if num1 > num2 and not approximately equal -pub fn real_gt(num1: f64, num2: f64) -> bool -{ - (num1 > num2) && !real_is_equal(num1, num2) -} - -/// Compare if num1 > num2 or approximately equal -pub fn real_gte(num1: f64, num2: f64) -> bool -{ - (num1 > num2) || real_is_equal(num1, num2) -} - -/// Compare if two points are approximately equal -pub fn point_is_equal(pt1: Point, pt2: Point, accuracy: f64) -> bool -{ - point_is_approx(pt1, pt2, accuracy, 1.) -} - -/// Compare if two points are approximately equal to a smaller degree than -/// point_is_equal -pub fn point_is_approx(pt1: Point, pt2: Point, accuracy: f64, scale: f64) -> bool -{ - (pt2 - pt1).hypot() <= get_epsilon((pt1.x).max(pt1.y).max(pt2.x).max(pt2.y), accuracy * scale) -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn test_real_comparisons() { - // Equals - assert!(real_is_equal(1., 1.)); - assert!(real_is_equal(1000000., 1000000.)); - assert!(real_is_equal(f64::INFINITY, f64::INFINITY)); - assert!(!real_is_equal(f64::INFINITY, f64::NEG_INFINITY)); - assert!(real_is_equal(f64::EPSILON, f64::EPSILON)); - assert!(!real_is_equal(f64::EPSILON, -f64::EPSILON)); - assert!(real_is_equal(f64::EPSILON, 0.)); - assert!(real_is_equal(0., f64::EPSILON)); - - // Less than - assert!(real_lte(1., 1.)); - assert!(real_lte(1000000., 1000000.)); - assert!(real_lte(f64::INFINITY, f64::INFINITY)); - assert!(!real_lt(1., 1.)); - assert!(!real_lt(1000000., 1000000.)); - assert!(!real_lt(f64::INFINITY, f64::INFINITY)); - assert!(real_lt(0., f64::INFINITY)); - assert!(real_lt(f64::NEG_INFINITY, 0.)); - assert!(real_lt(-1., 1.)); - assert!(real_lt(f64::NEG_INFINITY, f64::INFINITY)); - - // Greater than - assert!(real_gte(1., 1.)); - assert!(real_gte(1000000., 1000000.)); - assert!(real_gte(f64::INFINITY, f64::INFINITY)); - assert!(!real_gt(1., 1.)); - assert!(!real_gt(1000000., 1000000.)); - assert!(!real_gt(f64::INFINITY, f64::INFINITY)); - assert!(real_gt(f64::INFINITY, 0.)); - assert!(real_gt(0., f64::NEG_INFINITY)); - assert!(real_gt(1., -1.)); - assert!(real_gt(f64::INFINITY, f64::NEG_INFINITY)); - } - - #[test] - fn test_point_comparisons() { - assert!(point_is_equal(Point::new(0., 0.), Point::new(0., 0.), 0.)); - assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.), 0.)); - assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), 0.)); - assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1.), 0.)); - assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), 0.)); - - let epsilon = 0.1; - assert!(point_is_equal(Point::new(0., 0.), Point::new(0. + epsilon, 0.), epsilon)); - assert!(point_is_equal(Point::new(1000000., 1000000.), Point::new(1000000., 1000000. + epsilon), epsilon)); - assert!(point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), epsilon)); - assert!(!point_is_equal(Point::new(0., 0.), Point::new(1., 1. + epsilon), epsilon)); - assert!(!point_is_equal(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), epsilon)); - } -} \ No newline at end of file From cb5a6725b92ba96ae361f910afcc2ad08e2f1bcf Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Fri, 4 Mar 2022 10:16:25 -0700 Subject: [PATCH 10/20] Reduced arclen epsilon to improve performance --- src/curve_intersections.rs | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 5285400a..15a5ad79 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -225,7 +225,7 @@ fn add_curve_intersections accuracy { + if curve2_arclen > accuracy { call_count = add_curve_intersections( curve2, curve1, From fc0bdaf20cff8fb80283631cfb8fc817d296a5cd Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 09:48:44 -0600 Subject: [PATCH 11/20] Fixed merge with master --- src/common.rs | 2 +- src/cubicbez.rs | 4 ++-- src/curve_intersections.rs | 16 ++++++++-------- src/param_curve.rs | 4 ++-- src/param_curve_clipping.rs | 12 ++++++------ 5 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/common.rs b/src/common.rs index 3c7f586e..dc53cedd 100644 --- a/src/common.rs +++ b/src/common.rs @@ -148,7 +148,7 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { /// Find real roots of linear equation. /// /// Return values of x for which c0 + c1 x = 0. -pub fn solve_linear(c0: f64, c1: f64) -> ArrayVec<[f64; 1]> { +pub fn solve_linear(c0: f64, c1: f64) -> ArrayVec { let mut result = ArrayVec::new(); let root = -c0 / c1; if root.is_finite() { diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 2fbc2452..91f89436 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -608,7 +608,7 @@ mod tests { cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez, }; - use arrayvec::{Array, ArrayVec}; + use arrayvec::ArrayVec; #[test] fn cubicbez_deriv() { @@ -892,7 +892,7 @@ mod tests { use crate::param_curve::ParamCurveBezierClipping; #[test] fn solve_t_for_xy() { - fn verify>(mut roots: ArrayVec, expected: &[f64]) { + fn verify(mut roots: ArrayVec, expected: &[f64]) { assert_eq!(expected.len(), roots.len()); let epsilon = 1e-6; roots.sort_by(|a, b| a.partial_cmp(b).unwrap()); diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 15a5ad79..3b8accbe 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -27,7 +27,7 @@ pub fn curve_curve_intersections ArrayVec<[(f64, f64); 9]> { +) -> ArrayVec<(f64, f64), 9> { let mut av = ArrayVec::new(); // Make sure we don't use too small of an accuracy @@ -64,7 +64,7 @@ fn check_for_overlap( curve2: &U, flags: CurveIntersectionFlags, accuracy: f64, - intersections: &mut ArrayVec<[(f64, f64); 9]>, + intersections: &mut ArrayVec<(f64, f64), 9>, ) -> bool { // If the two curves overlap, there are two main cases we need to account for: // 1) One curve is completely along some portion of the other curve, thus @@ -210,7 +210,7 @@ fn add_curve_intersections, domain2: &Range, - intersections: &mut ArrayVec<[(f64, f64); 9]>, + intersections: &mut ArrayVec<(f64, f64), 9>, flip: bool, mut recursion_count: u32, mut call_count: u32, @@ -458,7 +458,7 @@ fn add_point_curve_intersection, curve_domain: &Range, - intersections: &mut ArrayVec<[(f64, f64); 9]>, + intersections: &mut ArrayVec<(f64, f64), 9>, flip: bool, orig_curve1: &T, orig_curve2: &U, @@ -497,7 +497,7 @@ pub fn t_along_curve_for_point( curve: &T, accuracy: f64, near_pts_only: bool, -) -> ArrayVec<[f64; 4]> { +) -> ArrayVec { let mut result = ArrayVec::new(); // If both endpoints are approximately close, we only return 0.0. @@ -554,7 +554,7 @@ pub fn t_along_curve_for_point( pt: Point, curve: &T, accuracy: f64, - result: &mut ArrayVec<[f64; 4]>, + result: &mut ArrayVec, ) -> bool { if Point::is_near(curve.eval(t), pt, accuracy) { result.push(t); @@ -574,7 +574,7 @@ fn point_is_on_curve( pt: Point, curve: &T, accuracy: f64, -) -> (bool, ArrayVec<[f64; 4]>) { +) -> (bool, ArrayVec) { let t_values = t_along_curve_for_point(pt, curve, accuracy, true); for t in &t_values { if Point::is_near(pt, curve.eval(*t), accuracy) { @@ -590,7 +590,7 @@ fn add_intersection( t2: f64, orig_curve2: &U, flip: bool, - intersections: &mut ArrayVec<[(f64, f64); 9]>, + intersections: &mut ArrayVec<(f64, f64), 9>, flags: CurveIntersectionFlags, ) { // Note: orig_curve1 and orig_curve2 are actually swapped with flip == true diff --git a/src/param_curve.rs b/src/param_curve.rs index a2533d8a..cc13e3d2 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -223,9 +223,9 @@ pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveDeriv + ParamCurveExt } /// Find the time `t` at which the curve has the given x value - fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]>; + fn solve_t_for_x(&self, x: f64) -> ArrayVec; /// Find the time `t` at which the curve has the given x value - fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]>; + fn solve_t_for_y(&self, y: f64) -> ArrayVec; /// Returns the upper and lower convex hull fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); /// Returns the minimum and maximum distances of the "fat line" enclosing this curve diff --git a/src/param_curve_clipping.rs b/src/param_curve_clipping.rs index 1d3c18a6..819ccca9 100644 --- a/src/param_curve_clipping.rs +++ b/src/param_curve_clipping.rs @@ -27,7 +27,7 @@ fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { } impl ParamCurveBezierClipping for Line { - fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { if (self.p0.x - self.p1.x).abs() < f64::EPSILON { return ArrayVec::new(); } @@ -38,7 +38,7 @@ impl ParamCurveBezierClipping for Line { .filter(|&t| (0.0..=1.0).contains(&t)) .collect() } - fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_y(&self, y: f64) -> ArrayVec { if (self.p0.y - self.p1.y).abs() < f64::EPSILON { return ArrayVec::new(); } @@ -76,7 +76,7 @@ impl ParamCurveBezierClipping for Line { } impl ParamCurveBezierClipping for QuadBez { - fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { if self.is_linear(f64::EPSILON) && (self.p0.x - self.p2.x).abs() < f64::EPSILON { return ArrayVec::new(); } @@ -87,7 +87,7 @@ impl ParamCurveBezierClipping for QuadBez { .filter(|&t| (0.0..=1.0).contains(&t)) .collect() } - fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_y(&self, y: f64) -> ArrayVec { if self.is_linear(f64::EPSILON) && (self.p0.y - self.p2.y).abs() < f64::EPSILON { return ArrayVec::new(); } @@ -146,7 +146,7 @@ impl ParamCurveBezierClipping for QuadBez { } impl ParamCurveBezierClipping for CubicBez { - fn solve_t_for_x(&self, x: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { return ArrayVec::new(); } @@ -157,7 +157,7 @@ impl ParamCurveBezierClipping for CubicBez { .filter(|&t| (0.0..=1.0).contains(&t)) .collect() } - fn solve_t_for_y(&self, y: f64) -> ArrayVec<[f64; 3]> { + fn solve_t_for_y(&self, y: f64) -> ArrayVec { if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { return ArrayVec::new(); } From 1a2cb70407443e32f5bf9a3ac85cbe2f91c758b9 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 11:54:10 -0600 Subject: [PATCH 12/20] Fixed formatting --- src/cubicbez.rs | 3 +- src/curve_intersections.rs | 244 ++++++++++++++++++++++-------------- src/param_curve.rs | 4 +- src/param_curve_clipping.rs | 9 +- src/point.rs | 57 +++++++-- src/quadbez.rs | 8 +- 6 files changed, 205 insertions(+), 120 deletions(-) diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 91f89436..d517a9f1 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -6,8 +6,7 @@ use crate::MAX_EXTREMA; use crate::{Line, QuadSpline, Vec2}; use arrayvec::ArrayVec; -use crate::common::GAUSS_LEGENDRE_COEFFS_9; -use crate::common::{solve_quadratic}; +use crate::common::{solve_quadratic, GAUSS_LEGENDRE_COEFFS_9}; use crate::{ Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 3b8accbe..a933b1ba 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,6 +1,5 @@ use crate::ParamCurveBezierClipping; -use crate::{ParamCurve, ParamCurveExtrema}; -use crate::{Point, Rect, Affine}; +use crate::{Affine, ParamCurve, ParamCurveExtrema, Point, Rect}; use arrayvec::ArrayVec; use std::ops::Range; @@ -91,7 +90,7 @@ fn check_for_overlap( assert!(t_values_other.1.len() >= 1); // This has to be true or we really messed up! if t_values_other.1.len() == 1 { - return (true, t_values_other.1[0]) + return (true, t_values_other.1[0]); } // Things get kind of complicated if we have more than one t value @@ -103,7 +102,7 @@ fn check_for_overlap( // Find the t-value with the closest derivative let mut deriv = Point::new(f64::MAX, f64::MAX); let mut index = 0; - for (i, t) in t_values_other.1.iter().enumerate() { + for (i, t) in t_values_other.1.iter().enumerate() { let curve_deriv_at_t_other = curve_deriv_other.eval(*t); let d = (curve_deriv_at_t_other - deriv).hypot2(); let d_compare = (curve_deriv_at_t_this - deriv).hypot2(); @@ -122,10 +121,10 @@ fn check_for_overlap( let end2_on_1 = t_value_on_other(1., curve2, curve1, accuracy); // At least 2 endpoints must lie on the other curve for us to proceed - let count = if start1_on_2.0 { 1 } else { 0 } + - if end1_on_2. 0 { 1 } else { 0 } + - if start2_on_1.0 { 1 } else { 0 } + - if end2_on_1.0 { 1 } else { 0 }; + let count = if start1_on_2.0 { 1 } else { 0 } + + if end1_on_2.0 { 1 } else { 0 } + + if start2_on_1.0 { 1 } else { 0 } + + if end2_on_1.0 { 1 } else { 0 }; if count < 2 { return false; // No overlap } @@ -139,13 +138,25 @@ fn check_for_overlap( // checked the endpoints, we just need to check 3 intermediate positions. fn get_projection(t_value: (bool, f64), default_value: f64) -> f64 { - if t_value.0 { t_value.1 } else { default_value } + if t_value.0 { + t_value.1 + } else { + default_value + } } // Determine the range of curve1 that lies along curve2 let overlap_along_curve1 = ( - if start1_on_2.0 { 0. } else { (get_projection(start2_on_1, 1.)).min(get_projection(end2_on_1, 1.)) }, - if end1_on_2.0 { 1. } else { (get_projection(start2_on_1, 0.)).max(get_projection(end2_on_1, 0.)) }, + if start1_on_2.0 { + 0. + } else { + (get_projection(start2_on_1, 1.)).min(get_projection(end2_on_1, 1.)) + }, + if end1_on_2.0 { + 1. + } else { + (get_projection(start2_on_1, 0.)).max(get_projection(end2_on_1, 0.)) + }, ); assert!(overlap_along_curve1.0 <= overlap_along_curve1.1); @@ -154,10 +165,28 @@ fn check_for_overlap( // overlap_along_curve1.0. // Note: It's not guaranteed that overlap_along_curve2.0 < overlap_along_curve2.1 let overlap_along_curve2 = ( - if start1_on_2.0 { start1_on_2.1 } - else { let p1 = get_projection(start2_on_1, 1.); let p2 = get_projection(end2_on_1, 1.); if p1 < p2 { 0. } else { 1. } }, - if end1_on_2.0 { end1_on_2.1 } - else { let p1 = get_projection(start2_on_1, 0.); let p2 = get_projection(end2_on_1, 0.); if p1 < p2 { 1. } else { 0. } }, + if start1_on_2.0 { + start1_on_2.1 + } else { + let p1 = get_projection(start2_on_1, 1.); + let p2 = get_projection(end2_on_1, 1.); + if p1 < p2 { + 0. + } else { + 1. + } + }, + if end1_on_2.0 { + end1_on_2.1 + } else { + let p1 = get_projection(start2_on_1, 0.); + let p2 = get_projection(end2_on_1, 0.); + if p1 < p2 { + 1. + } else { + 0. + } + }, ); // Check the 3 intermediate locations @@ -179,15 +208,22 @@ fn check_for_overlap( if flags.intersects(CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS) { // Add the endpoints as intersections if necessary if (!start1_on_2.0 && (!overlap_along_curve2_0_is_start && !overlap_along_curve2_0_is_end)) - || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && start1_on_2.0) - || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && overlap_along_curve2_0_is_start) - || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && overlap_along_curve2_0_is_end) { + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) + && start1_on_2.0) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) + && overlap_along_curve2_0_is_start) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) + && overlap_along_curve2_0_is_end) + { intersections.push((overlap_along_curve1.0, overlap_along_curve2.0)); } if (!end1_on_2.0 && (!overlap_along_curve2_1_is_start && !overlap_along_curve2_1_is_end)) || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && end1_on_2.0) - || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && overlap_along_curve2_1_is_start) - || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && overlap_along_curve2_1_is_end) { + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) + && overlap_along_curve2_1_is_start) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) + && overlap_along_curve2_1_is_end) + { intersections.push((overlap_along_curve1.1, overlap_along_curve2.1)); } } @@ -227,7 +263,11 @@ fn add_curve_intersections { - match clip_perp { - Some((min_perp, max_perp)) => (min_norm.max(min_perp), max_norm.min(max_perp)), - None => (min_norm, max_norm), - } - } - None => { - match clip_perp { - Some((min_perp, max_perp)) => (min_perp, max_perp), - None => return call_count, - } + Some((min_norm, max_norm)) => match clip_perp { + Some((min_perp, max_perp)) => (min_norm.max(min_perp), max_norm.min(max_perp)), + None => (min_norm, max_norm), + }, + None => match clip_perp { + Some((min_perp, max_perp)) => (min_perp, max_perp), + None => return call_count, }, }; // t_min_clip and t_max_clip are (0, 1)-based, so project them back to get the new restricted // range: - let new_domain1 = &(domain_value_at_t(domain1, t_min_clip)..domain_value_at_t(domain1, t_max_clip)); + let new_domain1 = + &(domain_value_at_t(domain1, t_min_clip)..domain_value_at_t(domain1, t_max_clip)); // Reduce curve1 to the part that might intersect curve2. let curve1 = &orig_curve1.subsegment(new_domain1.clone()); @@ -319,7 +356,7 @@ fn add_curve_intersections( @@ -520,7 +569,7 @@ pub fn t_along_curve_for_point( for t in params { let t = *t; - if near_pts_only && !Point::is_near(pt, curve.eval(t), 1e-3) { + if near_pts_only && !Point::is_near(pt, curve.eval(t), 1e-3) { continue; } @@ -596,10 +645,9 @@ fn add_intersection( // Note: orig_curve1 and orig_curve2 are actually swapped with flip == true let (t1_flipped, t2_flipped) = (t1, t2); let (t1, t2) = if flip { (t2, t1) } else { (t1, t2) }; - + // Discard endpoint/endpoint intersections if desired - let keep_intersection = - (t1 > 0.0 && t1 < 1.0 && t2 > 0.0 && t2 < 1.0) + let keep_intersection = (t1 > 0.0 && t1 < 1.0 && t2 > 0.0 && t2 < 1.0) || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION) && t1 <= 0.0) || (flags.contains(CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION) && t1 >= 1.0) || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION) && t2 <= 0.0) @@ -608,8 +656,7 @@ fn add_intersection( return; } - if !flags.contains(CurveIntersectionFlags::KEEP_DUPLICATE_INTERSECTIONS) - { + if !flags.contains(CurveIntersectionFlags::KEEP_DUPLICATE_INTERSECTIONS) { let (pt1, pt2) = (orig_curve1.eval(t1_flipped), orig_curve2.eval(t2_flipped)); // We can get repeated intersections when we split a curve at an intersection point, or when @@ -617,8 +664,15 @@ fn add_intersection( // line process breaks down. for i in 0..intersections.len() { let (old_t1, old_t2) = intersections[i]; - let (old_t1_flipped, old_t2_flipped) = if flip { (old_t2, old_t1) } else { (old_t1, old_t2) }; - let (old_pt1, old_pt2) = (orig_curve1.eval(old_t1_flipped), orig_curve2.eval(old_t2_flipped)); + let (old_t1_flipped, old_t2_flipped) = if flip { + (old_t2, old_t1) + } else { + (old_t1, old_t2) + }; + let (old_pt1, old_pt2) = ( + orig_curve1.eval(old_t1_flipped), + orig_curve2.eval(old_t2_flipped), + ); // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" // duplicate intersection, so always keep the one that minimizes sample distance. @@ -678,7 +732,9 @@ fn restrict_curve_to_perpendicular_fat_line< let baseline2 = curve2.baseline(); let center_baseline2 = baseline2.bounding_box().center().to_vec2(); - let affine_baseline2 = Affine::translate(center_baseline2) * Affine::rotate(PI / 2.) * Affine::translate(-center_baseline2); + let affine_baseline2 = Affine::translate(center_baseline2) + * Affine::rotate(PI / 2.) + * Affine::translate(-center_baseline2); let baseline2 = affine_baseline2 * baseline2; let (mut top, mut bottom) = curve1.convex_hull_from_line(&baseline2); let (d_min, d_max) = curve2.fat_line_min_max(&baseline2); @@ -765,7 +821,7 @@ fn rectangles_overlap(r1: &Rect, r2: &Rect) -> bool { #[cfg(test)] mod tests { use super::*; - use crate::{CubicBez, QuadBez, Line}; + use crate::{CubicBez, Line, QuadBez}; fn do_single_test( curve1: &T, @@ -791,11 +847,7 @@ mod tests { assert_eq!(arr2.len(), count); } - fn test_t_along_curve( - curve: &T, - t_test: f64, - accuracy: f64, - ) { + fn test_t_along_curve(curve: &T, t_test: f64, accuracy: f64) { let pt_test = curve.eval(t_test); let t_values = t_along_curve_for_point(pt_test, curve, accuracy, true); assert!(!t_values.is_empty()); @@ -804,8 +856,7 @@ mod tests { for t in &t_values { let pt = curve.eval(*t); - if (*t - t_test) <= 1e-9 - { + if (*t - t_test) <= 1e-9 { found_t = true; } @@ -814,21 +865,14 @@ mod tests { assert!(found_t) } - fn test_t( - curve: &T, - ) { + fn test_t(curve: &T) { let tenths = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; for &t in tenths.iter() { test_t_along_curve(curve, t, 0.); } } - fn test_overlapping( - curve: &T, - t1: f64, - t2: f64, - accuracy: f64, - ) { + fn test_overlapping(curve: &T, t1: f64, t2: f64, accuracy: f64) { assert!(t1 < t2); do_double_test( @@ -836,28 +880,28 @@ mod tests { &curve.subsegment(t1..1.0), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - accuracy + accuracy, ); do_double_test( &curve.subsegment(0.0..t2), &curve.subsegment(1.0..t1), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - accuracy + accuracy, ); do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t1..t2), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - accuracy + accuracy, ); do_double_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t2..t1), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - accuracy + accuracy, ); do_single_test( @@ -865,149 +909,159 @@ mod tests { &curve.subsegment(t1..t2), 1, CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection - accuracy + accuracy, ); do_single_test( &curve.subsegment(t1..t2), &curve.subsegment(0.0..1.0), 0, CurveIntersectionFlags::KEEP_CURVE2_T0_INTERSECTION, // Only include curve2 t=0 endpoint intersection - accuracy + accuracy, ); do_single_test( &curve.subsegment(0.0..1.0), &curve.subsegment(t1..t2), 1, CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection - accuracy + accuracy, ); do_single_test( &curve.subsegment(t1..t2), &curve.subsegment(0.0..1.0), 0, CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION, // Only include curve2 t=1 endpoint intersection - accuracy + accuracy, ); } #[test] fn test_cubic_cubic_intersections() { - test_t(&CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0))); - test_t(&CubicBez::new((0.0, 0.0), (-1.0, 0.0), (1.0, 0.0), (1.0, 0.0))); - + test_t(&CubicBez::new( + (0.0, 0.0), + (0.3, -1.0), + (0.7, -1.0), + (1.0, 0.0), + )); + test_t(&CubicBez::new( + (0.0, 0.0), + (-1.0, 0.0), + (1.0, 0.0), + (1.0, 0.0), + )); + do_double_test( &CubicBez::new((0.0, 0.0), (0.0, 1.0), (0.0, 1.0), (1.0, 1.0)), &CubicBez::new((0.0, 1.0), (1.0, 1.0), (1.0, 1.0), (1.0, 0.0)), 1, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((48.0, 84.0), (104.0, 176.0), (190.0, 37.0), (121.0, 75.0)), &CubicBez::new((68.0, 145.0), (74.0, 6.0), (143.0, 197.0), (138.0, 55.0)), 4, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.5, 1.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.0, 1.0), (0.5, 0.0), (0.5, 0.0), (1.0, 1.0)), 2, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &CubicBez::new((0.0, 0.0), (2.5, 0.5), (-1.5, 0.5), (1.0, 0.0)), 9, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), 3, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, 0.75), (1.0, 0.25)), 1, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &QuadBez::new((0.0, 0.25), (0.5, -0.25), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 2, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &QuadBez::new((0.0, 0.0), (0.5, 1.0), (1.0, 0.0)), &Line::new((0.0, 0.5), (1.0, 0.5)), 1, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.2, 0.0), (0.5, 3.0), (0.5, -2.0), (0.8, 1.0)), &Line::new((0.0, 0.5), (1.0, 0.25)), 3, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); - + // These tests intersect at the endpoints which is only sometimes useful to know about. do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((1.0, 0.0), (0.7, -1.0), (0.3, -1.0), (0.0, 0.0)), 1, CurveIntersectionFlags::KEEP_CURVE1_T0_INTERSECTION, // Include curve1 t=0 intersection - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 2, CurveIntersectionFlags::KEEP_ALL_ENDPOINT_INTERSECTIONS, // Include endpoint intersections - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 0, CurveIntersectionFlags::NONE, // Discard endpoint intersections - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); do_double_test( &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), &CubicBez::new((0.0, 0.0), (0.3, -0.5), (0.7, -0.5), (1.0, 0.0)), 1, CurveIntersectionFlags::KEEP_CURVE1_T1_INTERSECTION, // Include curve1 t=1 intersection - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); // Test curves that lie exactly along one-another @@ -1015,7 +1069,7 @@ mod tests { &CubicBez::new((0.0, 0.0), (0.3, -1.0), (0.7, -1.0), (1.0, 0.0)), 0.2, 0.8, - 0. + 0., ); // (A previous version of the code was returning two practically identical @@ -1035,7 +1089,7 @@ mod tests { ), 1, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); // On these curves the algorithm runs to a state at which the new clipped domain1 becomes a @@ -1057,7 +1111,7 @@ mod tests { ), 3, CurveIntersectionFlags::NONE, // Standard algorithm - 0., // Tightest possible accuracy + 0., // Tightest possible accuracy ); } } diff --git a/src/param_curve.rs b/src/param_curve.rs index cc13e3d2..7efb475c 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -216,7 +216,9 @@ pub trait ParamCurveExtrema: ParamCurve { } /// A parameterized curve that can be used in the Bezier clipping algorithm -pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen { +pub trait ParamCurveBezierClipping: + ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen +{ /// Returns a line from the curve's start point to its end point fn baseline(&self) -> Line { Line::new(self.start(), self.end()) diff --git a/src/param_curve_clipping.rs b/src/param_curve_clipping.rs index 819ccca9..ca0eb0f5 100644 --- a/src/param_curve_clipping.rs +++ b/src/param_curve_clipping.rs @@ -1,7 +1,7 @@ //! A trait for clipping parametrized curves. -use crate::common::{solve_cubic, solve_quadratic, solve_linear}; -use crate::{CubicBez, QuadBez, Line, Point, ParamCurve, ParamCurveBezierClipping}; +use crate::common::{solve_cubic, solve_linear, solve_quadratic}; +use crate::{CubicBez, Line, ParamCurve, ParamCurveBezierClipping, Point, QuadBez}; use arrayvec::ArrayVec; // Note that the line is unbounded here! @@ -16,8 +16,7 @@ fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { let b = b / len; let c = c / len; - if a.is_infinite() || b.is_infinite() || c.is_infinite() - { + if a.is_infinite() || b.is_infinite() || c.is_infinite() { // Can't compute distance from zero-length line, so return distance // from p0 instead return (p - l.p0).hypot(); @@ -66,7 +65,7 @@ impl ParamCurveBezierClipping for Line { fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { let d0 = signed_distance_from_ray_to_point(l, self.start()); let d1 = signed_distance_from_ray_to_point(l, self.end()); - + let p0 = Point::new(0.0, d0); let p1 = Point::new(1.0, d1); diff --git a/src/point.rs b/src/point.rs index e581ce29..8f84912a 100644 --- a/src/point.rs +++ b/src/point.rs @@ -184,21 +184,20 @@ impl Point { } /// Compare if two points are approximately equal - pub fn is_near(pt1: Point, pt2: Point, accuracy: f64) -> bool - { + pub fn is_near(pt1: Point, pt2: Point, accuracy: f64) -> bool { Point::is_approx(pt1, pt2, accuracy, 1.) } /// Compare if two points are approximately equal to some scaled accuracy - pub fn is_approx(pt1: Point, pt2: Point, accuracy: f64, scale: f64) -> bool - { + pub fn is_approx(pt1: Point, pt2: Point, accuracy: f64, scale: f64) -> bool { if pt1.x.is_infinite() || pt1.y.is_infinite() || pt2.x.is_infinite() || pt2.y.is_infinite() { return pt1 == pt2; } let accuracy_scaled = accuracy * scale; - let epsilon = Point::epsilon(pt1, accuracy_scaled).max(Point::epsilon(pt2, accuracy_scaled)); + let epsilon = + Point::epsilon(pt1, accuracy_scaled).max(Point::epsilon(pt2, accuracy_scaled)); (pt2 - pt1).hypot() <= epsilon } } @@ -345,17 +344,49 @@ mod tests { #[test] fn test_point_comparisons() { assert!(Point::is_near(Point::new(0., 0.), Point::new(0., 0.), 0.)); - assert!(Point::is_near(Point::new(1000000., 1000000.), Point::new(1000000., 1000000.), 0.)); - assert!(Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), 0.)); + assert!(Point::is_near( + Point::new(1000000., 1000000.), + Point::new(1000000., 1000000.), + 0., + )); + assert!(Point::is_near( + Point::new(f64::INFINITY, f64::INFINITY), + Point::new(f64::INFINITY, f64::INFINITY), + 0., + )); assert!(!Point::is_near(Point::new(0., 0.), Point::new(1., 1.), 0.)); - assert!(!Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), 0.)); + assert!(!Point::is_near( + Point::new(f64::INFINITY, f64::INFINITY), + Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), + 0., + )); let epsilon = 0.1; - assert!(Point::is_near(Point::new(0., 0.), Point::new(0. + epsilon, 0.), epsilon)); - assert!(Point::is_near(Point::new(1000000., 1000000.), Point::new(1000000., 1000000. + epsilon), epsilon)); - assert!(Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::INFINITY, f64::INFINITY), epsilon)); - assert!(!Point::is_near(Point::new(0., 0.), Point::new(1., 1. + epsilon), epsilon)); - assert!(!Point::is_near(Point::new(f64::INFINITY, f64::INFINITY), Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), epsilon)); + assert!(Point::is_near( + Point::new(0., 0.), + Point::new(0. + epsilon, 0.), + epsilon + )); + assert!(Point::is_near( + Point::new(1000000., 1000000.), + Point::new(1000000., 1000000. + epsilon), + epsilon + )); + assert!(Point::is_near( + Point::new(f64::INFINITY, f64::INFINITY), + Point::new(f64::INFINITY, f64::INFINITY), + epsilon + )); + assert!(!Point::is_near( + Point::new(0., 0.), + Point::new(1., 1. + epsilon), + epsilon + )); + assert!(!Point::is_near( + Point::new(f64::INFINITY, f64::INFINITY), + Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY), + epsilon + )); } } diff --git a/src/quadbez.rs b/src/quadbez.rs index a0ebe22f..dc54c096 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -7,9 +7,9 @@ use arrayvec::ArrayVec; use crate::common::solve_cubic; use crate::MAX_EXTREMA; use crate::{ - Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, - ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, - Rect, Shape, Vec2 + Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, + ParamCurveBezierClipping, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, + ParamCurveNearest, PathEl, Point, Rect, Shape, Vec2, }; /// A single quadratic Bézier segment. @@ -282,7 +282,7 @@ impl ParamCurveArclen for QuadBez { if a2.is_infinite() { // The arclength is zero - return 0. + return 0.; } let v0 = 0.25 * a2 * a2 * b * (2.0 * sabc - c2) + sabc; From c23c4ef0bd42a4b405ce0a0803803ca920821db3 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 14:12:43 -0600 Subject: [PATCH 13/20] Moved fn baseline to ParamCurve --- src/cubicbez.rs | 5 ++--- src/param_curve.rs | 10 +++++----- src/quadbez.rs | 4 ++-- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/cubicbez.rs b/src/cubicbez.rs index d517a9f1..303c93ca 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -8,9 +8,8 @@ use arrayvec::ArrayVec; use crate::common::{solve_quadratic, GAUSS_LEGENDRE_COEFFS_9}; use crate::{ - Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveBezierClipping, - ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, - QuadBez, Rect, Shape, + Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature, + ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, QuadBez, Rect, Shape, }; const MAX_SPLINE_SPLIT: usize = 100; diff --git a/src/param_curve.rs b/src/param_curve.rs index 7efb475c..e5f014ef 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -41,6 +41,11 @@ pub trait ParamCurve: Sized { fn end(&self) -> Point { self.eval(1.0) } + + /// Returns a line from the curve's start point to its end point + fn baseline(&self) -> Line { + Line::new(self.start(), self.end()) + } } // TODO: I might not want to have separate traits for all these. @@ -219,11 +224,6 @@ pub trait ParamCurveExtrema: ParamCurve { pub trait ParamCurveBezierClipping: ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen { - /// Returns a line from the curve's start point to its end point - fn baseline(&self) -> Line { - Line::new(self.start(), self.end()) - } - /// Find the time `t` at which the curve has the given x value fn solve_t_for_x(&self, x: f64) -> ArrayVec; /// Find the time `t` at which the curve has the given x value diff --git a/src/quadbez.rs b/src/quadbez.rs index dc54c096..0cfa09c0 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -8,8 +8,8 @@ use crate::common::solve_cubic; use crate::MAX_EXTREMA; use crate::{ Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, - ParamCurveBezierClipping, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, - ParamCurveNearest, PathEl, Point, Rect, Shape, Vec2, + ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, + Rect, Shape, Vec2, }; /// A single quadratic Bézier segment. From 5bfbf1e167950445c3e6856923933958f4ee8831 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 14:21:13 -0600 Subject: [PATCH 14/20] Organized code to reduce ParamCurveBezierClipping spaghetti --- src/cubicbez.rs | 27 +--- src/curve_intersections.rs | 277 +++++++++++++++++++++++++++++++++++- src/lib.rs | 2 - src/param_curve.rs | 14 -- src/param_curve_clipping.rs | 236 ------------------------------ 5 files changed, 276 insertions(+), 280 deletions(-) delete mode 100644 src/param_curve_clipping.rs diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 303c93ca..b176db76 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -1,10 +1,10 @@ //! Cubic Bézier segments. +use arrayvec::ArrayVec; use std::ops::{Mul, Range}; use crate::MAX_EXTREMA; use crate::{Line, QuadSpline, Vec2}; -use arrayvec::ArrayVec; use crate::common::{solve_quadratic, GAUSS_LEGENDRE_COEFFS_9}; use crate::{ @@ -606,7 +606,6 @@ mod tests { cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez, }; - use arrayvec::ArrayVec; #[test] fn cubicbez_deriv() { @@ -886,28 +885,4 @@ mod tests { converted[0].points()[2].distance(Point::new(88639.0 / 90.0, 52584.0 / 90.0)) < 0.0001 ); } - - use crate::param_curve::ParamCurveBezierClipping; - #[test] - fn solve_t_for_xy() { - fn verify(mut roots: ArrayVec, expected: &[f64]) { - assert_eq!(expected.len(), roots.len()); - let epsilon = 1e-6; - roots.sort_by(|a, b| a.partial_cmp(b).unwrap()); - - for i in 0..expected.len() { - assert!((roots[i] - expected[i]).abs() < epsilon); - } - } - - let curve = CubicBez::new((0.0, 0.0), (0.0, 8.0), (10.0, 8.0), (10.0, 0.0)); - verify(curve.solve_t_for_x(5.0), &[0.5]); - verify(curve.solve_t_for_y(6.0), &[0.5]); - - { - let curve = CubicBez::new((0.0, 10.0), (0.0, 10.0), (10.0, 10.0), (10.0, 10.0)); - - verify(curve.solve_t_for_y(10.0), &[]); - } - } } diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index a933b1ba..56de6a79 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,8 +1,258 @@ -use crate::ParamCurveBezierClipping; -use crate::{Affine, ParamCurve, ParamCurveExtrema, Point, Rect}; +//! Functions used for clipping parametrized curves and finding intersections between them + +use crate::common::{solve_cubic, solve_linear, solve_quadratic}; +use crate::{ + Affine, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveDeriv, ParamCurveExtrema, + Point, QuadBez, Rect, +}; use arrayvec::ArrayVec; use std::ops::Range; +/// A trait for descibing a parameterized curve that can be used in the Bezier clipping algorithm +pub trait ParamCurveBezierClipping: + ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen +{ + /// Find the time `t` at which the curve has the given x value + fn solve_t_for_x(&self, x: f64) -> ArrayVec; + /// Find the time `t` at which the curve has the given x value + fn solve_t_for_y(&self, y: f64) -> ArrayVec; + /// Returns the upper and lower convex hull + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); + /// Returns the minimum and maximum distances of the "fat line" enclosing this curve + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64); +} + +// Note that the line is unbounded here! +fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { + let vec2 = l.p1 - l.p0; + let a = -vec2.y; + let b = vec2.x; + let c = -(a * l.start().x + b * l.start().y); + + let len = (a * a + b * b).sqrt(); + let a = a / len; + let b = b / len; + let c = c / len; + + if a.is_infinite() || b.is_infinite() || c.is_infinite() { + // Can't compute distance from zero-length line, so return distance + // from p0 instead + return (p - l.p0).hypot(); + } + + a * p.x + b * p.y + c +} + +impl ParamCurveBezierClipping for Line { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { + if (self.p0.x - self.p1.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b) = self.parameters(); + solve_linear(b.x - x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec { + if (self.p0.y - self.p1.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b) = self.parameters(); + solve_linear(b.y - y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + + // Calculate min and max distances + let d_min = f64::min(d0.min(d1), 0.0); + let d_max = f64::max(d0.max(d1), 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0, d1); + + // The hull is simply the line itself + (vec![p0, p1], vec![p0, p1]) + } +} + +impl ParamCurveBezierClipping for QuadBez { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { + if self.is_linear(f64::EPSILON) && (self.p0.x - self.p2.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b, c) = self.parameters(); + solve_quadratic(c.x - x, b.x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec { + if self.is_linear(f64::EPSILON) && (self.p0.y - self.p2.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b, c) = self.parameters(); + solve_quadratic(c.y - y, b.y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); + + // Shrink d1 to a more accurate number, but test the distance compared + // to d0 and d3 in case we're looking at a perpendicular fat line + let factor = 1.0 / 2.0; + //let factor = 1.0; + let d1_1 = (d1 - d0) * factor + d0; + let d1_2 = (d1 - d2) * factor + d2; + // let d1_1 = d1; + // let d1_2 = d1; + + // Calculate min and max distances + let d_min = f64::min(d0.min(d1_1).min(d1_2).min(d2), 0.0); + let d_max = f64::max(d0.max(d1_1).max(d1_2).max(d2), 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.p1); + let d2 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0 / 2.0, d1); + let p2 = Point::new(1.0, d2); + // Compute the vertical signed distance of p1 and p2 from [p0, p3]. + let dist1 = d1 - (d0 + d2) / 2.0; + + // Compute the hull assuming p1 is on top - we'll switch later if needed. + let mut hull = (vec![p0, p1, p2], vec![p0, p2]); + + // Flip the hull if needed: + if dist1 < 0.0 { + hull = (hull.1, hull.0); + } + + hull + } +} + +impl ParamCurveBezierClipping for CubicBez { + fn solve_t_for_x(&self, x: f64) -> ArrayVec { + if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { + return ArrayVec::new(); + } + let (a, b, c, d) = self.parameters(); + solve_cubic(d.x - x, c.x, b.x, a.x) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + fn solve_t_for_y(&self, y: f64) -> ArrayVec { + if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { + return ArrayVec::new(); + } + + let (a, b, c, d) = self.parameters(); + solve_cubic(d.y - y, c.y, b.y, a.y) + .iter() + .copied() + .filter(|&t| (0.0..=1.0).contains(&t)) + .collect() + } + + fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { + // Get the signed distance to each point + let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); + let d3 = signed_distance_from_ray_to_point(&baseline, self.p3); + + // Shrink d1 and d2 to more accurate numbers + let factor: f64 = if ((d1 - d0) * (d2 - d3)) > 0.0 { + 3.0 / 4.0 + } else { + 4.0 / 9.0 + }; + let d1 = (d1 - d0) * factor + d0; + let d2 = (d2 - d3) * factor + d3; + + // Calculate min and max distances + let d_min = f64::min(d0.min(d1).min(d2).min(d3), 0.0); + let d_max = f64::max(d0.max(d1).max(d2).max(d3), 0.0); + + (d_min, d_max) + } + + fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { + let d0 = signed_distance_from_ray_to_point(l, self.start()); + let d1 = signed_distance_from_ray_to_point(l, self.p1); + let d2 = signed_distance_from_ray_to_point(l, self.p2); + let d3 = signed_distance_from_ray_to_point(l, self.end()); + + let p0 = Point::new(0.0, d0); + let p1 = Point::new(1.0 / 3.0, d1); + let p2 = Point::new(2.0 / 3.0, d2); + let p3 = Point::new(1.0, d3); + // Compute the vertical signed distance of p1 and p2 from [p0, p3]. + let dist1 = d1 - (2.0 * d0 + d3) / 3.0; + let dist2 = d2 - (d0 + 2.0 * d3) / 3.0; + + // Compute the hull assuming p1 is on top - we'll switch later if needed. + let mut hull = if dist1 * dist2 < 0.0 { + // p1 and p2 lie on opposite sides of [p0, p3], so the hull is a quadrilateral: + (vec![p0, p1, p3], vec![p0, p2, p3]) + } else { + // p1 and p2 lie on the same side of [p0, p3]. The hull can be a triangle or a + // quadrilateral, and [p0, p3] is part of the hull. The hull is a triangle if the vertical + // distance of one of the middle points p1, p2 is <= half the vertical distance of the + // other middle point. + let dist1 = dist1.abs(); + let dist2 = dist2.abs(); + if dist1 >= 2.0 * dist2 { + (vec![p0, p1, p3], vec![p0, p3]) + } else if dist2 >= 2.0 * dist1 { + (vec![p0, p2, p3], vec![p0, p3]) + } else { + (vec![p0, p1, p2, p3], vec![p0, p3]) + } + }; + + // Flip the hull if needed: + if dist1 < 0.0 || (dist1 == 0.0 && dist2 < 0.0) { + hull = (hull.1, hull.0); + } + + hull + } +} + bitflags::bitflags! { pub struct CurveIntersectionFlags: u32 { const NONE = 0b00000000; @@ -823,6 +1073,29 @@ mod tests { use super::*; use crate::{CubicBez, Line, QuadBez}; + #[test] + fn solve_t_for_xy() { + fn verify(mut roots: ArrayVec, expected: &[f64]) { + assert_eq!(expected.len(), roots.len()); + let epsilon = 1e-6; + roots.sort_by(|a, b| a.partial_cmp(b).unwrap()); + + for i in 0..expected.len() { + assert!((roots[i] - expected[i]).abs() < epsilon); + } + } + + let curve = CubicBez::new((0.0, 0.0), (0.0, 8.0), (10.0, 8.0), (10.0, 0.0)); + verify(curve.solve_t_for_x(5.0), &[0.5]); + verify(curve.solve_t_for_y(6.0), &[0.5]); + + { + let curve = CubicBez::new((0.0, 10.0), (0.0, 10.0), (10.0, 10.0), (10.0, 10.0)); + + verify(curve.solve_t_for_y(10.0), &[]); + } + } + fn do_single_test( curve1: &T, curve2: &U, diff --git a/src/lib.rs b/src/lib.rs index cb897da7..8432385f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -90,7 +90,6 @@ mod insets; mod line; mod mindist; mod param_curve; -mod param_curve_clipping; mod point; mod quadbez; mod quadspline; @@ -113,7 +112,6 @@ pub use crate::ellipse::*; pub use crate::insets::*; pub use crate::line::*; pub use crate::param_curve::*; -pub use crate::param_curve_clipping::*; pub use crate::point::*; pub use crate::quadbez::*; pub use crate::quadspline::*; diff --git a/src/param_curve.rs b/src/param_curve.rs index e5f014ef..2c4a7275 100644 --- a/src/param_curve.rs +++ b/src/param_curve.rs @@ -219,17 +219,3 @@ pub trait ParamCurveExtrema: ParamCurve { bbox } } - -/// A parameterized curve that can be used in the Bezier clipping algorithm -pub trait ParamCurveBezierClipping: - ParamCurve + ParamCurveDeriv + ParamCurveExtrema + ParamCurveArclen -{ - /// Find the time `t` at which the curve has the given x value - fn solve_t_for_x(&self, x: f64) -> ArrayVec; - /// Find the time `t` at which the curve has the given x value - fn solve_t_for_y(&self, y: f64) -> ArrayVec; - /// Returns the upper and lower convex hull - fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); - /// Returns the minimum and maximum distances of the "fat line" enclosing this curve - fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64); -} diff --git a/src/param_curve_clipping.rs b/src/param_curve_clipping.rs deleted file mode 100644 index ca0eb0f5..00000000 --- a/src/param_curve_clipping.rs +++ /dev/null @@ -1,236 +0,0 @@ -//! A trait for clipping parametrized curves. - -use crate::common::{solve_cubic, solve_linear, solve_quadratic}; -use crate::{CubicBez, Line, ParamCurve, ParamCurveBezierClipping, Point, QuadBez}; -use arrayvec::ArrayVec; - -// Note that the line is unbounded here! -fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { - let vec2 = l.p1 - l.p0; - let a = -vec2.y; - let b = vec2.x; - let c = -(a * l.start().x + b * l.start().y); - - let len = (a * a + b * b).sqrt(); - let a = a / len; - let b = b / len; - let c = c / len; - - if a.is_infinite() || b.is_infinite() || c.is_infinite() { - // Can't compute distance from zero-length line, so return distance - // from p0 instead - return (p - l.p0).hypot(); - } - - a * p.x + b * p.y + c -} - -impl ParamCurveBezierClipping for Line { - fn solve_t_for_x(&self, x: f64) -> ArrayVec { - if (self.p0.x - self.p1.x).abs() < f64::EPSILON { - return ArrayVec::new(); - } - let (a, b) = self.parameters(); - solve_linear(b.x - x, a.x) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - fn solve_t_for_y(&self, y: f64) -> ArrayVec { - if (self.p0.y - self.p1.y).abs() < f64::EPSILON { - return ArrayVec::new(); - } - - let (a, b) = self.parameters(); - solve_linear(b.y - y, a.y) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - - fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { - // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); - - // Calculate min and max distances - let d_min = f64::min(d0.min(d1), 0.0); - let d_max = f64::max(d0.max(d1), 0.0); - - (d_min, d_max) - } - - fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { - let d0 = signed_distance_from_ray_to_point(l, self.start()); - let d1 = signed_distance_from_ray_to_point(l, self.end()); - - let p0 = Point::new(0.0, d0); - let p1 = Point::new(1.0, d1); - - // The hull is simply the line itself - (vec![p0, p1], vec![p0, p1]) - } -} - -impl ParamCurveBezierClipping for QuadBez { - fn solve_t_for_x(&self, x: f64) -> ArrayVec { - if self.is_linear(f64::EPSILON) && (self.p0.x - self.p2.x).abs() < f64::EPSILON { - return ArrayVec::new(); - } - let (a, b, c) = self.parameters(); - solve_quadratic(c.x - x, b.x, a.x) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - fn solve_t_for_y(&self, y: f64) -> ArrayVec { - if self.is_linear(f64::EPSILON) && (self.p0.y - self.p2.y).abs() < f64::EPSILON { - return ArrayVec::new(); - } - - let (a, b, c) = self.parameters(); - solve_quadratic(c.y - y, b.y, a.y) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - - fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { - // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); - let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); - - // Shrink d1 to a more accurate number, but test the distance compared - // to d0 and d3 in case we're looking at a perpendicular fat line - let factor = 1.0 / 2.0; - //let factor = 1.0; - let d1_1 = (d1 - d0) * factor + d0; - let d1_2 = (d1 - d2) * factor + d2; - // let d1_1 = d1; - // let d1_2 = d1; - - // Calculate min and max distances - let d_min = f64::min(d0.min(d1_1).min(d1_2).min(d2), 0.0); - let d_max = f64::max(d0.max(d1_1).max(d1_2).max(d2), 0.0); - - (d_min, d_max) - } - - fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { - let d0 = signed_distance_from_ray_to_point(l, self.start()); - let d1 = signed_distance_from_ray_to_point(l, self.p1); - let d2 = signed_distance_from_ray_to_point(l, self.end()); - - let p0 = Point::new(0.0, d0); - let p1 = Point::new(1.0 / 2.0, d1); - let p2 = Point::new(1.0, d2); - // Compute the vertical signed distance of p1 and p2 from [p0, p3]. - let dist1 = d1 - (d0 + d2) / 2.0; - - // Compute the hull assuming p1 is on top - we'll switch later if needed. - let mut hull = (vec![p0, p1, p2], vec![p0, p2]); - - // Flip the hull if needed: - if dist1 < 0.0 { - hull = (hull.1, hull.0); - } - - hull - } -} - -impl ParamCurveBezierClipping for CubicBez { - fn solve_t_for_x(&self, x: f64) -> ArrayVec { - if self.is_linear(f64::EPSILON) && (self.p0.x - self.p3.x).abs() < f64::EPSILON { - return ArrayVec::new(); - } - let (a, b, c, d) = self.parameters(); - solve_cubic(d.x - x, c.x, b.x, a.x) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - fn solve_t_for_y(&self, y: f64) -> ArrayVec { - if self.is_linear(f64::EPSILON) && (self.p0.y - self.p3.y).abs() < f64::EPSILON { - return ArrayVec::new(); - } - - let (a, b, c, d) = self.parameters(); - solve_cubic(d.y - y, c.y, b.y, a.y) - .iter() - .copied() - .filter(|&t| (0.0..=1.0).contains(&t)) - .collect() - } - - fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { - // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); - let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); - let d3 = signed_distance_from_ray_to_point(&baseline, self.p3); - - // Shrink d1 and d2 to more accurate numbers - let factor: f64 = if ((d1 - d0) * (d2 - d3)) > 0.0 { - 3.0 / 4.0 - } else { - 4.0 / 9.0 - }; - let d1 = (d1 - d0) * factor + d0; - let d2 = (d2 - d3) * factor + d3; - - // Calculate min and max distances - let d_min = f64::min(d0.min(d1).min(d2).min(d3), 0.0); - let d_max = f64::max(d0.max(d1).max(d2).max(d3), 0.0); - - (d_min, d_max) - } - - fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec) { - let d0 = signed_distance_from_ray_to_point(l, self.start()); - let d1 = signed_distance_from_ray_to_point(l, self.p1); - let d2 = signed_distance_from_ray_to_point(l, self.p2); - let d3 = signed_distance_from_ray_to_point(l, self.end()); - - let p0 = Point::new(0.0, d0); - let p1 = Point::new(1.0 / 3.0, d1); - let p2 = Point::new(2.0 / 3.0, d2); - let p3 = Point::new(1.0, d3); - // Compute the vertical signed distance of p1 and p2 from [p0, p3]. - let dist1 = d1 - (2.0 * d0 + d3) / 3.0; - let dist2 = d2 - (d0 + 2.0 * d3) / 3.0; - - // Compute the hull assuming p1 is on top - we'll switch later if needed. - let mut hull = if dist1 * dist2 < 0.0 { - // p1 and p2 lie on opposite sides of [p0, p3], so the hull is a quadrilateral: - (vec![p0, p1, p3], vec![p0, p2, p3]) - } else { - // p1 and p2 lie on the same side of [p0, p3]. The hull can be a triangle or a - // quadrilateral, and [p0, p3] is part of the hull. The hull is a triangle if the vertical - // distance of one of the middle points p1, p2 is <= half the vertical distance of the - // other middle point. - let dist1 = dist1.abs(); - let dist2 = dist2.abs(); - if dist1 >= 2.0 * dist2 { - (vec![p0, p1, p3], vec![p0, p3]) - } else if dist2 >= 2.0 * dist1 { - (vec![p0, p2, p3], vec![p0, p3]) - } else { - (vec![p0, p1, p2, p3], vec![p0, p3]) - } - }; - - // Flip the hull if needed: - if dist1 < 0.0 || (dist1 == 0.0 && dist2 < 0.0) { - hull = (hull.1, hull.0); - } - - hull - } -} From 9f1c40f9880ac5e2c3cdf3699574a10a0bc27ac1 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 15:47:43 -0600 Subject: [PATCH 15/20] Added disclaimer for curve intersection functions --- src/curve_intersections.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 56de6a79..2effe112 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -1,4 +1,6 @@ //! Functions used for clipping parametrized curves and finding intersections between them +//! +//! Caution: All functions in this file are provisional and may change in the future! use crate::common::{solve_cubic, solve_linear, solve_quadratic}; use crate::{ From 531d5387a3de0fa615b140ed49ec2a1a4ec5f6ba Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Thu, 7 Apr 2022 15:57:05 -0600 Subject: [PATCH 16/20] Fixed running tests in release mode --- src/bezpath.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/bezpath.rs b/src/bezpath.rs index 3153be66..fe76a230 100644 --- a/src/bezpath.rs +++ b/src/bezpath.rs @@ -1280,6 +1280,7 @@ mod tests { #[test] #[should_panic(expected = "no open subpath")] + #[cfg(debug_assertions)] // Only provides proper panic message in debug mode fn test_elements_to_segments_starts_on_closepath() { let mut path = BezPath::new(); path.close_path(); @@ -1302,6 +1303,7 @@ mod tests { #[test] #[should_panic(expected = "no open subpath")] + #[cfg(debug_assertions)] // Only provides proper panic message in debug mode fn test_must_not_start_on_quad() { let mut path = BezPath::new(); path.quad_to((5.0, 5.0), (10.0, 10.0)); From 5a6ac30c574b0d7b8e60ae2d583faf2e34452eb7 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Tue, 19 Apr 2022 11:20:46 -0600 Subject: [PATCH 17/20] Fixed clippy issues --- src/curve_intersections.rs | 71 +++++++++++++++----------------------- 1 file changed, 28 insertions(+), 43 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 2effe112..d9aebc41 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -72,8 +72,8 @@ impl ParamCurveBezierClipping for Line { fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); + let d0 = signed_distance_from_ray_to_point(baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(baseline, self.p1); // Calculate min and max distances let d_min = f64::min(d0.min(d1), 0.0); @@ -121,9 +121,9 @@ impl ParamCurveBezierClipping for QuadBez { fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); - let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); + let d0 = signed_distance_from_ray_to_point(baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(baseline, self.p2); // Shrink d1 to a more accurate number, but test the distance compared // to d0 and d3 in case we're looking at a perpendicular fat line @@ -191,10 +191,10 @@ impl ParamCurveBezierClipping for CubicBez { fn fat_line_min_max(&self, baseline: &Line) -> (f64, f64) { // Get the signed distance to each point - let d0 = signed_distance_from_ray_to_point(&baseline, self.p0); - let d1 = signed_distance_from_ray_to_point(&baseline, self.p1); - let d2 = signed_distance_from_ray_to_point(&baseline, self.p2); - let d3 = signed_distance_from_ray_to_point(&baseline, self.p3); + let d0 = signed_distance_from_ray_to_point(baseline, self.p0); + let d1 = signed_distance_from_ray_to_point(baseline, self.p1); + let d2 = signed_distance_from_ray_to_point(baseline, self.p2); + let d3 = signed_distance_from_ray_to_point(baseline, self.p3); // Shrink d1 and d2 to more accurate numbers let factor: f64 = if ((d1 - d0) * (d2 - d3)) > 0.0 { @@ -340,7 +340,7 @@ fn check_for_overlap( return (false, 0.0); } - assert!(t_values_other.1.len() >= 1); // This has to be true or we really messed up! + assert!(!t_values_other.1.is_empty()); // This has to be true or we really messed up! if t_values_other.1.len() == 1 { return (true, t_values_other.1[0]); } @@ -521,15 +521,11 @@ fn add_curve_intersections( - pt_curve: &T, - pt_curve_is_curve1: bool, - curve: &U, - pt_domain: &Range, - curve_domain: &Range, + curves: (&T, &U), + domains: (&Range, &Range), intersections: &mut ArrayVec<(f64, f64), 9>, flip: bool, - orig_curve1: &T, - orig_curve2: &U, + orig_curves: (&T, &U), flags: CurveIntersectionFlags, accuracy: f64, ) { - // We assume pt is curve1 when we add intersections below. - let flip = if pt_curve_is_curve1 { flip } else { !flip }; - // Get the mid-point - let pt = pt_curve.eval(0.5); - let pt_t = (pt_domain.start + pt_domain.end) * 0.5; + let pt = curves.0.eval(0.5); + let pt_t = (domains.0.start + domains.0.end) * 0.5; // Generally speaking |curve| will be quite small at this point, so see if we can get away with // just finding one of the points along the curve (even though there may be multiple). - let results = t_along_curve_for_point(pt, curve, accuracy, false); + let results = t_along_curve_for_point(pt, curves.1, accuracy, false); let curve_t = results .iter() .fold((-1., f64::MAX), |(t, d_sq), result| { - let result_pos = curve.eval(*result); + let result_pos = curves.1.eval(*result); let result_d_sq = (pt - result_pos).hypot2(); if result_d_sq < d_sq { return (*result, result_d_sq); @@ -781,12 +766,12 @@ fn add_point_curve_intersection( return (true, t_values.clone()); // Point lies along curve } } - (false, t_values.clone()) // Point is not along curve + (false, t_values) // Point is not along curve } fn add_intersection( From 7ddb666329ffcc018acf7069204b3c4312719a2c Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Wed, 20 Apr 2022 12:34:15 -0600 Subject: [PATCH 18/20] Replaced vec object with splice --- src/curve_intersections.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index d9aebc41..d5e6a4d6 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -981,8 +981,8 @@ fn restrict_curve_to_perpendicular_fat_line< // Returns the min and max values at which the convex hull enters the fat line min/max offset lines. fn clip_convex_hull_to_fat_line( - hull_top: &mut Vec, - hull_bottom: &mut Vec, + hull_top: &mut [Point], + hull_bottom: &mut [Point], d_min: f64, d_max: f64, ) -> Option<(f64, f64)> { From 099da3e610cc5a50b5123a93129591f82aa3a716 Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Tue, 23 Aug 2022 10:14:41 -0600 Subject: [PATCH 19/20] Addressed pull request comments --- src/common.rs | 8 +- src/cubicbez.rs | 2 + src/curve_intersections.rs | 11 +- src/fat_line.rs | 111 +++++ src/image.rs | 805 +++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 +- src/pathops.rs | 429 ++++++++++++++++++++ src/point.rs | 6 + src/quadbez.rs | 11 +- 9 files changed, 1376 insertions(+), 9 deletions(-) create mode 100644 src/fat_line.rs create mode 100644 src/image.rs create mode 100644 src/pathops.rs diff --git a/src/common.rs b/src/common.rs index dc53cedd..5f786729 100644 --- a/src/common.rs +++ b/src/common.rs @@ -51,6 +51,8 @@ impl FloatExt for f32 { /// See: /// /// Return values of x for which c0 + c1 x + c2 x² + c3 x³ = 0. +/// +/// Note: The input arguments are in increasing exponent order pub fn solve_cubic(c0: f64, c1: f64, c2: f64, c3: f64) -> ArrayVec { let mut result = ArrayVec::new(); let c3_recip = c3.recip(); @@ -103,6 +105,8 @@ pub fn solve_cubic(c0: f64, c1: f64, c2: f64, c3: f64) -> ArrayVec { /// the other root might be out of representable range. In the degenerate /// case where all coefficients are zero, so that all values of x satisfy /// the equation, a single `0.0` is returned. +/// +/// Note: The input arguments are in increasing exponent order pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { let mut result = ArrayVec::new(); let sc0 = c0 * c2.recip(); @@ -148,13 +152,15 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { /// Find real roots of linear equation. /// /// Return values of x for which c0 + c1 x = 0. +/// +/// Note: The input arguments are in increasing exponent order pub fn solve_linear(c0: f64, c1: f64) -> ArrayVec { let mut result = ArrayVec::new(); let root = -c0 / c1; if root.is_finite() { result.push(root); } else if c0 == 0.0 && c1 == 0.0 { - // Degenerate case + // Degenerate case, this happens the line is horizontal, hence when c1 == 0. result.push(0.0); } result diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 46869efe..f1b9d925 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -239,6 +239,8 @@ impl CubicBez { /// Get the parameters such that the curve can be represented by the following formula: /// B(t) = a*t^3 + b*t^2 + c*t + d + /// + /// Note: Values returned are in decresing exponent order pub fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) { let c = (self.p1 - self.p0) * 3.0; let b = (self.p2 - self.p1) * 3.0 - c; diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index d5e6a4d6..19e89ed5 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -16,7 +16,7 @@ pub trait ParamCurveBezierClipping: { /// Find the time `t` at which the curve has the given x value fn solve_t_for_x(&self, x: f64) -> ArrayVec; - /// Find the time `t` at which the curve has the given x value + /// Find the time `t` at which the curve has the given y value fn solve_t_for_y(&self, y: f64) -> ArrayVec; /// Returns the upper and lower convex hull fn convex_hull_from_line(&self, l: &Line) -> (Vec, Vec); @@ -31,10 +31,10 @@ fn signed_distance_from_ray_to_point(l: &Line, p: Point) -> f64 { let b = vec2.x; let c = -(a * l.start().x + b * l.start().y); - let len = (a * a + b * b).sqrt(); - let a = a / len; - let b = b / len; - let c = c / len; + let len_inv = (a * a + b * b).sqrt().recip(); + let a = a * len_inv; + let b = b * len_inv; + let c = c * len_inv; if a.is_infinite() || b.is_infinite() || c.is_infinite() { // Can't compute distance from zero-length line, so return distance @@ -328,7 +328,6 @@ fn check_for_overlap( // curve 2: O-----------O // Use this function to get the t-value on another curve that corresponds to the t-value on this curve - #[inline] fn t_value_on_other( t_this: f64, curve_this: &T, diff --git a/src/fat_line.rs b/src/fat_line.rs new file mode 100644 index 00000000..284ffbdd --- /dev/null +++ b/src/fat_line.rs @@ -0,0 +1,111 @@ +// DO NOT COMMIT THIS FILE + +use std::{ops, fmt}; + +pub use crate::{Point, Vec2, Line}; + +/***************************************************************************** + * Struct + *****************************************************************************/ + +/// A struct that defines a fat line +#[derive(Debug, Copy, Clone)] +pub struct FatLine +{ + /// The line's start point. + pub p0: Point, + /// The line's end point. + pub p1: Point, + /// The min distance from the center-line (zero or less) + pub dmin: f64, + /// The max distance from the center-line (zero or more) + pub dmax: f64, +} + +/***************************************************************************** + * Implementation + *****************************************************************************/ + +impl FatLine +{ + /// Gives the distance from a point to the center line + fn distance_from_line(line: &Line, pt: &Point) -> f64 + { + // Calculate line parameters + let mut a = line.p1.y - line.p0.y; + let mut b = line.p0.x - line.p1.x; + let mut c = line.p0.x * -line.p1.y + line.p0.y * line.p1.x; + let length = (a * a + b * b).sqrt(); + + // Normalize parameters + a /= length; + b /= length; + c /= length; + + // Calculate distance + a * pt.x + b * pt.y + c + } + + /// Creates a fat line from a control polygon + pub fn from_control_poly(control_poly: &Vec) -> FatLine + { + // Get start/end points of fat line + let line = Line { + p0: control_poly.first().unwrap().clone(), + p1: control_poly.last().unwrap().clone(), + }; + + // Determine min/max distances + let mut mind: f64 = 0.; + let mut maxd: f64 = 0.; + println!("test poly: {}", control_poly.len()); + for pt in control_poly + { + let d = FatLine::distance_from_line(&line, pt); + mind = mind.min(d); + maxd = maxd.max(d); + + println!("test dist: {} -> {} {}", d, mind, maxd); + } + + FatLine { + p0: line.p0, + p1: line.p1, + dmin: mind, + dmax: maxd + } + } +} + +/***************************************************************************** + * Ops + *****************************************************************************/ + +impl ops::Add for FatLine +{ + type Output = Self; + + /// Adds a Vec2 to a fat line + fn add(self, other: Vec2) -> Self::Output + { + FatLine { + p0: self.p0 + other, + p1: self.p1 + other, + dmin: self.dmin, + dmax: self.dmax, + } + } +} + +/***************************************************************************** + * Format + *****************************************************************************/ + +impl fmt::Display for FatLine +{ + /// Defines how fat lines are printed to the output + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result + { + write!(f, "Fattie({}, {}, [{}, {}])", self.p0, self.p1, self.dmin, self.dmax) + } +} \ No newline at end of file diff --git a/src/image.rs b/src/image.rs new file mode 100644 index 00000000..324fa5a6 --- /dev/null +++ b/src/image.rs @@ -0,0 +1,805 @@ +// DO NOT COMMIT THIS FILE + +use std::{fmt}; +use std::any::Any; + +use raqote::*; +use crate::{ + Affine, CubicBez, QuadBez, Shape, ParamCurve, Point, Vec2, Size, Rect, + Line, FatLine +}; +use arrayvec::ArrayVec; + +/***************************************************************************** + * Struct + *****************************************************************************/ + +/// A struct that produces an image with the contained items drawn on it +pub struct Image +{ + /// The image scale + scale: f64, + /// The offset between the top-left corner of the item bounds and the image origin + offset: Vec2, + /// The white space around the edge of the image in pixels + margins: i32, + /// The fixed size of the image + fixed_size: Size, + /// Whether or not the image should use the fixed_size field (or the image scale) + use_fixed_size: bool, + /// The items to draw to the image + items: Vec>, +} + +/***************************************************************************** + * Enums + *****************************************************************************/ + + #[derive(Debug)] +enum CurveType { + CurveNone, + CurveCubic, + CurveQuad, + CurveLinear, +} + +fn get_curve_type_from_curve_param(c1: &dyn Any) -> CurveType +{ + if c1.is::() { return CurveType::CurveCubic } + if c1.is::() { return CurveType::CurveQuad } + if c1.is::() { return CurveType::CurveLinear } + CurveType::CurveNone +} + +/***************************************************************************** + * Trait - Drawable Element + *****************************************************************************/ + +trait Drawable +{ + fn bounds(&self) -> Rect; + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget); +} + +struct DrawableCubicBez(pub CubicBez, pub u32, pub bool); +impl Drawable for DrawableCubicBez { + #[inline(never)] + fn bounds(&self) -> Rect + { + let bbox = if self.2 { Rect::ZERO } else { self.0.bounding_box() }; + // let bbox = bbox.union_pt(self.0.p1); + // let bbox = bbox.union_pt(self.0.p2); + bbox + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + if self.2 + { + image.draw_cubic_bez_with_opacity(dt, &self.0, self.1, 0x40); + } + else + { + image.draw_cubic_bez(dt, &self.0, self.1); + image.draw_point(dt, &self.0.p0, 5., self.1, true); + image.draw_point_with_opacity(dt, &self.0.p1, 5., self.1, 0x80, true); + // image.draw_dashed_line_with_opacity(dt, &Line::new(self.0.p0, self.0.p1), self.1, 0x40); + image.draw_line_with_opacity(dt, &Line::new(self.0.p0, self.0.p1), self.1, 0x40); + image.draw_point_with_opacity(dt, &self.0.p2, 5., self.1, 0x80, true); + image.draw_point(dt, &self.0.p3, 5., self.1, true); + // image.draw_dashed_line_with_opacity(dt, &Line::new(self.0.p3, self.0.p2), self.1, 0x40); + image.draw_line_with_opacity(dt, &Line::new(self.0.p3, self.0.p2), self.1, 0x40); + } + } +} + +struct DrawableQuadBezier(pub QuadBez, pub u32, pub bool); +impl Drawable for DrawableQuadBezier { + #[inline(never)] + fn bounds(&self) -> Rect + { + let bbox = if self.2 { Rect::ZERO } else { self.0.bounding_box() }; + // let bbox = bbox.union_pt(self.0.p1); + bbox + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + if self.2 + { + image.draw_quad_bez_with_opacity(dt, &self.0, self.1, 0x40); + } + else + { + image.draw_quad_bez(dt, &self.0, self.1); + image.draw_point(dt, &self.0.p0, 5., self.1, true); + image.draw_point_with_opacity(dt, &self.0.p1, 5., self.1, 0x80, true); + image.draw_point(dt, &self.0.p2, 5., self.1, true); + // image.draw_dashed_line_with_opacity(dt, &Line::new(self.0.p0, self.0.p1), self.1, 0x40); + // image.draw_dashed_line_with_opacity(dt, &Line::new(self.0.p2, self.0.p1), self.1, 0x40); + image.draw_line_with_opacity(dt, &Line::new(self.0.p0, self.0.p1), self.1, 0x40); + image.draw_line_with_opacity(dt, &Line::new(self.0.p2, self.0.p1), self.1, 0x40); + } + } +} + +struct DrawableLine(pub Line, pub u32, pub bool); +impl Drawable for DrawableLine { + #[inline(never)] + fn bounds(&self) -> Rect + { + let bbox = if self.2 { Rect::ZERO } else { self.0.bounding_box() }; + bbox + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + if self.2 + { + image.draw_line_with_opacity(dt, &self.0, self.1, 0x40); + } + else + { + image.draw_line(dt, &self.0, self.1); + image.draw_point(dt, &self.0.p0, 5., self.1, true); + image.draw_point(dt, &self.0.p1, 5., self.1, true); + } + } +} + +struct DrawableFatLine(pub FatLine, pub u32); +impl Drawable for DrawableFatLine { + #[inline(never)] + fn bounds(&self) -> Rect + { + // Use empty rect so that it's not used in calculating total bounds + // Note: Fat lines are technically infinite + Rect::ZERO + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + image.draw_fat_line(dt, &self.0, self.1); + } +} + +struct DrawablePoint(pub Point, pub u32, pub f32, pub bool); +impl Drawable for DrawablePoint +{ + #[inline(never)] + fn bounds(&self) -> Rect + { + // Use empty rect so that it's not used in calculating total bounds + Rect::ZERO + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + image.draw_point(dt, &self.0, self.2, self.1, self.3); + } +} + +struct DrawableIntersections( + pub ArrayVec<(f64, f64), 9>, + pub T, + pub U, + pub u32, + pub f32, + pub bool + ); +impl Drawable for DrawableIntersections where + T: ParamCurve, + U: ParamCurve +{ + #[inline(never)] + fn bounds(&self) -> Rect + { + // Use empty rect so that it's not used in calculating total bounds + Rect::ZERO + } + + #[inline(never)] + fn draw_to_image(&self, image: &Image, dt: &mut DrawTarget) + { + for &i in &self.0 + { + image.draw_point(dt, &self.1.eval(i.0), self.4, self.3, self.5); + image.draw_point(dt, &self.2.eval(i.1), self.4 + 5., self.3, false); + } + } +} + +/***************************************************************************** + * Implementation + *****************************************************************************/ + +impl Image +{ + /// Create an image that scales all items by a constant amount + pub fn from_scale(scale_: f64, margins_: i32) -> Image + { + Image { + scale: scale_, + offset: Vec2::new(0., 0.), + margins: margins_, + fixed_size: Size::new(0., 0.), + use_fixed_size: false, + items: vec![], + } + } + + /// Create an image with a fixed size, all items are scaled to fit + pub fn from_fixed_size(size_: Size, margins_: i32) -> Image + { + Image { + scale: 1.0, + offset: Vec2::new(0., 0.), + margins: margins_, + fixed_size: size_, + use_fixed_size: true, + items: vec![], + } + } + + /// Adds a cubic bezier to the list of items + pub fn add_cubic_bez(&mut self, new_curve: &CubicBez, color: u32) + { + self.items.push(Box::new(DrawableCubicBez(new_curve.clone(), color, false))); + } + + /// Adds a quadratic bezier to the list of items + pub fn add_quad_bez(&mut self, new_curve: &QuadBez, color: u32) + { + self.items.push(Box::new(DrawableQuadBezier(new_curve.clone(), color, false))); + } + + /// Adds a line to the list of items + pub fn add_line(&mut self, new_curve: &Line, color: u32) + { + self.items.push(Box::new(DrawableLine(new_curve.clone(), color, false))); + } + + /// Adds a path segment to the list of items if possible + pub fn add_path_segment(&mut self, new_curve: &dyn Any, color: u32) + { + match get_curve_type_from_curve_param(new_curve) { + CurveType::CurveCubic => { self.add_cubic_bez(new_curve.downcast_ref::().unwrap(), color); }, + CurveType::CurveQuad => { self.add_quad_bez(new_curve.downcast_ref::().unwrap(), color); }, + CurveType::CurveLinear => { self.add_line(new_curve.downcast_ref::().unwrap(), color); }, + _ => { println!("Cannot add invalid cubic bezier to image!"); return }, + } + } + + /// Adds a phantom cubic bezier to the list of items + pub fn add_phantom_cubic_bez(&mut self, new_curve: &CubicBez, color: u32) + { + self.items.push(Box::new(DrawableCubicBez(new_curve.clone(), color, true))); + } + + /// Adds a phantom quadratic bezier to the list of items + pub fn add_phantom_quad_bez(&mut self, new_curve: &QuadBez, color: u32) + { + self.items.push(Box::new(DrawableQuadBezier(new_curve.clone(), color, true))); + } + + /// Adds a phantom line to the list of items + pub fn add_phantom_line(&mut self, new_curve: &Line, color: u32) + { + self.items.push(Box::new(DrawableLine(new_curve.clone(), color, true))); + } + + /// Adds a phantom path segment to the list of items if possible + pub fn add_phantom_path_segment(&mut self, new_curve: &dyn Any, color: u32) + { + match get_curve_type_from_curve_param(new_curve) { + CurveType::CurveCubic => { self.add_phantom_cubic_bez(new_curve.downcast_ref::().unwrap(), color); }, + CurveType::CurveQuad => { self.add_phantom_quad_bez(new_curve.downcast_ref::().unwrap(), color); }, + CurveType::CurveLinear => { self.add_phantom_line(new_curve.downcast_ref::().unwrap(), color); }, + _ => { println!("Cannot add invalid cubic bezier to image!"); return }, + } + } + + /// Adds a fat line to the list of items + pub fn add_fat_line(&mut self, fat_line: FatLine, color: u32) + { + self.items.push(Box::new(DrawableFatLine(fat_line, color))); + } + + /// Adds a point to the list of items + pub fn add_point(&mut self, pt: Point, color: u32, radius: f32, do_fill: bool) + { + self.items.push(Box::new(DrawablePoint(pt, color, radius, do_fill))); + } + + /// Adds a set of interections to the list of items + pub fn add_intersections( + &mut self, + intersections: &ArrayVec<(f64, f64), 9>, + c1: T, + c2: U, + color: u32, + radius: f32, + fill_first: bool) + { + self.items.push(Box::new(DrawableIntersections(intersections.clone(), c1, c2, color, radius, fill_first))); + } + + /// Adds a set of interections to the list of items if possible + pub fn add_intersections_cubic( + &mut self, + intersections: &ArrayVec<(f64, f64), 9>, + c1: &dyn Any, + c2: &dyn Any, + color: u32, + radius: f32, + fill_first: bool) + { + let t1 = get_curve_type_from_curve_param(c1); + let t2 = get_curve_type_from_curve_param(c2); + + match t1 { + CurveType::CurveCubic => { + match t2 { + CurveType::CurveCubic => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveQuad => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveLinear => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + _ => { println!("Cannot add invalid curve intersections to image!"); return }, + } + } + CurveType::CurveQuad => { + match t2 { + CurveType::CurveCubic => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveQuad => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveLinear => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + _ => { println!("Cannot add invalid curve intersections to image!"); return }, + } + } + CurveType::CurveLinear => { + match t2 { + CurveType::CurveCubic => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveQuad => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + CurveType::CurveLinear => { self.add_intersections(intersections, *c1.downcast_ref::().unwrap(), *c2.downcast_ref::().unwrap(), color, radius, fill_first); }, + _ => { println!("Cannot add invalid curve intersections to image!"); return }, + } + } + _ => { println!("Cannot add invalid curve intersections to image!"); return }, + } + } + + /// Gets the desired image size from the bounding rect of the items + fn image_size_from_bounds(&mut self, bounds: &Rect) -> Size + { + if self.use_fixed_size + { + let scale_x = (self.fixed_size.width - (self.margins as f64) * 2.) / bounds.width(); + let scale_y = (self.fixed_size.height - (self.margins as f64) * 2.) / bounds.height(); + self.scale = scale_x.min(scale_y); + self.offset = Vec2::new(self.fixed_size.width, self.fixed_size.height) * 0.5 - bounds.center().to_vec2() * self.scale; + return self.fixed_size; + } + else + { + self.offset = Vec2::new(self.margins as f64, self.margins as f64) - bounds.origin().to_vec2() * self.scale; + return bounds.size() * self.scale + Size::new(self.margins as f64, self.margins as f64) * 2.; + } + } + + /// Creates a color that is lighter than the given color + // fn lighter_color(color: u32) -> u32 + // { + // let r = (((color >> 16) as u8) as u32 + 0xFF) / 2; + // let g = (((color >> 8) as u8) as u32 + 0xFF) / 2; + // let b = (((color >> 0) as u8) as u32 + 0xFF) / 2; + // (r << 16) | (g << 8) | (b << 0) + // } + + /// Creates a raqote source from the given color + fn source_from_color(color: u32) -> Source<'static> + { + Source::Solid( + SolidSource { + r: (color >> 16) as u8, + g: (color >> 8) as u8, + b: (color >> 0) as u8, + a: 0xFF // full opacity + } + ) + } + + /// Creates a raqote source from the given color with a given opacity (0-255) + fn source_from_color_with_opacity(color: u32, alpha: u8) -> Source<'static> + { + Source::Solid( + SolidSource { + r: (color >> 16) as u8, + g: (color >> 8) as u8, + b: (color >> 0) as u8, + a: alpha + } + ) + } + + /// Strokes the given path onto the image + fn stroke_path_raw(dt: &mut DrawTarget, path: &Path, color: &Source<'static>) + { + dt.stroke( + &path, + &color, + &StrokeStyle { + cap: LineCap::Round, + join: LineJoin::Round, + width: 3., + miter_limit: 2., + dash_array: vec![], + dash_offset: 0., + }, + &DrawOptions::new() + ); + } + + /// Strokes the given path onto the image + fn stroke_path(dt: &mut DrawTarget, path: &Path, color: u32) + { + Image::stroke_path_raw(dt, path, &Image::source_from_color(color)); + } + + /// Strokes the given path onto the image + fn stroke_path_with_opacity(dt: &mut DrawTarget, path: &Path, color: u32, alpha: u8) + { + Image::stroke_path_raw(dt, path, &Image::source_from_color_with_opacity(color, alpha)); + } + + /// Strokes the given path onto the image with a dash pattern + fn stroke_dashed_path_raw(dt: &mut DrawTarget, path: &Path, color: &Source<'static>) + { + dt.stroke( + &path, + &color, + &StrokeStyle { + cap: LineCap::Round, + join: LineJoin::Round, + width: 3., + miter_limit: 2., + dash_array: vec![10., 5.], + dash_offset: 0., + }, + &DrawOptions::new() + ); + } + + /// Strokes the given path onto the image with a dash pattern + fn stroke_dashed_path(dt: &mut DrawTarget, path: &Path, color: u32) + { + Image::stroke_dashed_path_raw(dt, path, &Image::source_from_color(color)); + } + + /// Fills the given path onto the image + fn fill_path_raw(dt: &mut DrawTarget, path: &Path, color: &Source<'static>) + { + dt.fill( + &path, + &color, + &DrawOptions::new() + ); + } + + /// Fills the given path onto the image + // fn fill_path(dt: &mut DrawTarget, path: &Path, color: u32) + // { + // Image::fill_path_raw(dt, path, &Image::source_from_color(color)); + // } + + /// Fills the given path onto the image with an opacity (0-255) + fn fill_path_with_opacity(dt: &mut DrawTarget, path: &Path, color: u32, alpha: u8) + { + Image::fill_path_raw(dt, path, &Image::source_from_color_with_opacity(color, alpha)); + } + + /// Gets the affine transform necessary to draw an item onto the image + fn get_affine_transform(&self) -> Affine + { + Affine::translate(self.offset) * Affine::scale(self.scale) + } + + /// Gets the affine transform necessary to draw an item onto the image + fn get_transformed_point(&self, pt: &Point) -> Point + { + self.get_affine_transform() * *pt + } + + /// Draws a point onto the image + fn draw_point_raw(&self, dt: &mut DrawTarget, pt: &Point, radius: f32, color: &Source<'static>, do_fill: bool) + { + let transformed_point = self.get_transformed_point(pt); + + let mut pb = PathBuilder::new(); + pb.arc( + transformed_point.x as f32, transformed_point.y as f32, + radius, + 0., (360.0_f32).to_radians() + ); + let path = pb.finish(); + + Image::stroke_path_raw(dt, &path, &color); + if do_fill + { + Image::fill_path_raw(dt, &path, &color); + } + } + + /// Draws a point onto the image + fn draw_point(&self, dt: &mut DrawTarget, pt: &Point, radius: f32, color: u32, do_fill: bool) + { + self.draw_point_raw(dt, pt, radius, &Image::source_from_color(color), do_fill); + } + + /// Draws a point onto the image with an opacity (0-255) + fn draw_point_with_opacity(&self, dt: &mut DrawTarget, pt: &Point, radius: f32, color: u32, alpha: u8, do_fill: bool) + { + self.draw_point_raw(dt, pt, radius, &Image::source_from_color_with_opacity(color, alpha), do_fill); + } + + /// Draws a line onto the image + fn draw_line_raw(&self, dt: &mut DrawTarget, line: &Line, color: &Source<'static>) + { + let transformed_line = Line::new( + self.get_transformed_point(&line.p0), + self.get_transformed_point(&line.p1) + ); + + let mut pb = PathBuilder::new(); + pb.move_to(transformed_line.p0.x as f32, transformed_line.p0.y as f32); + pb.line_to(transformed_line.p1.x as f32, transformed_line.p1.y as f32); + let path = pb.finish(); + + Image::stroke_path_raw(dt, &path, &color); + } + + /// Draws a line onto the image + fn draw_line(&self, dt: &mut DrawTarget, line: &Line, color: u32) + { + self.draw_line_raw(dt, line, &Image::source_from_color(color)); + } + + /// Draws a line onto the image with an opacity (0-255) + fn draw_line_with_opacity(&self, dt: &mut DrawTarget, line: &Line, color: u32, alpha: u8) + { + self.draw_line_raw(dt, line, &Image::source_from_color_with_opacity(color, alpha)); + } + + /// Draws a dashed line onto the image + // fn draw_dashed_line_raw(&self, dt: &mut DrawTarget, line: &Line, color: &Source<'static>) + // { + // let transformed_line = Line::new( + // self.get_transformed_point(&line.p0), + // self.get_transformed_point(&line.p1) + // ); + + // let mut pb = PathBuilder::new(); + // pb.move_to(transformed_line.p0.x as f32, transformed_line.p0.y as f32); + // pb.line_to(transformed_line.p1.x as f32, transformed_line.p1.y as f32); + // let path = pb.finish(); + + // Image::stroke_dashed_path_raw(dt, &path, &color); + // } + + /// Draws a dashed line onto the image + // fn draw_dashed_line(&self, dt: &mut DrawTarget, line: &Line, color: u32) + // { + // self.draw_dashed_line_raw(dt, line, &Image::source_from_color(color)); + // } + + /// Draws a dashed line onto the image with an opacity (0-255) + // fn draw_dashed_line_with_opacity(&self, dt: &mut DrawTarget, line: &Line, color: u32, alpha: u8) + // { + // self.draw_dashed_line_raw(dt, line, &Image::source_from_color_with_opacity(color, alpha)); + // } + + /// Draws a cubic bezier onto the image + fn draw_cubic_bez(&self, dt: &mut DrawTarget, curve: &CubicBez, color: u32) + { + let transformed_curve = self.get_affine_transform() * *curve; + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_curve.p0.x as f32, transformed_curve.p0.y as f32 + ); + pb.cubic_to( + transformed_curve.p1.x as f32, transformed_curve.p1.y as f32, + transformed_curve.p2.x as f32, transformed_curve.p2.y as f32, + transformed_curve.p3.x as f32, transformed_curve.p3.y as f32 + ); + let path = pb.finish(); + + Image::stroke_path(dt, &path, color); + } + + /// Draws a cubic bezier onto the image + fn draw_cubic_bez_with_opacity(&self, dt: &mut DrawTarget, curve: &CubicBez, color: u32, alpha: u8) + { + let transformed_curve = self.get_affine_transform() * *curve; + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_curve.p0.x as f32, transformed_curve.p0.y as f32 + ); + pb.cubic_to( + transformed_curve.p1.x as f32, transformed_curve.p1.y as f32, + transformed_curve.p2.x as f32, transformed_curve.p2.y as f32, + transformed_curve.p3.x as f32, transformed_curve.p3.y as f32 + ); + let path = pb.finish(); + + Image::stroke_path_with_opacity(dt, &path, color, alpha); + } + + /// Draws a quadratic bezier onto the image + fn draw_quad_bez(&self, dt: &mut DrawTarget, curve: &QuadBez, color: u32) + { + let transformed_curve = self.get_affine_transform() * *curve; + //let transformed_curve = (*curve * self.scale) + self.offset; + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_curve.p0.x as f32, transformed_curve.p0.y as f32 + ); + pb.quad_to( + transformed_curve.p1.x as f32, transformed_curve.p1.y as f32, + transformed_curve.p2.x as f32, transformed_curve.p2.y as f32 + ); + let path = pb.finish(); + Image::stroke_path(dt, &path, color); + } + + /// Draws a quadratic bezier onto the image + fn draw_quad_bez_with_opacity(&self, dt: &mut DrawTarget, curve: &QuadBez, color: u32, alpha: u8) + { + let transformed_curve = self.get_affine_transform() * *curve; + //let transformed_curve = (*curve * self.scale) + self.offset; + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_curve.p0.x as f32, transformed_curve.p0.y as f32 + ); + pb.quad_to( + transformed_curve.p1.x as f32, transformed_curve.p1.y as f32, + transformed_curve.p2.x as f32, transformed_curve.p2.y as f32 + ); + let path = pb.finish(); + Image::stroke_path_with_opacity(dt, &path, color, alpha); + } + + /// Offsets a line by a signed distance + fn offset_line(line: &Line, distance: f64) -> Line + { + let vec = line.p1 - line.p0; + let vec = Vec2::new(vec.y, -vec.x); + let vec = (vec / vec.hypot()) * -distance; + + line.clone() + vec + } + + /// Extends a line by a signed distance + fn extend_line(line: &Line, distance: f64) -> Line + { + let t_extend = distance / (line.p1 - line.p0).hypot(); + Line::new(line.eval(-t_extend), line.eval(1. + t_extend)) + } + + /// Draws a fat line onto the image + fn draw_fat_line(&self, dt: &mut DrawTarget, fl: &FatLine, color: u32) + { + let img_size = Vec2::new(dt.width() as f64, dt.height() as f64).hypot(); + + let line = Line::new(fl.p0, fl.p1); + let line_dmin = Image::offset_line(&line, fl.dmin); + let line_dmax = Image::offset_line(&line, fl.dmax); + + let transform = self.get_affine_transform(); + let transformed_line = Image::extend_line(&(transform * line), img_size); + let transformed_line_dmin = Image::extend_line(&(transform * line_dmin), img_size); + let transformed_line_dmax = Image::extend_line(&(transform * line_dmax), img_size); + // let transformed_line = transform * line; + // let transformed_line_dmin = transform * line_dmin; + // let transformed_line_dmax = transform * line_dmax; + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_line.p0.x as f32, transformed_line.p0.y as f32, + ); + pb.line_to( + transformed_line.p1.x as f32, transformed_line.p1.y as f32, + ); + let path = pb.finish(); + Image::stroke_path(dt, &path, color); + + let mut pb = PathBuilder::new(); + pb.move_to( + transformed_line_dmin.p0.x as f32, transformed_line_dmin.p0.y as f32, + ); + pb.line_to( + transformed_line_dmin.p1.x as f32, transformed_line_dmin.p1.y as f32, + ); + pb.line_to( + transformed_line_dmax.p1.x as f32, transformed_line_dmax.p1.y as f32, + ); + pb.line_to( + transformed_line_dmax.p0.x as f32, transformed_line_dmax.p0.y as f32, + ); + pb.line_to( + transformed_line_dmin.p0.x as f32, transformed_line_dmin.p0.y as f32, + ); + let path = pb.finish(); + // Image::fill_path_with_opacity(dt, &path, Image::lighter_color(color), 0xB0); + Image::fill_path_with_opacity(dt, &path, color, 0xB0); + Image::stroke_dashed_path(dt, &path, color); + // Image::stroke_path(dt, &path, color); + } + + /// Saves the image to a file + pub fn save(&mut self, file_name: &str) + { + let mut bounds = Rect::ZERO; + for drawable in self.items.iter() + { + let dbounds = drawable.bounds(); + if bounds == Rect::ZERO + { + bounds = dbounds; + } + else if dbounds != Rect::ZERO + { + bounds = bounds.union(dbounds); + } + } + + if bounds.size() == Size::ZERO + { + println!("cannot create image from zero bounds: {}", bounds); + return; + } + + let image_size = self.image_size_from_bounds(&bounds); + let mut dt = DrawTarget::new(image_size.width as i32, image_size.height as i32); + + let mut background = PathBuilder::new(); + background.rect(0.0, 0.0, dt.width() as f32, dt.height() as f32); + background.close(); + let background = background.finish(); + dt.fill( + &background, + &Source::Solid(SolidSource { r: 0xFF, g: 0xFF, b: 0xFF, a: 0xFF }), + &DrawOptions::new() + ); + + for drawable in self.items.iter() + { + drawable.draw_to_image(self, &mut dt); + } + + let save_success = dt.write_png(file_name); + match save_success + { + Ok(_) => println!("Image saved! {}", file_name), + Err(e) => println!("Unable to save image: {:?}", e), + } + } +} + +/***************************************************************************** + * Format + *****************************************************************************/ + +impl fmt::Display for Image +{ + /// Defines how Image is printed to the output + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result + { + write!(f, "Img(scale: {:?}, margins: {}, items: {})", self.scale, self.margins, self.items.len()) + } +} \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 8432385f..38908032 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -84,7 +84,7 @@ mod bezpath; mod circle; pub mod common; mod cubicbez; -mod curve_intersections; +pub mod curve_intersections; mod ellipse; mod insets; mod line; diff --git a/src/pathops.rs b/src/pathops.rs new file mode 100644 index 00000000..6de2ca2a --- /dev/null +++ b/src/pathops.rs @@ -0,0 +1,429 @@ +/// Boolean operations on BezPaths +/// +/// The methodology approximately follows that which is described here: +/// https://losingfight.com/blog/2011/07/07/how-to-implement-boolean-operations-on-bezier-paths-part-1/ +/// https://losingfight.com/blog/2011/07/08/how-to-implement-boolean-operations-on-bezier-paths-part-2/ +/// https://losingfight.com/blog/2011/07/09/how-to-implement-boolean-operations-on-bezier-paths-part-3/ + +use crate::{Point, Rect, BezPath, PathEl, Shape}; +use crate::{Image, Size, Line}; // DO NOT COMMIT +use std::sync::atomic::{AtomicUsize, Ordering}; // DO NOT COMMIT + +/// An operation that can be performed on a path +#[derive(Clone, Debug, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum PathOperation<'a> { + /// Subtracts the other path from the current path + Subtract(&'a BezPath), + /// Like subtract but doesn't close current path along second path + Outside(&'a BezPath), + /// Gets the inersection of this path with another + Intersect(&'a BezPath), + /// Like intersection but doesn't close current path along second path + Inside(&'a BezPath), + /// Unites this path with another + UniteWith(&'a BezPath), + /// Gets the exclusive-or operation with another path + Xor(&'a BezPath), +} + +bitflags::bitflags! { + // private struct, used in do_path_op + struct ClosedPathOperationsFlags: u32 { + const DISCARD_SEGS = 0b00000000; + const KEEP_A_OUTSIDE_B = 0b00000001; + const KEEP_A_INSIDE_B = 0b00000010; + const KEEP_B_INSIDE_A = 0b00000100; + const KEEP_B_OUTSIDE_A = 0b00001000; + const KEEP_ALL_SEGS = Self::KEEP_A_OUTSIDE_B.bits + | Self::KEEP_A_INSIDE_B.bits + | Self::KEEP_B_INSIDE_A.bits + | Self::KEEP_B_OUTSIDE_A.bits; + } +} + +/// Compute an operation between two paths +pub fn path_path_operation( + path_a: &BezPath, + operation: PathOperation, +) -> BezPath { + match operation { + PathOperation::Subtract(path_b) => { + // Subtract b from a + do_path_op( + path_a, + path_b, + ClosedPathOperationsFlags::KEEP_A_OUTSIDE_B + | ClosedPathOperationsFlags::KEEP_B_INSIDE_A, + true, + 0., + ) + } + PathOperation::Outside(path_b) => { + // Get a outside of b + do_path_op( + path_a, + path_b, + ClosedPathOperationsFlags::KEEP_A_OUTSIDE_B, + false, + 0., + ) + } + PathOperation::Intersect(path_b) => { + // Get the intersections of a and b + do_path_op( + path_a, + path_b, + ClosedPathOperationsFlags::KEEP_A_INSIDE_B + | ClosedPathOperationsFlags::KEEP_B_INSIDE_A, + true, + 0., + ) + } + PathOperation::Inside(path_b) => { + // Get a inside of b + do_path_op( + path_a, + path_b, + ClosedPathOperationsFlags::KEEP_A_INSIDE_B, + false, + 0., + ) + } + PathOperation::UniteWith(path_b) => { + // Units a and b + do_path_op( + path_a, + path_b, + ClosedPathOperationsFlags::KEEP_A_OUTSIDE_B + | ClosedPathOperationsFlags::KEEP_B_OUTSIDE_A, + true, + 0., + ) + } + PathOperation::Xor(path_b) => { + // Get a XOR b + // We can actually do a simple shortcut with this one. All we have + // to do is append path_b to path_a. Since they are both even-odd + // fill, everything is taken care of just by the logic of the fill. + let mut path_final = path_a.clone(); + path_final.extend(path_b.iter()); + path_final + } + } +} + +#[derive(Clone, Debug)] +struct PathSectionInfo<'a> { + /// The index in the array of elements so that we can quickly and easily + /// modify the array + pub index: usize, + /// Path index, ex: 0 for path1, and 1 for path2 + pub path_index: usize, + /// The path for the section + pub section_path: &'a BezPath, + /// The bounding rect of the path section + pub section_bounds: Rect, + /// A flag representing whether we have processed this section + pub processed: bool, + /// A flag determining whether we will keep this section in the final + /// result, only valid if processed == true + pub keep: bool, +} + +impl PathSectionInfo<'_> { + pub fn new( + index: usize, + path_index: usize, + section_path: &BezPath, + ) -> PathSectionInfo { + PathSectionInfo { + index: index, + path_index: path_index, + section_path: section_path, + section_bounds: section_path.bounding_box(), + processed: false, + keep: false + } + } +} + +/// Populates +fn populate_elements_list<'a>( + elements: &'a Vec>, + path_index: usize, + sections: &'a Vec, +) -> Vec> { + let mut new_elements = elements.clone(); + let mut test = sections.iter() + .enumerate() + .map(|(index, path)| { + // Create an info struct + PathSectionInfo::new( + elements.len() + index, + path_index, + &path, + ) + }) + .collect::>(); + new_elements.append(&mut test); + new_elements +} + +/// Determines whihc +fn flag_path_sections_to_keep( + elements: &mut Vec, + process_winding_results: &dyn Fn(Vec<(f64, i32, usize, usize, bool)>, &mut Vec), + accuracy: f64, +) -> Vec { + while let Some(info1) = elements.iter().find(|info| !info.processed) { + // We've picked the first unprocessed element and now we need to draw a + // ray through it to test if we should keep it + + // Create a ray that travels through our path + let ray_is_horizontal = info1.section_bounds.width().abs() <= info1.section_bounds.height().abs(); + let ray_position = if ray_is_horizontal { (info1.section_bounds.min_y() + info1.section_bounds.max_y()) / 2. } else { (info1.section_bounds.min_x() + info1.section_bounds.max_x()) / 2. }; + let ray_point = if ray_is_horizontal { Point::new(info1.section_bounds.max_x() + 1., ray_position) } else { Point::new(ray_position, info1.section_bounds.max_y() + 1.) }; + let ray_fn = if ray_is_horizontal { BezPath::winding_x_results } else { BezPath::winding_y_results }; + + // Loop through all paths and determine which ones intersect our ray + let mut winding_results = elements.iter().by_ref() + .map(|info2| { + let winding_array = ray_fn(&*info2.section_path, ray_point).into_iter() + .map(|(pt, winding)| (pt, winding, info2.index.clone(), info2.path_index.clone(), info2.processed.clone())) + .collect::>(); + winding_array + }) + .flatten() + .collect::>(); + + // Sort by the position of the intersection with the ray + winding_results.sort_by(|(a, _, _, _, _), (b, _, _, _, _)| (a).partial_cmp(b).unwrap()); + + process_winding_results( + winding_results, + elements + ); + } + + // Remove flagged sections of the path, add a bool representing whether or + // not the path section should be connected to the previous one in the + // final result. + let mut elements = elements.into_iter() + .filter_map(|info| if info.keep { Some((info.section_path, false, false)) } else { None }) + .filter(|(path, _, _)| !path.elements().is_empty()) + .collect::>(); + + // Re-close paths + for index1 in 0..elements.len() { + if let Some(path_end1) = elements[index1].0.get_seg_end(elements[index1].0.elements().len() - 1) { + let accuracy2 = accuracy * accuracy; + let mut candidates = elements.iter() + .enumerate() + .filter(|(index2, _)| *index2 > index1) + .map(|(index2, (path2, _, _))| (index2, path2.get_seg_end(0).unwrap(), path2.get_seg_end(path2.elements().len() - 1).unwrap())) + .filter(|&(_, start, end)| (start - path_end1).hypot2() <= accuracy2 || (end - path_end1).hypot2() <= accuracy2) + .collect::>(); + candidates.sort_by(|&(_, start1, end1), &(_, start2, end2)| { + let d1 = (start1 - path_end1).hypot2().min((end1 - path_end1).hypot2()); + let d2 = (start2 - path_end1).hypot2().min((end2 - path_end1).hypot2()); + d1.partial_cmp(&d2).unwrap() + }); + + if let Some(&(connected_index, start, end)) = candidates.get(0) { + // Set a flag so that we connect to path1 later by removing path2's move-to + elements[connected_index].1 = true; + + // Determine if we need to reverse the segment + elements[connected_index].2 = (start - path_end1).hypot2() > (end - path_end1).hypot2(); + + // Move index to be next in line + elements.swap(index1 + 1, connected_index); + } + } + }; + + if !elements.is_empty() { + assert!(!elements[0].1); // Fist element can't conenct to the previous one + } + + elements.into_iter() + .map(|(path, connect_to_previous, reverse)| { + let p_rev = if reverse { path.reverse() } else { BezPath::new() }; + let p_correct = if reverse { &p_rev } else { path }; + BezPath::from_vec(p_correct.elements()[if connect_to_previous { 1 } else { 0 }..].to_vec()) + }) + .flatten() + .collect::>() +} + +/// This function converts a winding-fill path into an even-odd filled path +/// Note: Since there is no winding fill property of path, the path is assumed +/// to have a winding fill when this function is called +/// Note: Open paths are closed +pub fn convert_path_to_even_odd( + path: &BezPath, + accuracy: f64, +) -> BezPath { + let bounds = path.bounding_box(); + let accuracy = Point::epsilon(Point::new(bounds.max_x(), bounds.max_y()), accuracy); + + // Make a copy of the path + let mut path = path.clone(); + + // Close all subpaths because this path is treated as being closed anyways + path.close_subpaths(); + + // Break apart at self intersections + BezPath::break_at_self_intersections(&mut path, accuracy); + let paths = path.split_at_moves(); + + // Convert the list of elements into a vector of tuples to make tracking elements easier + let elements = Vec::::new(); + let elements = populate_elements_list(&elements, 0, &paths); + let mut elements = elements; // Make mutable + + fn process_winding_results( + winding_results: Vec<(f64, i32, usize, usize, bool)>, + elements: &mut Vec) { + winding_results + .into_iter() + .scan(0, |winding_total, (_, w, index, _, processed2)| { + let winding_before = *winding_total; + *winding_total = *winding_total + w; + Some((winding_before, *winding_total, index, processed2)) + }) + .for_each(|(w_before, w_current, index, _)| { + let ref mut info = elements[index]; + + // Discard if crossing the path didn't change whether or not we're inside the original path + info.keep = (w_before != 0) != (w_current != 0); + + // Mark as processed + info.processed = true; + }); + } + + let elements = flag_path_sections_to_keep( + &mut elements, + &process_winding_results, + accuracy + ); + BezPath::from_vec(elements) +} + +/// Performs an operation on the current path +/// Note: Paths are assumed to be even-odd filled +fn do_path_op( + path_a: &BezPath, + path_b: &BezPath, + closed_flags: ClosedPathOperationsFlags, + connect_nearby_paths: bool, + accuracy: f64, +) -> BezPath { + let bounds = path_a.bounding_box().union(path_b.bounding_box()); + let accuracy = Point::epsilon(Point::new(bounds.max_x(), bounds.max_y()), accuracy); + + // Make copies of the paths + let mut path_a = path_a.clone(); + let mut path_b = path_b.clone(); + + if connect_nearby_paths { + // We treat all paths as if they are closed, so we'll just go ahead and + // close the paths + path_a.close_subpaths(); + path_b.close_subpaths(); + } + + // Find all intersections between the two paths and break them apart at + // those locations + BezPath::break_at_intersections(&mut path_a, &mut path_b, accuracy); + let paths_a = path_a.split_at_moves(); + let paths_b = path_b.split_at_moves(); + + // Make a list of all elements in the paths + let elements = Vec::::new(); + let elements = populate_elements_list(&elements, 0, &paths_a); + let elements = populate_elements_list(&elements, 1, &paths_b); + let mut elements = elements; // Make mutable + + let process_winding_results = | + winding_results: Vec<(f64, i32, usize, usize, bool)>, + elements: &mut Vec, + | { + winding_results.into_iter() + .scan((0, 0), |winding_total, (_, w, index, path_index, processed2)| { + if path_index == 0 { winding_total.0 = winding_total.0 + w } else { winding_total.1 = winding_total.1 + w }; + Some((*winding_total, index, path_index, processed2)) + }) + .for_each(|(w_current, index, path_index, _)| { + let ref mut info = elements[index]; + + // Use closed_flags to determine if we need to remove or keep the path section + let is_path_a = path_index == 0; + let inside_path_a = (w_current.0 % 2) != 0; + let inside_path_b = (w_current.1 % 2) != 0; + info.keep = (is_path_a && closed_flags.contains(ClosedPathOperationsFlags::KEEP_A_OUTSIDE_B) && !inside_path_b) + || (is_path_a && closed_flags.contains(ClosedPathOperationsFlags::KEEP_A_INSIDE_B) && inside_path_b) + || (!is_path_a && closed_flags.contains(ClosedPathOperationsFlags::KEEP_B_INSIDE_A) && inside_path_a) + || (!is_path_a && closed_flags.contains(ClosedPathOperationsFlags::KEEP_B_OUTSIDE_A) && !inside_path_a); + + + // Mark as processed + info.processed = true; + }); + }; + + let elements = flag_path_sections_to_keep( + &mut elements, + &process_winding_results, + accuracy + ); + BezPath::from_vec(elements) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{Shape, Circle}; + use crate::{ParamCurve}; // DO NOT COMMIT + + #[test] + fn test_circle_unite_circle() { + let c1 = Circle::new((0., 0.), 10.).to_path(1e-3); + let c2 = Circle::new((5., 0.), 10.).to_path(1e-3); + let res = path_path_operation(&c1, PathOperation::UniteWith(&c2)); + + // DO NOT COMMIT + static COUNTER_PATHOP: AtomicUsize = AtomicUsize::new(0); + let mut img1 = Image::from_fixed_size(Size::new(1000., 1000.), 20); + let intersects = c1.intersections(&c2, 0.); + for i in &intersects { + let seg1 = c1.get_seg(i.0.0).unwrap(); + let seg2 = c2.get_seg(i.1.0).unwrap(); + + img1.add_point(seg1.eval(i.0.1), 0xFF00FF, 10., true); + img1.add_point(seg2.eval(i.1.1), 0xFF00FF, 15., false); + } + img1.add_path(&res, 0xFF00FF, true); + img1.add_path(&c1, 0x800000, false); + img1.add_path(&c2, 0x008000, false); + img1.save(&format!("_path_op_{}.png", COUNTER_PATHOP.fetch_add(1, Ordering::Relaxed))); + + assert_eq!(res.elements().len(), 9); + } + + #[test] + fn test_convert_to_even_odd() { + let mut path = BezPath::new(); + path.move_to((0., 0.)); + path.line_to((0., 1.)); + path.curve_to((2., 0.5), (-1., 0.5), (1., 1.)); + path.line_to((-0.5, 0.)); + path.close_path(); + + let path2 = convert_path_to_even_odd(&path, 0.); + + assert_eq!(path2.elements().len(), 10); + } +} diff --git a/src/point.rs b/src/point.rs index 8f84912a..873fe2b3 100644 --- a/src/point.rs +++ b/src/point.rs @@ -174,6 +174,12 @@ impl Point { /// If we're comparing distances between samples of curves, our epsilon should /// depend on how big the points we're comparing are. This function returns an /// epsilon appropriate for the size of pt. + /// + /// Note: In general, epsilon will be equal to accuracy, but in cases where + /// accuracy is very small, we want a larger epsilon value to help account for + /// floating point errors (after adding, subtracting, etc.). This larger epsilon + /// value is 0.0000000001 times the largest value in the point which was + /// determined experimentally and may change in the future. pub fn epsilon(pt: Point, accuracy: f64) -> f64 { // Get max number in the point let max = f64::max(f64::abs(pt.x), f64::abs(pt.y)); diff --git a/src/quadbez.rs b/src/quadbez.rs index 0cfa09c0..289ee01f 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -95,6 +95,8 @@ impl QuadBez { /// Get the parameters such that the curve can be represented by the following formula: /// B(t) = a*t^2 + b*t + c + /// + /// Note: Values returned are in decresing exponent order pub fn parameters(&self) -> (Vec2, Vec2, Vec2) { let c = self.p0.to_vec2(); let b = (self.p1 - self.p0) * 2.0; @@ -281,7 +283,7 @@ impl ParamCurveArclen for QuadBez { let ba_c2 = b * a2 + c2; if a2.is_infinite() { - // The arclength is zero + // The arclength is zero, this happens when all control points lie on top of each other return 0.; } @@ -451,6 +453,13 @@ mod tests { ); } + #[test] + fn quadbez_arclen_zero() { + let q = QuadBez::new((-1.0, 0.0), (-1.0, 0.0), (-1.0, 0.0)); + let accuracy = 1e-11; + assert_eq!(q.arclen(accuracy), 0.); + } + #[test] fn quadbez_subsegment() { let q = QuadBez::new((3.1, 4.1), (5.9, 2.6), (5.3, 5.8)); From ba180f7a6ed9aaec505728452be01872c0cb8b9d Mon Sep 17 00:00:00 2001 From: Tyler Scott Date: Tue, 23 Aug 2022 13:01:59 -0600 Subject: [PATCH 20/20] Cargo clippy and fmt fixes --- src/common.rs | 6 +++--- src/cubicbez.rs | 2 +- src/curve_intersections.rs | 12 ++++++++++++ src/quadbez.rs | 2 +- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/common.rs b/src/common.rs index 5f786729..506deb33 100644 --- a/src/common.rs +++ b/src/common.rs @@ -51,7 +51,7 @@ impl FloatExt for f32 { /// See: /// /// Return values of x for which c0 + c1 x + c2 x² + c3 x³ = 0. -/// +/// /// Note: The input arguments are in increasing exponent order pub fn solve_cubic(c0: f64, c1: f64, c2: f64, c3: f64) -> ArrayVec { let mut result = ArrayVec::new(); @@ -105,7 +105,7 @@ pub fn solve_cubic(c0: f64, c1: f64, c2: f64, c3: f64) -> ArrayVec { /// the other root might be out of representable range. In the degenerate /// case where all coefficients are zero, so that all values of x satisfy /// the equation, a single `0.0` is returned. -/// +/// /// Note: The input arguments are in increasing exponent order pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { let mut result = ArrayVec::new(); @@ -152,7 +152,7 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { /// Find real roots of linear equation. /// /// Return values of x for which c0 + c1 x = 0. -/// +/// /// Note: The input arguments are in increasing exponent order pub fn solve_linear(c0: f64, c1: f64) -> ArrayVec { let mut result = ArrayVec::new(); diff --git a/src/cubicbez.rs b/src/cubicbez.rs index f1b9d925..b9b7714b 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -239,7 +239,7 @@ impl CubicBez { /// Get the parameters such that the curve can be represented by the following formula: /// B(t) = a*t^3 + b*t^2 + c*t + d - /// + /// /// Note: Values returned are in decresing exponent order pub fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) { let c = (self.p1 - self.p0) * 3.0; diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs index 19e89ed5..3b5257b4 100644 --- a/src/curve_intersections.rs +++ b/src/curve_intersections.rs @@ -256,18 +256,29 @@ impl ParamCurveBezierClipping for CubicBez { } bitflags::bitflags! { + /// Flags that dictate how certain edge cases are handled when finding the intersections + /// between two path segments. pub struct CurveIntersectionFlags: u32 { + /// The default mode const NONE = 0b00000000; + /// Intersections at t==0 on curve 1 are kept const KEEP_CURVE1_T0_INTERSECTION = 0b00000001; + /// Intersections at t==1 on curve 1 are kept const KEEP_CURVE1_T1_INTERSECTION = 0b00000010; + /// Intersections at t==0 on curve 2 are kept const KEEP_CURVE2_T0_INTERSECTION = 0b00000100; + /// Intersections at t==1 on curve 2 are kept const KEEP_CURVE2_T1_INTERSECTION = 0b00001000; + /// Intersections at t==0 and t==1 on curve 1 are kept const KEEP_CURVE1_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE1_T0_INTERSECTION.bits | Self::KEEP_CURVE1_T1_INTERSECTION.bits; + /// Intersections at t==0 and t==1 on curve 2 are kept const KEEP_CURVE2_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE2_T0_INTERSECTION.bits | Self::KEEP_CURVE2_T1_INTERSECTION.bits; + /// Intersections at t==0 and t==1 on curve 1 and curve 2 are kept const KEEP_ALL_ENDPOINT_INTERSECTIONS = Self::KEEP_CURVE1_ENDPOINT_INTERSECTIONS.bits | Self::KEEP_CURVE2_ENDPOINT_INTERSECTIONS.bits; + /// All intersections that lie at the same location are kept const KEEP_DUPLICATE_INTERSECTIONS = 0b00010000; } } @@ -777,6 +788,7 @@ fn add_point_curve_intersection( pt: Point, curve: &T, diff --git a/src/quadbez.rs b/src/quadbez.rs index 289ee01f..61ae5967 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -95,7 +95,7 @@ impl QuadBez { /// Get the parameters such that the curve can be represented by the following formula: /// B(t) = a*t^2 + b*t + c - /// + /// /// Note: Values returned are in decresing exponent order pub fn parameters(&self) -> (Vec2, Vec2, Vec2) { let c = self.p0.to_vec2();