use super::collate::Collated;
use crate::mesh::Mesh;
use nalgebra::{Matrix4, Vector4};
const ABS_FLOOR_M: f64 = 1e-6;
const ULP_FACTOR: f64 = 9.5367431640625e-7;
fn tolerance(mag: f64) -> f64 {
(mag * ULP_FACTOR).max(ABS_FLOOR_M)
}
fn stored_magnitude(positions: &[f32], v: usize) -> f64 {
(positions[v * 3] as f64)
.abs()
.max((positions[v * 3 + 1] as f64).abs())
.max((positions[v * 3 + 2] as f64).abs())
}
fn reconstructed_vertex_error(
template_origin: [f64; 3],
template_positions: &[f32],
target_origin: [f64; 3],
target_positions: &[f32],
rel: &Matrix4<f64>,
v: usize,
) -> f64 {
let tx = template_origin[0] + template_positions[v * 3] as f64;
let ty = template_origin[1] + template_positions[v * 3 + 1] as f64;
let tz = template_origin[2] + template_positions[v * 3 + 2] as f64;
let w = rel * Vector4::new(tx, ty, tz, 1.0);
let (rx, ry, rz) = (w.x / w.w, w.y / w.w, w.z / w.w);
let gx = target_origin[0] + target_positions[v * 3] as f64;
let gy = target_origin[1] + target_positions[v * 3 + 1] as f64;
let gz = target_origin[2] + target_positions[v * 3 + 2] as f64;
((rx - gx).powi(2) + (ry - gy).powi(2) + (rz - gz).powi(2)).sqrt()
}
pub(super) fn verify_pairing(
template_origin: [f64; 3],
template_positions: &[f32],
target_origin: [f64; 3],
target_positions: &[f32],
rel: &Matrix4<f64>,
) -> bool {
let n = template_positions.len() / 3;
if target_positions.len() / 3 != n {
return false;
}
for v in 0..n {
let err = reconstructed_vertex_error(
template_origin,
template_positions,
target_origin,
target_positions,
rel,
v,
);
let mag = stored_magnitude(template_positions, v).max(stored_magnitude(target_positions, v));
if !mag.is_finite() {
return false;
}
if err.is_nan() || err > tolerance(mag) {
return false;
}
}
true
}
pub fn verify_recomposition(meshes: &[Mesh], collated: &Collated) -> f64 {
let mut max_err = 0.0f64;
for tmpl in &collated.templates {
let template = &meshes[tmpl.template_index];
for occ in &tmpl.occurrences {
let target = &meshes[occ.mesh_index];
let rel = Matrix4::from_row_slice(&occ.transform.map(|v| v as f64));
let n = template.positions.len() / 3;
if target.positions.len() / 3 != n {
max_err = f64::INFINITY;
continue;
}
for v in 0..n {
let err = reconstructed_vertex_error(
template.origin,
&template.positions,
target.origin,
&target.positions,
&rel,
v,
);
if err.is_nan() {
max_err = f64::INFINITY;
continue;
}
if err > max_err {
max_err = err;
}
}
}
}
max_err
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn tolerance_scales_with_magnitude_above_the_floor() {
assert_eq!(tolerance(0.0), ABS_FLOOR_M);
assert_eq!(tolerance(1.0), ABS_FLOOR_M);
let far = tolerance(220_000.0);
assert!(
far > ABS_FLOOR_M && far < 1.0,
"expected a sub-metre tolerance at 220km, got {far}"
);
}
#[test]
fn identical_geometry_through_identity_rel_passes() {
let positions: [f32; 6] = [0.0, 0.0, 0.0, 1.0, 2.0, 3.0];
let identity = Matrix4::identity();
assert!(verify_pairing(
[0.0, 0.0, 0.0],
&positions,
[0.0, 0.0, 0.0],
&positions,
&identity,
));
}
#[test]
fn a_non_finite_transform_is_rejected() {
let positions: [f32; 6] = [0.0, 0.0, 0.0, 1.0, 2.0, 3.0];
let mut rel = Matrix4::identity();
rel[(0, 3)] = f64::NAN;
assert!(
!verify_pairing([0.0; 3], &positions, [0.0; 3], &positions, &rel),
"a NaN transform must not pass verification"
);
}
#[test]
fn a_non_finite_position_is_rejected() {
let finite: [f32; 3] = [0.0, 0.0, 0.0];
let nan: [f32; 3] = [f32::NAN, 0.0, 0.0];
assert!(
!verify_pairing([0.0; 3], &finite, [0.0; 3], &nan, &Matrix4::identity()),
"a NaN baked target position must not pass verification"
);
assert!(
!verify_pairing([0.0; 3], &nan, [0.0; 3], &finite, &Matrix4::identity()),
"a NaN baked template position must not pass verification"
);
}
#[test]
fn an_infinite_target_position_is_rejected() {
let template: [f32; 3] = [0.0, 0.0, 0.0];
let target: [f32; 3] = [f32::INFINITY, 0.0, 0.0];
assert!(
!verify_pairing(
[0.0; 3],
&template,
[0.0; 3],
&target,
&Matrix4::identity()
),
"an infinite baked position must not pass verification"
);
}
#[test]
fn a_translated_target_fails_identity_rel() {
let template: [f32; 3] = [0.0, 0.0, 0.0];
let target: [f32; 3] = [2.0, 0.0, 0.0];
let identity = Matrix4::identity();
assert!(!verify_pairing(
[0.0, 0.0, 0.0],
&template,
[0.0, 0.0, 0.0],
&target,
&identity,
));
}
#[test]
fn a_genuine_match_still_passes_after_the_nan_guard() {
let template: [f32; 6] = [0.0, 0.0, 0.0, 1.0, 2.0, 3.0];
let rel = Matrix4::new_translation(&nalgebra::Vector3::new(5.0, -2.0, 0.5));
let target: [f32; 6] = [5.0, -2.0, 0.5, 6.0, 0.0, 3.5];
assert!(
verify_pairing(
[0.0, 0.0, 0.0],
&template,
[0.0, 0.0, 0.0],
&target,
&rel,
),
"a genuine matching pair must still pass verification"
);
}
}