fpm_rs/diagnostics/
coverage.rs1use serde::{Deserialize, Serialize};
2
3use crate::{experiment::KVector, model::ImagePlaneModel};
4
5#[derive(Clone, Debug, Serialize, Deserialize)]
7pub struct CropIndexDiagnostics {
8 pub illumination_index: usize,
10
11 pub x_start: usize,
13 pub x_end: usize,
15
16 pub y_start: usize,
18 pub y_end: usize,
20
21 pub center_x: f64,
23 pub center_y: f64,
25}
26
27#[derive(Clone, Debug, Serialize, Deserialize)]
29pub struct FourierCoverageDiagnostics {
30 pub synthetic_na: Option<f64>,
32
33 pub pupil_radius_px: f64,
35
36 pub pupil_centers_px: Vec<[f64; 2]>,
38
39 pub illumination_na: Vec<f64>,
41
42 pub crop_indices: Vec<CropIndexDiagnostics>,
44
45 pub overlap_shape: Option<[usize; 2]>,
47}
48
49pub 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}