fpm_rs/simulation/
camera.rs

1use rand::Rng;
2use rand_distr::{Distribution, Normal, Poisson};
3use serde::{Deserialize, Serialize};
4
5use crate::{Result, error::Error, model::ImagePlaneModel};
6
7/// Concrete detector pipeline from optical intensity to digital counts.
8#[derive(Clone, Debug, Serialize, Deserialize)]
9pub struct CameraModel {
10    /// Expected photoelectrons per pixel at unit optical intensity.
11    pub photons_per_pixel: f64,
12    pub gain_counts_per_electron: f64,
13    pub offset_counts: f64,
14    pub read_noise_electrons: f64,
15    pub dark_current_electrons: f64,
16    /// Whether to Poisson-sample photoelectrons and dark current.
17    pub shot_noise: bool,
18    /// Multiplicative detector sensitivity for each pixel.
19    pub pixel_sensitivity: Option<Vec<f64>>,
20    pub bit_depth: Option<u8>,
21    pub saturation_counts: Option<f64>,
22    pub quantize: bool,
23    pub bad_pixels: Vec<usize>,
24    pub bad_pixel_value_counts: Option<f64>,
25}
26
27impl Default for CameraModel {
28    fn default() -> Self {
29        Self {
30            photons_per_pixel: 1_000.0,
31            gain_counts_per_electron: 1.0,
32            offset_counts: 0.0,
33            read_noise_electrons: 0.0,
34            dark_current_electrons: 0.0,
35            shot_noise: false,
36            pixel_sensitivity: None,
37            bit_depth: Some(16),
38            saturation_counts: None,
39            quantize: true,
40            bad_pixels: Vec::new(),
41            bad_pixel_value_counts: None,
42        }
43    }
44}
45
46impl CameraModel {
47    pub fn new() -> Self {
48        Self::default()
49    }
50
51    /// Unit detector response without noise, clipping, or quantization.
52    pub fn ideal() -> Self {
53        Self {
54            photons_per_pixel: 1.0,
55            gain_counts_per_electron: 1.0,
56            offset_counts: 0.0,
57            read_noise_electrons: 0.0,
58            dark_current_electrons: 0.0,
59            shot_noise: false,
60            pixel_sensitivity: None,
61            bit_depth: None,
62            saturation_counts: None,
63            quantize: false,
64            bad_pixels: Vec::new(),
65            bad_pixel_value_counts: None,
66        }
67    }
68
69    pub fn photons_per_pixel(mut self, value: f64) -> Self {
70        self.photons_per_pixel = value;
71        self
72    }
73
74    pub fn gain(mut self, counts_per_electron: f64) -> Self {
75        self.gain_counts_per_electron = counts_per_electron;
76        self
77    }
78
79    pub fn offset_counts(mut self, value: f64) -> Self {
80        self.offset_counts = value;
81        self
82    }
83
84    pub fn read_noise_electrons(mut self, value: f64) -> Self {
85        self.read_noise_electrons = value;
86        self
87    }
88
89    pub fn dark_current_electrons(mut self, value: f64) -> Self {
90        self.dark_current_electrons = value;
91        self
92    }
93
94    pub fn shot_noise(mut self, enabled: bool) -> Self {
95        self.shot_noise = enabled;
96        self
97    }
98
99    pub fn pixel_sensitivity(mut self, values: Vec<f64>) -> Self {
100        self.pixel_sensitivity = Some(values);
101        self
102    }
103
104    pub fn bit_depth(mut self, bits: u8) -> Self {
105        self.bit_depth = Some(bits);
106        self
107    }
108
109    pub fn saturation(mut self, counts: f64) -> Self {
110        self.saturation_counts = Some(counts);
111        self
112    }
113
114    pub fn quantize(mut self, enabled: bool) -> Self {
115        self.quantize = enabled;
116        self
117    }
118
119    pub fn bad_pixels(mut self, indices: Vec<usize>, value_counts: f64) -> Self {
120        self.bad_pixels = indices;
121        self.bad_pixel_value_counts = Some(value_counts);
122        self
123    }
124
125    pub fn validate(&self) -> Result<()> {
126        for (name, value) in [
127            ("photons_per_pixel", self.photons_per_pixel),
128            ("gain_counts_per_electron", self.gain_counts_per_electron),
129        ] {
130            if !value.is_finite() || value <= 0.0 {
131                return Err(Error::InvalidParameter {
132                    name,
133                    reason: "must be finite and positive".into(),
134                });
135            }
136        }
137        if !self.offset_counts.is_finite()
138            || !self.read_noise_electrons.is_finite()
139            || self.read_noise_electrons < 0.0
140            || !self.dark_current_electrons.is_finite()
141            || self.dark_current_electrons < 0.0
142        {
143            return Err(Error::InvalidParameter {
144                name: "camera noise/offset",
145                reason: "offset must be finite; read noise and dark current must be non-negative"
146                    .into(),
147            });
148        }
149        if self.bit_depth.is_some_and(|bits| bits == 0 || bits > 32) {
150            return Err(Error::InvalidParameter {
151                name: "bit_depth",
152                reason: "must be between 1 and 32".into(),
153            });
154        }
155        if self
156            .saturation_counts
157            .is_some_and(|value| !value.is_finite() || value <= 0.0)
158        {
159            return Err(Error::InvalidParameter {
160                name: "saturation_counts",
161                reason: "must be finite and positive".into(),
162            });
163        }
164        if self
165            .bad_pixel_value_counts
166            .is_some_and(|value| !value.is_finite() || value < 0.0)
167        {
168            return Err(Error::InvalidParameter {
169                name: "bad_pixel_value_counts",
170                reason: "must be finite and non-negative".into(),
171            });
172        }
173        if !self.bad_pixels.is_empty()
174            && self.bad_pixel_value_counts.is_none()
175            && !self.maximum_count().is_finite()
176        {
177            return Err(Error::InvalidParameter {
178                name: "bad_pixels",
179                reason: "a finite bad-pixel value or camera maximum is required".into(),
180            });
181        }
182        if self.pixel_sensitivity.as_ref().is_some_and(|values| {
183            values
184                .iter()
185                .any(|value| !value.is_finite() || *value < 0.0)
186        }) {
187            return Err(Error::InvalidParameter {
188                name: "pixel_sensitivity",
189                reason: "values must be finite and non-negative".into(),
190            });
191        }
192        let mut bad_pixels = self.bad_pixels.clone();
193        bad_pixels.sort_unstable();
194        if bad_pixels.windows(2).any(|pair| pair[0] == pair[1]) {
195            return Err(Error::InvalidParameter {
196                name: "bad_pixels",
197                reason: "indices must be unique".into(),
198            });
199        }
200        Ok(())
201    }
202
203    pub fn validate_for_frame(&self, frame_len: usize) -> Result<()> {
204        self.validate()?;
205        if self
206            .pixel_sensitivity
207            .as_ref()
208            .is_some_and(|values| values.len() != frame_len)
209        {
210            return Err(Error::InvalidParameter {
211                name: "pixel_sensitivity",
212                reason: format!("must contain {frame_len} values"),
213            });
214        }
215        if self.bad_pixels.iter().any(|&pixel| pixel >= frame_len) {
216            return Err(Error::InvalidParameter {
217                name: "bad_pixels",
218                reason: format!("indices must be below the frame size {frame_len}"),
219            });
220        }
221        Ok(())
222    }
223
224    /// Converts one optical-intensity frame into detector counts in place.
225    pub fn measure_frame<R: Rng + ?Sized>(&self, intensity: &mut [f64], rng: &mut R) -> Result<()> {
226        self.validate_for_frame(intensity.len())?;
227        let read_noise = (self.read_noise_electrons > 0.0)
228            .then(|| Normal::new(0.0, self.read_noise_electrons))
229            .transpose()
230            .map_err(|error| Error::Numerical(error.to_string()))?;
231        for (pixel, value) in intensity.iter_mut().enumerate() {
232            let sensitivity = self
233                .pixel_sensitivity
234                .as_ref()
235                .map_or(1.0, |values| values[pixel]);
236            let expected_electrons =
237                value.max(0.0) * sensitivity * self.photons_per_pixel + self.dark_current_electrons;
238            let mut electrons = if self.shot_noise && expected_electrons > 0.0 {
239                Poisson::new(expected_electrons)
240                    .map_err(|error| Error::Numerical(error.to_string()))?
241                    .sample(rng)
242            } else {
243                expected_electrons
244            };
245            if let Some(distribution) = &read_noise {
246                electrons += distribution.sample(rng);
247            }
248            let mut counts = electrons * self.gain_counts_per_electron + self.offset_counts;
249            counts = counts.clamp(0.0, self.maximum_count());
250            if self.quantize {
251                counts = counts.round();
252            }
253            *value = counts;
254        }
255        for &pixel in &self.bad_pixels {
256            let mut value = self
257                .bad_pixel_value_counts
258                .unwrap_or_else(|| self.maximum_count())
259                .clamp(0.0, self.maximum_count());
260            if self.quantize {
261                value = value.round();
262            }
263            intensity[pixel] = value;
264        }
265        Ok(())
266    }
267
268    /// Compiles known uniform linear response into a reconstruction model.
269    ///
270    /// Pixel sensitivity, stochastic noise, clipping, quantization, and bad
271    /// pixels remain detector effects and are not compiled into the model.
272    pub fn compile_reconstruction_model(
273        &self,
274        mut model: ImagePlaneModel,
275    ) -> Result<ImagePlaneModel> {
276        self.validate_for_frame(model.image_shape.0 * model.image_shape.1)?;
277        let scale = self.photons_per_pixel * self.gain_counts_per_electron;
278        model.frame_gains = Some(match &model.frame_gains {
279            Some(gains) => gains.iter().map(|gain| gain * scale).collect(),
280            None => vec![scale; model.frame_count()],
281        });
282        let additive_counts =
283            self.dark_current_electrons * self.gain_counts_per_electron + self.offset_counts;
284        if let Some(background) = &mut model.background {
285            for value in background {
286                *value = *value * scale + additive_counts;
287            }
288        } else if additive_counts != 0.0 {
289            model.background = Some(vec![
290                additive_counts;
291                model.image_shape.0 * model.image_shape.1
292            ]);
293        }
294        model.validate()?;
295        Ok(model)
296    }
297
298    fn maximum_count(&self) -> f64 {
299        let digital_maximum = self
300            .bit_depth
301            .map_or(f64::INFINITY, |bits| (2_f64).powi(bits as i32) - 1.0);
302        self.saturation_counts
303            .unwrap_or(f64::INFINITY)
304            .min(digital_maximum)
305    }
306}
fpm_rs/simulation/camera.rs

fpm_rs/simulation/
camera.rs
