ifc_lite_processing/appearance/
registration.rs

1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at https://mozilla.org/MPL/2.0/.
4//! Bounded proper-rigid Kabsch solve. Held-out observations never enter fitting.
5use super::registration_types::*;
6use nalgebra::{linalg::SVD, DMatrix, Matrix3, Vector3};
7use sha2::{Digest, Sha256};
8use std::collections::HashSet;
9
10const ALGORITHM: &str = "ifclite-rigid-correspondence-v1";
11const MAX_POINTS: usize = 256;
12const MIN_NON_COLLINEARITY: f64 = 1e-10;
13
14fn label(value: &str) -> bool {
15    !value.trim().is_empty() && value.len() <= 256 && !value.chars().any(char::is_control)
16}
17fn validate(r: &ScanRegistrationRequest) -> Result<(), String> {
18    for frame in [&r.source_frame, &r.target_frame] {
19        if !label(&frame.frame_key)
20            || frame.asset_sha256.len() != 64
21            || !frame
22                .asset_sha256
23                .bytes()
24                .all(|b| b.is_ascii_digit() || (b'a'..=b'f').contains(&b))
25        {
26            return Err(
27                "Registration requires exact frame keys and lowercase SHA-256 asset identities"
28                    .into(),
29            );
30        }
31    }
32    if r.fit.len() < 3 || r.fit.len() > MAX_POINTS || r.held_out.len() > MAX_POINTS {
33        return Err("Registration needs 3–256 fit points and at most 256 held-out points".into());
34    }
35    let mut ids = HashSet::new();
36    let mut source_ids = HashSet::new();
37    let mut target_ids = HashSet::new();
38    let mut sources = HashSet::new();
39    let mut targets = HashSet::new();
40    // Canonicalize signed zero so renaming an identical coordinate cannot leak
41    // a fitting observation into the check set. Neighborhood overlap is the
42    // host's responsibility: this boundary only knows their stable identities.
43    let key = |p: [f64; 3]| p.map(|v| if v == 0. { 0 } else { v.to_bits() });
44    for p in r.fit.iter().chain(&r.held_out) {
45        if !label(&p.id)
46            || !label(&p.source_observation)
47            || !label(&p.target_feature)
48            || !p
49                .source
50                .iter()
51                .chain(&p.target)
52                .all(|v| v.is_finite() && v.abs() <= 1e12)
53        {
54            return Err("Registration needs bounded observation identities and finite coordinates within 1e12 metres".into());
55        }
56        if !ids.insert(&p.id)
57            || !source_ids.insert(&p.source_observation)
58            || !target_ids.insert(&p.target_feature)
59            || !sources.insert(key(p.source))
60            || !targets.insert(key(p.target))
61        {
62            return Err("Fit and held-out observations must be distinct: duplicate ID, feature or coordinate".into());
63        }
64    }
65    Ok(())
66}
67fn centroid(points: &[ScanCorrespondence], source: bool) -> Vector3<f64> {
68    let point = |p: &ScanCorrespondence| Vector3::from(if source { p.source } else { p.target });
69    let anchor = point(&points[0]);
70    anchor
71        + points
72            .iter()
73            .map(|p| point(p) - anchor)
74            .sum::<Vector3<f64>>()
75            / points.len() as f64
76}
77fn spread(points: &[Vector3<f64>]) -> Result<RegistrationSpread, String> {
78    // Avoid nalgebra's specialized Matrix3 SVD: it squares the matrix and loses
79    // small singular directions in thin, but accepted, correspondence sets.
80    // Dynamic dimensions select the direct bidiagonal algorithm instead.
81    let scatter = points
82        .iter()
83        .fold(Matrix3::zeros(), |sum, p| sum + p * p.transpose());
84    let values = SVD::try_new(
85        DMatrix::from_column_slice(3, 3, scatter.as_slice()),
86        false,
87        false,
88        f64::EPSILON,
89        128,
90    )
91    .ok_or("Registration scatter decomposition did not converge within its budget")?
92    .singular_values;
93    if !values.iter().all(|v| v.is_finite())
94        || values[0] <= 0.
95        || values[1] / values[0] < MIN_NON_COLLINEARITY
96    {
97        return Err(
98            "Degenerate registration: fitting points coincide or are effectively collinear".into(),
99        );
100    }
101    Ok(RegistrationSpread {
102        singular_values: std::array::from_fn(|i| values[i]),
103        non_collinearity_ratio: values[1] / values[0],
104        non_planarity_ratio: values[2] / values[0],
105    })
106}
107fn residuals(
108    points: &[ScanCorrespondence],
109    rotation: &Matrix3<f64>,
110    source: Vector3<f64>,
111    target: Vector3<f64>,
112) -> RegistrationResiduals {
113    let points: Vec<_> = points
114        .iter()
115        .map(|p| {
116            let error =
117                target + rotation * (Vector3::from(p.source) - source) - Vector3::from(p.target);
118            CorrespondenceResidual {
119                id: p.id.clone(),
120                vector_metres: error.into(),
121                distance_metres: error.norm(),
122            }
123        })
124        .collect();
125    let rms_metres = (!points.is_empty()).then(|| {
126        (points
127            .iter()
128            .map(|p| p.distance_metres.powi(2))
129            .sum::<f64>()
130            / points.len() as f64)
131            .sqrt()
132    });
133    let max_metres = points.iter().map(|p| p.distance_metres).reduce(f64::max);
134    RegistrationResiduals {
135        points,
136        rms_metres,
137        max_metres,
138    }
139}
140
141/// Pure fit in declared metre frames. The host must verify source bytes/frame
142/// revisions before using this report and must not promote it to an accuracy
143/// verdict without independently selected and spatially distributed checks.
144pub fn register_scan_correspondences(
145    r: &ScanRegistrationRequest,
146) -> Result<ScanRegistrationReport, String> {
147    validate(r)?;
148    let source_anchor = centroid(&r.fit, true);
149    let target_anchor = centroid(&r.fit, false);
150    let sources: Vec<_> = r
151        .fit
152        .iter()
153        .map(|p| Vector3::from(p.source) - source_anchor)
154        .collect();
155    let targets: Vec<_> = r
156        .fit
157        .iter()
158        .map(|p| Vector3::from(p.target) - target_anchor)
159        .collect();
160    let source_spread = spread(&sources)?;
161    let target_spread = spread(&targets)?;
162    let covariance = sources
163        .iter()
164        .zip(&targets)
165        .fold(Matrix3::zeros(), |sum, (s, t)| sum + s * t.transpose());
166    let svd = SVD::try_new(
167        DMatrix::from_column_slice(3, 3, covariance.as_slice()),
168        true,
169        true,
170        f64::EPSILON,
171        128,
172    )
173    .ok_or("Registration rotation decomposition did not converge within its budget")?;
174    if svd.singular_values[0] <= 0.
175        || svd.singular_values[1] / svd.singular_values[0] < MIN_NON_COLLINEARITY
176    {
177        return Err(
178            "Degenerate registration: correspondence covariance cannot constrain a unique rotation"
179                .into(),
180        );
181    }
182    let dynamic_u = svd.u.ok_or("Registration SVD omitted source basis")?;
183    let u = Matrix3::from_column_slice(dynamic_u.as_slice());
184    let dynamic_vt = svd
185        .v_t
186        .ok_or("Registration SVD omitted destination basis")?;
187    let v = Matrix3::from_column_slice(dynamic_vt.as_slice()).transpose();
188    let mut correction = Matrix3::identity();
189    correction[(2, 2)] = (v * u.transpose()).determinant().signum();
190    let rotation = v * correction * u.transpose();
191    if !rotation.iter().all(|v| v.is_finite()) || (rotation.determinant() - 1.).abs() > 1e-10 {
192        return Err("Registration could not produce a finite proper rotation".into());
193    }
194    let mut digest = Sha256::new();
195    digest.update(ALGORITHM.as_bytes());
196    digest.update([0]);
197    digest.update(
198        serde_json::to_vec(r).map_err(|e| format!("Cannot bind registration request: {e}"))?,
199    );
200    let mut diagnostics =
201        vec!["Mathematical rigid fit only; no automatic registration or accuracy approval".into()];
202    if r.held_out.is_empty() {
203        diagnostics.push("No held-out evidence supplied".into());
204    }
205    if r.fit.len() < 4 || r.held_out.len() < 4 {
206        diagnostics.push("F4 acceptance requires at least four fitting and four independently selected, spatially distributed checks".into());
207    }
208    if source_spread.non_planarity_ratio < 1e-10 || target_spread.non_planarity_ratio < 1e-10 {
209        diagnostics.push("Planar non-collinear fitting points are valid; check coverage at different heights separately".into());
210    }
211    if correction[(2, 2)] < 0. && svd.singular_values[2] / svd.singular_values[0] > 1e-10 {
212        diagnostics.push("Unconstrained solution would reflect; proper rotation retained and mismatch remains in residuals".into());
213    }
214    Ok(ScanRegistrationReport {
215        request_sha256: format!("{:x}", digest.finalize()),
216        algorithm: ALGORITHM.into(),
217        source_frame: r.source_frame.clone(),
218        target_frame: r.target_frame.clone(),
219        rotation: std::array::from_fn(|i| std::array::from_fn(|j| rotation[(i, j)])),
220        source_anchor: source_anchor.into(),
221        target_anchor: target_anchor.into(),
222        fit: residuals(&r.fit, &rotation, source_anchor, target_anchor),
223        held_out: residuals(&r.held_out, &rotation, source_anchor, target_anchor),
224        source_spread,
225        target_spread,
226        diagnostics,
227    })
228}
229
230#[cfg(test)]
231#[path = "registration_tests.rs"]
232mod tests;
ifc_lite_processing/appearance/registration.rs

ifc_lite_processing/appearance/
registration.rs
