axiolid_mesh_compile/station/
accuracy.rs

1//! How far a node placed at a station lies from its exact placement (#405,
2//! ADR 0082 amendment).
3//!
4//! [`placement`](super::placement) puts a source's local `x`, `y`, `z` on
5//! the oriented tangent, lateral and up of the frame the kernel reads at
6//! the station, its origin at the offsets read in the base frame. That
7//! frame is exact on a line ([`ResolvedPlacement::exact`]); on any other
8//! basis it is the exact frame law's at a distance read by an estimate, and
9//! [`placement_accuracy`] states how far the placement lies from the exact
10//! one: [`PlacementAccuracy::Exact`], [`PlacementAccuracy::Certified`] with
11//! a [`PlacementBound`], or [`PlacementAccuracy::Unbounded`] by name.
12//!
13//! # The exact placement
14//!
15//! The source's local `x`, `y`, `z` on the oriented axes of the exact
16//! frame law (`axiolid_reference::station::FrameLaw`, #316, #361) at the
17//! station's distance, snapped by its window (#294) and read from its seam
18//! side (#263), its origin at the offsets read in the exact base frame, a
19//! section's or the plan's. The kernel's frame is compared with it, never
20//! trusted: `FrameLaw::error_at` gives the error `e_p` of the read frame's
21//! point and `e_a` of each of its axes.
22//!
23//! # The bound
24//!
25//! With `L1 = |lateral| + |vertical| + |longitudinal|` of the offsets, the
26//! station's point `P = O + lateral L + vertical U + longitudinal T` is
27//! within `e_p + e_a L1` of the exact one, plus the rounding of computing
28//! it; the placement's axes are the base frame's, or its axes turned by
29//! the orientation, each a combination of the three with unit components
30//! `w`, `|w|_1 <= sqrt 3`: within `sqrt 3 e_a` of the exact ones, plus the
31//! rounding of turning them. A local point `q` then lands within
32//! `point_error + axis_error |q|_1` of its exact place
33//! ([`PlacementBound::displacement`]); the bound is the same for every
34//! point of a mesh, since `|q|_1` is convex. The rounding of applying the
35//! transform in `f64` is not part of it: the mesh compiler adds it as it
36//! adds it to every instance (#295).
37//!
38//! # Unbounded
39//!
40//! By name ([`UNCERTIFIED_STATION_FRAME`]), where `FrameLaw::error_at`
41//! certifies nothing: frames read by quadrature (an ellipse, a B-spline, an
42//! offset), chains, banked curves over intrinsic or chain plans, intrinsic
43//! space curves with torsion or tilted, elevated curves with an intrinsic
44//! profile or another plan, a skewed circle, a steep tangent, a basis
45//! placed at a station whose own frame is not exact, and a station within
46//! its distance's resolution of a seam it is not read on.
47
48use axiolid_contracts::{GeomError, GeomResult};
49use axiolid_core::{Point3, Scalar};
50use axiolid_model::{GeometryGraph, GeometryNode, NodeId, StationOffsets};
51use axiolid_reference::station::POINT_ROUNDING;
52
53use super::basis::Basis;
54use super::bound::reading as frame_reading;
55use super::resolve_parts;
56
57/// Why a placement whose frame has no certified error is unbounded.
58pub(crate) const UNCERTIFIED_STATION_FRAME: &str =
59    "a node placed at a station whose frame has no certified error (a frame read by quadrature \
60     on an ellipse, a B-spline or an offset, a chain, a banked curve over a chain \
61     plan, an intrinsic space curve with torsion or tilted, an elevated curve with an intrinsic \
62     profile or another plan, a skewed circle, a steep tangent, a basis placed at a station \
63     whose frame is not exact, or a station within its distance's resolution of a seam it is \
64     not read on)";
65
66/// The rounding of turning the base frame by an orientation, per axis, for
67/// an axis and reference direction at a right angle: turning is conditioned
68/// by `1 / |axis x reference|` (Gram-Schmidt of two nearly parallel
69/// vectors).
70const TURN_ROUNDING: Scalar = 1e-14;
71
72/// The margin on every composed term, over the rounding of composing it.
73const MARGIN: Scalar = 1e-12;
74
75/// How far a placement at a station lies from the exact one, see
76/// [placing at a station](super#placing-at-a-station-264).
77#[non_exhaustive]
78#[derive(Debug, Clone, Copy, PartialEq)]
79pub struct PlacementBound {
80    /// How far the placement's origin, the station's point, lies from the
81    /// exact one.
82    pub point_error: Scalar,
83    /// How far each of the placement's axes lies from the exact one.
84    pub axis_error: Scalar,
85}
86
87impl PlacementBound {
88    /// How far the image of the source's local point `local` may lie from
89    /// its exact image: `point_error + axis_error (|x| + |y| + |z|)`.
90    #[must_use]
91    pub fn displacement(&self, local: Point3) -> Scalar {
92        self.point_error + self.axis_error * (local.x.abs() + local.y.abs() + local.z.abs())
93    }
94}
95
96/// How far a node placed at a station lies from its exact placement.
97#[non_exhaustive]
98#[derive(Debug, Clone, Copy, PartialEq)]
99pub enum PlacementAccuracy {
100    /// The frame is exact, rounding aside (`ResolvedPlacement::exact`): a
101    /// line basis, itself placed, if at all, in exact frames.
102    Exact,
103    /// The placement lies within the bound of the exact one.
104    Certified(PlacementBound),
105    /// No bound is certified, for the reason named.
106    Unbounded(&'static str),
107}
108
109/// How far the placement [`placement`](super::placement) resolves for the
110/// [`InstanceAtStation`] node `id` lies from the exact one, its frame law's
111/// quadratures certified to `tolerance` (see
112/// [placing at a station](super#placing-at-a-station-264)).
113///
114/// [`PlacementAccuracy::Exact`] exactly when the placement is
115/// ([`ResolvedPlacement::exact`](super::ResolvedPlacement::exact)); a
116/// certified bound is finite and not negative. Both read the frame the
117/// placement reads, bit for bit.
118///
119/// # Errors
120///
121/// A node that is not an [`InstanceAtStation`], a `tolerance` that is not
122/// finite and positive, and every refusal of [`placement`](super::placement),
123/// each by name.
124///
125/// [`InstanceAtStation`]: axiolid_model::InstanceAtStation
126pub fn placement_accuracy(
127    graph: &GeometryGraph,
128    id: NodeId,
129    tolerance: Scalar,
130) -> GeomResult<PlacementAccuracy> {
131    let Some(GeometryNode::InstanceAtStation(placed)) = graph.get(id) else {
132        return Err(GeomError::InvalidInput(format!(
133            "node {id:?} is not an instance at a station"
134        )));
135    };
136    if !(tolerance.is_finite() && tolerance > 0.0) {
137        return Err(GeomError::InvalidInput(format!(
138            "station: a certified placement needs a finite, positive tolerance, not {tolerance}"
139        )));
140    }
141    let station = &placed.station;
142    let parts = resolve_parts(graph, station, 0)?;
143    if parts.exact {
144        return Ok(PlacementAccuracy::Exact);
145    }
146    let unbounded = PlacementAccuracy::Unbounded(UNCERTIFIED_STATION_FRAME);
147    // A curve placed at a station whose frame is not exact has no exact
148    // frame law to compare with.
149    if matches!(
150        parts.basis,
151        Basis::Atomic {
152            placed_exactly: false,
153            ..
154        }
155    ) {
156        return Ok(unbounded);
157    }
158    let Some(reading) = frame_reading(station.station.frame) else {
159        return Ok(unbounded);
160    };
161    let law = parts.basis.frame_law(tolerance)?;
162    let Some(error) = law.error_at(parts.distance, station.seam, reading) else {
163        return Ok(unbounded);
164    };
165    // `None` for the base frame, else the sine of the angle between the
166    // orientation's axis and reference direction.
167    let turn = if station.orientation.is_base() {
168        None
169    } else {
170        match station.orientation.unit_axes() {
171            Ok((axis, reference)) => Some(axis.cross(reference).length()),
172            Err(_) => return Ok(unbounded),
173        }
174    };
175    if turn.is_some_and(|sin| !(sin.is_finite() && sin > 0.0)) {
176        return Ok(unbounded);
177    }
178    let bound = compose(
179        (error.point, error.axis),
180        &station.station.station.offsets,
181        turn,
182        parts.resolved.point,
183    );
184    Ok(
185        if bound.point_error.is_finite() && bound.axis_error.is_finite() {
186            PlacementAccuracy::Certified(bound)
187        } else {
188            unbounded
189        },
190    )
191}
192
193/// The placement's bound from the errors `(read_point, read_axis)` of the
194/// frame read: `turn` the sine of the angle between the orientation's axis
195/// and reference direction, `None` when it leaves the base frame,
196/// `offsets` read in it and `point` the station's point.
197fn compose(
198    (read_point, read_axis): (Scalar, Scalar),
199    offsets: &StationOffsets,
200    turn: Option<Scalar>,
201    point: Point3,
202) -> PlacementBound {
203    let reach = offsets.lateral.abs() + offsets.vertical.abs() + offsets.longitudinal.abs();
204    let scale = (point.abs().max_element() + reach).max(1.0);
205    let point_error = (read_point + read_axis * reach + POINT_ROUNDING * scale) * (1.0 + MARGIN);
206    let axis_error = match turn {
207        None => read_axis,
208        // Each turned axis combines the base axes with unit components.
209        Some(sin) => 3.0_f64.sqrt() * read_axis + TURN_ROUNDING / sin,
210    } * (1.0 + MARGIN);
211    PlacementBound {
212        point_error,
213        axis_error,
214    }
215}
216
217#[cfg(test)]
218mod tests {
219    use super::*;
220    use axiolid_core::Vec3;
221
222    /// A frame read: its point and its tangent, lateral and up.
223    #[derive(Clone, Copy)]
224    struct Read {
225        point: Vec3,
226        axes: [Vec3; 3],
227    }
228
229    fn base() -> Read {
230        // Orthonormal, not axis-aligned.
231        let t = Vec3::new(0.6, 0.8, 0.0);
232        let l = Vec3::new(-0.8, 0.6, 0.0);
233        Read {
234            point: Vec3::new(1000.0, 2000.0, 30.0),
235            axes: [t, l, Vec3::Z],
236        }
237    }
238
239    /// The read frame off by `a` in its point and `b` in each axis, all in
240    /// the direction `n`: the worst case the bound allows.
241    fn perturbed(read: Read, a: Scalar, b: Scalar, n: Vec3) -> Read {
242        Read {
243            point: read.point + a * n,
244            axes: read.axes.map(|axis| axis + b * n),
245        }
246    }
247
248    fn offsets() -> StationOffsets {
249        // Positive, so the axes' errors add.
250        StationOffsets::new(2.0, 1.0, 0.5)
251    }
252
253    /// The station's point: the origin moved by the offsets in the frame.
254    fn station(read: Read, offsets: &StationOffsets) -> Vec3 {
255        read.point
256            + offsets.longitudinal * read.axes[0]
257            + offsets.lateral * read.axes[1]
258            + offsets.vertical * read.axes[2]
259    }
260
261    /// A local point placed in the axes `axes`, at `origin`.
262    fn placed(origin: Vec3, axes: [Vec3; 3], local: Point3) -> Vec3 {
263        origin + local.x * axes[0] + local.y * axes[1] + local.z * axes[2]
264    }
265
266    /// Axes turned by the unit components `w` of each of the base axes.
267    fn turned(axes: [Vec3; 3], w: Vec3) -> [Vec3; 3] {
268        let mix = |u: Vec3| u.x * axes[0] + u.y * axes[1] + u.z * axes[2];
269        [
270            mix(w),
271            mix(Vec3::new(w.y, w.z, w.x)),
272            mix(Vec3::new(w.z, w.x, w.y)),
273        ]
274    }
275
276    fn at() -> Point3 {
277        let p = station(base(), &offsets());
278        Point3::new(p.x, p.y, p.z)
279    }
280
281    const A: Scalar = 3e-7;
282    const B: Scalar = 2e-9;
283
284    fn n() -> Vec3 {
285        Vec3::new(1.0, -2.0, 2.0) / 3.0
286    }
287
288    #[test]
289    fn the_point_error_is_the_station_points_worst_displacement() {
290        let error = (A, B);
291        let offsets = offsets();
292        let exact = base();
293        let read = perturbed(exact, A, B, n());
294        let moved = (station(read, &offsets) - station(exact, &offsets)).length();
295        let reach = 3.5; // 2 + 1 + 0.5
296                         // The worst case reaches `a + b L1` exactly.
297        assert!((moved - (A + B * reach)).abs() <= 1e-12, "{moved}");
298        let bound = compose(error, &offsets, None, at());
299        assert!(bound.point_error >= moved, "{bound:?} < {moved}");
300        assert!(
301            bound.point_error <= moved * (1.0 + 1e-9) + 1e-9,
302            "{bound:?}"
303        );
304        // No offsets, no lever on the axes.
305        let flat = compose(error, &StationOffsets::ZERO, None, at());
306        assert!(flat.point_error < A + 1e-9 && flat.point_error >= A);
307    }
308
309    #[test]
310    fn a_turned_axis_moves_by_sqrt_three_of_the_axis_error() {
311        let w = Vec3::ONE / 3.0_f64.sqrt();
312        let exact = base();
313        let read = perturbed(exact, A, B, n());
314        let (e, r) = (turned(exact.axes, w), turned(read.axes, w));
315        let moved = (0..3)
316            .map(|i| (e[i] - r[i]).length())
317            .fold(0.0, Scalar::max);
318        assert!((moved - 3.0_f64.sqrt() * B).abs() <= 1e-15, "{moved}");
319        let error = (A, B);
320        let off = offsets();
321        let turned_bound = compose(error, &off, Some(1.0), at());
322        assert!(
323            turned_bound.axis_error >= moved,
324            "{turned_bound:?} < {moved}"
325        );
326        assert!(turned_bound.axis_error <= moved * (1.0 + 1e-6) + 1e-13);
327        // The base frame's axes move by `b` and the bound is not looser.
328        let base_bound = compose(error, &off, None, at());
329        assert!(base_bound.axis_error >= B && base_bound.axis_error <= B * (1.0 + 1e-9));
330        assert!(turned_bound.axis_error > base_bound.axis_error);
331    }
332
333    #[test]
334    fn a_body_vertex_moves_by_the_one_norm_of_its_local_point() {
335        let off = offsets();
336        let exact = base();
337        let read = perturbed(exact, A, B, n());
338        let error = (A, B);
339        let w = Vec3::ONE / 3.0_f64.sqrt();
340        for k in [0.0, 1.0, -3.0, 6.0] {
341            let local = Point3::new(k, k, k);
342            // The base frame.
343            let bound = compose(error, &off, None, at());
344            let moved = (placed(station(read, &off), read.axes, local)
345                - placed(station(exact, &off), exact.axes, local))
346            .length();
347            let shown = bound.displacement(local);
348            assert!(shown >= moved, "k {k}: {shown} < {moved}");
349            // `|q|_1` times the axes' error, and no looser than the margin
350            // where the lever's signs agree with the offsets'.
351            assert!(
352                k < 0.0 || shown <= moved * (1.0 + 1e-9) + 1e-9,
353                "k {k}: {shown} {moved}"
354            );
355            // A turned frame, components `(1, 1, 1) / sqrt 3`.
356            let bound = compose(error, &off, Some(1.0), at());
357            let moved = (placed(station(read, &off), turned(read.axes, w), local)
358                - placed(station(exact, &off), turned(exact.axes, w), local))
359            .length();
360            let shown = bound.displacement(local);
361            assert!(shown >= moved, "turned, k {k}: {shown} < {moved}");
362            assert!(
363                k < 0.0 || shown <= moved * (1.0 + 1e-6) + 1e-9,
364                "turned, k {k}: {shown} {moved}"
365            );
366        }
367    }
368}
axiolid_mesh_compile/station/accuracy.rs

axiolid_mesh_compile/station/
accuracy.rs
