geom.rs
1 //! Geometric primitives, like points and lines. 2 3 use arrayvec::ArrayVec; 4 use kurbo::{CubicBez, ParamCurve, ParamCurveExtrema}; 5 6 use crate::curve::{monic_quadratic_roots, solve_t_for_y, solve_x_for_y}; 7 use crate::num::CheapOrderedFloat; 8 9 /// A two-dimensional point. 10 /// 11 /// Points are sorted by `y` and then by `x`, for the convenience of our sweep-line 12 /// algorithm (which moves in increasing `y`). 13 #[cfg_attr(test, derive(serde::Serialize, serde::Deserialize))] 14 #[derive(Clone, Copy, PartialEq)] 15 pub struct Point { 16 /// Vertical coordinate. 17 /// 18 /// Although it isn't important for functionality, the documentation and method naming 19 /// assumes that larger values are down. 20 pub y: f64, 21 /// Horizontal component. 22 /// 23 /// Although it isn't important for functionality, the documentation and method naming 24 /// assumes that larger values are to the right. 25 pub x: f64, 26 } 27 28 impl Ord for Point { 29 fn cmp(&self, other: &Self) -> std::cmp::Ordering { 30 ( 31 CheapOrderedFloat::from(self.y), 32 CheapOrderedFloat::from(self.x), 33 ) 34 .cmp(&( 35 CheapOrderedFloat::from(other.y), 36 CheapOrderedFloat::from(other.x), 37 )) 38 } 39 } 40 41 impl PartialOrd for Point { 42 #[inline(always)] 43 fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> { 44 Some(self.cmp(other)) 45 } 46 } 47 48 impl Eq for Point {} 49 50 impl std::fmt::Debug for Point { 51 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { 52 write!(f, "({:?}, {:?})", self.x, self.y) 53 } 54 } 55 56 impl Point { 57 /// Create a new point. 58 /// 59 /// Note that the `x` coordinate comes first. This might be a tiny bit 60 /// confusing because we're sorting by `y` coordinate first, but `(x, y)` is 61 /// the only sane order (prove me wrong). 62 pub fn new(x: f64, y: f64) -> Self { 63 debug_assert!(x.is_finite()); 64 debug_assert!(y.is_finite()); 65 Point { x, y } 66 } 67 68 /// Compute an affine combination between `self` and `other`; that is, `(1 - t) * self + t * other`. 69 /// 70 /// Panics if a NaN is encountered, which might happen if some input is infinite. 71 pub fn affine(&self, other: &Self, t: f64) -> Self { 72 Point { 73 x: (1.0 - t) * self.x + t * other.x, 74 y: (1.0 - t) * self.y + t * other.y, 75 } 76 } 77 78 /// Compatibility with `kurbo` points. 79 pub fn to_kurbo(self) -> kurbo::Point { 80 kurbo::Point::new(self.x, self.y) 81 } 82 83 /// Compatibility with `kurbo` points. 84 pub fn from_kurbo(c: kurbo::Point) -> Self { 85 Self::new(c.x, c.y) 86 } 87 } 88 89 impl From<(f64, f64)> for Point { 90 fn from((x, y): (f64, f64)) -> Self { 91 Self { x, y } 92 } 93 } 94 95 impl From<Point> for kurbo::Point { 96 fn from(p: Point) -> Self { 97 kurbo::Point::new(p.x, p.y) 98 } 99 } 100 101 impl From<kurbo::Point> for Point { 102 fn from(p: kurbo::Point) -> Self { 103 Point::new(p.x, p.y) 104 } 105 } 106 107 /// A contour segment. 108 #[derive(Clone, PartialEq, Eq)] 109 pub struct Segment { 110 inner: SegmentInner, 111 } 112 113 #[derive(Clone, PartialEq, Eq)] 114 enum SegmentInner { 115 Line { 116 p0: Point, 117 p1: Point, 118 }, 119 Cubic { 120 p0: Point, 121 p1: Point, 122 p2: Point, 123 p3: Point, 124 }, 125 } 126 127 impl std::fmt::Debug for Segment { 128 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { 129 match self.inner { 130 SegmentInner::Line { p0, p1 } => write!(f, "{p0:?} -- {p1:?}"), 131 SegmentInner::Cubic { p0, p1, p2, p3 } => { 132 write!(f, "{p0:?} -- {p1:?} -- {p2:?} -- {p3:?}") 133 } 134 } 135 } 136 } 137 138 // Imagine that `cub` is basically a monotonic cubic, in that we didn't find any 139 // roots of its derivative. It could still have some starting and ending 140 // tangents that aren't quite right, so fix them up if necessary. 141 // 142 // For numerical reasons, the output of this may not actually be strictly 143 // monotonic. But it should be close. 144 fn force_monotonic(mut cub: CubicBez) -> Option<CubicBez> { 145 if cub.p0.y < cub.p3.y { 146 cub.p1.y = cub.p0.y.max(cub.p1.y); 147 cub.p2.y = cub.p3.y.min(cub.p2.y); 148 Some(cub) 149 } else if cub.p3.y < cub.p0.y { 150 cub.p1.y = cub.p0.y.min(cub.p1.y); 151 cub.p2.y = cub.p3.y.max(cub.p2.y); 152 Some(cub) 153 } else if cub.p0 != cub.p3 { 154 // It's a horizontal segment (or very close to one). Replace it with 155 // a true horizontal segment. 156 Some(CubicBez { 157 p0: cub.p0, 158 p1: cub.p0 + (cub.p3 - cub.p0) * (1.0 / 3.0), 159 p2: cub.p0 + (cub.p3 - cub.p0) * (2.0 / 3.0), 160 p3: cub.p3, 161 }) 162 } else { 163 None 164 } 165 } 166 167 pub(crate) fn monotonic_pieces(cub: CubicBez) -> ArrayVec<CubicBez, 3> { 168 let mut ret = ArrayVec::new(); 169 let q0 = cub.p1.y - cub.p0.y; 170 let q1 = cub.p2.y - cub.p1.y; 171 let q2 = cub.p3.y - cub.p2.y; 172 173 // Convert to the representation a t^2 + b t + c. 174 let c = q0; 175 let b = (q1 - q0) * 2.0; 176 let a = q0 - q1 * 2.0 + q2; 177 178 // Convert to the monic representation t^2 + b t + c. 179 let scaled_c = c * a.recip(); 180 let scaled_b = b * a.recip(); 181 182 let mut roots: ArrayVec<f64, 3> = ArrayVec::new(); 183 if scaled_c.is_infinite() || scaled_b.is_infinite() { 184 // A was small; treat it as a linear equation. 185 // We could use scaled_c here, but the originals should be more accurate. 186 roots.push(-c / b) 187 } else { 188 let (r0, r1) = monic_quadratic_roots(scaled_b, scaled_c); 189 roots.push(r0); 190 if r1 != r0 { 191 roots.push(r1); 192 } 193 } 194 195 let mut last_r = 0.0; 196 //dbg!(&roots); 197 // TODO: better handling for roots that are very close to 0.0 or 1.0 198 for r in roots { 199 if r > 0.0 && r < 1.0 { 200 let piece_before = cub.subsegment(last_r..r); 201 if let Some(c) = force_monotonic(piece_before) { 202 ret.push(c) 203 } 204 last_r = r; 205 } 206 } 207 208 let piece_before = cub.subsegment(last_r..1.0); 209 if let Some(c) = force_monotonic(piece_before) { 210 ret.push(c) 211 } 212 213 ret 214 } 215 216 impl Segment { 217 /// Create a new cubic segment that must be increasing in `y`. 218 pub fn monotonic_cubic(p0: Point, p1: Point, p2: Point, p3: Point) -> Self { 219 if p3.y == p0.y { 220 // Ensure that horizontal segments are just represented by straight lines. 221 // TODO: maybe we should do something about degenerate S-shaped curves? 222 Self::straight(p0, p3) 223 } else { 224 Self { 225 inner: SegmentInner::Cubic { p0, p1, p2, p3 }, 226 } 227 } 228 } 229 230 /// Create a new segment that's just a straight line. 231 /// 232 /// `start` must be less than `end`. 233 pub fn straight(start: Point, end: Point) -> Self { 234 debug_assert!(start <= end); 235 Self { 236 inner: SegmentInner::Line { p0: start, p1: end }, 237 } 238 } 239 240 /// The starting point (smallest in the sweep-line order) of this segment. 241 pub fn start(&self) -> Point { 242 match self.inner { 243 SegmentInner::Line { p0, .. } | SegmentInner::Cubic { p0, .. } => p0, 244 } 245 } 246 247 /// The ending point (largest in the sweep-line order) of this segment. 248 pub fn end(&self) -> Point { 249 match self.inner { 250 SegmentInner::Line { p1, .. } => p1, 251 SegmentInner::Cubic { p3, .. } => p3, 252 } 253 } 254 255 /// Compatibility with `kurbo`'s cubics. 256 pub fn to_kurbo_cubic(&self) -> kurbo::CubicBez { 257 match self.inner { 258 SegmentInner::Line { p0, p1 } => { 259 let p0 = p0.to_kurbo(); 260 let p3 = p1.to_kurbo(); 261 let p1 = p0 + (1. / 3.) * (p3 - p0); 262 let p2 = p0 + (2. / 3.) * (p3 - p0); 263 kurbo::CubicBez { p0, p1, p2, p3 } 264 } 265 SegmentInner::Cubic { p0, p1, p2, p3 } => kurbo::CubicBez { 266 p0: p0.to_kurbo(), 267 p1: p1.to_kurbo(), 268 p2: p2.to_kurbo(), 269 p3: p3.to_kurbo(), 270 }, 271 } 272 } 273 274 /// TODO: maybe From instead? 275 pub fn from_kurbo_cubic(c: kurbo::CubicBez) -> Self { 276 Self::monotonic_cubic( 277 Point::from_kurbo(c.p0), 278 Point::from_kurbo(c.p1), 279 Point::from_kurbo(c.p2), 280 Point::from_kurbo(c.p3), 281 ) 282 } 283 284 /// Is this segment just a straight line? 285 pub fn is_line(&self) -> bool { 286 matches!(self.inner, SegmentInner::Line { .. }) 287 } 288 289 /// A crude lower bound on our minimum horizontal position. 290 pub fn min_x(&self) -> f64 { 291 match self.inner { 292 SegmentInner::Line { p0, p1 } => p0.x.min(p1.x), 293 SegmentInner::Cubic { p0, p1, p2, p3 } => p0.x.min(p1.x).min(p2.x).min(p3.x), 294 } 295 } 296 297 /// A crude upper bound on our maximum horizontal position. 298 pub fn max_x(&self) -> f64 { 299 match self.inner { 300 SegmentInner::Line { p0, p1 } => p0.x.max(p1.x), 301 SegmentInner::Cubic { p0, p1, p2, p3 } => p0.x.max(p1.x).max(p2.x).max(p3.x), 302 } 303 } 304 305 /// Our `x` coordinate at the given `y` coordinate. 306 /// 307 /// Horizontal segments will return their largest `x` coordinate. 308 /// 309 /// # Panics 310 /// 311 /// Panics if `y` is outside the `y` range of this segment. 312 pub fn at_y(&self, y: f64) -> f64 { 313 debug_assert!( 314 (self.start().y..=self.end().y).contains(&y), 315 "segment {self:?}, y={y:?}" 316 ); 317 318 match self.inner { 319 SegmentInner::Line { p0, p1 } => { 320 if self.is_horizontal() { 321 p1.x 322 } else { 323 let t = (y - p0.y) / (p1.y - p0.y); 324 p0.x + t * (p1.x - p0.x) 325 } 326 } 327 SegmentInner::Cubic { .. } => solve_x_for_y(self.to_kurbo_cubic(), y), 328 } 329 } 330 331 fn local_bbox(&self, y: f64, eps: f64) -> kurbo::Rect { 332 let start_y = (y - eps).max(self.start().y); 333 let end_y = (y + eps).min(self.end().y); 334 335 match self.inner { 336 SegmentInner::Line { .. } => { 337 let start_x = self.at_y(start_y); 338 let end_x = self.at_y(end_y); 339 kurbo::Rect::from_points((start_x, start_y), (end_x, end_y)) 340 } 341 SegmentInner::Cubic { .. } => { 342 let c = self.to_kurbo_cubic(); 343 let t_min = solve_t_for_y(c, start_y); 344 let t_max = solve_t_for_y(c, end_y); 345 346 c.subsegment(t_min..t_max).bounding_box() 347 } 348 } 349 } 350 351 /// Returns a lower bound on the `x` coordinate near `y`. 352 /// 353 /// More precisely, take a Minkowski sum between this curve 354 /// and a small square. This returns a lower bound on the `x` 355 /// coordinate of the resulting set at height `y`. 356 pub fn lower(&self, y: f64, eps: f64) -> f64 { 357 self.local_bbox(y, eps).min_x() - eps 358 } 359 360 /// Returns an upper bound on the `x` coordinate near `y`. 361 /// 362 /// More precisely, take a Minkowski sum between this curve 363 /// and a small square. This returns a upper bound on the `x` 364 /// coordinate of the resulting set at height `y`. 365 pub fn upper(&self, y: f64, eps: f64) -> f64 { 366 self.local_bbox(y, eps).max_x() + eps 367 } 368 369 /// Returns true if this segment is exactly horizontal. 370 pub fn is_horizontal(&self) -> bool { 371 match self.inner { 372 SegmentInner::Line { p0, p1 } => p0.y == p1.y, 373 SegmentInner::Cubic { p0, p3, .. } => { 374 debug_assert_ne!( 375 p0.y, p3.y, 376 "horizontal segments should be represented by lines" 377 ); 378 false 379 } 380 } 381 } 382 } 383 384 #[cfg(test)] 385 pub(crate) mod tests { 386 use super::*; 387 use crate::num::tests::Reasonable; 388 use proptest::prelude::*; 389 390 impl Reasonable for Point { 391 type Strategy = BoxedStrategy<Point>; 392 393 fn reasonable() -> Self::Strategy { 394 (f64::reasonable(), f64::reasonable()) 395 .prop_map(|(x, y)| Point::new(x, y)) 396 .boxed() 397 } 398 } 399 400 // impl Reasonable for Segment { 401 // type Strategy = BoxedStrategy<Segment>; 402 403 // fn reasonable() -> Self::Strategy { 404 // (Point::reasonable(), Point::reasonable()) 405 // .prop_map(|(start, end)| { 406 // if start <= end { 407 // Segment::new(start, end) 408 // } else { 409 // Segment::new(end, start) 410 // } 411 // }) 412 // .boxed() 413 // } 414 // } 415 416 // fn segment_and_y() -> BoxedStrategy<(Segment, f64)> { 417 // Segment::reasonable() 418 // .prop_flat_map(|s| { 419 // let y0 = s.start.y; 420 // let y1 = s.end.y; 421 422 // // proptest's float sampler doesn't like a range like x..=x 423 // // https://github.com/proptest-rs/proptest/issues/343 424 // if y0 < y1 { 425 // (Just(s), (y0..=y1).boxed()) 426 // } else { 427 // (Just(s), Just(y0).boxed()) 428 // } 429 // }) 430 // .boxed() 431 // } 432 }