fpm_rs/diagnostics/
coverage.rs1use 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}