Skip to main content

fpm_rs/diagnostics/
coverage.rs

1use serde::{Deserialize, Serialize};
2
3use crate::{experiment::KVector, model::ImagePlaneModel};
4
5/// Pixel bounds and center of one source crop in the high-resolution Fourier grid.
6#[derive(Clone, Debug, Serialize, Deserialize)]
7pub struct CropIndexDiagnostics {
8    /// Zero-based individual illumination-source index.
9    pub illumination_index: usize,
10
11    /// Inclusive first Fourier column (`x`).
12    pub x_start: usize,
13    /// Exclusive last Fourier column (`x`).
14    pub x_end: usize,
15
16    /// Inclusive first Fourier row (`y`).
17    pub y_start: usize,
18    /// Exclusive last Fourier row (`y`).
19    pub y_end: usize,
20
21    /// Crop center along the Fourier column axis, in grid pixels.
22    pub center_x: f64,
23    /// Crop center along the Fourier row axis, in grid pixels.
24    pub center_y: f64,
25}
26
27/// Summary of individual-source coverage on the high-resolution Fourier grid.
28#[derive(Clone, Debug, Serialize, Deserialize)]
29pub struct FourierCoverageDiagnostics {
30    /// Dimensionless synthetic numerical aperture, when wavelength metadata is available.
31    pub synthetic_na: Option<f64>,
32
33    /// Approximate circular-pupil radius in low-resolution Fourier-grid pixels.
34    pub pupil_radius_px: f64,
35
36    /// Source crop centers as `[x_column, y_row]` in high-resolution grid pixels.
37    pub pupil_centers_px: Vec<[f64; 2]>,
38
39    /// Dimensionless illumination NA for each individual source.
40    pub illumination_na: Vec<f64>,
41
42    /// Integer crop diagnostics in individual source order.
43    pub crop_indices: Vec<CropIndexDiagnostics>,
44
45    /// High-resolution `[height, width]` of pixels covered by at least two pupils.
46    pub overlap_shape: Option<[usize; 2]>,
47}
48
49/// Validates `model` and derives source centers, crop bounds, NA, and overlap coverage.
50pub fn compute_fourier_coverage(
51    model: &ImagePlaneModel,
52) -> crate::Result<FourierCoverageDiagnostics> {
53    model.validate()?;
54    let wavelength = model.sampling.wavelength.unwrap_or(0.0);
55    let illumination_na = model
56        .k_vectors
57        .iter()
58        .map(|vector| k_vector_na(vector, wavelength))
59        .collect();
60    let crop_indices = model
61        .crop_indices
62        .crops
63        .iter()
64        .enumerate()
65        .map(|(illumination_index, crop)| {
66            let offset = model.source_offset(illumination_index)?;
67            Ok(CropIndexDiagnostics {
68                illumination_index,
69                x_start: crop.start_col,
70                x_end: crop.start_col + crop.width,
71                y_start: crop.start_row,
72                y_end: crop.start_row + crop.height,
73                center_x: crop.start_col as f64 + crop.width as f64 / 2.0 + offset.column,
74                center_y: crop.start_row as f64 + crop.height as f64 / 2.0 + offset.row,
75            })
76        })
77        .collect::<crate::Result<Vec<_>>>()?;
78    let pupil_radius_px = pupil_radius_px(&model.pupil);
79    let pupil_centers_px = crop_indices
80        .iter()
81        .map(|crop| [crop.center_x, crop.center_y])
82        .collect();
83    Ok(FourierCoverageDiagnostics {
84        synthetic_na: model.sampling.synthetic_na,
85        pupil_radius_px,
86        pupil_centers_px,
87        illumination_na,
88        crop_indices,
89        overlap_shape: None,
90    })
91}
92
93fn k_vector_na(vector: &KVector, wavelength: f64) -> f64 {
94    if wavelength <= 0.0 {
95        0.0
96    } else {
97        vector.kx.hypot(vector.ky) * wavelength / std::f64::consts::TAU
98    }
99}
100
101fn pupil_radius_px(pupil: &crate::model::Pupil) -> f64 {
102    let shape = pupil.shape();
103    let center_y = shape.0 as f64 / 2.0;
104    let center_x = shape.1 as f64 / 2.0;
105    let mut max_radius: f64 = 0.0;
106    for row in 0..shape.0 {
107        for col in 0..shape.1 {
108            let index = row * shape.1 + col;
109            if pupil.support.as_slice()[index] != 0 {
110                let radius = (row as f64 - center_y).hypot(col as f64 - center_x);
111                max_radius = max_radius.max(radius);
112            }
113        }
114    }
115    max_radius
116}