axiolid_evaluate/bound/
positions.rs

1//! Certified positions of intrinsic curves (#317).
2//!
3//! An intrinsic curve's point is the integral of its unit tangent `f =
4//! (cos theta, sin theta)` in the start frame, `theta` the law's closed-form
5//! heading. [`IntrinsicPositions`] integrates it by the composite midpoint
6//! rule and bounds the error of every point it returns, so a polyline
7//! through such points can carry a bound rather than the evaluator's
8//! quadrature estimate (`crate::arc_length`).
9//!
10//! # The remainder
11//!
12//! On a panel of width `w` with midpoint `m`, `|int f - w f(m)| <= w^3 / 24
13//! sup |f''|` (the vector Taylor remainder). `f'' = theta'' f^perp -
14//! theta'^2 f`, so `|f''| <= |k'| + k^2`, bounded per seam-free piece of
15//! the law from its closed-form suprema (`crate::chain::curvature_bound`),
16//! as the certified station lengths bound their speeds (#294). Panels never
17//! straddle a seam, where `k` may jump. The panel width is chosen so the
18//! remainders over the whole curve sum to half the requested error; each
19//! point reports the remainders actually summed up to it, plus:
20//!
21//! - the heading's rounding, at most `64 eps (1 + K L) (1 + seams)` radians
22//!   (`K` bounding `|k|`, `L` the length), which turns each tangent by as
23//!   much, so a point at `s` moves by at most `s` times it;
24//! - the sums' rounding, `16 eps s` (compensated sums);
25//! - the start frame's stretch `sigma`, applied to all of these, and
26//!   `8 eps (|origin| + sigma s)` for placing the point.
27//!
28//! A planar space curve (torsion identically zero) is the same integral in
29//! its start frame's `x`, `y` plane: its tangent is `cos theta x + sin
30//! theta y`.
31//!
32//! Points are marched forward from the last one read, so reading them in
33//! increasing arc length costs one pass; reading backwards starts again.
34
35use axiolid_core::{Frame2, Frame3, Point2, Point3, Scalar, Vec2};
36use axiolid_curve::{Intrinsic2, Intrinsic3};
37
38use super::{frame_stretch2, frame_stretch3, inflate};
39use crate::chain::{curvature_bound, harmonic_rounding};
40
41/// Most midpoint panels one curve's positions may use.
42pub const MAX_POSITION_PANELS: usize = 1 << 22;
43
44/// Where a marched point is placed: by a 2D curve's start frame, or in the
45/// plane of a planar space curve's start frame.
46#[derive(Debug, Clone, Copy)]
47pub(super) enum Placement {
48    Plane(Frame2),
49    Space(Frame3),
50}
51
52/// Certified positions of an intrinsic curve: the unit tangent of its
53/// closed-form heading integrated by the composite midpoint rule, each
54/// point with a certified bound on its error (the remainder `w^3 / 24
55/// (|k'| + k^2)` per panel, panels split at seams, and the heading's, the
56/// sums' and the placement's rounding).
57#[derive(Debug, Clone)]
58pub struct IntrinsicPositions {
59    /// The curve in its own start frame.
60    local: Intrinsic2,
61    placement: Placement,
62    stretch: Scalar,
63    origin: Scalar,
64    /// Seam-free pieces `(lo, hi, sup |f''|)`.
65    pieces: Vec<(Scalar, Scalar, Scalar)>,
66    step: Scalar,
67    heading_rounding: Scalar,
68    at: Scalar,
69    sum: [(Scalar, Scalar); 2],
70    remainder: Scalar,
71}
72
73/// Neumaier's compensated addition into `(sum, carry)`.
74pub(super) fn add(acc: &mut (Scalar, Scalar), value: Scalar) {
75    let total = acc.0 + value;
76    acc.1 += if acc.0.abs() >= value.abs() {
77        (acc.0 - total) + value
78    } else {
79        (value - total) + acc.0
80    };
81    acc.0 = total;
82}
83
84pub(super) fn unit_frame() -> Frame2 {
85    Frame2 {
86        origin: Point2::ZERO,
87        x: Vec2::X,
88        y: Vec2::Y,
89    }
90}
91
92impl IntrinsicPositions {
93    /// Positions of a 2D intrinsic curve within about `error`, half of it
94    /// spent on the midpoint remainders; `None` for a non-positive or
95    /// non-finite error or length, a law with no closed-form bound, or more
96    /// than [`MAX_POSITION_PANELS`] panels.
97    #[must_use]
98    pub fn new2(curve: &Intrinsic2, error: Scalar) -> Option<Self> {
99        let local = Intrinsic2::new(unit_frame(), curve.curvature.clone(), curve.length);
100        Self::build(
101            local,
102            Placement::Plane(curve.start),
103            frame_stretch2(&curve.start),
104            curve.start.origin.abs().max_element(),
105            error,
106        )
107    }
108
109    /// Positions of a planar space curve (its torsion law identically
110    /// zero, [`Intrinsic3::is_planar`]); `None` as for [`Self::new2`], or
111    /// for a curve with torsion.
112    #[must_use]
113    pub fn new3(curve: &Intrinsic3, error: Scalar) -> Option<Self> {
114        if !curve.is_planar() {
115            return None;
116        }
117        let local = Intrinsic2::new(unit_frame(), curve.curvature.clone(), curve.length);
118        Self::build(
119            local,
120            Placement::Space(curve.start),
121            frame_stretch3(&curve.start),
122            curve.start.origin.abs().max_element(),
123            error,
124        )
125    }
126
127    fn build(
128        local: Intrinsic2,
129        placement: Placement,
130        stretch: Scalar,
131        origin: Scalar,
132        error: Scalar,
133    ) -> Option<Self> {
134        let length = local.length;
135        if !(error.is_finite() && error > 0.0 && length.is_finite() && length > 0.0) {
136            return None;
137        }
138        let law = &local.curvature;
139        let rate = law.derivative();
140        let seams = law.seams_within(length);
141        let mut ends = vec![0.0];
142        ends.extend(seams.iter().copied());
143        ends.push(length);
144        let mut pieces = Vec::with_capacity(ends.len() - 1);
145        let mut bend: Scalar = 0.0;
146        for pair in ends.windows(2) {
147            let (lo, hi) = (pair[0], pair[1]);
148            let k = curvature_bound(law, lo, hi)?;
149            let m = inflate(curvature_bound(&rate, lo, hi)? + k * k);
150            bend = bend.max(m);
151            pieces.push((lo, hi, m));
152        }
153        let budget = 0.5 * error / stretch;
154        let step = if bend > 0.0 {
155            (24.0 * budget / (length * bend)).sqrt().min(length)
156        } else {
157            length
158        };
159        if !(step.is_finite() && step > 0.0) || length / step > MAX_POSITION_PANELS as Scalar {
160            return None;
161        }
162        let turning = curvature_bound(law, 0.0, length)? * length;
163        let heading_rounding =
164            64.0 * Scalar::EPSILON * (1.0 + turning) * (1.0 + seams.len() as Scalar)
165                + harmonic_rounding(law, length)?;
166        Some(Self {
167            local,
168            placement,
169            stretch,
170            origin,
171            pieces,
172            step,
173            heading_rounding,
174            at: 0.0,
175            sum: [(0.0, 0.0); 2],
176            remainder: 0.0,
177        })
178    }
179
180    /// March the local integral to `s`.
181    fn march(&mut self, s: Scalar) -> Option<()> {
182        if !(s.is_finite() && (0.0..=self.local.length).contains(&s)) {
183            return None;
184        }
185        if s < self.at {
186            self.at = 0.0;
187            self.sum = [(0.0, 0.0); 2];
188            self.remainder = 0.0;
189        }
190        for &(lo, hi, bend) in &self.pieces {
191            let (a, b) = (lo.max(self.at), hi.min(s));
192            if b <= a {
193                continue;
194            }
195            let n = ((b - a) / self.step).ceil().max(1.0);
196            if n > MAX_POSITION_PANELS as Scalar {
197                return None;
198            }
199            let w = (b - a) / n;
200            for i in 0..n as usize {
201                let theta = self.local.heading_at(a + (i as Scalar + 0.5) * w)?;
202                add(&mut self.sum[0], w * theta.cos());
203                add(&mut self.sum[1], w * theta.sin());
204            }
205            self.remainder += n * w * w * w / 24.0 * bend;
206        }
207        self.at = s;
208        Some(())
209    }
210
211    /// The marched point in the start frame's coordinates, and its error
212    /// there.
213    fn local_point(&self) -> (Point2, Scalar) {
214        let p = Point2::new(self.sum[0].0 + self.sum[0].1, self.sum[1].0 + self.sum[1].1);
215        let error = self.remainder + self.at * (self.heading_rounding + 16.0 * Scalar::EPSILON);
216        (p, error)
217    }
218
219    /// A local error through the start frame, with the placement's rounding.
220    fn placed_error(&self, local: Scalar) -> Scalar {
221        inflate(
222            self.stretch * local + 8.0 * Scalar::EPSILON * (self.origin + self.stretch * self.at),
223        )
224    }
225
226    /// The point at arc length `s` of a 2D curve and a certified bound on
227    /// its distance from the exact point; `None` outside `[0, length]`, for
228    /// a law that does not integrate there, or for a space curve.
229    pub fn point2(&mut self, s: Scalar) -> Option<(Point2, Scalar)> {
230        let Placement::Plane(frame) = self.placement else {
231            return None;
232        };
233        self.march(s)?;
234        let (p, local) = self.local_point();
235        let out = Point2::new(
236            frame.origin.x + frame.x.x * p.x + frame.y.x * p.y,
237            frame.origin.y + frame.x.y * p.x + frame.y.y * p.y,
238        );
239        Some((out, self.placed_error(local)))
240    }
241
242    /// [`Self::point2`] for a planar space curve.
243    pub fn point3(&mut self, s: Scalar) -> Option<(Point3, Scalar)> {
244        let Placement::Space(frame) = self.placement else {
245            return None;
246        };
247        self.march(s)?;
248        let (p, local) = self.local_point();
249        Some((
250            frame.origin + frame.x * p.x + frame.y * p.y,
251            self.placed_error(local),
252        ))
253    }
254
255    /// The point at arc length `s` of the 2D curve's offset by `distance`
256    /// to its left (the curve of [`super::offset_point2`]) and a certified
257    /// bound on its distance from the exact one: the point's, plus `|d|`
258    /// times the turn of the left normal by the heading's rounding through
259    /// the start frame (its stretch squared over its determinant bounds its
260    /// condition) and the normalisation's rounding.
261    pub fn offset_point2(&mut self, s: Scalar, distance: Scalar) -> Option<(Point2, Scalar)> {
262        let Placement::Plane(frame) = self.placement else {
263            return None;
264        };
265        let (p, error) = self.point2(s)?;
266        let theta = self.local.heading_at(s)?;
267        let (c, n) = (theta.cos(), theta.sin());
268        let t = Vec2::new(frame.x.x * c + frame.y.x * n, frame.x.y * c + frame.y.y * n);
269        let size = t.length();
270        let det = (frame.x.x * frame.y.y - frame.x.y * frame.y.x).abs();
271        if !(size.is_finite() && size > 0.0 && det > 0.0 && distance.is_finite()) {
272            return None;
273        }
274        let left = Vec2::new(-t.y, t.x) / size;
275        let condition = self.stretch * self.stretch / det;
276        let turn = condition * (self.heading_rounding + 16.0 * Scalar::EPSILON);
277        Some((
278            p + left * distance,
279            inflate(error + distance.abs() * (turn + 8.0 * Scalar::EPSILON)),
280        ))
281    }
282}
axiolid_evaluate/bound/positions.rs

axiolid_evaluate/bound/
positions.rs
