axiolid_evaluate/bound/
position_table.rs

1//! A certified, random-access position table of intrinsic curves (#361).
2//!
3//! [`IntrinsicPositions`](super::IntrinsicPositions) (#317) marches the
4//! midpoint rule forward from the last point it read: its answers depend on
5//! the order of the queries, a backward read starts again from zero, and
6//! its `O(h^2)` rule needs about a million panels for a 10 km alignment at
7//! a thousandth of a millimetre. [`IntrinsicPositionTable`] reads any point
8//! of the same curve in any order, in at most `CHECKPOINT_PANELS` (64)
9//! panels of a fourth-order rule, and bounds the error of every point it
10//! returns.
11//!
12//! # The rule
13//!
14//! A point is the integral of the unit tangent `f = e^(i theta)` in the
15//! curve's start frame, `theta' = k`, `theta'' = k'`, `theta''' = k''`,
16//! `theta'''' = k'''` the closed-form derivatives of the curvature law. By
17//! the chain rule
18//!
19//! `f'''' = (i theta'''' - 4 theta' theta''' - 3 theta''^2 - 6 i theta'^2 theta'' + theta'^4) f`,
20//!
21//! so on a seam-free piece of the law
22//!
23//! `|f''''| <= B4 = sup|k'''| + 4 sup|k| sup|k''| + 3 sup k'^2 + 6 sup k^2 sup|k'| + sup k^4`,
24//!
25//! each supremum from the law's closed-form bound term by term
26//! (`crate::chain::curvature_bound`). The two-point Gauss-Legendre rule on
27//! a panel of width `w` has the error `w^5 / 4320 f''''(xi)` for a real
28//! integrand. Its Peano kernel keeps one sign, so the vector form
29//! `|E| <= w^5 / 4320 B4` holds without a factor `sqrt 2`.
30//!
31//! # The grid
32//!
33//! Each seam-free piece of `[0, end]` (the law may jump at a seam, so no
34//! panel straddles one) is cut into equal panels no wider than `step`,
35//! chosen so the remainders over the whole curve sum to half the requested
36//! error. Every 64 panels the table stores the running sums and the
37//! remainder so far. A point at `s` starts from the last checkpoint at or
38//! before the panel holding `s`, sums the full panels up to it, and adds
39//! one partial panel from the grid line to `s` by the same rule, with its
40//! own remainder `(s - grid)^5 / 4320 B4`. A query is a pure function of
41//! `s`.
42//!
43//! The rounding terms are those of the midpoint positions: the heading's,
44//! at most `64 eps (1 + K L) (1 + seams)` radians (`K` bounding `|k|`, `L`
45//! the end) plus the harmonics' own, `8 eps |A/w| (2 + |p| + |w| L)` each
46//! (`crate::chain::harmonic_rounding`; the cosine difference of a slow
47//! harmonic cancels, so its rounding is not proportional to the turning),
48//! which turns each tangent by as much so a point at `s` moves by
49//! at most `s` times it; the compensated sums' `16 eps s`; and the start
50//! frame's stretch, applied to all of these, with `8 eps (|origin| +
51//! stretch s)` for placing the point.
52//!
53//! [`IntrinsicPositionTable::max_error`] is one uniform bound for the whole
54//! table: all the remainders plus the roundings at the end (#393).
55
56use axiolid_core::{Point2, Point3, Scalar};
57use axiolid_curve::{Intrinsic2, Intrinsic3};
58
59use super::positions::{add, unit_frame, Placement};
60use super::{frame_stretch2, frame_stretch3, inflate};
61use crate::chain::{curvature_bound, harmonic_rounding};
62
63/// Most panels one table may use.
64pub const MAX_TABLE_PANELS: usize = 1 << 22;
65
66/// Panels between two stored checkpoints, and so most full panels a query
67/// sums.
68const CHECKPOINT_PANELS: usize = 64;
69
70/// How far past the curve's length an end may lie, relative to `max(1,
71/// length)`, and still be read as the length.
72const END_SLACK: Scalar = 1e-9;
73
74/// The running sums at a panel boundary.
75#[derive(Debug, Clone, Copy)]
76struct Checkpoint {
77    sum: [(Scalar, Scalar); 2],
78    remainder: Scalar,
79}
80
81/// A seam-free piece of `[0, end]`, cut into equal panels.
82#[derive(Debug, Clone, Copy)]
83struct Span {
84    lo: Scalar,
85    hi: Scalar,
86    width: Scalar,
87    panels: usize,
88    /// The index of its first panel in the whole grid.
89    base: usize,
90    /// `sup |f''''|` over it.
91    b4: Scalar,
92}
93
94/// Certified positions of an intrinsic curve over `[0, end]`, read in any
95/// order: the unit tangent of its closed-form heading integrated by the
96/// composite two-point Gauss-Legendre rule, each point with a certified
97/// bound on its error (see the module documentation).
98#[derive(Debug, Clone)]
99pub struct IntrinsicPositionTable {
100    /// The curve in its own start frame.
101    local: Intrinsic2,
102    placement: Placement,
103    stretch: Scalar,
104    origin: Scalar,
105    end: Scalar,
106    spans: Vec<Span>,
107    checkpoints: Vec<Checkpoint>,
108    heading_rounding: Scalar,
109    panels: usize,
110    /// The remainders of every panel together.
111    total_remainder: Scalar,
112}
113
114/// The remainder of one panel of width `w`.
115fn remainder(width: Scalar, b4: Scalar) -> Scalar {
116    width.powi(5) / 4320.0 * b4
117}
118
119impl IntrinsicPositionTable {
120    /// The table of a 2D intrinsic curve over `[0, end]`, its points
121    /// within about `error`, half of it spent on the remainders; `None`
122    /// for an error or an end that is not finite and positive, an end past
123    /// the curve's length, a law with no closed-form bound of `k` to `k'''`
124    /// on a seam-free piece, or more than [`MAX_TABLE_PANELS`] panels.
125    #[must_use]
126    pub fn new2(curve: &Intrinsic2, error: Scalar, end: Scalar) -> Option<Self> {
127        let local = Intrinsic2::new(unit_frame(), curve.curvature.clone(), curve.length);
128        Self::build(
129            local,
130            Placement::Plane(curve.start),
131            frame_stretch2(&curve.start),
132            curve.start.origin.abs().max_element(),
133            error,
134            end,
135        )
136    }
137
138    /// The table of a planar space curve (its torsion law identically zero,
139    /// [`Intrinsic3::is_planar`]); `None` as for [`Self::new2`], or for a
140    /// curve with torsion.
141    #[must_use]
142    pub fn new3(curve: &Intrinsic3, error: Scalar, end: Scalar) -> Option<Self> {
143        if !curve.is_planar() {
144            return None;
145        }
146        let local = Intrinsic2::new(unit_frame(), curve.curvature.clone(), curve.length);
147        Self::build(
148            local,
149            Placement::Space(curve.start),
150            frame_stretch3(&curve.start),
151            curve.start.origin.abs().max_element(),
152            error,
153            end,
154        )
155    }
156
157    fn build(
158        local: Intrinsic2,
159        placement: Placement,
160        stretch: Scalar,
161        origin: Scalar,
162        error: Scalar,
163        end: Scalar,
164    ) -> Option<Self> {
165        let length = local.length;
166        if !(error.is_finite() && error > 0.0 && length.is_finite() && length > 0.0) {
167            return None;
168        }
169        if !(end.is_finite() && end > 0.0 && end <= length + END_SLACK * length.max(1.0)) {
170            return None;
171        }
172        let end = end.min(length);
173        let law = &local.curvature;
174        let rate = law.derivative();
175        let second = rate.derivative();
176        let third = second.derivative();
177        let seams: Vec<Scalar> = law
178            .seams_within(end)
179            .into_iter()
180            .filter(|&s| s > 0.0 && s < end)
181            .collect();
182        let mut ends = vec![0.0];
183        ends.extend(seams.iter().copied());
184        ends.push(end);
185        let mut pieces = Vec::with_capacity(ends.len() - 1);
186        let mut worst: Scalar = 0.0;
187        for pair in ends.windows(2) {
188            let (lo, hi) = (pair[0], pair[1]);
189            let k = curvature_bound(law, lo, hi)?;
190            let k1 = curvature_bound(&rate, lo, hi)?;
191            let k2 = curvature_bound(&second, lo, hi)?;
192            let k3 = curvature_bound(&third, lo, hi)?;
193            let b4 = inflate(k3 + 4.0 * k * k2 + 3.0 * k1 * k1 + 6.0 * k * k * k1 + k.powi(4));
194            if !b4.is_finite() {
195                return None;
196            }
197            worst = worst.max(b4);
198            pieces.push((lo, hi, b4));
199        }
200        let budget = 0.5 * error / stretch;
201        let step = if worst > 0.0 {
202            (4320.0 * budget / (end * worst)).powf(0.25).min(end)
203        } else {
204            end
205        };
206        if !(step.is_finite() && step > 0.0) {
207            return None;
208        }
209        let mut spans = Vec::with_capacity(pieces.len());
210        let mut total = 0usize;
211        for (lo, hi, b4) in pieces {
212            let count = ((hi - lo) / step).ceil().max(1.0);
213            if count > MAX_TABLE_PANELS as Scalar || total + count as usize > MAX_TABLE_PANELS {
214                return None;
215            }
216            let panels = count as usize;
217            spans.push(Span {
218                lo,
219                hi,
220                width: (hi - lo) / count,
221                panels,
222                base: total,
223                b4,
224            });
225            total += panels;
226        }
227        let turning = curvature_bound(law, 0.0, end)? * end;
228        let heading_rounding =
229            64.0 * Scalar::EPSILON * (1.0 + turning) * (1.0 + seams.len() as Scalar)
230                + harmonic_rounding(law, end)?;
231        let mut table = Self {
232            local,
233            placement,
234            stretch,
235            origin,
236            end,
237            spans,
238            checkpoints: Vec::with_capacity(total / CHECKPOINT_PANELS + 1),
239            heading_rounding,
240            panels: total,
241            total_remainder: 0.0,
242        };
243        let mut state = Checkpoint {
244            sum: [(0.0, 0.0); 2],
245            remainder: 0.0,
246        };
247        for index in 0..total {
248            if index % CHECKPOINT_PANELS == 0 {
249                table.checkpoints.push(state);
250            }
251            table.add_panel(index, &mut state)?;
252        }
253        table.total_remainder = state.remainder;
254        Some(table)
255    }
256
257    /// The number of panels the table is cut into.
258    #[must_use]
259    pub fn panels(&self) -> usize {
260        self.panels
261    }
262
263    /// The span holding grid panel `index`.
264    fn span_of(&self, index: usize) -> &Span {
265        let at = self.spans.partition_point(|s| s.base + s.panels <= index);
266        &self.spans[at.min(self.spans.len() - 1)]
267    }
268
269    /// The start of panel `i` of `span`.
270    fn grid(span: &Span, i: usize) -> Scalar {
271        if i == 0 {
272            span.lo
273        } else {
274            span.lo + i as Scalar * span.width
275        }
276    }
277
278    /// Add the integral of `f` over `[a, a + w]` by two Gauss points to
279    /// `sum`.
280    fn rule(&self, a: Scalar, w: Scalar, sum: &mut [(Scalar, Scalar); 2]) -> Option<()> {
281        let half = 0.5 * w;
282        let offset = half / (3.0 as Scalar).sqrt();
283        for x in [a + half - offset, a + half + offset] {
284            let theta = self.local.heading_at(x)?;
285            add(&mut sum[0], half * theta.cos());
286            add(&mut sum[1], half * theta.sin());
287        }
288        Some(())
289    }
290
291    /// Add grid panel `index` to `state`.
292    fn add_panel(&self, index: usize, state: &mut Checkpoint) -> Option<()> {
293        let span = self.span_of(index);
294        let a = Self::grid(span, index - span.base);
295        self.rule(a, span.width, &mut state.sum)?;
296        state.remainder += remainder(span.width, span.b4);
297        Some(())
298    }
299
300    /// The point at `s` in the start frame's coordinates, and its error
301    /// there.
302    fn local_point(&self, s: Scalar) -> Option<(Point2, Scalar)> {
303        if !(s.is_finite() && (0.0..=self.end).contains(&s)) {
304            return None;
305        }
306        let at = self.spans.partition_point(|span| span.hi < s);
307        let span = self.spans.get(at)?;
308        let i = (((s - span.lo) / span.width).floor().max(0.0) as usize).min(span.panels - 1);
309        let index = span.base + i;
310        let first = index - index % CHECKPOINT_PANELS;
311        let mut state = self.checkpoints[index / CHECKPOINT_PANELS];
312        for panel in first..index {
313            self.add_panel(panel, &mut state)?;
314        }
315        let from = Self::grid(span, i);
316        let rest = (s - from).max(0.0);
317        if rest > 0.0 {
318            self.rule(from, rest, &mut state.sum)?;
319            state.remainder += remainder(rest, span.b4);
320        }
321        let p = Point2::new(
322            state.sum[0].0 + state.sum[0].1,
323            state.sum[1].0 + state.sum[1].1,
324        );
325        let error = state.remainder + s * (self.heading_rounding + 16.0 * Scalar::EPSILON);
326        Some((p, error))
327    }
328
329    /// A local error through the start frame, with the placement's
330    /// rounding.
331    fn placed_error(&self, local: Scalar, s: Scalar) -> Scalar {
332        inflate(self.stretch * local + 8.0 * Scalar::EPSILON * (self.origin + self.stretch * s))
333    }
334
335    /// The point at arc length `s` of a 2D curve and a certified bound on
336    /// its distance from the exact point: the same bits whatever was read
337    /// before. `None` outside `[0, end]`, for a law that does not integrate
338    /// there, or for a space curve.
339    #[must_use]
340    pub fn point2(&self, s: Scalar) -> Option<(Point2, Scalar)> {
341        let Placement::Plane(frame) = self.placement else {
342            return None;
343        };
344        let (p, local) = self.local_point(s)?;
345        let out = Point2::new(
346            frame.origin.x + frame.x.x * p.x + frame.y.x * p.y,
347            frame.origin.y + frame.x.y * p.x + frame.y.y * p.y,
348        );
349        Some((out, self.placed_error(local, s)))
350    }
351
352    /// [`Self::point2`] for a planar space curve.
353    #[must_use]
354    pub fn point3(&self, s: Scalar) -> Option<(Point3, Scalar)> {
355        let Placement::Space(frame) = self.placement else {
356            return None;
357        };
358        let (p, local) = self.local_point(s)?;
359        Some((
360            frame.origin + frame.x * p.x + frame.y * p.y,
361            self.placed_error(local, s),
362        ))
363    }
364
365    /// A certified bound on the error of every point the table returns
366    /// (the error of [`Self::point2`] and [`Self::point3`] grows with `s`
367    /// and is never above it): the remainders of all panels and the
368    /// roundings of the end, through the start frame.
369    #[must_use]
370    pub fn max_error(&self) -> Scalar {
371        let local = inflate(
372            self.total_remainder + self.end * (self.heading_rounding + 16.0 * Scalar::EPSILON),
373        );
374        self.placed_error(local, self.end)
375    }
376
377    /// A certified bound on the rounding of the heading, in radians: what a
378    /// tangent read from it may be turned by.
379    #[must_use]
380    pub fn heading_rounding(&self) -> Scalar {
381        self.heading_rounding
382    }
383}
axiolid_evaluate/bound/position_table.rs

axiolid_evaluate/bound/
position_table.rs
