rusty_psf/core/utils.rs
1use ndarray::{Array1, Array2, Array3, Axis};
2use crate::types::OpticalParameters;
3use num_traits::Float;
4
5/// Generate a centered coordinate grid.
6///
7/// Creates a 1D array of coordinates centered around zero with specified size and pixel spacing.
8///
9/// # Arguments
10/// * `size` - Number of points in the grid
11/// * `pixel_size` - Physical size of each pixel
12///
13/// # Returns
14/// A 1D array of evenly spaced coordinates
15pub fn make_grid<T: Float>(size: usize, pixel_size: T) -> Array1<T> {
16 let half = T::from(size).unwrap() / (T::from(2).unwrap());
17 let range = Array1::linspace(-half * pixel_size, (half - T::one()) * pixel_size, size);
18 range
19}
20
21/// Generate a radial coordinate grid.
22///
23/// Creates a 2D array of radial coordinates (distance from center) for each point.
24///
25/// # Arguments
26/// * `size` - Size of the grid (will be size × size)
27/// * `pixel_size` - Physical size of each pixel
28///
29/// # Returns
30/// A 2D array containing the radial distance from the center for each point
31pub fn make_radial_grid(size: usize, pixel_size: f64) -> Array2<f64> {
32 let coords = make_grid(size, pixel_size);
33 let mut radial = Array2::zeros((size, size));
34
35 for i in 0..size {
36 for j in 0..size {
37 radial[[i, j]] = (coords[i].powi(2) + coords[j].powi(2)).sqrt();
38 }
39 }
40 radial
41}
42
43/// Normalize PSF to unit sum.
44///
45/// Normalizes a PSF array so that the sum of all elements equals 1.0.
46/// This ensures energy conservation in the PSF.
47///
48/// # Arguments
49/// * `psf` - 3D array containing PSF values
50///
51/// # Returns
52/// The normalized PSF array
53pub fn normalize_psf(mut psf: Array3<f64>) -> Array3<f64> {
54 let sum = psf.sum();
55 if sum != 0.0 {
56 psf /= sum;
57 }
58 psf
59}
60
61/// Calculate the optimal grid size based on optical parameters.
62///
63/// Determines the optimal grid size for PSF calculations based on the optical
64/// parameters and desired oversampling factor. The size is rounded up to the
65/// next power of 2 for efficient FFT calculations.
66///
67/// # Arguments
68/// * `params` - Optical parameters of the system
69/// * `oversampling` - Oversampling factor for the grid
70///
71/// # Returns
72/// The optimal grid size (power of 2)
73pub fn optimal_grid_size(params: &OpticalParameters, oversampling: f64) -> usize {
74 let min_size = (2.0 * params.na / params.wavelength * oversampling).ceil() as usize;
75 // Round up to next power of 2 for FFT efficiency
76 let mut size = 1;
77 while size < min_size {
78 size *= 2;
79 }
80 size
81}
82
83/// Calculate the center of mass for a 3D array.
84///
85/// Computes the center of mass coordinates for each dimension of a 3D array.
86/// This is useful for finding the central position of a PSF.
87///
88/// # Arguments
89/// * `arr` - 3D array to analyze
90///
91/// # Returns
92/// Array of [x, y, z] coordinates of the center of mass
93pub fn center_of_mass(arr: &Array3<f64>) -> [f64; 3] {
94 let total_mass = arr.sum();
95 if total_mass == 0.0 {
96 return [0.0, 0.0, 0.0];
97 }
98
99 let dims = arr.shape();
100 let mut com = [0.0; 3];
101
102 for i in 0..3 {
103 let axis_coords: Array1<f64> = Array1::linspace(0.0, (dims[i] - 1) as f64, dims[i]);
104 let sum = arr.sum_axis(Axis(i));
105 com[i] = (axis_coords * sum).sum() / total_mass;
106 }
107
108 com
109}