transverse.rs
1 //! Conservative checks for transversal intersections. 2 //! 3 //! We go to some effort to output curves with the correct horizontal 4 //! order. When curves are close to one another, our general solution 5 //! involves approximation with axis-aligned quadratics. This works, 6 //! but can lead to lots of subdivision. This module is about 7 //! improving a common case: when two curves intersect each other 8 //! transversally, they are close to one another near the intersection 9 //! point, but we can still tell that they are ordered. Hence, we 10 //! can avoid subdividing. 11 12 use kurbo::{CubicBez, ParamCurve, ParamCurveDeriv, QuadBez}; 13 14 use crate::curve::solve_t_for_y; 15 16 // Finds the minimum of `a t^2 + 2 b t (1 - t) + c (1 - t)^2` in the given range of t. 17 fn quad_max_between(a: f64, b: f64, c: f64, range: std::ops::RangeInclusive<f64>) -> f64 { 18 // We don't check for inf or nan here, but our subsequent comparisons will filter 19 // those out. (The order of the comparisons is important, because nan always compares 20 // false.) 21 let critical_t = (c - b) / (a - b + (c - b)); 22 23 let t0 = *range.start(); 24 let t1 = *range.end(); 25 26 let eval = |t: f64| a * t * t + 2.0 * b * t * (1.0 - t) + c * (1.0 - t) * (1.0 - t); 27 let mut max = eval(t0).max(eval(t1)); 28 29 if critical_t > t0 && critical_t < t1 { 30 max = max.max(eval(critical_t)); 31 } 32 max 33 } 34 35 fn quad_min_between(a: f64, b: f64, c: f64, range: std::ops::RangeInclusive<f64>) -> f64 { 36 -quad_max_between(-a, -b, -c, range) 37 } 38 39 pub fn transversal_after(left: CubicBez, right: CubicBez, y1: f64) -> bool { 40 debug_assert_eq!(left.start().y, right.start().y); 41 debug_assert!(left.start().x <= right.start().x); 42 debug_assert!(left.start().y < y1); 43 44 // First, check if they're transversal at the initial point. We exclude 45 // curves whose tangents are close to degenerate here, just to make things 46 // nicer numerically. 47 let left_tangent = left.p1 - left.p0; 48 let right_tangent = right.p1 - right.p0; 49 if left_tangent.hypot2() < 1e-20 || right_tangent.hypot2() < 1e-20 { 50 return false; 51 } 52 let cross = left_tangent.cross(right_tangent); 53 if left_tangent.dot(right_tangent) > 0.0 54 && cross * cross < 1e-2 * left_tangent.hypot2() * right_tangent.hypot2() 55 { 56 return false; 57 } 58 59 let t1_left = solve_t_for_y(left, y1); 60 let QuadBez { 61 p0: l0, 62 p1: l1, 63 p2: l2, 64 } = left.deriv(); 65 let left_max_dx_dt = quad_max_between(l2.x, l1.x, l0.x, 0.0..=t1_left); 66 67 let t1_right = solve_t_for_y(right, y1); 68 let QuadBez { 69 p0: r0, 70 p1: r1, 71 p2: r2, 72 } = right.deriv(); 73 let right_min_dx_dt = quad_min_between(r2.x, r1.x, r0.x, 0.0..=t1_right); 74 75 if left_max_dx_dt < 0.0 && right_min_dx_dt > 0.0 { 76 return true; 77 } 78 79 let left_dy_dt = if left_max_dx_dt < 0.0 { 80 quad_max_between(l2.y, l1.y, l0.y, 0.0..=t1_left) 81 } else { 82 quad_min_between(l2.y, l1.y, l0.y, 0.0..=t1_left) 83 }; 84 85 let right_dy_dt = if right_min_dx_dt < 0.0 { 86 quad_min_between(r2.y, r1.y, r0.y, 0.0..=t1_right) 87 } else { 88 quad_max_between(r2.y, r1.y, r0.y, 0.0..=t1_right) 89 }; 90 91 left_max_dx_dt * right_dy_dt < right_min_dx_dt * left_dy_dt 92 } 93 94 pub fn transversal_before(left: CubicBez, right: CubicBez, y0: f64) -> bool { 95 // Reflect everything in y. 96 let p = |q: kurbo::Point| kurbo::Point { x: q.x, y: -q.y }; 97 let c = |d: CubicBez| CubicBez { 98 p0: p(d.p3), 99 p1: p(d.p2), 100 p2: p(d.p1), 101 p3: p(d.p0), 102 }; 103 transversal_after(c(left), c(right), -y0) 104 } 105 106 #[cfg(test)] 107 mod tests { 108 use super::*; 109 110 #[test] 111 fn basic_transversal() { 112 let left = CubicBez { 113 p0: (0., 0.).into(), 114 p1: (-1., 1.).into(), 115 p2: (-1., 1.).into(), 116 p3: (0., 2.).into(), 117 }; 118 let right = CubicBez { 119 p0: (0., 0.).into(), 120 p1: (1., 1.).into(), 121 p2: (1., 1.).into(), 122 p3: (0., 2.).into(), 123 }; 124 125 assert!(transversal_after(left, right, 0.5)); 126 127 // It isn't really necessary for the implementation to pass 128 // these two tests; they just validate that our particular 129 // implementation does what we think. Our implementation looks 130 // only at derivatives, and the derivatives change direction 131 // at y = 1.0. 132 assert!(transversal_after(left, right, 0.9)); 133 assert!(!transversal_after(left, right, 1.1)); 134 } 135 136 #[test] 137 fn parallel() { 138 let left = CubicBez { 139 p0: (0., 0.).into(), 140 p1: (-1., 1.).into(), 141 p2: (-1., 1.).into(), 142 p3: (0., 2.).into(), 143 }; 144 let right = CubicBez { 145 p0: (0., 0.).into(), 146 p1: (-1., 1.).into(), 147 p2: (1., 1.).into(), 148 p3: (0., 2.).into(), 149 }; 150 151 assert!(!transversal_after(left, right, 0.01)); 152 } 153 154 #[test] 155 fn parallel_and_horizontal() { 156 let left = CubicBez { 157 p0: (0., 0.).into(), 158 p1: (-1., 0.).into(), 159 p2: (-1., 1.).into(), 160 p3: (0., 2.).into(), 161 }; 162 let right = CubicBez { 163 p0: (0., 0.).into(), 164 p1: (-1., 0.).into(), 165 p2: (1., 1.).into(), 166 p3: (0., 2.).into(), 167 }; 168 169 assert!(!transversal_after(left, right, 0.01)); 170 } 171 172 #[test] 173 fn antiparallel_and_horizontal() { 174 let left = CubicBez { 175 p0: (0., 0.).into(), 176 p1: (-1., 0.).into(), 177 p2: (-1., 1.).into(), 178 p3: (0., 2.).into(), 179 }; 180 let right = CubicBez { 181 p0: (0., 0.).into(), 182 p1: (1., 0.).into(), 183 p2: (1., 1.).into(), 184 p3: (0., 2.).into(), 185 }; 186 187 assert!(transversal_after(left, right, 0.01)); 188 } 189 }