fpm_rs/metrics/complex_field/
statistics.rs

1//! Metrics calculated from one complex-valued field.
2
3use ndarray::ArrayView2;
4use num_complex::Complex64;
5
6use crate::{
7    Result,
8    backend::{Backend, CpuBackend, FftDirection},
9};
10
11/// Azimuthally averaged power in annular bins of a two-dimensional Fourier transform.
12///
13/// [`radial_fourier_spectrum`] uses integer-radius bins centered on the DC
14/// component. `radius_px[b]` is the radius, in Fourier-grid pixels, of bin
15/// `b`; `power[b]` is the mean normalized power in that bin; and
16/// `sample_count[b]` is the number of Fourier samples contributing to it.
17///
18/// Each underlying Fourier-pixel power is normalized so the sum before radial
19/// averaging equals the input field's squared norm (Parseval's identity). No
20/// physical-frequency calibration is applied: supplying that requires the
21/// image sampling pitch.
22#[derive(Clone, Debug, PartialEq)]
23pub struct RadialFourierSpectrum {
24    /// Mean annular radius in Fourier-grid pixels for each non-empty bin.
25    pub radius_px: Vec<f64>,
26    /// Mean Fourier power per bin, normalized by total field element count squared.
27    pub power: Vec<f64>,
28    /// Number of Fourier pixels accumulated into each bin.
29    pub sample_count: Vec<usize>,
30}
31
32/// Calculate the radial Fourier power spectrum of a two-dimensional complex field.
33///
34/// This applies a forward FFT, treats DC as the centre of the Fourier grid, and
35/// averages normalized power over annular bins with integer Fourier-pixel radius.
36pub fn radial_fourier_spectrum(field: ArrayView2<'_, Complex64>) -> Result<RadialFourierSpectrum> {
37    let shape = field.dim();
38    let mut spectrum: Vec<_> = field.iter().copied().collect();
39    let backend = CpuBackend::new(shape, shape)?;
40    let mut column_scratch = vec![Complex64::default(); shape.0];
41    backend.fft2(
42        &mut spectrum,
43        shape,
44        FftDirection::Forward,
45        &mut column_scratch,
46    )?;
47
48    let center_row = shape.0 / 2;
49    let center_column = shape.1 / 2;
50    let max_radius = ((center_row as f64).hypot(center_column as f64)).floor() as usize;
51    let mut power = vec![0.0; max_radius + 1];
52    let mut sample_count = vec![0_usize; max_radius + 1];
53    let normalization = field.len() as f64;
54
55    for row in 0..shape.0 {
56        let frequency_row = if row <= shape.0 / 2 {
57            row as isize
58        } else {
59            row as isize - shape.0 as isize
60        };
61        for column in 0..shape.1 {
62            let frequency_column = if column <= shape.1 / 2 {
63                column as isize
64            } else {
65                column as isize - shape.1 as isize
66            };
67            let radius = ((frequency_row * frequency_row + frequency_column * frequency_column)
68                as f64)
69                .sqrt()
70                .floor() as usize;
71            power[radius] += spectrum[row * shape.1 + column].norm_sqr() * normalization;
72            sample_count[radius] += 1;
73        }
74    }
75
76    for (value, count) in power.iter_mut().zip(&sample_count) {
77        *value /= *count as f64;
78    }
79
80    Ok(RadialFourierSpectrum {
81        radius_px: (0..power.len()).map(|bin| bin as f64).collect(),
82        power,
83        sample_count,
84    })
85}
86
87#[cfg(test)]
88mod tests {
89    use super::*;
90    use ndarray::Array2;
91
92    #[test]
93    fn constant_field_has_power_only_at_dc() {
94        let field = Array2::from_elem((4, 4), Complex64::new(1.0, 0.0));
95
96        let spectrum = radial_fourier_spectrum(field.view()).unwrap();
97
98        assert_eq!(spectrum.radius_px, vec![0.0, 1.0, 2.0]);
99        assert_eq!(spectrum.sample_count, vec![1, 8, 7]);
100        assert!((spectrum.power[0] - 16.0).abs() < 1e-12);
101        assert!(spectrum.power[1..].iter().all(|value| value.abs() < 1e-12));
102    }
103}
fpm_rs/metrics/complex_field/statistics.rs

fpm_rs/metrics/complex_field/
statistics.rs
