1// Copyright 2018 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! Cubic Bézier segments.
5
6use alloc::vec;
7use alloc::vec::Vec;
8use core::ops::{Mul, Range};
9
10use crate::MAX_EXTREMA;
11use crate::{Line, QuadSpline, Vec2};
12use arrayvec::ArrayVec;
13
14use crate::common::{
15 solve_quadratic, GAUSS_LEGENDRE_COEFFS_16_HALF, GAUSS_LEGENDRE_COEFFS_24_HALF,
16 GAUSS_LEGENDRE_COEFFS_8, GAUSS_LEGENDRE_COEFFS_8_HALF,
17};
18use crate::{
19 Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature,
20 ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, QuadBez, Rect, Shape,
21};
22
23#[cfg(not(feature = "std"))]
24use crate::common::FloatFuncs;
25
26const MAX_SPLINE_SPLIT: usize = 100;
27
28/// A single cubic Bézier segment.
29#[derive(Clone, Copy, Debug, PartialEq)]
30#[cfg_attr(feature = "schemars", derive(schemars::JsonSchema))]
31#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
32#[allow(missing_docs)]
33pub struct CubicBez {
34 pub p0: Point,
35 pub p1: Point,
36 pub p2: Point,
37 pub p3: Point,
38}
39
40/// An iterator which produces quadratic Bézier segments.
41struct ToQuads {
42 c: CubicBez,
43 i: usize,
44 n: usize,
45}
46
47impl CubicBez {
48 /// Create a new cubic Bézier segment.
49 #[inline]
50 pub fn new<P: Into<Point>>(p0: P, p1: P, p2: P, p3: P) -> CubicBez {
51 CubicBez {
52 p0: p0.into(),
53 p1: p1.into(),
54 p2: p2.into(),
55 p3: p3.into(),
56 }
57 }
58
59 /// Convert to quadratic Béziers.
60 ///
61 /// The iterator returns the start and end parameter in the cubic of each quadratic
62 /// segment, along with the quadratic.
63 ///
64 /// Note that the resulting quadratic Béziers are not in general G1 continuous;
65 /// they are optimized for minimizing distance error.
66 ///
67 /// This iterator will always produce at least one `QuadBez`.
68 #[inline]
69 pub fn to_quads(&self, accuracy: f64) -> impl Iterator<Item = (f64, f64, QuadBez)> {
70 // The maximum error, as a vector from the cubic to the best approximating
71 // quadratic, is proportional to the third derivative, which is constant
72 // across the segment. Thus, the error scales down as the third power of
73 // the number of subdivisions. Our strategy then is to subdivide `t` evenly.
74 //
75 // This is an overestimate of the error because only the component
76 // perpendicular to the first derivative is important. But the simplicity is
77 // appealing.
78
79 // This magic number is the square of 36 / sqrt(3).
80 // See: http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
81 let max_hypot2 = 432.0 * accuracy * accuracy;
82 let p1x2 = 3.0 * self.p1.to_vec2() - self.p0.to_vec2();
83 let p2x2 = 3.0 * self.p2.to_vec2() - self.p3.to_vec2();
84 let err = (p2x2 - p1x2).hypot2();
85 let n = ((err / max_hypot2).powf(1. / 6.0).ceil() as usize).max(1);
86
87 ToQuads { c: *self, n, i: 0 }
88 }
89
90 /// Return a [`QuadSpline`] approximating this cubic Bézier.
91 ///
92 /// Returns `None` if no suitable approximation is found within the given
93 /// tolerance.
94 pub fn approx_spline(&self, accuracy: f64) -> Option<QuadSpline> {
95 (1..=MAX_SPLINE_SPLIT).find_map(|n| self.approx_spline_n(n, accuracy))
96 }
97
98 // Approximate a cubic curve with a quadratic spline of `n` curves
99 fn approx_spline_n(&self, n: usize, accuracy: f64) -> Option<QuadSpline> {
100 if n == 1 {
101 return self
102 .try_approx_quadratic(accuracy)
103 .map(|quad| QuadSpline::new(vec![quad.p0, quad.p1, quad.p2]));
104 }
105 let mut cubics = self.split_into_n(n);
106
107 // The above function guarantees that the iterator returns n items,
108 // which is why we're unwrapping things with wild abandon.
109 let mut next_cubic = cubics.next().unwrap();
110 let mut next_q1: Point = next_cubic.approx_quad_control(0.0);
111 let mut q2 = self.p0;
112 let mut d1 = Vec2::ZERO;
113 let mut spline = vec![self.p0, next_q1];
114 for i in 1..=n {
115 let current_cubic: CubicBez = next_cubic;
116 let q0 = q2;
117 let q1 = next_q1;
118 q2 = if i < n {
119 next_cubic = cubics.next().unwrap();
120 next_q1 = next_cubic.approx_quad_control(i as f64 / (n - 1) as f64);
121
122 spline.push(next_q1);
123 q1.midpoint(next_q1)
124 } else {
125 current_cubic.p3
126 };
127 let d0 = d1;
128 d1 = q2.to_vec2() - current_cubic.p3.to_vec2();
129
130 if d1.hypot() > accuracy
131 || !CubicBez::new(
132 d0.to_point(),
133 q0.lerp(q1, 2.0 / 3.0) - current_cubic.p1.to_vec2(),
134 q2.lerp(q1, 2.0 / 3.0) - current_cubic.p2.to_vec2(),
135 d1.to_point(),
136 )
137 .fit_inside(accuracy)
138 {
139 return None;
140 }
141 }
142 spline.push(self.p3);
143 Some(QuadSpline::new(spline))
144 }
145
146 fn approx_quad_control(&self, t: f64) -> Point {
147 let p1 = self.p0 + (self.p1 - self.p0) * 1.5;
148 let p2 = self.p3 + (self.p2 - self.p3) * 1.5;
149 p1.lerp(p2, t)
150 }
151
152 /// Approximate a cubic with a single quadratic
153 ///
154 /// Returns a quadratic approximating the given cubic that maintains
155 /// endpoint tangents if that is within tolerance, or None otherwise.
156 fn try_approx_quadratic(&self, accuracy: f64) -> Option<QuadBez> {
157 if let Some(q1) = Line::new(self.p0, self.p1).crossing_point(Line::new(self.p2, self.p3)) {
158 let c1 = self.p0.lerp(q1, 2.0 / 3.0);
159 let c2 = self.p3.lerp(q1, 2.0 / 3.0);
160 if !CubicBez::new(
161 Point::ZERO,
162 c1 - self.p1.to_vec2(),
163 c2 - self.p2.to_vec2(),
164 Point::ZERO,
165 )
166 .fit_inside(accuracy)
167 {
168 return None;
169 }
170 return Some(QuadBez::new(self.p0, q1, self.p3));
171 }
172 None
173 }
174
175 fn split_into_n(&self, n: usize) -> impl Iterator<Item = CubicBez> {
176 // for certain small values of `n` we precompute all our values.
177 // if we have six or fewer items we precompute them.
178 let mut storage = ArrayVec::<_, 6>::new();
179
180 match n {
181 1 => storage.push(*self),
182 2 => {
183 let (l, r) = self.subdivide();
184 storage.try_extend_from_slice(&[r, l]).unwrap();
185 }
186 3 => {
187 let (left, mid, right) = self.subdivide_3();
188 storage.try_extend_from_slice(&[right, mid, left]).unwrap();
189 }
190 4 => {
191 let (l, r) = self.subdivide();
192 let (ll, lr) = l.subdivide();
193 let (rl, rr) = r.subdivide();
194 storage.try_extend_from_slice(&[rr, rl, lr, ll]).unwrap();
195 }
196 6 => {
197 let (l, r) = self.subdivide();
198 let (l1, l2, l3) = l.subdivide_3();
199 let (r1, r2, r3) = r.subdivide_3();
200 storage
201 .try_extend_from_slice(&[r3, r2, r1, l3, l2, l1])
202 .unwrap();
203 }
204 _ => (),
205 }
206
207 // a limitation of returning 'impl Trait' is that the implementation
208 // can only return a single concrete type; that is you cannot return
209 // Vec::into_iter() from one branch, and then HashSet::into_iter from
210 // another branch.
211 //
212 // This means we have to get a bit fancy, and have a single concrete
213 // type that represents both of our possible cases.
214
215 let mut storage = if storage.is_empty() {
216 None
217 } else {
218 Some(storage)
219 };
220
221 // used in the fallback case
222 let mut i = 0;
223 let (a, b, c, d) = self.parameters();
224 let dt = 1.0 / n as f64;
225 let delta_2 = dt * dt;
226 let delta_3 = dt * delta_2;
227
228 core::iter::from_fn(move || {
229 // if storage exists, we use it exclusively
230 if let Some(storage) = storage.as_mut() {
231 return storage.pop();
232 }
233
234 // if storage does not exist, we are exclusively working down here.
235 if i >= n {
236 return None;
237 }
238
239 let t1 = i as f64 * dt;
240 let t1_2 = t1 * t1;
241 let a1 = a * delta_3;
242 let b1 = (3.0 * a * t1 + b) * delta_2;
243 let c1 = (2.0 * b * t1 + c + 3.0 * a * t1_2) * dt;
244 let d1 = a * t1 * t1_2 + b * t1_2 + c * t1 + d;
245 let result = CubicBez::from_parameters(a1, b1, c1, d1);
246 i += 1;
247 Some(result)
248 })
249 }
250
251 fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) {
252 let c = (self.p1 - self.p0) * 3.0;
253 let b = (self.p2 - self.p1) * 3.0 - c;
254 let d = self.p0.to_vec2();
255 let a = self.p3.to_vec2() - d - c - b;
256 (a, b, c, d)
257 }
258
259 /// Rust port of cu2qu [calc_cubic_points](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L63-L68)
260 fn from_parameters(a: Vec2, b: Vec2, c: Vec2, d: Vec2) -> Self {
261 let p0 = d.to_point();
262 let p1 = c.div_exact(3.0).to_point() + d;
263 let p2 = (b + c).div_exact(3.0).to_point() + p1.to_vec2();
264 let p3 = (a + d + c + b).to_point();
265 CubicBez::new(p0, p1, p2, p3)
266 }
267
268 fn subdivide_3(&self) -> (CubicBez, CubicBez, CubicBez) {
269 let (p0, p1, p2, p3) = (
270 self.p0.to_vec2(),
271 self.p1.to_vec2(),
272 self.p2.to_vec2(),
273 self.p3.to_vec2(),
274 );
275 // The original Python cu2qu code here does not use division operator to divide by 27 but
276 // instead uses multiplication by the reciprocal 1 / 27. We want to match it exactly
277 // to avoid any floating point differences, hence in this particular case we do not use div_exact.
278 // I could directly use the Vec2 Div trait (also implemented as multiplication by reciprocal)
279 // but I prefer to be explicit here.
280 // Source: https://github.com/fonttools/fonttools/blob/85c80be/Lib/fontTools/cu2qu/cu2qu.py#L215-L218
281 // See also: https://github.com/linebender/kurbo/issues/272
282 let one_27th = 27.0_f64.recip();
283 let mid1 = ((8.0 * p0 + 12.0 * p1 + 6.0 * p2 + p3) * one_27th).to_point();
284 let deriv1 = (p3 + 3.0 * p2 - 4.0 * p0) * one_27th;
285 let mid2 = ((p0 + 6.0 * p1 + 12.0 * p2 + 8.0 * p3) * one_27th).to_point();
286 let deriv2 = (4.0 * p3 - 3.0 * p1 - p0) * one_27th;
287
288 let left = CubicBez::new(
289 self.p0,
290 (2.0 * p0 + p1).div_exact(3.0).to_point(),
291 mid1 - deriv1,
292 mid1,
293 );
294 let mid = CubicBez::new(mid1, mid1 + deriv1, mid2 - deriv2, mid2);
295 let right = CubicBez::new(
296 mid2,
297 mid2 + deriv2,
298 (p2 + 2.0 * p3).div_exact(3.0).to_point(),
299 self.p3,
300 );
301 (left, mid, right)
302 }
303
304 /// Does this curve fit inside the given distance from the origin?
305 ///
306 /// Rust port of cu2qu [cubic_farthest_fit_inside](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L281)
307 fn fit_inside(&self, distance: f64) -> bool {
308 if self.p2.to_vec2().hypot() <= distance && self.p1.to_vec2().hypot() <= distance {
309 return true;
310 }
311 let mid =
312 (self.p0.to_vec2() + 3.0 * (self.p1.to_vec2() + self.p2.to_vec2()) + self.p3.to_vec2())
313 * 0.125;
314 if mid.hypot() > distance {
315 return false;
316 }
317 // Split in two. Note that cu2qu here uses a 3/8 subdivision. I don't know why.
318 let (left, right) = self.subdivide();
319 left.fit_inside(distance) && right.fit_inside(distance)
320 }
321
322 /// Is this cubic Bezier curve finite?
323 #[inline]
324 pub fn is_finite(&self) -> bool {
325 self.p0.is_finite() && self.p1.is_finite() && self.p2.is_finite() && self.p3.is_finite()
326 }
327
328 /// Is this cubic Bezier curve NaN?
329 #[inline]
330 pub fn is_nan(&self) -> bool {
331 self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan() || self.p3.is_nan()
332 }
333
334 /// Determine the inflection points.
335 ///
336 /// Return value is t parameter for the inflection points of the curve segment.
337 /// There are a maximum of two for a cubic Bézier.
338 ///
339 /// See <https://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html>
340 /// for the theory.
341 pub fn inflections(&self) -> ArrayVec<f64, 2> {
342 let a = self.p1 - self.p0;
343 let b = (self.p2 - self.p1) - a;
344 let c = (self.p3 - self.p0) - 3. * (self.p2 - self.p1);
345 solve_quadratic(a.cross(b), a.cross(c), b.cross(c))
346 .iter()
347 .copied()
348 .filter(|t| *t >= 0.0 && *t <= 1.0)
349 .collect()
350 }
351}
352
353/// An iterator for cubic beziers.
354pub struct CubicBezIter {
355 cubic: CubicBez,
356 ix: usize,
357}
358
359impl Shape for CubicBez {
360 type PathElementsIter<'iter> = CubicBezIter;
361
362 #[inline]
363 fn path_elements(&self, _tolerance: f64) -> CubicBezIter {
364 CubicBezIter {
365 cubic: *self,
366 ix: 0,
367 }
368 }
369
370 fn area(&self) -> f64 {
371 0.0
372 }
373
374 #[inline]
375 fn perimeter(&self, accuracy: f64) -> f64 {
376 self.arclen(accuracy)
377 }
378
379 fn winding(&self, _pt: Point) -> i32 {
380 0
381 }
382
383 #[inline]
384 fn bounding_box(&self) -> Rect {
385 ParamCurveExtrema::bounding_box(self)
386 }
387}
388
389impl Iterator for CubicBezIter {
390 type Item = PathEl;
391
392 fn next(&mut self) -> Option<PathEl> {
393 self.ix += 1;
394 match self.ix {
395 1 => Some(PathEl::MoveTo(self.cubic.p0)),
396 2 => Some(PathEl::CurveTo(self.cubic.p1, self.cubic.p2, self.cubic.p3)),
397 _ => None,
398 }
399 }
400}
401
402impl ParamCurve for CubicBez {
403 #[inline]
404 fn eval(&self, t: f64) -> Point {
405 let mt = 1.0 - t;
406 let v = self.p0.to_vec2() * (mt * mt * mt)
407 + (self.p1.to_vec2() * (mt * mt * 3.0)
408 + (self.p2.to_vec2() * (mt * 3.0) + self.p3.to_vec2() * t) * t)
409 * t;
410 v.to_point()
411 }
412
413 fn subsegment(&self, range: Range<f64>) -> CubicBez {
414 let (t0, t1) = (range.start, range.end);
415 let p0 = self.eval(t0);
416 let p3 = self.eval(t1);
417 let d = self.deriv();
418 let scale = (t1 - t0) * (1.0 / 3.0);
419 let p1 = p0 + scale * d.eval(t0).to_vec2();
420 let p2 = p3 - scale * d.eval(t1).to_vec2();
421 CubicBez { p0, p1, p2, p3 }
422 }
423
424 /// Subdivide into halves, using de Casteljau.
425 #[inline]
426 fn subdivide(&self) -> (CubicBez, CubicBez) {
427 let pm = self.eval(0.5);
428 (
429 CubicBez::new(
430 self.p0,
431 self.p0.midpoint(self.p1),
432 ((self.p0.to_vec2() + self.p1.to_vec2() * 2.0 + self.p2.to_vec2()) * 0.25)
433 .to_point(),
434 pm,
435 ),
436 CubicBez::new(
437 pm,
438 ((self.p1.to_vec2() + self.p2.to_vec2() * 2.0 + self.p3.to_vec2()) * 0.25)
439 .to_point(),
440 self.p2.midpoint(self.p3),
441 self.p3,
442 ),
443 )
444 }
445
446 #[inline]
447 fn start(&self) -> Point {
448 self.p0
449 }
450
451 #[inline]
452 fn end(&self) -> Point {
453 self.p3
454 }
455}
456
457impl ParamCurveDeriv for CubicBez {
458 type DerivResult = QuadBez;
459
460 #[inline]
461 fn deriv(&self) -> QuadBez {
462 QuadBez::new(
463 (3.0 * (self.p1 - self.p0)).to_point(),
464 (3.0 * (self.p2 - self.p1)).to_point(),
465 (3.0 * (self.p3 - self.p2)).to_point(),
466 )
467 }
468}
469
470fn arclen_quadrature_core(coeffs: &[(f64, f64)], dm: Vec2, dm1: Vec2, dm2: Vec2) -> f64 {
471 coeffsimpl Iterator
472 .iter()
473 .map(|&(wi: f64, xi: f64)| {
474 let d: Vec2 = dm + dm2 * (xi * xi);
475 let dpx: f64 = (d + dm1 * xi).hypot();
476 let dmx: f64 = (d - dm1 * xi).hypot();
477 (2.25f64.sqrt() * wi) * (dpx + dmx)
478 })
479 .sum::<f64>()
480}
481
482fn arclen_rec(c: &CubicBez, accuracy: f64, depth: usize) -> f64 {
483 let d03 = c.p3 - c.p0;
484 let d01 = c.p1 - c.p0;
485 let d12 = c.p2 - c.p1;
486 let d23 = c.p3 - c.p2;
487 let lp_lc = d01.hypot() + d12.hypot() + d23.hypot() - d03.hypot();
488 let dd1 = d12 - d01;
489 let dd2 = d23 - d12;
490 // It might be faster to do direct multiplies, the data dependencies would be shorter.
491 // The following values don't have the factor of 3 for first deriv
492 let dm = 0.25 * (d01 + d23) + 0.5 * d12; // first derivative at midpoint
493 let dm1 = 0.5 * (dd2 + dd1); // second derivative at midpoint
494 let dm2 = 0.25 * (dd2 - dd1); // 0.5 * (third derivative at midpoint)
495
496 let est = GAUSS_LEGENDRE_COEFFS_8
497 .iter()
498 .map(|&(wi, xi)| {
499 wi * {
500 let d_norm2 = (dm + dm1 * xi + dm2 * (xi * xi)).hypot2();
501 let dd_norm2 = (dm1 + dm2 * (2.0 * xi)).hypot2();
502 dd_norm2 / d_norm2
503 }
504 })
505 .sum::<f64>();
506 let est_gauss8_error = (est.powi(3) * 2.5e-6).min(3e-2) * lp_lc;
507 if est_gauss8_error < accuracy {
508 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_8_HALF, dm, dm1, dm2);
509 }
510 let est_gauss16_error = (est.powi(6) * 1.5e-11).min(9e-3) * lp_lc;
511 if est_gauss16_error < accuracy {
512 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_16_HALF, dm, dm1, dm2);
513 }
514 let est_gauss24_error = (est.powi(9) * 3.5e-16).min(3.5e-3) * lp_lc;
515 if est_gauss24_error < accuracy || depth >= 20 {
516 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_24_HALF, dm, dm1, dm2);
517 }
518 let (c0, c1) = c.subdivide();
519 arclen_rec(&c0, accuracy * 0.5, depth + 1) + arclen_rec(&c1, accuracy * 0.5, depth + 1)
520}
521
522impl ParamCurveArclen for CubicBez {
523 /// Arclength of a cubic Bézier segment.
524 ///
525 /// This is an adaptive subdivision approach using Legendre-Gauss quadrature
526 /// in the base case, and an error estimate to decide when to subdivide.
527 fn arclen(&self, accuracy: f64) -> f64 {
528 arclen_rec(self, accuracy, depth:0)
529 }
530}
531
532impl ParamCurveArea for CubicBez {
533 #[inline]
534 fn signed_area(&self) -> f64 {
535 (self.p0.x * (6.0 * self.p1.y + 3.0 * self.p2.y + self.p3.y)
536 + 3.0
537 * (self.p1.x * (-2.0 * self.p0.y + self.p2.y + self.p3.y)
538 - self.p2.x * (self.p0.y + self.p1.y - 2.0 * self.p3.y))
539 - self.p3.x * (self.p0.y + 3.0 * self.p1.y + 6.0 * self.p2.y))
540 * (1.0 / 20.0)
541 }
542}
543
544impl ParamCurveNearest for CubicBez {
545 /// Find the nearest point, using subdivision.
546 fn nearest(&self, p: Point, accuracy: f64) -> Nearest {
547 let mut best_r: Option = None;
548 let mut best_t: f64 = 0.0;
549 for (t0: f64, t1: f64, q: QuadBez) in self.to_quads(accuracy) {
550 let nearest: Nearest = q.nearest(p, accuracy);
551 if best_r
552 .map(|best_r| nearest.distance_sq < best_r)
553 .unwrap_or(default:true)
554 {
555 best_t = t0 + nearest.t * (t1 - t0);
556 best_r = Some(nearest.distance_sq);
557 }
558 }
559 Nearest {
560 t: best_t,
561 distance_sq: best_r.unwrap(),
562 }
563 }
564}
565
566impl ParamCurveCurvature for CubicBez {}
567
568impl ParamCurveExtrema for CubicBez {
569 fn extrema(&self) -> ArrayVec<f64, MAX_EXTREMA> {
570 fn one_coord(result: &mut ArrayVec<f64, MAX_EXTREMA>, d0: f64, d1: f64, d2: f64) {
571 let a: f64 = d0 - 2.0 * d1 + d2;
572 let b: f64 = 2.0 * (d1 - d0);
573 let c: f64 = d0;
574 let roots: ArrayVec = solve_quadratic(c0:c, c1:b, c2:a);
575 for &t: f64 in &roots {
576 if t > 0.0 && t < 1.0 {
577 result.push(element:t);
578 }
579 }
580 }
581 let mut result: ArrayVec = ArrayVec::new();
582 let d0: Vec2 = self.p1 - self.p0;
583 let d1: Vec2 = self.p2 - self.p1;
584 let d2: Vec2 = self.p3 - self.p2;
585 one_coord(&mut result, d0:d0.x, d1:d1.x, d2:d2.x);
586 one_coord(&mut result, d0:d0.y, d1:d1.y, d2:d2.y);
587 result.sort_by(|a: &f64, b: &f64| a.partial_cmp(b).unwrap());
588 result
589 }
590}
591
592impl Mul<CubicBez> for Affine {
593 type Output = CubicBez;
594
595 #[inline]
596 fn mul(self, c: CubicBez) -> CubicBez {
597 CubicBez {
598 p0: self * c.p0,
599 p1: self * c.p1,
600 p2: self * c.p2,
601 p3: self * c.p3,
602 }
603 }
604}
605
606impl Iterator for ToQuads {
607 type Item = (f64, f64, QuadBez);
608
609 fn next(&mut self) -> Option<(f64, f64, QuadBez)> {
610 if self.i == self.n {
611 return None;
612 }
613 let t0: f64 = self.i as f64 / self.n as f64;
614 let t1: f64 = (self.i + 1) as f64 / self.n as f64;
615 let seg: CubicBez = self.c.subsegment(range:t0..t1);
616 let p1x2: Vec2 = 3.0 * seg.p1.to_vec2() - seg.p0.to_vec2();
617 let p2x2: Vec2 = 3.0 * seg.p2.to_vec2() - seg.p3.to_vec2();
618 let result: QuadBez = QuadBez::new(seg.p0, ((p1x2 + p2x2) / 4.0).to_point(), p2:seg.p3);
619 self.i += 1;
620 Some((t0, t1, result))
621 }
622
623 fn size_hint(&self) -> (usize, Option<usize>) {
624 let remaining: usize = self.n - self.i;
625 (remaining, Some(remaining))
626 }
627}
628
629/// Convert multiple cubic Bézier curves to quadratic splines.
630///
631/// Ensures that the resulting splines have the same number of control points.
632///
633/// Rust port of cu2qu [cubic_approx_quadratic](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L322)
634pub fn cubics_to_quadratic_splines(curves: &[CubicBez], accuracy: f64) -> Option<Vec<QuadSpline>> {
635 let mut result: Vec = Vec::new();
636 let mut split_order: usize = 0;
637
638 while split_order <= MAX_SPLINE_SPLIT {
639 split_order += 1;
640 result.clear();
641
642 for curve: &CubicBez in curves {
643 match curve.approx_spline_n(n:split_order, accuracy) {
644 Some(spline: QuadSpline) => result.push(spline),
645 None => break,
646 }
647 }
648
649 if result.len() == curves.len() {
650 return Some(result);
651 }
652 }
653 None
654}
655#[cfg(test)]
656mod tests {
657 use crate::{
658 cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen,
659 ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez,
660 QuadSpline,
661 };
662
663 #[test]
664 fn cubicbez_deriv() {
665 // y = x^2
666 let c = CubicBez::new(
667 (0.0, 0.0),
668 (1.0 / 3.0, 0.0),
669 (2.0 / 3.0, 1.0 / 3.0),
670 (1.0, 1.0),
671 );
672 let deriv = c.deriv();
673
674 let n = 10;
675 for i in 0..=n {
676 let t = (i as f64) * (n as f64).recip();
677 let delta = 1e-6;
678 let p = c.eval(t);
679 let p1 = c.eval(t + delta);
680 let d_approx = (p1 - p) * delta.recip();
681 let d = deriv.eval(t).to_vec2();
682 assert!((d - d_approx).hypot() < delta * 2.0);
683 }
684 }
685
686 #[test]
687 fn cubicbez_arclen() {
688 // y = x^2
689 let c = CubicBez::new(
690 (0.0, 0.0),
691 (1.0 / 3.0, 0.0),
692 (2.0 / 3.0, 1.0 / 3.0),
693 (1.0, 1.0),
694 );
695 let true_arclen = 0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln();
696 for i in 0..12 {
697 let accuracy = 0.1f64.powi(i);
698 let error = c.arclen(accuracy) - true_arclen;
699 assert!(error.abs() < accuracy);
700 }
701 }
702
703 #[test]
704 fn cubicbez_inv_arclen() {
705 // y = x^2 / 100
706 let c = CubicBez::new(
707 (0.0, 0.0),
708 (100.0 / 3.0, 0.0),
709 (200.0 / 3.0, 100.0 / 3.0),
710 (100.0, 100.0),
711 );
712 let true_arclen = 100.0 * (0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln());
713 for i in 0..12 {
714 let accuracy = 0.1f64.powi(i);
715 let n = 10;
716 for j in 0..=n {
717 let arc = (j as f64) * ((n as f64).recip() * true_arclen);
718 let t = c.inv_arclen(arc, accuracy * 0.5);
719 let actual_arc = c.subsegment(0.0..t).arclen(accuracy * 0.5);
720 assert!(
721 (arc - actual_arc).abs() < accuracy,
722 "at accuracy {accuracy:e}, wanted {actual_arc} got {arc}"
723 );
724 }
725 }
726 // corner case: user passes accuracy larger than total arc length
727 let accuracy = true_arclen * 1.1;
728 let arc = true_arclen * 0.5;
729 let t = c.inv_arclen(arc, accuracy);
730 let actual_arc = c.subsegment(0.0..t).arclen(accuracy);
731 assert!(
732 (arc - actual_arc).abs() < 2.0 * accuracy,
733 "at accuracy {accuracy:e}, want {actual_arc} got {arc}"
734 );
735 }
736
737 #[test]
738 fn cubicbez_inv_arclen_accuracy() {
739 let c = CubicBez::new((0.2, 0.73), (0.35, 1.08), (0.85, 1.08), (1.0, 0.73));
740 let true_t = c.inv_arclen(0.5, 1e-12);
741 for i in 1..12 {
742 let accuracy = (0.1f64).powi(i);
743 let approx_t = c.inv_arclen(0.5, accuracy);
744 assert!((approx_t - true_t).abs() <= accuracy);
745 }
746 }
747
748 #[test]
749 #[allow(clippy::float_cmp)]
750 fn cubicbez_signed_area_linear() {
751 // y = 1 - x
752 let c = CubicBez::new(
753 (1.0, 0.0),
754 (2.0 / 3.0, 1.0 / 3.0),
755 (1.0 / 3.0, 2.0 / 3.0),
756 (0.0, 1.0),
757 );
758 let epsilon = 1e-12;
759 assert_eq!((Affine::rotate(0.5) * c).signed_area(), 0.5);
760 assert!(((Affine::rotate(0.5) * c).signed_area() - 0.5).abs() < epsilon);
761 assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.0).abs() < epsilon);
762 assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.0).abs() < epsilon);
763 }
764
765 #[test]
766 fn cubicbez_signed_area() {
767 // y = 1 - x^3
768 let c = CubicBez::new((1.0, 0.0), (2.0 / 3.0, 1.0), (1.0 / 3.0, 1.0), (0.0, 1.0));
769 let epsilon = 1e-12;
770 assert!((c.signed_area() - 0.75).abs() < epsilon);
771 assert!(((Affine::rotate(0.5) * c).signed_area() - 0.75).abs() < epsilon);
772 assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.25).abs() < epsilon);
773 assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.25).abs() < epsilon);
774 }
775
776 #[test]
777 fn cubicbez_nearest() {
778 fn verify(result: Nearest, expected: f64) {
779 assert!(
780 (result.t - expected).abs() < 1e-6,
781 "got {result:?} expected {expected}"
782 );
783 }
784 // y = x^3
785 let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
786 verify(c.nearest((0.1, 0.001).into(), 1e-6), 0.1);
787 verify(c.nearest((0.2, 0.008).into(), 1e-6), 0.2);
788 verify(c.nearest((0.3, 0.027).into(), 1e-6), 0.3);
789 verify(c.nearest((0.4, 0.064).into(), 1e-6), 0.4);
790 verify(c.nearest((0.5, 0.125).into(), 1e-6), 0.5);
791 verify(c.nearest((0.6, 0.216).into(), 1e-6), 0.6);
792 verify(c.nearest((0.7, 0.343).into(), 1e-6), 0.7);
793 verify(c.nearest((0.8, 0.512).into(), 1e-6), 0.8);
794 verify(c.nearest((0.9, 0.729).into(), 1e-6), 0.9);
795 verify(c.nearest((1.0, 1.0).into(), 1e-6), 1.0);
796 verify(c.nearest((1.1, 1.1).into(), 1e-6), 1.0);
797 verify(c.nearest((-0.1, 0.0).into(), 1e-6), 0.0);
798 let a = Affine::rotate(0.5);
799 verify((a * c).nearest(a * Point::new(0.1, 0.001), 1e-6), 0.1);
800 }
801
802 // ensure to_quads returns something given colinear points
803 #[test]
804 fn degenerate_to_quads() {
805 let c = CubicBez::new((0., 9.), (6., 6.), (12., 3.0), (18., 0.0));
806 let quads = c.to_quads(1e-6).collect::<Vec<_>>();
807 assert_eq!(quads.len(), 1, "{:?}", &quads);
808 }
809
810 #[test]
811 fn cubicbez_extrema() {
812 // y = x^2
813 let q = CubicBez::new((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0));
814 let extrema = q.extrema();
815 assert_eq!(extrema.len(), 1);
816 assert!((extrema[0] - 0.5).abs() < 1e-6);
817
818 let q = CubicBez::new((0.4, 0.5), (0.0, 1.0), (1.0, 0.0), (0.5, 0.4));
819 let extrema = q.extrema();
820 assert_eq!(extrema.len(), 4);
821 }
822
823 #[test]
824 fn cubicbez_toquads() {
825 // y = x^3
826 let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
827 for i in 0..10 {
828 let accuracy = 0.1f64.powi(i);
829 let mut worst: f64 = 0.0;
830 for (_count, (t0, t1, q)) in c.to_quads(accuracy).enumerate() {
831 let epsilon = 1e-12;
832 assert!((q.start() - c.eval(t0)).hypot() < epsilon);
833 assert!((q.end() - c.eval(t1)).hypot() < epsilon);
834 let n = 4;
835 for j in 0..=n {
836 let t = (j as f64) * (n as f64).recip();
837 let p = q.eval(t);
838 let err = (p.y - p.x.powi(3)).abs();
839 worst = worst.max(err);
840 assert!(err < accuracy, "got {err} wanted {accuracy}");
841 }
842 }
843 }
844 }
845
846 #[test]
847 fn cubicbez_approx_spline() {
848 let c1 = CubicBez::new(
849 (550.0, 258.0),
850 (1044.0, 482.0),
851 (2029.0, 1841.0),
852 (1934.0, 1554.0),
853 );
854
855 let quad = c1.try_approx_quadratic(344.0);
856 let expected = QuadBez::new(
857 Point::new(550.0, 258.0),
858 Point::new(1673.665720592873, 767.5164401068898),
859 Point::new(1934.0, 1554.0),
860 );
861 assert!(quad.is_some());
862 assert_eq!(quad.unwrap(), expected);
863
864 let quad = c1.try_approx_quadratic(343.0);
865 assert!(quad.is_none());
866
867 let spline = c1.approx_spline_n(2, 343.0);
868 assert!(spline.is_some());
869 let spline = spline.unwrap();
870 let expected = vec![
871 Point::new(550.0, 258.0),
872 Point::new(920.5, 426.0),
873 Point::new(2005.25, 1769.25),
874 Point::new(1934.0, 1554.0),
875 ];
876 assert_eq!(spline.points().len(), expected.len());
877 for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
878 assert!(got.distance(wanted) < 5.0)
879 }
880
881 let spline = c1.approx_spline(5.0);
882 let expected = vec![
883 Point::new(550.0, 258.0),
884 Point::new(673.5, 314.0),
885 Point::new(984.8777777777776, 584.2666666666667),
886 Point::new(1312.6305555555557, 927.825),
887 Point::new(1613.1194444444443, 1267.425),
888 Point::new(1842.7055555555555, 1525.8166666666666),
889 Point::new(1957.75, 1625.75),
890 Point::new(1934.0, 1554.0),
891 ];
892 assert!(spline.is_some());
893 let spline = spline.unwrap();
894 assert_eq!(spline.points().len(), expected.len());
895 for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
896 assert!(got.distance(wanted) < 5.0)
897 }
898 }
899
900 #[test]
901 fn cubicbez_cubics_to_quadratic_splines() {
902 let curves = vec![
903 CubicBez::new(
904 (550.0, 258.0),
905 (1044.0, 482.0),
906 (2029.0, 1841.0),
907 (1934.0, 1554.0),
908 ),
909 CubicBez::new(
910 (859.0, 384.0),
911 (1998.0, 116.0),
912 (1596.0, 1772.0),
913 (8.0, 1824.0),
914 ),
915 CubicBez::new(
916 (1090.0, 937.0),
917 (418.0, 1300.0),
918 (125.0, 91.0),
919 (104.0, 37.0),
920 ),
921 ];
922 let converted = cubics_to_quadratic_splines(&curves, 5.0);
923 assert!(converted.is_some());
924 let converted = converted.unwrap();
925 assert_eq!(converted[0].points().len(), 8);
926 assert_eq!(converted[1].points().len(), 8);
927 assert_eq!(converted[2].points().len(), 8);
928 assert!(converted[0].points()[1].distance(Point::new(673.5, 314.0)) < 0.0001);
929 assert!(
930 converted[0].points()[2].distance(Point::new(88639.0 / 90.0, 52584.0 / 90.0)) < 0.0001
931 );
932 }
933
934 #[test]
935 fn cubicbez_approx_spline_div_exact() {
936 // Ensure rounding behavior for division matches fonttools
937 // cu2qu.
938 // See <https://github.com/linebender/kurbo/issues/272>
939 let cubic = CubicBez::new(
940 Point::new(408.0, 321.0),
941 Point::new(408.0, 452.0),
942 Point::new(342.0, 560.0),
943 Point::new(260.0, 560.0),
944 );
945 let spline = cubic.approx_spline(1.0).unwrap();
946 assert_eq!(
947 spline.points(),
948 &[
949 Point::new(408.0, 321.0),
950 // Previous behavior produced 386.49999999999994 for the
951 // y coordinate leading to inconsistent rounding.
952 Point::new(408.0, 386.5),
953 Point::new(368.16666666666663, 495.0833333333333),
954 Point::new(301.0, 560.0),
955 Point::new(260.0, 560.0)
956 ]
957 )
958 }
959
960 #[test]
961 fn cubicbez_inflections() {
962 let c = CubicBez::new((0., 0.), (0.8, 1.), (0.2, 1.), (1., 0.));
963 let inflections = c.inflections();
964 assert_eq!(inflections.len(), 2);
965 assert!((inflections[0] - 0.311018).abs() < 1e-6);
966 assert!((inflections[1] - 0.688982).abs() < 1e-6);
967 let c = CubicBez::new((0., 0.), (1., 1.), (2., -1.), (3., 0.));
968 let inflections = c.inflections();
969 assert_eq!(inflections.len(), 1);
970 assert!((inflections[0] - 0.5).abs() < 1e-6);
971 let c = CubicBez::new((0., 0.), (1., 1.), (2., 1.), (3., 0.));
972 let inflections = c.inflections();
973 assert_eq!(inflections.len(), 0);
974 }
975
976 #[test]
977 fn cubic_to_quadratic_matches_python() {
978 // from https://github.com/googlefonts/fontmake-rs/issues/217
979 let cubic = CubicBez {
980 p0: (796.0, 319.0).into(),
981 p1: (727.0, 314.0).into(),
982 p2: (242.0, 303.0).into(),
983 p3: (106.0, 303.0).into(),
984 };
985
986 // FontTools can approximate this curve successfully in 7 splits, we can too
987 assert!(cubic.approx_spline_n(7, 1.0).is_some());
988
989 // FontTools can solve this with accuracy 0.001, we can too
990 assert!(cubics_to_quadratic_splines(&[cubic], 0.001).is_some());
991 }
992
993 #[test]
994 fn cubics_to_quadratic_splines_matches_python() {
995 // https://github.com/linebender/kurbo/pull/273
996 let light = CubicBez::new((378., 608.), (378., 524.), (355., 455.), (266., 455.));
997 let regular = CubicBez::new((367., 607.), (367., 511.), (338., 472.), (243., 472.));
998 let bold = CubicBez::new(
999 (372.425, 593.05),
1000 (372.425, 524.95),
1001 (355.05, 485.95),
1002 (274., 485.95),
1003 );
1004 let qsplines = cubics_to_quadratic_splines(&[light, regular, bold], 1.0).unwrap();
1005 assert_eq!(
1006 qsplines,
1007 [
1008 QuadSpline::new(vec![
1009 (378.0, 608.0).into(),
1010 (378.0, 566.0).into(),
1011 (359.0833333333333, 496.5).into(),
1012 (310.5, 455.0).into(),
1013 (266.0, 455.0).into(),
1014 ]),
1015 QuadSpline::new(vec![
1016 (367.0, 607.0).into(),
1017 (367.0, 559.0).into(),
1018 // Previous behavior produced 496.5 for the y coordinate
1019 (344.5833333333333, 499.49999999999994).into(),
1020 (290.5, 472.0).into(),
1021 (243.0, 472.0).into(),
1022 ]),
1023 QuadSpline::new(vec![
1024 (372.425, 593.05).into(),
1025 (372.425, 559.0).into(),
1026 (356.98333333333335, 511.125).into(),
1027 (314.525, 485.95).into(),
1028 (274.0, 485.95).into(),
1029 ]),
1030 ]
1031 )
1032 }
1033}
1034