Skip to main content

fpm_rs/diagnostics/
coverage.rs

1use serde::{Deserialize, Serialize};
2
3use crate::{experiment::KVector, model::ImagePlaneModel};
4
5#[derive(Clone, Debug, Serialize, Deserialize)]
6pub struct CropIndexDiagnostics {
7    pub illumination_index: usize,
8
9    pub x_start: usize,
10    pub x_end: usize,
11
12    pub y_start: usize,
13    pub y_end: usize,
14
15    pub center_x: f64,
16    pub center_y: f64,
17}
18
19#[derive(Clone, Debug, Serialize, Deserialize)]
20pub struct FourierCoverageDiagnostics {
21    pub synthetic_na: Option<f64>,
22
23    pub pupil_radius_px: f64,
24
25    pub pupil_centers_px: Vec<[f64; 2]>,
26
27    pub illumination_na: Vec<f64>,
28
29    pub crop_indices: Vec<CropIndexDiagnostics>,
30
31    pub overlap_shape: Option<[usize; 2]>,
32}
33
34pub fn compute_fourier_coverage(
35    model: &ImagePlaneModel,
36) -> crate::Result<FourierCoverageDiagnostics> {
37    model.validate()?;
38    let wavelength = model.sampling.wavelength.unwrap_or(0.0);
39    let illumination_na = model
40        .k_vectors
41        .iter()
42        .map(|vector| k_vector_na(vector, wavelength))
43        .collect();
44    let crop_indices = model
45        .crop_indices
46        .crops
47        .iter()
48        .enumerate()
49        .map(|(illumination_index, crop)| {
50            let offset = model.source_offset(illumination_index)?;
51            Ok(CropIndexDiagnostics {
52                illumination_index,
53                x_start: crop.start_col,
54                x_end: crop.start_col + crop.width,
55                y_start: crop.start_row,
56                y_end: crop.start_row + crop.height,
57                center_x: crop.start_col as f64 + crop.width as f64 / 2.0 + offset.column,
58                center_y: crop.start_row as f64 + crop.height as f64 / 2.0 + offset.row,
59            })
60        })
61        .collect::<crate::Result<Vec<_>>>()?;
62    let pupil_radius_px = pupil_radius_px(&model.pupil);
63    let pupil_centers_px = crop_indices
64        .iter()
65        .map(|crop| [crop.center_x, crop.center_y])
66        .collect();
67    Ok(FourierCoverageDiagnostics {
68        synthetic_na: model.sampling.synthetic_na,
69        pupil_radius_px,
70        pupil_centers_px,
71        illumination_na,
72        crop_indices,
73        overlap_shape: None,
74    })
75}
76
77fn k_vector_na(vector: &KVector, wavelength: f64) -> f64 {
78    if wavelength <= 0.0 {
79        0.0
80    } else {
81        vector.kx.hypot(vector.ky) * wavelength / std::f64::consts::TAU
82    }
83}
84
85fn pupil_radius_px(pupil: &crate::model::Pupil) -> f64 {
86    let shape = pupil.values.shape();
87    let center_y = shape.0 as f64 / 2.0;
88    let center_x = shape.1 as f64 / 2.0;
89    let mut max_radius: f64 = 0.0;
90    for row in 0..shape.0 {
91        for col in 0..shape.1 {
92            let index = row * shape.1 + col;
93            if pupil.support[index] {
94                let radius = (row as f64 - center_y).hypot(col as f64 - center_x);
95                max_radius = max_radius.max(radius);
96            }
97        }
98    }
99    max_radius
100}