fpm_rs/simulation/
synthetic_object.rs

1use num_complex::Complex64;
2use rand::{Rng, SeedableRng, rngs::StdRng};
3use rand_distr::{Distribution, Normal};
4use std::path::Path;
5
6use crate::{
7    Array2, Result, complex,
8    error::Error,
9    image_io::{GrayscaleScaling, load_grayscale},
10};
11
12#[derive(Clone, Debug)]
13pub struct SyntheticObject {
14    pub field: Array2<Complex64>,
15    pub label: Option<String>,
16}
17
18impl SyntheticObject {
19    pub fn new(field: Array2<Complex64>) -> Self {
20        Self { field, label: None }
21    }
22
23    pub fn from_amplitude_phase(amplitude: &Array2<f64>, phase: &Array2<f64>) -> Result<Self> {
24        Ok(Self::new(complex::from_amplitude_phase(amplitude, phase)?))
25    }
26
27    pub fn amplitude_only(amplitude: Array2<f64>) -> Result<Self> {
28        if amplitude
29            .as_slice()
30            .iter()
31            .any(|value| !value.is_finite() || *value < 0.0)
32        {
33            return Err(Error::InvalidParameter {
34                name: "amplitude",
35                reason: "values must be finite and non-negative".into(),
36            });
37        }
38        Ok(Self::new(
39            amplitude.map(|&value| Complex64::new(value, 0.0)),
40        ))
41    }
42
43    pub fn phase_only(phase: Array2<f64>) -> Result<Self> {
44        if phase.as_slice().iter().any(|value| !value.is_finite()) {
45            return Err(Error::InvalidParameter {
46                name: "phase",
47                reason: "values must be finite".into(),
48            });
49        }
50        Ok(Self::new(
51            phase.map(|&value| Complex64::from_polar(1.0, value)),
52        ))
53    }
54
55    pub fn from_amplitude_image(path: impl AsRef<Path>) -> Result<Self> {
56        let amplitude = load_grayscale(path, GrayscaleScaling::Unit)?;
57        Self::amplitude_only(amplitude)
58    }
59
60    /// Loads normalized grayscale amplitude and phase images. Phase values map
61    /// linearly from black/white to `-phase_extent`/`+phase_extent` radians.
62    pub fn from_amplitude_phase_images(
63        amplitude_path: impl AsRef<Path>,
64        phase_path: impl AsRef<Path>,
65        phase_extent: f64,
66    ) -> Result<Self> {
67        if !phase_extent.is_finite() || phase_extent <= 0.0 {
68            return Err(Error::InvalidParameter {
69                name: "phase_extent",
70                reason: "must be finite and positive".into(),
71            });
72        }
73        let amplitude = load_grayscale(amplitude_path, GrayscaleScaling::Unit)?;
74        let normalized_phase = load_grayscale(phase_path, GrayscaleScaling::Unit)?;
75        if amplitude.shape() != normalized_phase.shape() {
76            return Err(Error::InvalidShape(format!(
77                "amplitude image shape {:?} differs from phase image shape {:?}",
78                amplitude.shape(),
79                normalized_phase.shape()
80            )));
81        }
82        let phase = normalized_phase.map(|&value| (2.0 * value - 1.0) * phase_extent);
83        Self::from_amplitude_phase(&amplitude, &phase)
84    }
85
86    pub fn constant(shape: (usize, usize), amplitude: f64, phase: f64) -> Result<Self> {
87        validate_shape(shape)?;
88        validate_amplitude(amplitude)?;
89        Ok(Self::new(Array2::filled(
90            shape,
91            Complex64::from_polar(amplitude, phase),
92        )?))
93    }
94
95    pub fn phase_disk(shape: (usize, usize), radius_pixels: f64, phase_shift: f64) -> Result<Self> {
96        validate_shape(shape)?;
97        if !radius_pixels.is_finite() || radius_pixels <= 0.0 || !phase_shift.is_finite() {
98            return Err(Error::InvalidParameter {
99                name: "phase disk",
100                reason: "radius must be positive and phase must be finite".into(),
101            });
102        }
103        let center = ((shape.0 - 1) as f64 / 2.0, (shape.1 - 1) as f64 / 2.0);
104        let mut values = Vec::with_capacity(shape.0 * shape.1);
105        for row in 0..shape.0 {
106            for column in 0..shape.1 {
107                let radius = (row as f64 - center.0).hypot(column as f64 - center.1);
108                values.push(Complex64::from_polar(
109                    1.0,
110                    if radius <= radius_pixels {
111                        phase_shift
112                    } else {
113                        0.0
114                    },
115                ));
116            }
117        }
118        Ok(Self::new(Array2::from_vec(shape, values)?))
119    }
120
121    pub fn siemens_star(shape: (usize, usize), spokes: usize) -> Result<Self> {
122        validate_shape(shape)?;
123        if spokes < 2 {
124            return Err(Error::InvalidParameter {
125                name: "spokes",
126                reason: "must be at least 2".into(),
127            });
128        }
129        let center = ((shape.0 - 1) as f64 / 2.0, (shape.1 - 1) as f64 / 2.0);
130        let maximum_radius = shape.0.min(shape.1) as f64 * 0.46;
131        let mut values = Vec::with_capacity(shape.0 * shape.1);
132        for row in 0..shape.0 {
133            for column in 0..shape.1 {
134                let y = row as f64 - center.0;
135                let x = column as f64 - center.1;
136                let radius = x.hypot(y);
137                let amplitude =
138                    if radius <= maximum_radius && ((x.atan2(y) * spokes as f64).sin() >= 0.0) {
139                        0.25
140                    } else {
141                        1.0
142                    };
143                values.push(Complex64::new(amplitude, 0.0));
144            }
145        }
146        Ok(Self::new(Array2::from_vec(shape, values)?))
147    }
148
149    pub fn resolution_target(shape: (usize, usize)) -> Result<Self> {
150        validate_shape(shape)?;
151        let mut values = vec![Complex64::new(1.0, 0.0); shape.0 * shape.1];
152        let groups = [2_usize, 3, 4, 6, 8];
153        for (group, &period) in groups.iter().enumerate() {
154            let top = group * shape.0 / groups.len();
155            let bottom = (group + 1) * shape.0 / groups.len();
156            for row in top..bottom {
157                for column in 0..shape.1 {
158                    let dark = if group.is_multiple_of(2) {
159                        (column / period).is_multiple_of(2)
160                    } else {
161                        ((row - top) / period).is_multiple_of(2)
162                    };
163                    if dark {
164                        values[row * shape.1 + column] = Complex64::new(0.2, 0.0);
165                    }
166                }
167            }
168        }
169        Ok(Self::new(Array2::from_vec(shape, values)?))
170    }
171
172    pub fn random_phase(shape: (usize, usize), standard_deviation: f64, seed: u64) -> Result<Self> {
173        validate_shape(shape)?;
174        if !standard_deviation.is_finite() || standard_deviation < 0.0 {
175            return Err(Error::InvalidParameter {
176                name: "standard_deviation",
177                reason: "must be finite and non-negative".into(),
178            });
179        }
180        let mut rng = StdRng::seed_from_u64(seed);
181        let values = if standard_deviation == 0.0 {
182            vec![Complex64::new(1.0, 0.0); shape.0 * shape.1]
183        } else {
184            let distribution =
185                Normal::new(0.0, standard_deviation).map_err(|error| Error::InvalidParameter {
186                    name: "standard_deviation",
187                    reason: error.to_string(),
188                })?;
189            (0..shape.0 * shape.1)
190                .map(|_| Complex64::from_polar(1.0, distribution.sample(&mut rng)))
191                .collect()
192        };
193        Ok(Self::new(Array2::from_vec(shape, values)?))
194    }
195
196    pub fn particle_field(shape: (usize, usize), particles: usize, seed: u64) -> Result<Self> {
197        validate_shape(shape)?;
198        let mut values = vec![Complex64::new(1.0, 0.0); shape.0 * shape.1];
199        let mut rng = StdRng::seed_from_u64(seed);
200        for _ in 0..particles {
201            let row = rng.random_range(0..shape.0);
202            let column = rng.random_range(0..shape.1);
203            values[row * shape.1 + column] = Complex64::new(0.1, 0.0);
204        }
205        Ok(Self::new(Array2::from_vec(shape, values)?))
206    }
207
208    /// A deterministic mixed amplitude/phase target useful for smoke tests.
209    pub fn mixed_test_pattern(shape: (usize, usize)) -> Result<Self> {
210        validate_shape(shape)?;
211        let center = ((shape.0 - 1) as f64 / 2.0, (shape.1 - 1) as f64 / 2.0);
212        let scale = shape.0.min(shape.1) as f64;
213        let mut values = Vec::with_capacity(shape.0 * shape.1);
214        for row in 0..shape.0 {
215            for column in 0..shape.1 {
216                let y = row as f64 - center.0;
217                let x = column as f64 - center.1;
218                let radius = x.hypot(y);
219                let amplitude = if (x.abs() < 0.12 * scale && y.abs() < 0.35 * scale)
220                    || radius < 0.12 * scale
221                {
222                    0.45
223                } else {
224                    1.0
225                };
226                let phase = 0.8 * (-radius * radius / (0.08 * scale * scale)).exp()
227                    + 0.25 * (std::f64::consts::TAU * x / scale).sin();
228                values.push(Complex64::from_polar(amplitude, phase));
229            }
230        }
231        Ok(Self::new(Array2::from_vec(shape, values)?))
232    }
233
234    /// Smooth random phase and absorption blobs approximating a weak biological
235    /// specimen. This is a test object, not a tissue-specific physical model.
236    pub fn biological_like(shape: (usize, usize), features: usize, seed: u64) -> Result<Self> {
237        validate_shape(shape)?;
238        if features == 0 {
239            return Err(Error::InvalidParameter {
240                name: "features",
241                reason: "must be greater than zero".into(),
242            });
243        }
244        let mut rng = StdRng::seed_from_u64(seed);
245        let minimum_size = shape.0.min(shape.1) as f64;
246        let lower_sigma = (0.02 * minimum_size).max(1.0);
247        let upper_sigma = (0.12 * minimum_size).max(lower_sigma + 0.1);
248        let blobs: Vec<_> = (0..features)
249            .map(|_| {
250                (
251                    rng.random_range(0.0..shape.0 as f64),
252                    rng.random_range(0.0..shape.1 as f64),
253                    rng.random_range(lower_sigma..upper_sigma),
254                    rng.random_range(-0.4..1.0),
255                    rng.random_range(0.0..0.25),
256                )
257            })
258            .collect();
259        let mut values = Vec::with_capacity(shape.0 * shape.1);
260        for row in 0..shape.0 {
261            for column in 0..shape.1 {
262                let mut phase = 0.0;
263                let mut absorption = 0.0;
264                for &(center_row, center_column, sigma, phase_strength, absorption_strength) in
265                    &blobs
266                {
267                    let squared_radius =
268                        (row as f64 - center_row).powi(2) + (column as f64 - center_column).powi(2);
269                    let profile = (-squared_radius / (2.0 * sigma * sigma)).exp();
270                    phase += phase_strength * profile;
271                    absorption += absorption_strength * profile;
272                }
273                values.push(Complex64::from_polar(
274                    (1.0_f64 - absorption).max(0.1),
275                    phase,
276                ));
277            }
278        }
279        Ok(Self::new(Array2::from_vec(shape, values)?))
280    }
281
282    pub fn shape(&self) -> (usize, usize) {
283        self.field.shape()
284    }
285}
286
287impl From<Array2<Complex64>> for SyntheticObject {
288    fn from(field: Array2<Complex64>) -> Self {
289        Self::new(field)
290    }
291}
292
293fn validate_amplitude(amplitude: f64) -> Result<()> {
294    if !amplitude.is_finite() || amplitude < 0.0 {
295        Err(Error::InvalidParameter {
296            name: "amplitude",
297            reason: "must be finite and non-negative".into(),
298        })
299    } else {
300        Ok(())
301    }
302}
303
304fn validate_shape(shape: (usize, usize)) -> Result<()> {
305    if shape.0 == 0 || shape.1 == 0 {
306        Err(Error::InvalidShape(format!(
307            "synthetic object dimensions must be non-zero, got {shape:?}"
308        )))
309    } else {
310        Ok(())
311    }
312}
fpm_rs/simulation/synthetic_object.rs

fpm_rs/simulation/
synthetic_object.rs
