Coverage for pyhiperta/waveform_indexing.py: 100%
43 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-07 20:49 +0000
1"""Implements the indexing/mapping of waveforms from 1D array (hexagonal structure) to 2D square array.
3The 1D array representation is suited for vectorization of "per-element" operations such as calibration,
4Hillas computation, etc. while the 2D square representation is faster for operations requiring neighbors
5look-up, such as the tail-cut cleaning. It is faster to copy the data from one representation to the other
6than to perform the cleaning with the 1D representation.
7"""
9from typing import Tuple
11import numpy as np
12from scipy.spatial import cKDTree
15def _get_neighbors_indices(pixel_positions: np.ndarray) -> np.ndarray:
16 """Get the indices of the pixel neighbors in the array, for each pixel.
18 Parameters
19 ----------
20 pixel_positions : np.ndarray
21 Positions of the pixels in the camera. 1st dimension is the "x" axis, 2nd dimension is "y".
22 Shape: (2, number_of_pixels)
24 Notes
25 -----
26 The returned array has shape (Number_of_pixels, 6) even though some pixels have less than 6 neighbors. When that
27 is the case, the "extra-neighbors" have indices greater or equal than the number of pixels (invalid indices!)
29 Warnings
30 --------
31 This function is only intended to work for hexagonal camera pixel positions. It assumes a pixel has at most
32 6 neighbors and the neigbors are placed on a hexagonal grid.
34 Returns
35 -------
36 np.ndarray
37 Indices of the pixels neighbors, per pixel. Shape: (Number_of_pixels, 6)
38 """
40 dist_kdtree = cKDTree(pixel_positions.transpose())
42 # Query the tree for the 1st neighbours of each pixel to get the distance between 1st neighbours
43 # We query the tree for neighbour 2, because we query the tree for points identical to the tree points,
44 # so the first neighbour will have distance 0 and be the pixels themselves.
45 pixel_distances, _ = dist_kdtree.query(dist_kdtree.data, k=list(range(2, 3)), p=2, workers=-1)
46 dist_1st_neighbour = pixel_distances.min()
48 # Now query the tree for all neighbours at most (1 + sqrt(3)) / 2 to get all 1st neighbours
49 # We query for 6 neighbours, but not all pixels have 6 neighbours. When that is the case the
50 # neighbours indices of the are set to a value > the number of pixels.
51 _, neighbors_indices = dist_kdtree.query(
52 pixel_positions.transpose(),
53 k=list(range(2, 8)), # get 2nd to 7th neighbours
54 distance_upper_bound=(1.0 + np.sqrt(3)) / 2.0 * dist_1st_neighbour,
55 p=2, # 2 norm
56 workers=-1, # use all cpus
57 )
58 return neighbors_indices
61def waveform_1D_to_2D_maps(pixel_positions: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
62 """Create the 2D waveform array with "squared geometry" and 1D to 2D indices maps.
64 To perform fast "neighbour looking" operations such as convolution we need to store the waveform
65 as a 2D array. To fit the hexagonal geometry of the camera into a 2D square grid (the 2D array)
66 we perform a rotation (to align an axis with the 2D array axis) and then deform the other axis to
67 align the hexagonal neighbour on the square grid.
68 If we represent the pixels by their index in the 1D array:
69 0 1 1 4
70 2 3 4 -> 0 3 6
71 5 6 2 5
72 The algorithm is:
73 - find 3 neighbours forming a triangle in the hexagonal representation.
74 The pixels should be such that: the left-most pixel is in between the other
75 two pixels on the y coordinate, eg:
76 x
77 x
78 x
79 Such a combination is always possible as the angle in the triangles is 60° and
80 we search a space of 90°. The goal of this is to use the left-most pixel as our
81 2D grid origin, and the other pixels coordinates will give us the rotation and
82 deformation.
83 In the example, the chosen are pixels are 0, 1 and 3. (1's y is equal or higher than 0's y)
84 - Rotate all pixels coordinates such that the "down" neighbour is aligned with the 2D
85 grid rows
86 x x
87 x ->
88 x x x
89 - Normalize the distances between pixels, so that a pixel is at distance 1 from another
90 on row and column axis (makes next step easier).
91 - "Deform" along the row axis to align the pixels on the 2D columns. This is a linear
92 transformation of the x coordinate with respect to the y coordinate:
93 ___x x ___ is the deformation to apply for a distance of +1 y
94 ->
95 x x x x
96 - Bin the pixel positions to get their index in the 2D array. This gives the 1D to 2D
97 indices maps.
99 Notes
100 -----
101 This function is slow and runs only on cpu (scipy.cKDTree). It should be used only once to
102 get the mappings for a given camera geometry.
104 Parameters
105 ----------
106 pixel_positions : np.array
107 2D array (shape: (2, N_pixels)) containing the pixels x and y coordinates. This should
108 be the actual pixels positions in the camera. It is important that the referential in
109 which the position are measured has equal metrics on both axis (1m in the real world
110 is 1 on x _and_ y axis)
112 Returns
113 -------
114 Tuple[np.ndarray, np.ndarray]
115 Row and column indices (shape (N_pixels,) each) to use as 1D to 2D indexing maps.
116 """
118 neighbors_indices = _get_neighbors_indices(pixel_positions)
120 # find three points that are:
121 # - all neighbours of each other
122 # - the 1 point is the left most point (minimal x)
123 # - the 2nd point is the higher point (max y)
124 # - the 3rd point is the lower (min y)
125 def find_neighbours_trio():
126 # Try all pixel as the left-most pixel
127 for idx_left_most_p, left_most_p_neighbours in enumerate(neighbors_indices):
128 # for all neighbours to left-most pixel
129 for idx_up_neighbour in left_most_p_neighbours:
130 if (
131 idx_up_neighbour < pixel_positions.shape[1] # valid neighbour
132 and pixel_positions[0, idx_left_most_p] <= pixel_positions[0, idx_up_neighbour]
133 and pixel_positions[1, idx_left_most_p] <= pixel_positions[1, idx_up_neighbour]
134 ): # We have a good higher point. Now search for a most down point.
135 for idx_down_neighbour in neighbors_indices[idx_up_neighbour]:
136 if (
137 idx_down_neighbour < pixel_positions.shape[1] # valid neighbour
138 and idx_down_neighbour in left_most_p_neighbours # also a neighbour of left-most
139 and pixel_positions[0, idx_left_most_p] <= pixel_positions[0, idx_down_neighbour]
140 and pixel_positions[1, idx_left_most_p] >= pixel_positions[1, idx_down_neighbour]
141 ):
142 return (idx_left_most_p, idx_up_neighbour, idx_down_neighbour)
143 raise ValueError("Could not find trio of points matching conditions.")
145 idx_left_most, idx_up_neighbour, idx_down_neighbour = find_neighbours_trio()
147 # Get the angle required to align the left-most - down-most points line to the x axis
148 # That is the angle between the points line and the x-axis (not the reverse ! orientation is important)
149 #
150 # x left-most - x down-most
151 # .___________________________
152 # |
153 # | y left-most - y down-most
154 # |
155 # .
156 psi = np.arctan2(
157 pixel_positions[1, idx_left_most] - pixel_positions[1, idx_down_neighbour],
158 pixel_positions[0, idx_down_neighbour] - pixel_positions[0, idx_left_most],
159 )
160 # Apply rotation
161 pixel_positions = np.array([[np.cos(psi), -np.sin(psi)], [np.sin(psi), np.cos(psi)]]) @ pixel_positions
163 # Normalize the neighbours distance for both axis (make neighbours distance 1 on x and y)
164 pixel_positions[0, :] *= 1.0 / (pixel_positions[0, idx_down_neighbour] - pixel_positions[0, idx_left_most])
165 pixel_positions[1, :] *= 1.0 / (pixel_positions[1, idx_up_neighbour] - pixel_positions[1, idx_left_most])
166 # Apply linear deformation to the x coordinate based on y coordinate
167 pixel_positions[0, :] -= (pixel_positions[1, :] - pixel_positions[1, idx_left_most]) * (
168 pixel_positions[0, idx_up_neighbour] - pixel_positions[0, idx_left_most]
169 )
171 # Bin the pixel coordinates to find their index in 2D array
172 # After the normalization and deformation the pixels are placed at integer values on the grid
173 # To make sure the points are not near the bins edge, make the bins shifted by 0.5
174 # The values are then floored and converted to int by the astype
175 # Note: - what binning x would do: (x - (min(x) - shift)) / bin_size (with integer division)
176 # - bin size is 1.0, so no need to divide.
177 # - row indices are reversed, because the position in the matrix start in the top left, not
178 # bottom left as the y coordinates do.
179 row_indices = (pixel_positions[1] - (pixel_positions[1].min() - 0.5)).astype(np.int32)
180 row_indices = row_indices.max() - row_indices
181 col_indices = (pixel_positions[0] - (pixel_positions[0].min() - 0.5)).astype(np.int32)
183 return (row_indices, col_indices)
186def waveform1Dto2D(
187 maps1D_to_2D: Tuple[np.ndarray, np.ndarray], waveform: np.ndarray, waveform2D: np.ndarray = None
188) -> np.ndarray:
189 """Get the 2D array representation of waveform
191 If waveform2D is passed, the values of waveform2D will be updated with `waveform`'s.
192 Otherwise a new 2D array will be created.
194 Parameters
195 ----------
196 maps1D_to_2D : tuple of np.ndarray
197 The tuple holding the indices to go from 1D to 2D representation
198 waveform : np.ndarray
199 Waveform values in 1D shape.
200 waveform2D : np.ndarray, optionnal
201 An already existing 2D representation of the waveform, to update with `waveform`'s values.
203 Returns
204 -------
205 np.ndarray
206 A 2D array representing the waveform. (The array is independent of waveform, not a view)
207 """
208 # Initialize waveform2D it with 0s if needed.
209 if waveform2D is None:
210 waveform2D = np.zeros(
211 waveform.shape[:-1] + (maps1D_to_2D[0].max() + 1, maps1D_to_2D[1].max() + 1), dtype=waveform.dtype
212 )
214 waveform2D[..., *maps1D_to_2D] = waveform
215 return waveform2D
218def waveform2Dto1D(maps1D_to_2D: np.ndarray, waveform2D: np.ndarray) -> np.ndarray:
219 """Get the 1D array representing the waveform from the 2D representation `waveform2D`.
221 Notes
222 -----
223 The result is a view in `waveform2D`! Modifying the view will modify waveform2D!
225 Parameters
226 ----------
227 maps1D_to_2D : tuple of np.ndarray
228 The tuple holding the indices to go from 2D to 1D representation
229 waveform2D : np.ndarray
230 The 2D array holding the waveform values.
232 Returns
233 -------
234 np.ndarray
235 1D array holding the waveform values.
236 """
237 return np.ascontiguousarray(waveform2D[..., *maps1D_to_2D])
240def neighbors_only_stencil() -> np.ndarray:
241 """Return the stencil selecting __only the neighbors__ of a pixel in 2D array.
243 The pixel itself is not included: center element of the stencil is 0.
245 Notes
246 -----
247 The stencil values depend on the way the 1D to 2D conversion is performed. Depending on how the axis are deformed,
248 it could be transposed or reversed. The current implementation does
249 0 1 1 4
250 2 3 4 -> 0 3 6
251 5 6 2 5
252 so the stencil is 1 where the neighbors are and 0 otherwise (and 0 at the center to exclude current pixel).
254 """
255 return np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]], dtype=bool)
258def leakage_pixel_masks(pixel_positions: np.ndarray) -> np.ndarray:
259 """Returns a mask (boolean array) that multiplied with waveforms (1D) puts all pixels except the outer rings to 0.
261 Parameters
262 ----------
263 pixel_positions : np.ndarray
264 2D array (shape: (2, N_pixels)) containing the pixels x and y coordinates. This should
265 be the actual pixels positions in the camera. It is important that the referential in
266 which the position are measured has equal metrics on both axis (1m in the real world
267 is 1 on x _and_ y axis)
269 Returns
270 -------
271 leakage_mask: np.ndarray
272 Mask array with values True at the indices of pixels on the outermost rings, False otherwise.
273 The shape is (2, nb_pixels). leakage_mask[0, :] selects the pixels of the outermost ring (leakage 1)
274 while leakage[1:] selects the pixels of the 1st and 2nd outermost ring (leakage 2).
275 """
276 neighbors_indices = _get_neighbors_indices(pixel_positions)
277 # leakage 1 pixels are the one of the borders, so they have less than 6 neighbors.
278 # we use the fact that the neighbors_indices array has invalid indices for non-existing neighbors.
279 leakage_1_mask = (neighbors_indices >= pixel_positions.shape[1]).any(axis=1)
280 # leakage 2 pixels are the ones that are neighbors of leakage 1 pixels, but are not in leakage 1 themselves.
281 # find them by getting the indices of leakage 1 pixels with nonzero, then use isin to see if there are neighbors
282 # finally exclude the leakage 1 pixels
283 leakage_2_mask = np.isin(neighbors_indices, np.nonzero(leakage_1_mask)[0]).any(axis=1) & np.bitwise_not(
284 leakage_1_mask
285 )
286 return np.stack([leakage_1_mask, leakage_1_mask | leakage_2_mask], axis=0)