alpha2_geometry/
arclen.rs

1//! Arc lengths of path segments, from `+ - * /` and `sqrt` alone.
2//!
3//! kurbo's arc length calls `powi` and `powf`, which are the platform's
4//! math library when kurbo's `std` is on; these integrate the speed with a
5//! fixed 16-point Gauss-Legendre rule, halving a segment while the halves
6//! disagree with the whole, so the result is the same bits everywhere.
7
8use kurbo::{ParamCurve, ParamCurveDeriv, PathSeg, Vec2};
9
10/// How close an arc length is computed, in pixels: the halving stops once
11/// the halves agree with the whole within this.
12pub const ACCURACY: f64 = 1e-3;
13
14/// Halvings at most, so a pathological segment still ends.
15const MAX_DEPTH: u32 = 12;
16
17/// The 16-point Gauss-Legendre rule on `[-1, 1]`: positive abscissas and
18/// their weights (the rule is symmetric).
19const GAUSS_16: [(f64, f64); 8] = [
20    (0.095_012_509_837_637_44, 0.189_450_610_455_068_5),
21    (0.281_603_550_779_258_9, 0.182_603_415_044_923_6),
22    (0.458_016_777_657_227_4, 0.169_156_519_395_002_5),
23    (0.617_876_244_402_643_8, 0.149_595_988_816_576_7),
24    (0.755_404_408_355_003, 0.124_628_971_255_533_9),
25    (0.865_631_202_387_831_7, 0.095_158_511_682_492_78),
26    (0.944_575_023_073_232_6, 0.062_253_523_938_647_89),
27    (0.989_400_934_991_649_9, 0.027_152_459_411_754_09),
28];
29
30fn speed(derivative: Vec2) -> f64 {
31    crate::affine::length(derivative)
32}
33
34/// The speed of `segment` at parameter `t`.
35fn segment_speed(segment: &PathSeg, t: f64) -> f64 {
36    match segment {
37        PathSeg::Line(line) => speed(line.p1 - line.p0),
38        PathSeg::Quad(quad) => speed(quad.deriv().eval(t).to_vec2()),
39        PathSeg::Cubic(cubic) => speed(cubic.deriv().eval(t).to_vec2()),
40    }
41}
42
43fn gauss(segment: &PathSeg, t0: f64, t1: f64) -> f64 {
44    let half = 0.5 * (t1 - t0);
45    let middle = f64::midpoint(t0, t1);
46    let mut sum = 0.0;
47    for (x, w) in GAUSS_16 {
48        sum += w
49            * (segment_speed(segment, middle - half * x)
50                + segment_speed(segment, middle + half * x));
51    }
52    sum * half
53}
54
55fn adaptive(segment: &PathSeg, t0: f64, t1: f64, whole: f64, depth: u32) -> f64 {
56    let middle = f64::midpoint(t0, t1);
57    let left = gauss(segment, t0, middle);
58    let right = gauss(segment, middle, t1);
59    if depth >= MAX_DEPTH || (left + right - whole).abs() <= ACCURACY * 0.5 {
60        return left + right;
61    }
62    adaptive(segment, t0, middle, left, depth + 1) + adaptive(segment, middle, t1, right, depth + 1)
63}
64
65/// The length of `segment` between parameters `t0` and `t1`.
66#[must_use]
67pub fn length_between(segment: &PathSeg, t0: f64, t1: f64) -> f64 {
68    if t1 <= t0 {
69        return 0.0;
70    }
71    if let PathSeg::Line(line) = segment {
72        return speed(line.p1 - line.p0) * (t1 - t0);
73    }
74    let whole = gauss(segment, t0, t1);
75    adaptive(segment, t0, t1, whole, 0)
76}
77
78/// The length of the whole of `segment`.
79#[must_use]
80pub fn length(segment: &PathSeg) -> f64 {
81    length_between(segment, 0.0, 1.0)
82}
83
84/// The parameter at which `segment` has run `distance` of its `total`
85/// length (which [`length`] gave), within [`ACCURACY`]; 0 or 1 past its
86/// ends.
87#[must_use]
88pub fn parameter_at(segment: &PathSeg, distance: f64, total: f64) -> f64 {
89    if distance <= 0.0 || total <= 0.0 {
90        return 0.0;
91    }
92    if distance >= total {
93        return 1.0;
94    }
95    if let PathSeg::Line(_) = segment {
96        return distance / total;
97    }
98    // Bisection on the length so far: monotonic, so it cannot leave the bracket.
99    let (mut low, mut high) = (0.0_f64, 1.0_f64);
100    let mut t = distance / total;
101    for _ in 0..60 {
102        let run = length_between(segment, 0.0, t);
103        let error = run - distance;
104        if error.abs() <= ACCURACY * 0.5 {
105            break;
106        }
107        if error < 0.0 {
108            low = t;
109        } else {
110            high = t;
111        }
112        t = f64::midpoint(low, high);
113    }
114    t
115}
116
117/// The point at parameter `t`, as kurbo evaluates it (Bernstein form, `+ - *` only).
118#[must_use]
119pub fn point_at(segment: &PathSeg, t: f64) -> kurbo::Point {
120    segment.eval(t)
121}
alpha2_geometry/arclen.rs

alpha2_geometry/
arclen.rs
