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;
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)]
47enum 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)`.
74fn 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
84fn 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        Some(Self {
166            local,
167            placement,
168            stretch,
169            origin,
170            pieces,
171            step,
172            heading_rounding,
173            at: 0.0,
174            sum: [(0.0, 0.0); 2],
175            remainder: 0.0,
176        })
177    }
178
179    /// March the local integral to `s`.
180    fn march(&mut self, s: Scalar) -> Option<()> {
181        if !(s.is_finite() && (0.0..=self.local.length).contains(&s)) {
182            return None;
183        }
184        if s < self.at {
185            self.at = 0.0;
186            self.sum = [(0.0, 0.0); 2];
187            self.remainder = 0.0;
188        }
189        for &(lo, hi, bend) in &self.pieces {
190            let (a, b) = (lo.max(self.at), hi.min(s));
191            if b <= a {
192                continue;
193            }
194            let n = ((b - a) / self.step).ceil().max(1.0);
195            if n > MAX_POSITION_PANELS as Scalar {
196                return None;
197            }
198            let w = (b - a) / n;
199            for i in 0..n as usize {
200                let theta = self.local.heading_at(a + (i as Scalar + 0.5) * w)?;
201                add(&mut self.sum[0], w * theta.cos());
202                add(&mut self.sum[1], w * theta.sin());
203            }
204            self.remainder += n * w * w * w / 24.0 * bend;
205        }
206        self.at = s;
207        Some(())
208    }
209
210    /// The marched point in the start frame's coordinates, and its error
211    /// there.
212    fn local_point(&self) -> (Point2, Scalar) {
213        let p = Point2::new(self.sum[0].0 + self.sum[0].1, self.sum[1].0 + self.sum[1].1);
214        let error = self.remainder + self.at * (self.heading_rounding + 16.0 * Scalar::EPSILON);
215        (p, error)
216    }
217
218    /// A local error through the start frame, with the placement's rounding.
219    fn placed_error(&self, local: Scalar) -> Scalar {
220        inflate(
221            self.stretch * local + 8.0 * Scalar::EPSILON * (self.origin + self.stretch * self.at),
222        )
223    }
224
225    /// The point at arc length `s` of a 2D curve and a certified bound on
226    /// its distance from the exact point; `None` outside `[0, length]`, for
227    /// a law that does not integrate there, or for a space curve.
228    pub fn point2(&mut self, s: Scalar) -> Option<(Point2, Scalar)> {
229        let Placement::Plane(frame) = self.placement else {
230            return None;
231        };
232        self.march(s)?;
233        let (p, local) = self.local_point();
234        let out = Point2::new(
235            frame.origin.x + frame.x.x * p.x + frame.y.x * p.y,
236            frame.origin.y + frame.x.y * p.x + frame.y.y * p.y,
237        );
238        Some((out, self.placed_error(local)))
239    }
240
241    /// [`Self::point2`] for a planar space curve.
242    pub fn point3(&mut self, s: Scalar) -> Option<(Point3, Scalar)> {
243        let Placement::Space(frame) = self.placement else {
244            return None;
245        };
246        self.march(s)?;
247        let (p, local) = self.local_point();
248        Some((
249            frame.origin + frame.x * p.x + frame.y * p.y,
250            self.placed_error(local),
251        ))
252    }
253
254    /// The point at arc length `s` of the 2D curve's offset by `distance`
255    /// to its left (the curve of [`super::offset_point2`]) and a certified
256    /// bound on its distance from the exact one: the point's, plus `|d|`
257    /// times the turn of the left normal by the heading's rounding through
258    /// the start frame (its stretch squared over its determinant bounds its
259    /// condition) and the normalisation's rounding.
260    pub fn offset_point2(&mut self, s: Scalar, distance: Scalar) -> Option<(Point2, Scalar)> {
261        let Placement::Plane(frame) = self.placement else {
262            return None;
263        };
264        let (p, error) = self.point2(s)?;
265        let theta = self.local.heading_at(s)?;
266        let (c, n) = (theta.cos(), theta.sin());
267        let t = Vec2::new(frame.x.x * c + frame.y.x * n, frame.x.y * c + frame.y.y * n);
268        let size = t.length();
269        let det = (frame.x.x * frame.y.y - frame.x.y * frame.y.x).abs();
270        if !(size.is_finite() && size > 0.0 && det > 0.0 && distance.is_finite()) {
271            return None;
272        }
273        let left = Vec2::new(-t.y, t.x) / size;
274        let condition = self.stretch * self.stretch / det;
275        let turn = condition * (self.heading_rounding + 16.0 * Scalar::EPSILON);
276        Some((
277            p + left * distance,
278            inflate(error + distance.abs() * (turn + 8.0 * Scalar::EPSILON)),
279        ))
280    }
281}
axiolid_evaluate/bound/positions.rs

axiolid_evaluate/bound/
positions.rs
