gmt_dos_clients_crseo/pyramid/
processing.rs

1use std::ops::{Add, DivAssign, Mul, Sub};
2
3use crate::Processing;
4
5use super::{PyramidData, PyramidProcessor};
6
7impl<T> Processing for PyramidProcessor<T>
8where
9    T: std::fmt::Debug
10        + TryFrom<f32>
11        + Copy
12        + Add<Output = T>
13        + Sub<Output = T>
14        + DivAssign
15        + Mul<Output = T>
16        + PartialOrd
17        + 'static,
18    for<'a> &'a T: Sub<Output = T> + Add<Output = T>,
19    <T as TryFrom<f32>>::Error: std::fmt::Debug,
20{
21    type ProcessorData = PyramidData<T>;
22
23    fn processing(&self) -> Self::ProcessorData {
24        let (n, m) = self.frame.resolution;
25        assert!(n > 0 && m > 0, "the detector frame resolution is null");
26
27        let crseo::wavefrontsensor::LensletArray { n_side_lenslet, .. } = self.lenslet_array;
28        let n0 = n_side_lenslet / 2;
29        let n1 = n0 + n / 2;
30
31        /*
32        The detector frame of the pyramid is divided in 4 quadrants:
33        | I1 I2 |
34        | I3 I4 |
35         */
36
37        /* left =
38        | I1 |
39        | I3 |
40         */
41        let left = self.frame.value.chunks(n).skip(n0).take(n_side_lenslet);
42        /* right =
43        | I2 |
44        | I4 |
45         */
46        let right = self.frame.value.chunks(n).skip(n1).take(n_side_lenslet);
47        // (left - right, left + right)
48        let row_diff_sum: Vec<_> = left
49            .zip(right)
50            .flat_map(|(left, right)| {
51                left.iter()
52                    .zip(right)
53                    .map(|(left, right)| (left - right, left + right))
54                    .collect::<Vec<_>>()
55            })
56            .collect();
57        // top = ( I1 - I2 , I1 + I2 )
58        let top = (0..n)
59            .map(|i| row_diff_sum.iter().skip(i).step_by(m).collect::<Vec<_>>())
60            .skip(n0)
61            .take(n_side_lenslet);
62        // bottom = ( I3 -I4 , I3 + I4)
63        let bottom = (0..n)
64            .map(|i| row_diff_sum.iter().skip(i).step_by(m).collect::<Vec<_>>())
65            .skip(n1)
66            .take(n_side_lenslet);
67        /*
68        top - bottom = I1 - I2 + I3 - I4
69        top + bottom = I1 + I2 + I3 + I4
70         */
71        let (row_col_data, flux): (Vec<_>, Vec<_>) = top
72            .zip(bottom)
73            .flat_map(|(top, bottom)| {
74                top.iter()
75                    .zip(bottom)
76                    .map(|((top_diff, top_sum), (bottom_diff, bottom_sum))| {
77                        (*top_diff + *bottom_diff, *top_sum + *bottom_sum)
78                    })
79                    .collect::<Vec<_>>()
80            })
81            .unzip();
82
83        // top = | I1 I2 |
84        let top = (0..n)
85            .map(|i| {
86                self.frame
87                    .value
88                    .iter()
89                    .skip(i)
90                    .step_by(m)
91                    .collect::<Vec<_>>()
92            })
93            .skip(n0)
94            .take(n_side_lenslet);
95        // bottom = | I3 I4 |
96        let bottom = (0..n)
97            .map(|i| {
98                self.frame
99                    .value
100                    .iter()
101                    .skip(i)
102                    .step_by(m)
103                    .collect::<Vec<_>>()
104            })
105            .skip(n1)
106            .take(n_side_lenslet);
107        // top - bottom
108        let col_diff: Vec<_> = top
109            .zip(bottom)
110            .flat_map(|(top, bottom)| {
111                top.iter()
112                    .zip(bottom)
113                    .map(|(&&top, &bottom)| top - bottom)
114                    .collect::<Vec<_>>()
115            })
116            .collect();
117        // I1 - I3 + I2 - I4
118        let left = col_diff.chunks(n);
119        let right = col_diff.chunks(n);
120        let col_row_data: Vec<_> = left
121            .zip(right)
122            .flat_map(|(left, right)| {
123                left.iter()
124                    .skip(n0)
125                    .take(n_side_lenslet)
126                    .zip(right.iter().skip(n1).take(n_side_lenslet))
127                    .map(|(&left, &right)| left + right)
128                    .collect::<Vec<_>>()
129            })
130            .collect();
131
132        let mut _flux = flux.clone();
133        _flux.sort_by(|a, b| a.partial_cmp(b).unwrap());
134        let med_flux = match _flux.len() {
135            n if n % 2 == 0 => T::try_from(0.5).unwrap() * (_flux[n / 2] + _flux[1 + n / 2]),
136            n => _flux[(n + 1) / 2],
137        };
138
139        let (sx, sy): (Vec<_>, Vec<_>) = row_col_data
140            .into_iter()
141            .zip(col_row_data.into_iter())
142            .map(|(mut sx, mut sy)| {
143                sx /= med_flux;
144                sy /= med_flux;
145                (sx, sy)
146            })
147            .unzip();
148
149        /*         let (sx, sy): (Vec<_>, Vec<_>) = row_col_data
150        .into_iter()
151        .zip(col_row_data.into_iter())
152        .zip(flux.iter())
153        .map(|((mut sx, mut sy), flux)| {
154            sx /= *flux;
155            sy /= *flux;
156            (sx, sy)
157        })
158        .unzip(); */
159
160        PyramidData { sx, sy, flux }
161    }
162}
gmt_dos_clients_crseo/pyramid/processing.rs

gmt_dos_clients_crseo/pyramid/
processing.rs
