Coverage for src/pyhiperta/waveform_indexing.py: 100%
42 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-07 14:41 +0000
1# Copyright 2026 CNRS
2# This software is distributed under the terms of the CeCILL-C free software license.
4"""Implements the indexing/mapping of waveforms from 1D array (hexagonal structure) to 2D square array.
6The 1D array representation is suited for vectorization of "per-element" operations such as calibration,
7Hillas computation, etc. while the 2D square representation is faster for operations requiring neighbors
8look-up, such as the tail-cut cleaning. It is faster to copy the data from one representation to the other
9than to perform the cleaning with the 1D representation.
10"""
13import numpy as np
14from scipy.spatial import cKDTree
17def _get_neighbors_indices(pixel_positions: np.ndarray) -> np.ndarray:
18 """Get the indices of the pixel neighbors in the array, for each pixel.
20 Parameters
21 ----------
22 pixel_positions : np.ndarray
23 Positions of the pixels in the camera. 1st dimension is the "x" axis, 2nd dimension is "y".
24 Shape: (2, number_of_pixels)
26 Notes
27 -----
28 The returned array has shape (Number_of_pixels, 6) even though some pixels have less than 6 neighbors. When that
29 is the case, the "extra-neighbors" have indices greater or equal than the number of pixels (invalid indices!)
31 Warnings
32 --------
33 This function is only intended to work for hexagonal camera pixel positions. It assumes a pixel has at most
34 6 neighbors and the neigbors are placed on a hexagonal grid.
36 Returns
37 -------
38 np.ndarray
39 Indices of the pixels neighbors, per pixel. Shape: (Number_of_pixels, 6)
40 """
41 dist_kdtree = cKDTree(pixel_positions.transpose())
43 # Query the tree for the 1st neighbours of each pixel to get the distance between 1st neighbours
44 # We query the tree for neighbour 2, because we query the tree for points identical to the tree points,
45 # so the first neighbour will have distance 0 and be the pixels themselves.
46 pixel_distances, _ = dist_kdtree.query(dist_kdtree.data, k=list(range(2, 3)), p=2, workers=-1)
47 dist_1st_neighbour = pixel_distances.min()
49 # Now query the tree for all neighbours at most (1 + sqrt(3)) / 2 to get all 1st neighbours
50 # We query for 6 neighbours, but not all pixels have 6 neighbours. When that is the case the
51 # neighbours indices of the are set to a value > the number of pixels.
52 _, neighbors_indices = dist_kdtree.query(
53 pixel_positions.transpose(),
54 k=list(range(2, 8)), # get 2nd to 7th neighbours
55 distance_upper_bound=(1.0 + np.sqrt(3)) / 2.0 * dist_1st_neighbour,
56 p=2, # 2 norm
57 workers=-1, # use all cpus
58 )
59 return neighbors_indices
62def waveform_1D_to_2D_maps(pixel_positions: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
63 """Create the 2D waveform array with "squared geometry" and 1D to 2D indices maps.
65 To perform fast "neighbour looking" operations such as convolution we need to store the waveform
66 as a 2D array. To fit the hexagonal geometry of the camera into a 2D square grid (the 2D array)
67 we perform a rotation (to align an axis with the 2D array axis) and then deform the other axis to
68 align the hexagonal neighbour on the square grid.
69 If we represent the pixels by their index in the 1D array:
70 0 1 1 4
71 2 3 4 -> 0 3 6
72 5 6 2 5
73 The algorithm is:
74 - find 3 neighbours forming a triangle in the hexagonal representation.
75 The pixels should be such that: the left-most pixel is in between the other
76 two pixels on the y coordinate, eg:
77 x
78 x
79 x
80 Such a combination is always possible as the angle in the triangles is 60° and
81 we search a space of 90°. The goal of this is to use the left-most pixel as our
82 2D grid origin, and the other pixels coordinates will give us the rotation and
83 deformation.
84 In the example, the chosen are pixels are 0, 1 and 3. (1's y is equal or higher than 0's y)
85 - Rotate all pixels coordinates such that the "down" neighbour is aligned with the 2D
86 grid rows
87 x x
88 x ->
89 x x x
90 - Normalize the distances between pixels, so that a pixel is at distance 1 from another
91 on row and column axis (makes next step easier).
92 - "Deform" along the row axis to align the pixels on the 2D columns. This is a linear
93 transformation of the x coordinate with respect to the y coordinate:
94 ___x x ___ is the deformation to apply for a distance of +1 y
95 ->
96 x x x x
97 - Bin the pixel positions to get their index in the 2D array. This gives the 1D to 2D
98 indices maps.
100 Notes
101 -----
102 This function is slow and runs only on cpu (scipy.cKDTree). It should be used only once to
103 get the mappings for a given camera geometry.
105 Parameters
106 ----------
107 pixel_positions : np.array
108 2D array (shape: (2, N_pixels)) containing the pixels x and y coordinates. This should
109 be the actual pixels positions in the camera. It is important that the referential in
110 which the position are measured has equal metrics on both axis (1m in the real world
111 is 1 on x _and_ y axis)
113 Returns
114 -------
115 Tuple[np.ndarray, np.ndarray]
116 Row and column indices (shape (N_pixels,) each) to use as 1D to 2D indexing maps.
117 """
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)