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/Cargo.toml b/Cargo.toml index 3e214f17..83682bfb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,6 +15,7 @@ features = ["mint", "schemars", "serde"] [dependencies] arrayvec = "0.7.1" +bitflags = "1.3" [dependencies.mint] version = "0.5.1" diff --git a/src/bezpath.rs b/src/bezpath.rs index 7ba3db0f..1af37e81 100644 --- a/src/bezpath.rs +++ b/src/bezpath.rs @@ -1285,6 +1285,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(); @@ -1307,6 +1308,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)); diff --git a/src/common.rs b/src/common.rs index 323b96a0..506deb33 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,18 +105,16 @@ 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(); 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; } @@ -149,6 +149,23 @@ pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> ArrayVec { result } +/// 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, this happens the line is horizontal, hence when c1 == 0. + 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 af02009a..b9b7714b 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -1,13 +1,12 @@ //! 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; -use crate::common::GAUSS_LEGENDRE_COEFFS_9; +use crate::common::{solve_quadratic, GAUSS_LEGENDRE_COEFFS_9}; use crate::{ Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, QuadBez, Rect, Shape, @@ -238,7 +237,11 @@ 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 + /// + /// 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; let d = self.p0.to_vec2(); @@ -310,6 +313,13 @@ impl CubicBez { self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan() || self.p3.is_nan() } + /// Is this cubic Bezier curve a line? + #[inline] + 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 + } + /// Determine the inflection points. /// /// Return value is t parameter for the inflection points of the curve segment. diff --git a/src/curve_intersections.rs b/src/curve_intersections.rs new file mode 100644 index 00000000..3b5257b4 --- /dev/null +++ b/src/curve_intersections.rs @@ -0,0 +1,1388 @@ +//! 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::{ + 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 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); + /// 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_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 + // 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! { + /// 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; + } +} + +/// Compute the intersections between two Bézier curves +pub fn curve_curve_intersections( + curve1: &T, + curve2: &U, + flags: CurveIntersectionFlags, + accuracy: f64, +) -> ArrayVec<(f64, f64), 9> { + let mut av = ArrayVec::new(); + + // Make sure we don't use too small of an accuracy + let bounds = curve1.bounding_box().union(curve2.bounding_box()); + let accuracy = Point::epsilon(Point::new(bounds.max_x(), bounds.max_y()), accuracy); + + // 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, + &(0.0..1.0), + &(0.0..1.0), + &mut av, + false, + 0, + 0, + curve1, + curve2, + flags, + accuracy, + ); + av +} + +/// This function tests if one curve lies along another curve for some range +fn check_for_overlap( + 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: + // 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 + fn t_value_on_other( + 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, accuracy); + if !t_values_other.0 { + return (false, 0.0); + } + + 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]); + } + + // 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, 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 } + + 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 { + 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_near(curve1_pt, curve2_pt, accuracy) { + return false; + } + } + + 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! +} + +// 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, + flags: CurveIntersectionFlags, + accuracy: f64, +) -> u32 { + call_count += 1; + recursion_count += 1; + if call_count >= 4096 || recursion_count >= 128 { + return call_count; + } + + let arclen_epsilon = 1e-3; + + if orig_curve2 + .subsegment(domain2.start..domain2.end) + .arclen(arclen_epsilon) + <= accuracy + { + add_point_curve_intersection( + (curve2, curve1), + (domain2, domain1), + intersections, + !flip, // Invert because curve order is inverted + (orig_curve2, orig_curve1), + flags, + accuracy, + ); + 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, + flags, + accuracy, + ); + 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, + flags, + accuracy, + ); + 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; + } + + // 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)); + + // Reduce curve1 to the part that might intersect curve2. + let curve1 = &orig_curve1.subsegment(new_domain1.clone()); + + // Check if both curves are very small + let curve1_arclen = curve1.arclen(arclen_epsilon); + let curve2_arclen = curve2.arclen(arclen_epsilon); + if curve1_arclen <= accuracy && curve2_arclen <= accuracy { + let t1 = (new_domain1.start + new_domain1.end) * 0.5; + let t2 = (domain2.start + domain2.end) * 0.5; + + let curve1_pt = orig_curve1.eval(t1); + let curve2_pt = orig_curve2.eval(t2); + + if Point::is_near(curve1_pt, curve2_pt, accuracy) { + // Note: add_intersection tests if the intersection is an end-point + add_intersection(t1, orig_curve1, t2, orig_curve2, flip, intersections, flags); + return call_count; + } + } + + // (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 curve1.arclen(arclen_epsilon * 0.001) <= accuracy * 2. { + add_point_curve_intersection( + (curve1, curve2), + (new_domain1, domain2), + intersections, + flip, + (orig_curve1, orig_curve2), + flags, + accuracy, + ); + 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, + flags, + accuracy, + ); + 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, + flags, + accuracy, + ); + } 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, + flags, + accuracy, + ); + 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, + flags, + accuracy, + ); + } + } else { + // Iterate. + if curve2_arclen > accuracy { + call_count = add_curve_intersections( + curve2, + curve1, + domain2, + new_domain1, + intersections, + !flip, + recursion_count, + call_count, + orig_curve2, + orig_curve1, + flags, + accuracy, + ); + } 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, + flags, + accuracy, + ); + } + } + + call_count +} + +fn add_point_curve_intersection( + curves: (&T, &U), + domains: (&Range, &Range), + intersections: &mut ArrayVec<(f64, f64), 9>, + flip: bool, + orig_curves: (&T, &U), + flags: CurveIntersectionFlags, + accuracy: f64, +) { + // Get the mid-point + 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, curves.1, accuracy, false); + let curve_t = results + .iter() + .fold((-1., f64::MAX), |(t, d_sq), 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); + } + (t, d_sq) + }) + .0; + + if curve_t == -1. { + // No intersection found + return; + } + + let curve_t = domain_value_at_t(domains.1, curve_t); // convert to proper domain + add_intersection( + pt_t, + orig_curves.0, + curve_t, + orig_curves.1, + flip, + intersections, + flags, + ); +} + +/// Return all t values where the specified point lies on the curve +pub fn t_along_curve_for_point( + pt: Point, + curve: &T, + accuracy: f64, + near_pts_only: bool, +) -> ArrayVec { + let mut result = ArrayVec::new(); + + // If both endpoints are approximately close, we only return 0.0. + if Point::is_near(pt, curve.start(), accuracy) { + result.push(0.0); + return result; + } + if Point::is_near(pt, curve.end(), accuracy) { + result.push(1.0); + return result; + } + + // 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; + + 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 t.eq(u) { + 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 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, + accuracy: f64, + result: &mut ArrayVec, + ) -> bool { + if Point::is_near(curve.eval(t), pt, accuracy) { + result.push(t); + return true; + } + false + } + + for ex in curve.extrema() { + maybe_add(ex, pt, curve, accuracy, &mut result); + } + + result +} + +fn point_is_on_curve( + pt: Point, + curve: &T, + accuracy: f64, +) -> (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) { + return (true, t_values.clone()); // Point lies along curve + } + } + (false, t_values) // Point is not along curve +} + +fn add_intersection( + t1: f64, + orig_curve1: &T, + t2: f64, + orig_curve2: &U, + flip: bool, + 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 + 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) + || (flags.contains(CurveIntersectionFlags::KEEP_CURVE2_T1_INTERSECTION) && t2 >= 1.0); + if !keep_intersection { + return; + } + + 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 + // 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 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); + } + 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(&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) +} + +// 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 [Point], + hull_bottom: &mut [Point], + 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, 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, + count: usize, + flags: CurveIntersectionFlags, + accuracy: f64, + ) { + let arr1 = curve_curve_intersections(curve1, curve2, flags, accuracy); + assert_eq!(arr1.len(), count); + } + + fn do_double_test( + curve1: &T, + curve2: &U, + count: usize, + flags: CurveIntersectionFlags, + accuracy: f64, + ) { + 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); + } + + 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()); + + let mut found_t = false; + for t in &t_values { + let pt = curve.eval(*t); + + if (*t - t_test) <= 1e-9 { + found_t = true; + } + + assert!(Point::is_near(pt, pt_test, 1e-3)) + } + 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, 0.); + } + } + + fn test_overlapping(curve: &T, t1: f64, t2: f64, accuracy: f64) { + assert!(t1 < t2); + + do_double_test( + &curve.subsegment(0.0..t2), + &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( + &curve.subsegment(0.0..1.0), + &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, + ); + } + + #[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_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 + ); + 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. + 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 + ); + 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, + 0., + ); + + // (A previous version of the code was returning two practically identical + // intersection points here.) + do_double_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, + 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 + // 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_double_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, + CurveIntersectionFlags::NONE, // Standard algorithm + 0., // Tightest possible accuracy + ); + } +} 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 a45ded73..38908032 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -84,6 +84,7 @@ mod bezpath; mod circle; pub mod common; mod cubicbez; +pub mod curve_intersections; mod ellipse; mod insets; mod line; @@ -106,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::*; diff --git a/src/line.rs b/src/line.rs index 29d73912..c6a398d7 100644 --- a/src/line.rs +++ b/src/line.rs @@ -49,6 +49,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.rs b/src/param_curve.rs index df30721c..2c4a7275 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; @@ -40,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. 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 6abf7cfd..873fe2b3 100644 --- a/src/point.rs +++ b/src/point.rs @@ -170,6 +170,42 @@ 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. + /// + /// 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)); + + // 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 { @@ -310,6 +346,54 @@ 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/quadbez.rs b/src/quadbez.rs index 41ccfde2..61ae5967 100644 --- a/src/quadbez.rs +++ b/src/quadbez.rs @@ -9,7 +9,7 @@ use crate::MAX_EXTREMA; use crate::{ Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, - Rect, Shape, + Rect, Shape, Vec2, }; /// A single quadratic Bézier segment. @@ -93,6 +93,17 @@ 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 + /// + /// 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; + 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 { @@ -104,6 +115,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. @@ -265,6 +282,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, this happens when all control points lie on top of each other + 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 { @@ -431,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));