Skip to main content

fpm_rs/metrics/intensity/
comparison.rs

1//! Metrics comparing a reference and candidate intensity image.
2
3use serde::{Deserialize, Serialize};
4
5#[derive(Clone, Debug, Serialize, Deserialize)]
6pub struct IntensityComparisonMetrics {
7    /// Serialized as `measured_sum` inside diagnostic records for wire compatibility.
8    #[serde(rename = "measured_sum", alias = "reference_sum")]
9    pub reference_sum: f64,
10    /// Serialized as `predicted_sum` inside diagnostic records for wire compatibility.
11    #[serde(rename = "predicted_sum", alias = "candidate_sum")]
12    pub candidate_sum: f64,
13    pub residual_l1: f64,
14    pub residual_l2: f64,
15    pub residual_mean: f64,
16    pub residual_std: f64,
17    pub residual_max_abs: f64,
18    pub normalized_l2: f64,
19    pub saturated_pixels: Option<usize>,
20}
21
22pub fn compare_intensity(
23    reference: &[f64],
24    candidate: &[f64],
25    saturation_value: Option<f64>,
26) -> crate::Result<IntensityComparisonMetrics> {
27    compare_intensity_masked(reference, candidate, None, saturation_value)
28}
29
30pub fn compare_intensity_masked(
31    reference: &[f64],
32    candidate: &[f64],
33    mask: Option<&[u8]>,
34    saturation_value: Option<f64>,
35) -> crate::Result<IntensityComparisonMetrics> {
36    if reference.len() != candidate.len() || reference.is_empty() {
37        return Err(crate::Error::InvalidShape(format!(
38            "intensity comparison inputs have lengths {} and {}",
39            reference.len(),
40            candidate.len()
41        )));
42    }
43    if mask.is_some_and(|mask| mask.len() != reference.len()) {
44        return Err(crate::Error::LengthMismatch {
45            actual: mask.map_or(0, <[u8]>::len),
46            expected: reference.len(),
47            shape: (1, reference.len()),
48        });
49    }
50    let mut reference_sum = 0.0;
51    let mut candidate_sum = 0.0;
52    let mut residual_l1 = 0.0;
53    let mut residual_squared = 0.0;
54    let mut residual_sum = 0.0;
55    let mut residual_max_abs: f64 = 0.0;
56    let mut reference_squared = 0.0;
57    let mut saturated_pixels = 0;
58    let mut count = 0usize;
59    for (index, (&reference, &candidate)) in reference.iter().zip(candidate).enumerate() {
60        if mask.is_some_and(|mask| mask[index] == 0) {
61            continue;
62        }
63        if !reference.is_finite() || !candidate.is_finite() {
64            return Err(crate::Error::Numerical(
65                "intensity comparison contains a non-finite value".into(),
66            ));
67        }
68        let residual = candidate - reference;
69        reference_sum += reference;
70        candidate_sum += candidate;
71        residual_l1 += residual.abs();
72        residual_squared += residual * residual;
73        residual_sum += residual;
74        residual_max_abs = residual_max_abs.max(residual.abs());
75        reference_squared += reference * reference;
76        saturated_pixels += usize::from(saturation_value.is_some_and(|limit| reference >= limit));
77        count += 1;
78    }
79    let count_f64 = count as f64;
80    let residual_mean = (count > 0)
81        .then_some(residual_sum / count_f64)
82        .unwrap_or(0.0);
83    Ok(IntensityComparisonMetrics {
84        reference_sum,
85        candidate_sum,
86        residual_l1,
87        residual_l2: residual_squared.sqrt(),
88        residual_mean,
89        residual_std: if count == 0 {
90            0.0
91        } else {
92            (residual_squared / count_f64 - residual_mean * residual_mean)
93                .max(0.0)
94                .sqrt()
95        },
96        residual_max_abs,
97        normalized_l2: if count == 0 {
98            0.0
99        } else {
100            residual_squared.sqrt() / (reference_squared.sqrt() + f64::EPSILON)
101        },
102        saturated_pixels: saturation_value.map(|_| saturated_pixels),
103    })
104}