1// Copyright 2012 Google Inc.
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
7use super::point64::{Point64, SearchAxis};
8use super::quad64;
9use super::Scalar64;
10
11#[cfg(all(not(feature = "std"), feature = "no-std-float"))]
12use tiny_skia_path::NoStdFloat;
13
14pub const POINT_COUNT: usize = 4;
15const PI: f64 = 3.141592653589793;
16
17pub struct Cubic64Pair {
18 pub points: [Point64; 7],
19}
20
21pub struct Cubic64 {
22 pub points: [Point64; POINT_COUNT],
23}
24
25impl Cubic64 {
26 pub fn new(points: [Point64; POINT_COUNT]) -> Self {
27 Cubic64 { points }
28 }
29
30 pub fn as_f64_slice(&self) -> [f64; POINT_COUNT * 2] {
31 [
32 self.points[0].x,
33 self.points[0].y,
34 self.points[1].x,
35 self.points[1].y,
36 self.points[2].x,
37 self.points[2].y,
38 self.points[3].x,
39 self.points[3].y,
40 ]
41 }
42
43 pub fn point_at_t(&self, t: f64) -> Point64 {
44 if t == 0.0 {
45 return self.points[0];
46 }
47
48 if t == 1.0 {
49 return self.points[3];
50 }
51
52 let one_t = 1.0 - t;
53 let one_t2 = one_t * one_t;
54 let a = one_t2 * one_t;
55 let b = 3.0 * one_t2 * t;
56 let t2 = t * t;
57 let c = 3.0 * one_t * t2;
58 let d = t2 * t;
59 Point64::from_xy(
60 a * self.points[0].x
61 + b * self.points[1].x
62 + c * self.points[2].x
63 + d * self.points[3].x,
64 a * self.points[0].y
65 + b * self.points[1].y
66 + c * self.points[2].y
67 + d * self.points[3].y,
68 )
69 }
70
71 pub fn search_roots(
72 &self,
73 mut extrema: usize,
74 axis_intercept: f64,
75 x_axis: SearchAxis,
76 extreme_ts: &mut [f64; 6],
77 valid_roots: &mut [f64],
78 ) -> usize {
79 extrema += self.find_inflections(&mut extreme_ts[extrema..]);
80 extreme_ts[extrema] = 0.0;
81 extrema += 1;
82 extreme_ts[extrema] = 1.0;
83 debug_assert!(extrema < 6);
84 extreme_ts[0..extrema].sort_by(cmp_f64);
85 let mut valid_count = 0;
86 let mut index = 0;
87 while index < extrema {
88 let min = extreme_ts[index];
89 index += 1;
90 let max = extreme_ts[index];
91 if min == max {
92 continue;
93 }
94
95 let new_t = self.binary_search(min, max, axis_intercept, x_axis);
96 if new_t >= 0.0 {
97 if valid_count >= 3 {
98 return 0;
99 }
100
101 valid_roots[valid_count] = new_t;
102 valid_count += 1;
103 }
104 }
105
106 valid_count
107 }
108
109 fn find_inflections(&self, t_values: &mut [f64]) -> usize {
110 let ax = self.points[1].x - self.points[0].x;
111 let ay = self.points[1].y - self.points[0].y;
112 let bx = self.points[2].x - 2.0 * self.points[1].x + self.points[0].x;
113 let by = self.points[2].y - 2.0 * self.points[1].y + self.points[0].y;
114 let cx = self.points[3].x + 3.0 * (self.points[1].x - self.points[2].x) - self.points[0].x;
115 let cy = self.points[3].y + 3.0 * (self.points[1].y - self.points[2].y) - self.points[0].y;
116 quad64::roots_valid_t(
117 bx * cy - by * cx,
118 ax * cy - ay * cx,
119 ax * by - ay * bx,
120 t_values,
121 )
122 }
123
124 // give up when changing t no longer moves point
125 // also, copy point rather than recompute it when it does change
126 fn binary_search(&self, min: f64, max: f64, axis_intercept: f64, x_axis: SearchAxis) -> f64 {
127 let mut t = (min + max) / 2.0;
128 let mut step = (t - min) / 2.0;
129 let mut cubic_at_t = self.point_at_t(t);
130 let mut calc_pos = cubic_at_t.axis_coord(x_axis);
131 let mut calc_dist = calc_pos - axis_intercept;
132 loop {
133 let prior_t = min.max(t - step);
134 let less_pt = self.point_at_t(prior_t);
135 if less_pt.x.approximately_equal_half(cubic_at_t.x)
136 && less_pt.y.approximately_equal_half(cubic_at_t.y)
137 {
138 return -1.0; // binary search found no point at this axis intercept
139 }
140
141 let less_dist = less_pt.axis_coord(x_axis) - axis_intercept;
142 let last_step = step;
143 step /= 2.0;
144 let ok = if calc_dist > 0.0 {
145 calc_dist > less_dist
146 } else {
147 calc_dist < less_dist
148 };
149 if ok {
150 t = prior_t;
151 } else {
152 let next_t = t + last_step;
153 if next_t > max {
154 return -1.0;
155 }
156
157 let more_pt = self.point_at_t(next_t);
158 if more_pt.x.approximately_equal_half(cubic_at_t.x)
159 && more_pt.y.approximately_equal_half(cubic_at_t.y)
160 {
161 return -1.0; // binary search found no point at this axis intercept
162 }
163
164 let more_dist = more_pt.axis_coord(x_axis) - axis_intercept;
165 let ok = if calc_dist > 0.0 {
166 calc_dist <= more_dist
167 } else {
168 calc_dist >= more_dist
169 };
170 if ok {
171 continue;
172 }
173
174 t = next_t;
175 }
176
177 let test_at_t = self.point_at_t(t);
178 cubic_at_t = test_at_t;
179 calc_pos = cubic_at_t.axis_coord(x_axis);
180 calc_dist = calc_pos - axis_intercept;
181
182 if calc_pos.approximately_equal(axis_intercept) {
183 break;
184 }
185 }
186
187 t
188 }
189
190 pub fn chop_at(&self, t: f64) -> Cubic64Pair {
191 let mut dst = [Point64::zero(); 7];
192 if t == 0.5 {
193 dst[0] = self.points[0];
194 dst[1].x = (self.points[0].x + self.points[1].x) / 2.0;
195 dst[1].y = (self.points[0].y + self.points[1].y) / 2.0;
196 dst[2].x = (self.points[0].x + 2.0 * self.points[1].x + self.points[2].x) / 4.0;
197 dst[2].y = (self.points[0].y + 2.0 * self.points[1].y + self.points[2].y) / 4.0;
198 dst[3].x =
199 (self.points[0].x + 3.0 * (self.points[1].x + self.points[2].x) + self.points[3].x)
200 / 8.0;
201 dst[3].y =
202 (self.points[0].y + 3.0 * (self.points[1].y + self.points[2].y) + self.points[3].y)
203 / 8.0;
204 dst[4].x = (self.points[1].x + 2.0 * self.points[2].x + self.points[3].x) / 4.0;
205 dst[4].y = (self.points[1].y + 2.0 * self.points[2].y + self.points[3].y) / 4.0;
206 dst[5].x = (self.points[2].x + self.points[3].x) / 2.0;
207 dst[5].y = (self.points[2].y + self.points[3].y) / 2.0;
208 dst[6] = self.points[3];
209
210 Cubic64Pair { points: dst }
211 } else {
212 interp_cubic_coords_x(&self.points, t, &mut dst);
213 interp_cubic_coords_y(&self.points, t, &mut dst);
214 Cubic64Pair { points: dst }
215 }
216 }
217}
218
219pub fn coefficients(src: &[f64]) -> (f64, f64, f64, f64) {
220 let mut a: f64 = src[6]; // d
221 let mut b: f64 = src[4] * 3.0; // 3*c
222 let mut c: f64 = src[2] * 3.0; // 3*b
223 let d: f64 = src[0]; // a
224 a -= d - c + b; // A = -a + 3*b - 3*c + d
225 b += 3.0 * d - 2.0 * c; // B = 3*a - 6*b + 3*c
226 c -= 3.0 * d; // C = -3*a + 3*b
227
228 (a, b, c, d)
229}
230
231// from SkGeometry.cpp (and Numeric Solutions, 5.6)
232pub fn roots_valid_t(a: f64, b: f64, c: f64, d: f64, t: &mut [f64; 3]) -> usize {
233 let mut s = [0.0; 3];
234 let real_roots = roots_real(a, b, c, d, &mut s);
235 let mut found_roots = quad64::push_valid_ts(&s, real_roots, t);
236 'outer: for index in 0..real_roots {
237 let t_value = s[index];
238 if !t_value.approximately_one_or_less() && t_value.between(1.0, 1.00005) {
239 for idx2 in 0..found_roots {
240 if t[idx2].approximately_equal(1.0) {
241 continue 'outer;
242 }
243 }
244
245 debug_assert!(found_roots < 3);
246 t[found_roots] = 1.0;
247 found_roots += 1;
248 } else if !t_value.approximately_zero_or_more() && t_value.between(-0.00005, 0.0) {
249 for idx2 in 0..found_roots {
250 if t[idx2].approximately_equal(0.0) {
251 continue 'outer;
252 }
253 }
254
255 debug_assert!(found_roots < 3);
256 t[found_roots] = 0.0;
257 found_roots += 1;
258 }
259 }
260
261 found_roots
262}
263
264fn roots_real(a: f64, b: f64, c: f64, d: f64, s: &mut [f64; 3]) -> usize {
265 if a.approximately_zero()
266 && a.approximately_zero_when_compared_to(b)
267 && a.approximately_zero_when_compared_to(c)
268 && a.approximately_zero_when_compared_to(d)
269 {
270 // we're just a quadratic
271 return quad64::roots_real(b, c, d, s);
272 }
273
274 if d.approximately_zero_when_compared_to(a)
275 && d.approximately_zero_when_compared_to(b)
276 && d.approximately_zero_when_compared_to(c)
277 {
278 // 0 is one root
279 let mut num = quad64::roots_real(a, b, c, s);
280 for i in 0..num {
281 if s[i].approximately_zero() {
282 return num;
283 }
284 }
285
286 s[num] = 0.0;
287 num += 1;
288
289 return num;
290 }
291
292 if (a + b + c + d).approximately_zero() {
293 // 1 is one root
294 let mut num = quad64::roots_real(a, a + b, -d, s);
295 for i in 0..num {
296 if s[i].almost_dequal_ulps(1.0) {
297 return num;
298 }
299 }
300 s[num] = 1.0;
301 num += 1;
302 return num;
303 }
304
305 let (a, b, c) = {
306 let inv_a = 1.0 / a;
307 let a = b * inv_a;
308 let b = c * inv_a;
309 let c = d * inv_a;
310 (a, b, c)
311 };
312
313 let a2 = a * a;
314 let q = (a2 - b * 3.0) / 9.0;
315 let r = (2.0 * a2 * a - 9.0 * a * b + 27.0 * c) / 54.0;
316 let r2 = r * r;
317 let q3 = q * q * q;
318 let r2_minus_q3 = r2 - q3;
319 let adiv3 = a / 3.0;
320 let mut offset = 0;
321 if r2_minus_q3 < 0.0 {
322 // we have 3 real roots
323
324 // the divide/root can, due to finite precisions, be slightly outside of -1...1
325 let theta = (r / q3.sqrt()).bound(-1.0, 1.0).acos();
326 let neg2_root_q = -2.0 * q.sqrt();
327
328 let mut rr = neg2_root_q * (theta / 3.0).cos() - adiv3;
329 s[offset] = rr;
330 offset += 1;
331
332 rr = neg2_root_q * ((theta + 2.0 * PI) / 3.0).cos() - adiv3;
333 if !s[0].almost_dequal_ulps(rr) {
334 s[offset] = rr;
335 offset += 1;
336 }
337
338 rr = neg2_root_q * ((theta - 2.0 * PI) / 3.0).cos() - adiv3;
339 if !s[0].almost_dequal_ulps(rr) && (offset == 1 || !s[1].almost_dequal_ulps(rr)) {
340 s[offset] = rr;
341 offset += 1;
342 }
343 } else {
344 // we have 1 real root
345 let sqrt_r2_minus_q3 = r2_minus_q3.sqrt();
346 let mut a = r.abs() + sqrt_r2_minus_q3;
347 a = super::cube_root(a);
348 if r > 0.0 {
349 a = -a;
350 }
351
352 if a != 0.0 {
353 a += q / a;
354 }
355
356 let mut r2 = a - adiv3;
357 s[offset] = r2;
358 offset += 1;
359 if r2.almost_dequal_ulps(q3) {
360 r2 = -a / 2.0 - adiv3;
361 if !s[0].almost_dequal_ulps(r2) {
362 s[offset] = r2;
363 offset += 1;
364 }
365 }
366 }
367
368 offset
369}
370
371// Cubic64'(t) = At^2 + Bt + C, where
372// A = 3(-a + 3(b - c) + d)
373// B = 6(a - 2b + c)
374// C = 3(b - a)
375// Solve for t, keeping only those that fit between 0 < t < 1
376pub fn find_extrema(src: &[f64], t_values: &mut [f64]) -> usize {
377 // we divide A,B,C by 3 to simplify
378 let a: f64 = src[0];
379 let b: f64 = src[2];
380 let c: f64 = src[4];
381 let d: f64 = src[6];
382 let a2: f64 = d - a + 3.0 * (b - c);
383 let b2: f64 = 2.0 * (a - b - b + c);
384 let c2: f64 = b - a;
385
386 quad64::roots_valid_t(a:a2, b:b2, c:c2, t_values)
387}
388
389// Skia doesn't seems to care about NaN/inf during sorting, so we don't too.
390fn cmp_f64(a: &f64, b: &f64) -> core::cmp::Ordering {
391 if a < b {
392 core::cmp::Ordering::Less
393 } else if a > b {
394 core::cmp::Ordering::Greater
395 } else {
396 core::cmp::Ordering::Equal
397 }
398}
399
400// classic one t subdivision
401fn interp_cubic_coords_x(src: &[Point64; 4], t: f64, dst: &mut [Point64; 7]) {
402 use super::interp;
403
404 let ab: f64 = interp(a:src[0].x, b:src[1].x, t);
405 let bc: f64 = interp(a:src[1].x, b:src[2].x, t);
406 let cd: f64 = interp(a:src[2].x, b:src[3].x, t);
407 let abc: f64 = interp(a:ab, b:bc, t);
408 let bcd: f64 = interp(a:bc, b:cd, t);
409 let abcd: f64 = interp(a:abc, b:bcd, t);
410
411 dst[0].x = src[0].x;
412 dst[1].x = ab;
413 dst[2].x = abc;
414 dst[3].x = abcd;
415 dst[4].x = bcd;
416 dst[5].x = cd;
417 dst[6].x = src[3].x;
418}
419
420fn interp_cubic_coords_y(src: &[Point64; 4], t: f64, dst: &mut [Point64; 7]) {
421 use super::interp;
422
423 let ab: f64 = interp(a:src[0].y, b:src[1].y, t);
424 let bc: f64 = interp(a:src[1].y, b:src[2].y, t);
425 let cd: f64 = interp(a:src[2].y, b:src[3].y, t);
426 let abc: f64 = interp(a:ab, b:bc, t);
427 let bcd: f64 = interp(a:bc, b:cd, t);
428 let abcd: f64 = interp(a:abc, b:bcd, t);
429
430 dst[0].y = src[0].y;
431 dst[1].y = ab;
432 dst[2].y = abc;
433 dst[3].y = abcd;
434 dst[4].y = bcd;
435 dst[5].y = cd;
436 dst[6].y = src[3].y;
437}
438