1// Copyright 2006 The Android Open Source Project
2// Copyright 2020 Yevhenii Reizner
3//
4// Use of this source code is governed by a BSD-style license that can be
5// found in the LICENSE file.
6
7//! A collection of functions to work with Bezier paths.
8//!
9//! Mainly for internal use. Do not rely on it!
10
11#![allow(missing_docs)]
12
13use crate::{Point, Transform};
14
15use crate::f32x2_t::f32x2;
16use crate::floating_point::FLOAT_PI;
17use crate::scalar::{Scalar, SCALAR_NEARLY_ZERO, SCALAR_ROOT_2_OVER_2};
18
19use crate::floating_point::{NormalizedF32, NormalizedF32Exclusive};
20use crate::path_builder::PathDirection;
21
22#[cfg(all(not(feature = "std"), feature = "no-std-float"))]
23use crate::NoStdFloat;
24
25// use for : eval(t) == A * t^2 + B * t + C
26#[derive(Clone, Copy, Default, Debug)]
27pub struct QuadCoeff {
28 pub a: f32x2,
29 pub b: f32x2,
30 pub c: f32x2,
31}
32
33impl QuadCoeff {
34 pub fn from_points(points: &[Point; 3]) -> Self {
35 let c: f32x2 = points[0].to_f32x2();
36 let p1: f32x2 = points[1].to_f32x2();
37 let p2: f32x2 = points[2].to_f32x2();
38 let b: f32x2 = times_2(p1 - c);
39 let a: f32x2 = p2 - times_2(p1) + c;
40
41 QuadCoeff { a, b, c }
42 }
43
44 pub fn eval(&self, t: f32x2) -> f32x2 {
45 (self.a * t + self.b) * t + self.c
46 }
47}
48
49#[derive(Clone, Copy, Default, Debug)]
50pub struct CubicCoeff {
51 pub a: f32x2,
52 pub b: f32x2,
53 pub c: f32x2,
54 pub d: f32x2,
55}
56
57impl CubicCoeff {
58 pub fn from_points(points: &[Point; 4]) -> Self {
59 let p0: f32x2 = points[0].to_f32x2();
60 let p1: f32x2 = points[1].to_f32x2();
61 let p2: f32x2 = points[2].to_f32x2();
62 let p3: f32x2 = points[3].to_f32x2();
63 let three: f32x2 = f32x2::splat(3.0);
64
65 CubicCoeff {
66 a: p3 + three * (p1 - p2) - p0,
67 b: three * (p2 - times_2(p1) + p0),
68 c: three * (p1 - p0),
69 d: p0,
70 }
71 }
72
73 pub fn eval(&self, t: f32x2) -> f32x2 {
74 ((self.a * t + self.b) * t + self.c) * t + self.d
75 }
76}
77
78// TODO: to a custom type?
79pub fn new_t_values() -> [NormalizedF32Exclusive; 3] {
80 [NormalizedF32Exclusive::ANY; 3]
81}
82
83pub fn chop_quad_at(src: &[Point], t: NormalizedF32Exclusive, dst: &mut [Point; 5]) {
84 let p0: f32x2 = src[0].to_f32x2();
85 let p1: f32x2 = src[1].to_f32x2();
86 let p2: f32x2 = src[2].to_f32x2();
87 let tt: f32x2 = f32x2::splat(t.get());
88
89 let p01: f32x2 = interp(v0:p0, v1:p1, t:tt);
90 let p12: f32x2 = interp(v0:p1, v1:p2, t:tt);
91
92 dst[0] = Point::from_f32x2(p0);
93 dst[1] = Point::from_f32x2(p01);
94 dst[2] = Point::from_f32x2(interp(v0:p01, v1:p12, t:tt));
95 dst[3] = Point::from_f32x2(p12);
96 dst[4] = Point::from_f32x2(p2);
97}
98
99// From Numerical Recipes in C.
100//
101// Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
102// x1 = Q / A
103// x2 = C / Q
104pub fn find_unit_quad_roots(
105 a: f32,
106 b: f32,
107 c: f32,
108 roots: &mut [NormalizedF32Exclusive; 3],
109) -> usize {
110 if a == 0.0 {
111 if let Some(r) = valid_unit_divide(-c, b) {
112 roots[0] = r;
113 return 1;
114 } else {
115 return 0;
116 }
117 }
118
119 // use doubles so we don't overflow temporarily trying to compute R
120 let mut dr = f64::from(b) * f64::from(b) - 4.0 * f64::from(a) * f64::from(c);
121 if dr < 0.0 {
122 return 0;
123 }
124 dr = dr.sqrt();
125 let r = dr as f32;
126 if !r.is_finite() {
127 return 0;
128 }
129
130 let q = if b < 0.0 {
131 -(b - r) / 2.0
132 } else {
133 -(b + r) / 2.0
134 };
135
136 let mut roots_offset = 0;
137 if let Some(r) = valid_unit_divide(q, a) {
138 roots[roots_offset] = r;
139 roots_offset += 1;
140 }
141
142 if let Some(r) = valid_unit_divide(c, q) {
143 roots[roots_offset] = r;
144 roots_offset += 1;
145 }
146
147 if roots_offset == 2 {
148 if roots[0].get() > roots[1].get() {
149 roots.swap(0, 1);
150 } else if roots[0] == roots[1] {
151 // nearly-equal?
152 roots_offset -= 1; // skip the double root
153 }
154 }
155
156 roots_offset
157}
158
159pub fn chop_cubic_at2(src: &[Point; 4], t: NormalizedF32Exclusive, dst: &mut [Point]) {
160 let p0: f32x2 = src[0].to_f32x2();
161 let p1: f32x2 = src[1].to_f32x2();
162 let p2: f32x2 = src[2].to_f32x2();
163 let p3: f32x2 = src[3].to_f32x2();
164 let tt: f32x2 = f32x2::splat(t.get());
165
166 let ab: f32x2 = interp(v0:p0, v1:p1, t:tt);
167 let bc: f32x2 = interp(v0:p1, v1:p2, t:tt);
168 let cd: f32x2 = interp(v0:p2, v1:p3, t:tt);
169 let abc: f32x2 = interp(v0:ab, v1:bc, t:tt);
170 let bcd: f32x2 = interp(v0:bc, v1:cd, t:tt);
171 let abcd: f32x2 = interp(v0:abc, v1:bcd, t:tt);
172
173 dst[0] = Point::from_f32x2(p0);
174 dst[1] = Point::from_f32x2(ab);
175 dst[2] = Point::from_f32x2(abc);
176 dst[3] = Point::from_f32x2(abcd);
177 dst[4] = Point::from_f32x2(bcd);
178 dst[5] = Point::from_f32x2(cd);
179 dst[6] = Point::from_f32x2(p3);
180}
181
182// Quad'(t) = At + B, where
183// A = 2(a - 2b + c)
184// B = 2(b - a)
185// Solve for t, only if it fits between 0 < t < 1
186pub(crate) fn find_quad_extrema(a: f32, b: f32, c: f32) -> Option<NormalizedF32Exclusive> {
187 // At + B == 0
188 // t = -B / A
189 valid_unit_divide(numer:a - b, denom:a - b - b + c)
190}
191
192pub fn valid_unit_divide(mut numer: f32, mut denom: f32) -> Option<NormalizedF32Exclusive> {
193 if numer < 0.0 {
194 numer = -numer;
195 denom = -denom;
196 }
197
198 if denom == 0.0 || numer == 0.0 || numer >= denom {
199 return None;
200 }
201
202 let r: f32 = numer / denom;
203 NormalizedF32Exclusive::new(r)
204}
205
206fn interp(v0: f32x2, v1: f32x2, t: f32x2) -> f32x2 {
207 v0 + (v1 - v0) * t
208}
209
210fn times_2(value: f32x2) -> f32x2 {
211 value + value
212}
213
214// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2
215// F'(t) = 2 (b - a) + 2 (a - 2b + c) t
216// F''(t) = 2 (a - 2b + c)
217//
218// A = 2 (b - a)
219// B = 2 (a - 2b + c)
220//
221// Maximum curvature for a quadratic means solving
222// Fx' Fx'' + Fy' Fy'' = 0
223//
224// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2)
225pub(crate) fn find_quad_max_curvature(src: &[Point; 3]) -> NormalizedF32 {
226 let ax = src[1].x - src[0].x;
227 let ay = src[1].y - src[0].y;
228 let bx = src[0].x - src[1].x - src[1].x + src[2].x;
229 let by = src[0].y - src[1].y - src[1].y + src[2].y;
230
231 let mut numer = -(ax * bx + ay * by);
232 let mut denom = bx * bx + by * by;
233 if denom < 0.0 {
234 numer = -numer;
235 denom = -denom;
236 }
237
238 if numer <= 0.0 {
239 return NormalizedF32::ZERO;
240 }
241
242 if numer >= denom {
243 // Also catches denom=0
244 return NormalizedF32::ONE;
245 }
246
247 let t = numer / denom;
248 NormalizedF32::new(t).unwrap()
249}
250
251pub(crate) fn eval_quad_at(src: &[Point; 3], t: NormalizedF32) -> Point {
252 Point::from_f32x2(QuadCoeff::from_points(src).eval(f32x2::splat(t.get())))
253}
254
255pub(crate) fn eval_quad_tangent_at(src: &[Point; 3], tol: NormalizedF32) -> Point {
256 // The derivative equation is 2(b - a +(a - 2b +c)t). This returns a
257 // zero tangent vector when t is 0 or 1, and the control point is equal
258 // to the end point. In this case, use the quad end points to compute the tangent.
259 if (tol == NormalizedF32::ZERO && src[0] == src[1])
260 || (tol == NormalizedF32::ONE && src[1] == src[2])
261 {
262 return src[2] - src[0];
263 }
264
265 let p0: f32x2 = src[0].to_f32x2();
266 let p1: f32x2 = src[1].to_f32x2();
267 let p2: f32x2 = src[2].to_f32x2();
268
269 let b: f32x2 = p1 - p0;
270 let a: f32x2 = p2 - p1 - b;
271 let t: f32x2 = a * f32x2::splat(tol.get()) + b;
272
273 Point::from_f32x2(t + t)
274}
275
276// Looking for F' dot F'' == 0
277//
278// A = b - a
279// B = c - 2b + a
280// C = d - 3c + 3b - a
281//
282// F' = 3Ct^2 + 6Bt + 3A
283// F'' = 6Ct + 6B
284//
285// F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
286pub fn find_cubic_max_curvature<'a>(
287 src: &[Point; 4],
288 t_values: &'a mut [NormalizedF32; 3],
289) -> &'a [NormalizedF32] {
290 let mut coeff_x: [f32; 4] = formulate_f1_dot_f2(&[src[0].x, src[1].x, src[2].x, src[3].x]);
291 let coeff_y: [f32; 4] = formulate_f1_dot_f2(&[src[0].y, src[1].y, src[2].y, src[3].y]);
292
293 for i: usize in 0..4 {
294 coeff_x[i] += coeff_y[i];
295 }
296
297 let len: usize = solve_cubic_poly(&coeff_x, t_values);
298 &t_values[0..len]
299}
300
301// Looking for F' dot F'' == 0
302//
303// A = b - a
304// B = c - 2b + a
305// C = d - 3c + 3b - a
306//
307// F' = 3Ct^2 + 6Bt + 3A
308// F'' = 6Ct + 6B
309//
310// F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
311fn formulate_f1_dot_f2(src: &[f32; 4]) -> [f32; 4] {
312 let a: f32 = src[1] - src[0];
313 let b: f32 = src[2] - 2.0 * src[1] + src[0];
314 let c: f32 = src[3] + 3.0 * (src[1] - src[2]) - src[0];
315
316 [c * c, 3.0 * b * c, 2.0 * b * b + c * a, a * b]
317}
318
319/// Solve coeff(t) == 0, returning the number of roots that lie withing 0 < t < 1.
320/// coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
321///
322/// Eliminates repeated roots (so that all t_values are distinct, and are always
323/// in increasing order.
324fn solve_cubic_poly(coeff: &[f32; 4], t_values: &mut [NormalizedF32; 3]) -> usize {
325 if coeff[0].is_nearly_zero() {
326 // we're just a quadratic
327 let mut tmp_t = new_t_values();
328 let count = find_unit_quad_roots(coeff[1], coeff[2], coeff[3], &mut tmp_t);
329 for i in 0..count {
330 t_values[i] = tmp_t[i].to_normalized();
331 }
332
333 return count;
334 }
335
336 debug_assert!(coeff[0] != 0.0);
337
338 let inva = coeff[0].invert();
339 let a = coeff[1] * inva;
340 let b = coeff[2] * inva;
341 let c = coeff[3] * inva;
342
343 let q = (a * a - b * 3.0) / 9.0;
344 let r = (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) / 54.0;
345
346 let q3 = q * q * q;
347 let r2_minus_q3 = r * r - q3;
348 let adiv3 = a / 3.0;
349
350 if r2_minus_q3 < 0.0 {
351 // we have 3 real roots
352 // the divide/root can, due to finite precisions, be slightly outside of -1...1
353 let theta = (r / q3.sqrt()).bound(-1.0, 1.0).acos();
354 let neg2_root_q = -2.0 * q.sqrt();
355
356 t_values[0] = NormalizedF32::new_clamped(neg2_root_q * (theta / 3.0).cos() - adiv3);
357 t_values[1] = NormalizedF32::new_clamped(
358 neg2_root_q * ((theta + 2.0 * FLOAT_PI) / 3.0).cos() - adiv3,
359 );
360 t_values[2] = NormalizedF32::new_clamped(
361 neg2_root_q * ((theta - 2.0 * FLOAT_PI) / 3.0).cos() - adiv3,
362 );
363
364 // now sort the roots
365 sort_array3(t_values);
366 collapse_duplicates3(t_values)
367 } else {
368 // we have 1 real root
369 let mut a = r.abs() + r2_minus_q3.sqrt();
370 a = scalar_cube_root(a);
371 if r > 0.0 {
372 a = -a;
373 }
374
375 if a != 0.0 {
376 a += q / a;
377 }
378
379 t_values[0] = NormalizedF32::new_clamped(a - adiv3);
380 1
381 }
382}
383
384fn sort_array3(array: &mut [NormalizedF32; 3]) {
385 if array[0] > array[1] {
386 array.swap(a:0, b:1);
387 }
388
389 if array[1] > array[2] {
390 array.swap(a:1, b:2);
391 }
392
393 if array[0] > array[1] {
394 array.swap(a:0, b:1);
395 }
396}
397
398fn collapse_duplicates3(array: &mut [NormalizedF32; 3]) -> usize {
399 let mut len: usize = 3;
400
401 if array[1] == array[2] {
402 len = 2;
403 }
404
405 if array[0] == array[1] {
406 len = 1;
407 }
408
409 len
410}
411
412fn scalar_cube_root(x: f32) -> f32 {
413 x.powf(0.3333333)
414}
415
416// This is SkEvalCubicAt split into three functions.
417pub(crate) fn eval_cubic_pos_at(src: &[Point; 4], t: NormalizedF32) -> Point {
418 Point::from_f32x2(CubicCoeff::from_points(src).eval(f32x2::splat(t.get())))
419}
420
421// This is SkEvalCubicAt split into three functions.
422pub(crate) fn eval_cubic_tangent_at(src: &[Point; 4], t: NormalizedF32) -> Point {
423 // The derivative equation returns a zero tangent vector when t is 0 or 1, and the
424 // adjacent control point is equal to the end point. In this case, use the
425 // next control point or the end points to compute the tangent.
426 if (t.get() == 0.0 && src[0] == src[1]) || (t.get() == 1.0 && src[2] == src[3]) {
427 let mut tangent: Point = if t.get() == 0.0 {
428 src[2] - src[0]
429 } else {
430 src[3] - src[1]
431 };
432
433 if tangent.x == 0.0 && tangent.y == 0.0 {
434 tangent = src[3] - src[0];
435 }
436
437 tangent
438 } else {
439 eval_cubic_derivative(src, t)
440 }
441}
442
443fn eval_cubic_derivative(src: &[Point; 4], t: NormalizedF32) -> Point {
444 let p0: f32x2 = src[0].to_f32x2();
445 let p1: f32x2 = src[1].to_f32x2();
446 let p2: f32x2 = src[2].to_f32x2();
447 let p3: f32x2 = src[3].to_f32x2();
448
449 let coeff: QuadCoeff = QuadCoeff {
450 a: p3 + f32x2::splat(3.0) * (p1 - p2) - p0,
451 b: times_2(p2 - times_2(p1) + p0),
452 c: p1 - p0,
453 };
454
455 Point::from_f32x2(coeff.eval(f32x2::splat(t.get())))
456}
457
458// Cubic'(t) = At^2 + Bt + C, where
459// A = 3(-a + 3(b - c) + d)
460// B = 6(a - 2b + c)
461// C = 3(b - a)
462// Solve for t, keeping only those that fit between 0 < t < 1
463pub(crate) fn find_cubic_extrema(
464 a: f32,
465 b: f32,
466 c: f32,
467 d: f32,
468 t_values: &mut [NormalizedF32Exclusive; 3],
469) -> usize {
470 // we divide A,B,C by 3 to simplify
471 let aa: f32 = d - a + 3.0 * (b - c);
472 let bb: f32 = 2.0 * (a - b - b + c);
473 let cc: f32 = b - a;
474
475 find_unit_quad_roots(a:aa, b:bb, c:cc, roots:t_values)
476}
477
478// http://www.faculty.idc.ac.il/arik/quality/appendixA.html
479//
480// Inflection means that curvature is zero.
481// Curvature is [F' x F''] / [F'^3]
482// So we solve F'x X F''y - F'y X F''y == 0
483// After some canceling of the cubic term, we get
484// A = b - a
485// B = c - 2b + a
486// C = d - 3c + 3b - a
487// (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
488pub(crate) fn find_cubic_inflections<'a>(
489 src: &[Point; 4],
490 t_values: &'a mut [NormalizedF32Exclusive; 3],
491) -> &'a [NormalizedF32Exclusive] {
492 let ax: f32 = src[1].x - src[0].x;
493 let ay: f32 = src[1].y - src[0].y;
494 let bx: f32 = src[2].x - 2.0 * src[1].x + src[0].x;
495 let by: f32 = src[2].y - 2.0 * src[1].y + src[0].y;
496 let cx: f32 = src[3].x + 3.0 * (src[1].x - src[2].x) - src[0].x;
497 let cy: f32 = src[3].y + 3.0 * (src[1].y - src[2].y) - src[0].y;
498
499 let len: usize = find_unit_quad_roots(
500 a:bx * cy - by * cx,
501 b:ax * cy - ay * cx,
502 c:ax * by - ay * bx,
503 roots:t_values,
504 );
505
506 &t_values[0..len]
507}
508
509// Return location (in t) of cubic cusp, if there is one.
510// Note that classify cubic code does not reliably return all cusp'd cubics, so
511// it is not called here.
512pub(crate) fn find_cubic_cusp(src: &[Point; 4]) -> Option<NormalizedF32Exclusive> {
513 // When the adjacent control point matches the end point, it behaves as if
514 // the cubic has a cusp: there's a point of max curvature where the derivative
515 // goes to zero. Ideally, this would be where t is zero or one, but math
516 // error makes not so. It is not uncommon to create cubics this way; skip them.
517 if src[0] == src[1] {
518 return None;
519 }
520
521 if src[2] == src[3] {
522 return None;
523 }
524
525 // Cubics only have a cusp if the line segments formed by the control and end points cross.
526 // Detect crossing if line ends are on opposite sides of plane formed by the other line.
527 if on_same_side(src, 0, 2) || on_same_side(src, 2, 0) {
528 return None;
529 }
530
531 // Cubics may have multiple points of maximum curvature, although at most only
532 // one is a cusp.
533 let mut t_values = [NormalizedF32::ZERO; 3];
534 let max_curvature = find_cubic_max_curvature(src, &mut t_values);
535 for test_t in max_curvature {
536 if 0.0 >= test_t.get() || test_t.get() >= 1.0 {
537 // no need to consider max curvature on the end
538 continue;
539 }
540
541 // A cusp is at the max curvature, and also has a derivative close to zero.
542 // Choose the 'close to zero' meaning by comparing the derivative length
543 // with the overall cubic size.
544 let d_pt = eval_cubic_derivative(src, *test_t);
545 let d_pt_magnitude = d_pt.length_sqd();
546 let precision = calc_cubic_precision(src);
547 if d_pt_magnitude < precision {
548 // All three max curvature t values may be close to the cusp;
549 // return the first one.
550 return Some(NormalizedF32Exclusive::new_bounded(test_t.get()));
551 }
552 }
553
554 None
555}
556
557// Returns true if both points src[testIndex], src[testIndex+1] are in the same half plane defined
558// by the line segment src[lineIndex], src[lineIndex+1].
559fn on_same_side(src: &[Point; 4], test_index: usize, line_index: usize) -> bool {
560 let origin: Point = src[line_index];
561 let line: Point = src[line_index + 1] - origin;
562 let mut crosses: [f32; 2] = [0.0, 0.0];
563 for index: usize in 0..2 {
564 let test_line: Point = src[test_index + index] - origin;
565 crosses[index] = line.cross(test_line);
566 }
567
568 crosses[0] * crosses[1] >= 0.0
569}
570
571// Returns a constant proportional to the dimensions of the cubic.
572// Constant found through experimentation -- maybe there's a better way....
573fn calc_cubic_precision(src: &[Point; 4]) -> f32 {
574 (src[1].distance_to_sqd(pt:src[0])
575 + src[2].distance_to_sqd(pt:src[1])
576 + src[3].distance_to_sqd(pt:src[2]))
577 * 1e-8
578}
579
580#[derive(Copy, Clone, Default, Debug)]
581pub(crate) struct Conic {
582 pub points: [Point; 3],
583 pub weight: f32,
584}
585
586impl Conic {
587 pub fn new(pt0: Point, pt1: Point, pt2: Point, weight: f32) -> Self {
588 Conic {
589 points: [pt0, pt1, pt2],
590 weight,
591 }
592 }
593
594 pub fn from_points(points: &[Point], weight: f32) -> Self {
595 Conic {
596 points: [points[0], points[1], points[2]],
597 weight,
598 }
599 }
600
601 fn compute_quad_pow2(&self, tolerance: f32) -> Option<u8> {
602 if tolerance < 0.0 || !tolerance.is_finite() {
603 return None;
604 }
605
606 if !self.points[0].is_finite() || !self.points[1].is_finite() || !self.points[2].is_finite()
607 {
608 return None;
609 }
610
611 // Limit the number of suggested quads to approximate a conic
612 const MAX_CONIC_TO_QUAD_POW2: usize = 4;
613
614 // "High order approximation of conic sections by quadratic splines"
615 // by Michael Floater, 1993
616 let a = self.weight - 1.0;
617 let k = a / (4.0 * (2.0 + a));
618 let x = k * (self.points[0].x - 2.0 * self.points[1].x + self.points[2].x);
619 let y = k * (self.points[0].y - 2.0 * self.points[1].y + self.points[2].y);
620
621 let mut error = (x * x + y * y).sqrt();
622 let mut pow2 = 0;
623 for _ in 0..MAX_CONIC_TO_QUAD_POW2 {
624 if error <= tolerance {
625 break;
626 }
627
628 error *= 0.25;
629 pow2 += 1;
630 }
631
632 // Unlike Skia, we always expect `pow2` to be at least 1.
633 // Otherwise it produces ugly results.
634 Some(pow2.max(1))
635 }
636
637 // Chop this conic into N quads, stored continuously in pts[], where
638 // N = 1 << pow2. The amount of storage needed is (1 + 2 * N)
639 pub fn chop_into_quads_pow2(&self, pow2: u8, points: &mut [Point]) -> u8 {
640 debug_assert!(pow2 < 5);
641
642 points[0] = self.points[0];
643 subdivide(self, &mut points[1..], pow2);
644
645 let quad_count = 1 << pow2;
646 let pt_count = 2 * quad_count + 1;
647 if points.iter().take(pt_count).any(|n| !n.is_finite()) {
648 // if we generated a non-finite, pin ourselves to the middle of the hull,
649 // as our first and last are already on the first/last pts of the hull.
650 for p in points.iter_mut().take(pt_count - 1).skip(1) {
651 *p = self.points[1];
652 }
653 }
654
655 1 << pow2
656 }
657
658 fn chop(&self) -> (Conic, Conic) {
659 let scale = f32x2::splat((1.0 + self.weight).invert());
660 let new_w = subdivide_weight_value(self.weight);
661
662 let p0 = self.points[0].to_f32x2();
663 let p1 = self.points[1].to_f32x2();
664 let p2 = self.points[2].to_f32x2();
665 let ww = f32x2::splat(self.weight);
666
667 let wp1 = ww * p1;
668 let m = (p0 + times_2(wp1) + p2) * scale * f32x2::splat(0.5);
669 let mut m_pt = Point::from_f32x2(m);
670 if !m_pt.is_finite() {
671 let w_d = self.weight as f64;
672 let w_2 = w_d * 2.0;
673 let scale_half = 1.0 / (1.0 + w_d) * 0.5;
674 m_pt.x = ((self.points[0].x as f64
675 + w_2 * self.points[1].x as f64
676 + self.points[2].x as f64)
677 * scale_half) as f32;
678
679 m_pt.y = ((self.points[0].y as f64
680 + w_2 * self.points[1].y as f64
681 + self.points[2].y as f64)
682 * scale_half) as f32;
683 }
684
685 (
686 Conic {
687 points: [self.points[0], Point::from_f32x2((p0 + wp1) * scale), m_pt],
688 weight: new_w,
689 },
690 Conic {
691 points: [m_pt, Point::from_f32x2((wp1 + p2) * scale), self.points[2]],
692 weight: new_w,
693 },
694 )
695 }
696
697 pub fn build_unit_arc(
698 u_start: Point,
699 u_stop: Point,
700 dir: PathDirection,
701 user_transform: Transform,
702 dst: &mut [Conic; 5],
703 ) -> Option<&[Conic]> {
704 // rotate by x,y so that u_start is (1.0)
705 let x = u_start.dot(u_stop);
706 let mut y = u_start.cross(u_stop);
707
708 let abs_y = y.abs();
709
710 // check for (effectively) coincident vectors
711 // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
712 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
713 if abs_y <= SCALAR_NEARLY_ZERO
714 && x > 0.0
715 && ((y >= 0.0 && dir == PathDirection::CW) || (y <= 0.0 && dir == PathDirection::CCW))
716 {
717 return None;
718 }
719
720 if dir == PathDirection::CCW {
721 y = -y;
722 }
723
724 // We decide to use 1-conic per quadrant of a circle. What quadrant does [xy] lie in?
725 // 0 == [0 .. 90)
726 // 1 == [90 ..180)
727 // 2 == [180..270)
728 // 3 == [270..360)
729 //
730 let mut quadrant = 0;
731 if y == 0.0 {
732 quadrant = 2; // 180
733 debug_assert!((x + 1.0) <= SCALAR_NEARLY_ZERO);
734 } else if x == 0.0 {
735 debug_assert!(abs_y - 1.0 <= SCALAR_NEARLY_ZERO);
736 quadrant = if y > 0.0 { 1 } else { 3 }; // 90 / 270
737 } else {
738 if y < 0.0 {
739 quadrant += 2;
740 }
741
742 if (x < 0.0) != (y < 0.0) {
743 quadrant += 1;
744 }
745 }
746
747 let quadrant_points = [
748 Point::from_xy(1.0, 0.0),
749 Point::from_xy(1.0, 1.0),
750 Point::from_xy(0.0, 1.0),
751 Point::from_xy(-1.0, 1.0),
752 Point::from_xy(-1.0, 0.0),
753 Point::from_xy(-1.0, -1.0),
754 Point::from_xy(0.0, -1.0),
755 Point::from_xy(1.0, -1.0),
756 ];
757
758 const QUADRANT_WEIGHT: f32 = SCALAR_ROOT_2_OVER_2;
759
760 let mut conic_count = quadrant;
761 for i in 0..conic_count {
762 dst[i] = Conic::from_points(&quadrant_points[i * 2..], QUADRANT_WEIGHT);
763 }
764
765 // Now compute any remaing (sub-90-degree) arc for the last conic
766 let final_pt = Point::from_xy(x, y);
767 let last_q = quadrant_points[quadrant * 2]; // will already be a unit-vector
768 let dot = last_q.dot(final_pt);
769 debug_assert!(0.0 <= dot && dot <= 1.0 + SCALAR_NEARLY_ZERO);
770
771 if dot < 1.0 {
772 let mut off_curve = Point::from_xy(last_q.x + x, last_q.y + y);
773 // compute the bisector vector, and then rescale to be the off-curve point.
774 // we compute its length from cos(theta/2) = length / 1, using half-angle identity we get
775 // length = sqrt(2 / (1 + cos(theta)). We already have cos() when to computed the dot.
776 // This is nice, since our computed weight is cos(theta/2) as well!
777 let cos_theta_over_2 = ((1.0 + dot) / 2.0).sqrt();
778 off_curve.set_length(cos_theta_over_2.invert());
779 if !last_q.almost_equal(off_curve) {
780 dst[conic_count] = Conic::new(last_q, off_curve, final_pt, cos_theta_over_2);
781 conic_count += 1;
782 }
783 }
784
785 // now handle counter-clockwise and the initial unitStart rotation
786 let mut transform = Transform::from_sin_cos(u_start.y, u_start.x);
787 if dir == PathDirection::CCW {
788 transform = transform.pre_scale(1.0, -1.0);
789 }
790
791 transform = transform.post_concat(user_transform);
792
793 for conic in dst.iter_mut().take(conic_count) {
794 transform.map_points(&mut conic.points);
795 }
796
797 if conic_count == 0 {
798 None
799 } else {
800 Some(&dst[0..conic_count])
801 }
802 }
803}
804
805fn subdivide_weight_value(w: f32) -> f32 {
806 (0.5 + w * 0.5).sqrt()
807}
808
809fn subdivide<'a>(src: &Conic, mut points: &'a mut [Point], mut level: u8) -> &'a mut [Point] {
810 if level == 0 {
811 points[0] = src.points[1];
812 points[1] = src.points[2];
813 &mut points[2..]
814 } else {
815 let mut dst = src.chop();
816
817 let start_y = src.points[0].y;
818 let end_y = src.points[2].y;
819 if between(start_y, src.points[1].y, end_y) {
820 // If the input is monotonic and the output is not, the scan converter hangs.
821 // Ensure that the chopped conics maintain their y-order.
822 let mid_y = dst.0.points[2].y;
823 if !between(start_y, mid_y, end_y) {
824 // If the computed midpoint is outside the ends, move it to the closer one.
825 let closer_y = if (mid_y - start_y).abs() < (mid_y - end_y).abs() {
826 start_y
827 } else {
828 end_y
829 };
830 dst.0.points[2].y = closer_y;
831 dst.1.points[0].y = closer_y;
832 }
833
834 if !between(start_y, dst.0.points[1].y, dst.0.points[2].y) {
835 // If the 1st control is not between the start and end, put it at the start.
836 // This also reduces the quad to a line.
837 dst.0.points[1].y = start_y;
838 }
839
840 if !between(dst.1.points[0].y, dst.1.points[1].y, end_y) {
841 // If the 2nd control is not between the start and end, put it at the end.
842 // This also reduces the quad to a line.
843 dst.1.points[1].y = end_y;
844 }
845
846 // Verify that all five points are in order.
847 debug_assert!(between(start_y, dst.0.points[1].y, dst.0.points[2].y));
848 debug_assert!(between(
849 dst.0.points[1].y,
850 dst.0.points[2].y,
851 dst.1.points[1].y
852 ));
853 debug_assert!(between(dst.0.points[2].y, dst.1.points[1].y, end_y));
854 }
855
856 level -= 1;
857 points = subdivide(&dst.0, points, level);
858 subdivide(&dst.1, points, level)
859 }
860}
861
862// This was originally developed and tested for pathops: see SkOpTypes.h
863// returns true if (a <= b <= c) || (a >= b >= c)
864fn between(a: f32, b: f32, c: f32) -> bool {
865 (a - b) * (c - b) <= 0.0
866}
867
868pub(crate) struct AutoConicToQuads {
869 pub points: [Point; 64],
870 pub len: u8, // the number of quads
871}
872
873impl AutoConicToQuads {
874 pub fn compute(pt0: Point, pt1: Point, pt2: Point, weight: f32) -> Option<Self> {
875 let conic: Conic = Conic::new(pt0, pt1, pt2, weight);
876 let pow2: u8 = conic.compute_quad_pow2(tolerance:0.25)?;
877 let mut points: [Point; 64] = [Point::zero(); 64];
878 let len: u8 = conic.chop_into_quads_pow2(pow2, &mut points);
879 Some(AutoConicToQuads { points, len })
880 }
881}
882
883#[cfg(test)]
884mod tests {
885 use super::*;
886
887 #[test]
888 fn eval_cubic_at_1() {
889 let src = [
890 Point::from_xy(30.0, 40.0),
891 Point::from_xy(30.0, 40.0),
892 Point::from_xy(171.0, 45.0),
893 Point::from_xy(180.0, 155.0),
894 ];
895
896 assert_eq!(
897 eval_cubic_pos_at(&src, NormalizedF32::ZERO),
898 Point::from_xy(30.0, 40.0)
899 );
900 assert_eq!(
901 eval_cubic_tangent_at(&src, NormalizedF32::ZERO),
902 Point::from_xy(141.0, 5.0)
903 );
904 }
905
906 #[test]
907 fn find_cubic_max_curvature_1() {
908 let src = [
909 Point::from_xy(20.0, 160.0),
910 Point::from_xy(20.0001, 160.0),
911 Point::from_xy(160.0, 20.0),
912 Point::from_xy(160.0001, 20.0),
913 ];
914
915 let mut t_values = [NormalizedF32::ZERO; 3];
916 let t_values = find_cubic_max_curvature(&src, &mut t_values);
917
918 assert_eq!(
919 &t_values,
920 &[
921 NormalizedF32::ZERO,
922 NormalizedF32::new_clamped(0.5),
923 NormalizedF32::ONE,
924 ]
925 );
926 }
927}
928