spatialrust_metrics/
distance.rs

1//! Chamfer and Hausdorff distances between two point clouds.
2
3use spatialrust_core::{HasPositions3, PointCloud, SpatialError, SpatialResult};
4use spatialrust_search::{KdTree, NearestNeighborIndex};
5
6/// A bundle of directed and symmetric cloud-to-cloud distance statistics.
7///
8/// "Directed" means nearest-neighbor distances from one cloud to the other; the
9/// metrics are not symmetric on their own, so both directions are reported.
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub struct CloudDistances {
12    /// Mean squared NN distance from `a` to `b`.
13    pub mean_sq_a_to_b: f64,
14    /// Mean squared NN distance from `b` to `a`.
15    pub mean_sq_b_to_a: f64,
16    /// Largest NN distance from `a` to `b` (directed Hausdorff).
17    pub max_a_to_b: f64,
18    /// Largest NN distance from `b` to `a` (directed Hausdorff).
19    pub max_b_to_a: f64,
20}
21
22impl CloudDistances {
23    /// Symmetric Chamfer distance: the sum of the two mean squared NN distances.
24    #[must_use]
25    pub fn chamfer(&self) -> f64 {
26        self.mean_sq_a_to_b + self.mean_sq_b_to_a
27    }
28
29    /// Symmetric Hausdorff distance: the larger of the two directed maxima.
30    #[must_use]
31    pub fn hausdorff(&self) -> f64 {
32        self.max_a_to_b.max(self.max_b_to_a)
33    }
34}
35
36/// Computes directed NN distance statistics between `a` and `b` in both
37/// directions with a single KD-tree per side.
38pub fn cloud_distances(a: &PointCloud, b: &PointCloud) -> SpatialResult<CloudDistances> {
39    if a.is_empty() || b.is_empty() {
40        return Err(SpatialError::InvalidArgument("both clouds must be non-empty".to_owned()));
41    }
42    let (ax, ay, az) = a.positions3()?;
43    let (bx, by, bz) = b.positions3()?;
44
45    let tree_b = KdTree::from_slices(bx, by, bz);
46    let tree_a = KdTree::from_slices(ax, ay, az);
47
48    let (mean_sq_a_to_b, max_a_to_b) = directed(ax, ay, az, &tree_b);
49    let (mean_sq_b_to_a, max_b_to_a) = directed(bx, by, bz, &tree_a);
50
51    Ok(CloudDistances { mean_sq_a_to_b, mean_sq_b_to_a, max_a_to_b, max_b_to_a })
52}
53
54/// Mean squared and max NN distance from each query point to `target`.
55fn directed(x: &[f32], y: &[f32], z: &[f32], target: &KdTree) -> (f64, f64) {
56    let mut sum_sq = 0.0_f64;
57    let mut max = 0.0_f64;
58    for i in 0..x.len() {
59        if let Some(neighbor) = target.nearest_one(x[i], y[i], z[i]) {
60            let d_sq = f64::from(neighbor.distance_squared);
61            sum_sq += d_sq;
62            let d = d_sq.sqrt();
63            if d > max {
64                max = d;
65            }
66        }
67    }
68    (sum_sq / x.len() as f64, max)
69}
70
71/// Symmetric Chamfer distance between two clouds (sum of mean squared NN
72/// distances in both directions). Lower is better; zero for identical clouds.
73pub fn chamfer_distance(a: &PointCloud, b: &PointCloud) -> SpatialResult<f64> {
74    Ok(cloud_distances(a, b)?.chamfer())
75}
76
77/// Symmetric Hausdorff distance between two clouds (the largest NN distance in
78/// either direction). Captures the worst-case discrepancy.
79pub fn hausdorff_distance(a: &PointCloud, b: &PointCloud) -> SpatialResult<f64> {
80    Ok(cloud_distances(a, b)?.hausdorff())
81}
82
83#[cfg(test)]
84mod tests {
85    use super::*;
86    use spatialrust_core::{PointCloudBuilder, StandardSchemas};
87
88    fn cloud(points: &[[f32; 3]]) -> PointCloud {
89        let mut builder = PointCloudBuilder::new(StandardSchemas::point_xyz());
90        for p in points {
91            builder.push_point(*p).unwrap();
92        }
93        builder.build().unwrap()
94    }
95
96    #[test]
97    fn identical_clouds_have_zero_distance() {
98        let c = cloud(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]);
99        assert!(chamfer_distance(&c, &c).unwrap() < 1e-9);
100        assert!(hausdorff_distance(&c, &c).unwrap() < 1e-9);
101    }
102
103    #[test]
104    fn translation_shows_up_in_distances() {
105        let a = cloud(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]);
106        let b = cloud(&[[0.0, 0.5, 0.0], [1.0, 0.5, 0.0], [2.0, 0.5, 0.0]]);
107        // Every point's nearest neighbor is exactly 0.5 away.
108        let chamfer = chamfer_distance(&a, &b).unwrap();
109        assert!((chamfer - (0.25 + 0.25)).abs() < 1e-6);
110        let hausdorff = hausdorff_distance(&a, &b).unwrap();
111        assert!((hausdorff - 0.5).abs() < 1e-6);
112    }
113
114    #[test]
115    fn hausdorff_catches_a_single_outlier() {
116        let a = cloud(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]);
117        let b = cloud(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [10.0, 0.0, 0.0]]);
118        // The lone far point in b is 9.0 from its nearest neighbor in a.
119        assert!((hausdorff_distance(&a, &b).unwrap() - 9.0).abs() < 1e-6);
120    }
121}
spatialrust_metrics/distance.rs

spatialrust_metrics/
distance.rs
