fpm_rs/simulation/
camera.rs

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

fpm_rs/simulation/
camera.rs
