Skip to main content

fpm_rs/
evaluation.rs

1//! Composed evaluation of a reconstruction against reference data.
2
3use num_complex::Complex64;
4use serde::{Deserialize, Serialize};
5
6use crate::{
7    Array2, Result,
8    measurements::MeasurementRead,
9    metrics::{
10        complex_field::{ComplexFieldComparisonMetrics, compare_complex_fields_masked},
11        intensity::{IntensityComparisonMetrics, compare_intensity, compare_intensity_masked},
12        model::{
13            FrameGainComparisonMetrics, IlluminationPositionMetrics, PupilComparisonMetrics,
14            compare_frame_gains, compare_illumination_positions, compare_pupils,
15        },
16    },
17    model::{ForwardModel, FourierOffset, ImagePlaneModel},
18    reconstruction::{ReconstructionProblem, ReconstructionResult},
19};
20
21#[derive(Clone, Debug, Serialize, Deserialize)]
22pub struct FrameIntensityEvaluation {
23    pub per_frame: Vec<IntensityComparisonMetrics>,
24}
25
26#[derive(Clone, Debug, Serialize, Deserialize)]
27pub struct ReconstructionEvaluation {
28    pub object: ComplexFieldComparisonMetrics,
29    pub pupil: Option<PupilComparisonMetrics>,
30    pub illumination: Option<IlluminationPositionMetrics>,
31    pub frame_gains: Option<FrameGainComparisonMetrics>,
32    pub intensity: Option<FrameIntensityEvaluation>,
33}
34
35pub fn evaluate_reconstruction(
36    result: &ReconstructionResult,
37    reference_object: &Array2<Complex64>,
38    reference_model: Option<&ImagePlaneModel>,
39    valid_object_mask: Option<&Array2<u8>>,
40) -> Result<ReconstructionEvaluation> {
41    let object =
42        compare_complex_fields_masked(reference_object, &result.object, valid_object_mask)?;
43    let (pupil, illumination, frame_gains) = if let Some(reference_model) = reference_model {
44        let pupil = Some(compare_pupils(
45            &reference_model.pupil.values,
46            &result.recovered_pupil.values,
47            &reference_model.pupil.support,
48        )?);
49        let illumination = None;
50        let frame_gains = result
51            .recovered_frame_gains
52            .as_ref()
53            .map(|candidate| {
54                let reference = reference_model
55                    .frame_gains
56                    .as_deref()
57                    .unwrap_or_else(|| &[]);
58                if reference.is_empty() {
59                    compare_frame_gains(&vec![1.0; reference_model.frame_count()], candidate)
60                } else {
61                    compare_frame_gains(reference, candidate)
62                }
63            })
64            .transpose()?;
65        (pupil, illumination, frame_gains)
66    } else {
67        (None, None, None)
68    };
69    Ok(ReconstructionEvaluation {
70        object,
71        pupil,
72        illumination,
73        frame_gains,
74        intensity: None,
75    })
76}
77
78pub fn evaluate_reconstruction_with_problem<M: MeasurementRead>(
79    result: &ReconstructionResult,
80    problem: &ReconstructionProblem<M>,
81    reference_object: &Array2<Complex64>,
82    reference_model: Option<&ImagePlaneModel>,
83    valid_object_mask: Option<&Array2<u8>>,
84) -> Result<ReconstructionEvaluation> {
85    let mut evaluation =
86        evaluate_reconstruction(result, reference_object, reference_model, valid_object_mask)?;
87    problem.validate()?;
88    let intensity = evaluate_frame_intensity(result, problem)?;
89    if let Some(reference_model) = reference_model {
90        evaluation.illumination = Some(compare_illumination_positions(
91            reference_model,
92            &problem.model,
93            result.calibrated_illumination.as_deref(),
94        )?);
95    }
96    evaluation.intensity = Some(intensity);
97    Ok(evaluation)
98}
99
100pub fn evaluate_frame_intensity<M: MeasurementRead>(
101    result: &ReconstructionResult,
102    problem: &ReconstructionProblem<M>,
103) -> Result<FrameIntensityEvaluation> {
104    problem.validate()?;
105    let model = model_with_result_calibration(problem, result)?;
106    let forward = ForwardModel::new(&model)?;
107    let mut workspace = forward.workspace();
108    let mut candidate = vec![0.0; model.image_shape.0 * model.image_shape.1];
109    let mut per_frame = Vec::with_capacity(model.frame_count());
110    for frame in 0..model.frame_count() {
111        if problem.measurements.frame_weight(frame)? == 0.0 {
112            per_frame.push(compare_intensity(&candidate, &candidate, None)?);
113            continue;
114        }
115        forward.forward_intensity_into(
116            &result.object_spectrum,
117            &result.recovered_pupil,
118            frame,
119            &mut workspace,
120            &mut candidate,
121        )?;
122        let reference = problem.measurements.frame(frame)?;
123        per_frame.push(compare_intensity_masked(
124            &reference,
125            &candidate,
126            problem.measurements.frame_mask(frame)?,
127            None,
128        )?);
129    }
130    Ok(FrameIntensityEvaluation { per_frame })
131}
132
133fn model_with_result_calibration<M: MeasurementRead>(
134    problem: &ReconstructionProblem<M>,
135    result: &ReconstructionResult,
136) -> Result<ImagePlaneModel> {
137    let mut model = problem.model.clone();
138    model.frame_gains = result.recovered_frame_gains.clone();
139    model.background = result.recovered_background.clone();
140    if let Some(corrections) = &result.calibrated_illumination {
141        model.subpixel_offsets = Some(
142            corrections
143                .iter()
144                .enumerate()
145                .map(|(source, &(row, column))| {
146                    let base = problem.model.source_offset(source)?;
147                    Ok(FourierOffset::new(base.row + row, base.column + column))
148                })
149                .collect::<Result<Vec<_>>>()?,
150        );
151    }
152    model.validate()?;
153    Ok(model)
154}