tiny_skia/path64/
cubic64.rs

Help
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 = src[6]; // d
221    let mut b = src[4] * 3.0; // 3*c
222    let mut c = src[2] * 3.0; // 3*b
223    let d = 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 = src[0];
379    let b = src[2];
380    let c = src[4];
381    let d = src[6];
382    let a2 = d - a + 3.0 * (b - c);
383    let b2 = 2.0 * (a - b - b + c);
384    let c2 = b - a;
385
386    quad64::roots_valid_t(a2, b2, 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 = interp(src[0].x, src[1].x, t);
405    let bc = interp(src[1].x, src[2].x, t);
406    let cd = interp(src[2].x, src[3].x, t);
407    let abc = interp(ab, bc, t);
408    let bcd = interp(bc, cd, t);
409    let abcd = interp(abc, 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 = interp(src[0].y, src[1].y, t);
424    let bc = interp(src[1].y, src[2].y, t);
425    let cd = interp(src[2].y, src[3].y, t);
426    let abc = interp(ab, bc, t);
427    let bcd = interp(bc, cd, t);
428    let abcd = interp(abc, 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}
tiny_skia/path64/cubic64.rs

tiny_skia/path64/
cubic64.rs
