fpm_rs/metrics/complex_field/
statistics.rs

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

fpm_rs/metrics/complex_field/
statistics.rs
