Source code for pyhiperta.waveform_indexing
"""Implements the indexing/mapping of waveforms from 1D array (hexagonal structure) to 2D square array.
The 1D array representation is suited for vectorization of "per-element" operations such as calibration,
Hillas computation, etc. while the 2D square representation is faster for operations requiring neighbors
look-up, such as the tail-cut cleaning. It is faster to copy the data from one representation to the other
than to perform the cleaning with the 1D representation.
"""
from typing import Tuple
import numpy as np
from scipy.spatial import cKDTree
def _get_neighbors_indices(pixel_positions: np.ndarray) -> np.ndarray:
"""Get the indices of the pixel neighbors in the array, for each pixel.
Parameters
----------
pixel_positions : np.ndarray
Positions of the pixels in the camera. 1st dimension is the "x" axis, 2nd dimension is "y".
Shape: (2, number_of_pixels)
Notes
-----
The returned array has shape (Number_of_pixels, 6) even though some pixels have less than 6 neighbors. When that
is the case, the "extra-neighbors" have indices greater or equal than the number of pixels (invalid indices!)
Warnings
--------
This function is only intended to work for hexagonal camera pixel positions. It assumes a pixel has at most
6 neighbors and the neigbors are placed on a hexagonal grid.
Returns
-------
np.ndarray
Indices of the pixels neighbors, per pixel. Shape: (Number_of_pixels, 6)
"""
dist_kdtree = cKDTree(pixel_positions.transpose())
# Query the tree for the 1st neighbours of each pixel to get the distance between 1st neighbours
# We query the tree for neighbour 2, because we query the tree for points identical to the tree points,
# so the first neighbour will have distance 0 and be the pixels themselves.
pixel_distances, _ = dist_kdtree.query(dist_kdtree.data, k=list(range(2, 3)), p=2, workers=-1)
dist_1st_neighbour = pixel_distances.min()
# Now query the tree for all neighbours at most (1 + sqrt(3)) / 2 to get all 1st neighbours
# We query for 6 neighbours, but not all pixels have 6 neighbours. When that is the case the
# neighbours indices of the are set to a value > the number of pixels.
_, neighbors_indices = dist_kdtree.query(
pixel_positions.transpose(),
k=list(range(2, 8)), # get 2nd to 7th neighbours
distance_upper_bound=(1.0 + np.sqrt(3)) / 2.0 * dist_1st_neighbour,
p=2, # 2 norm
workers=-1, # use all cpus
)
return neighbors_indices
[docs]
def waveform_1D_to_2D_maps(pixel_positions: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
"""Create the 2D waveform array with "squared geometry" and 1D to 2D indices maps.
To perform fast "neighbour looking" operations such as convolution we need to store the waveform
as a 2D array. To fit the hexagonal geometry of the camera into a 2D square grid (the 2D array)
we perform a rotation (to align an axis with the 2D array axis) and then deform the other axis to
align the hexagonal neighbour on the square grid.
If we represent the pixels by their index in the 1D array:
0 1 1 4
2 3 4 -> 0 3 6
5 6 2 5
The algorithm is:
- find 3 neighbours forming a triangle in the hexagonal representation.
The pixels should be such that: the left-most pixel is in between the other
two pixels on the y coordinate, eg:
x
x
x
Such a combination is always possible as the angle in the triangles is 60° and
we search a space of 90°. The goal of this is to use the left-most pixel as our
2D grid origin, and the other pixels coordinates will give us the rotation and
deformation.
In the example, the chosen are pixels are 0, 1 and 3. (1's y is equal or higher than 0's y)
- Rotate all pixels coordinates such that the "down" neighbour is aligned with the 2D
grid rows
x x
x ->
x x x
- Normalize the distances between pixels, so that a pixel is at distance 1 from another
on row and column axis (makes next step easier).
- "Deform" along the row axis to align the pixels on the 2D columns. This is a linear
transformation of the x coordinate with respect to the y coordinate:
___x x ___ is the deformation to apply for a distance of +1 y
->
x x x x
- Bin the pixel positions to get their index in the 2D array. This gives the 1D to 2D
indices maps.
Notes
-----
This function is slow and runs only on cpu (scipy.cKDTree). It should be used only once to
get the mappings for a given camera geometry.
Parameters
----------
pixel_positions : np.array
2D array (shape: (2, N_pixels)) containing the pixels x and y coordinates. This should
be the actual pixels positions in the camera. It is important that the referential in
which the position are measured has equal metrics on both axis (1m in the real world
is 1 on x _and_ y axis)
Returns
-------
Tuple[np.ndarray, np.ndarray]
Row and column indices (shape (N_pixels,) each) to use as 1D to 2D indexing maps.
"""
neighbors_indices = _get_neighbors_indices(pixel_positions)
# find three points that are:
# - all neighbours of each other
# - the 1 point is the left most point (minimal x)
# - the 2nd point is the higher point (max y)
# - the 3rd point is the lower (min y)
def find_neighbours_trio():
# Try all pixel as the left-most pixel
for idx_left_most_p, left_most_p_neighbours in enumerate(neighbors_indices):
# for all neighbours to left-most pixel
for idx_up_neighbour in left_most_p_neighbours:
if (
idx_up_neighbour < pixel_positions.shape[1] # valid neighbour
and pixel_positions[0, idx_left_most_p] <= pixel_positions[0, idx_up_neighbour]
and pixel_positions[1, idx_left_most_p] <= pixel_positions[1, idx_up_neighbour]
): # We have a good higher point. Now search for a most down point.
for idx_down_neighbour in neighbors_indices[idx_up_neighbour]:
if (
idx_down_neighbour < pixel_positions.shape[1] # valid neighbour
and idx_down_neighbour in left_most_p_neighbours # also a neighbour of left-most
and pixel_positions[0, idx_left_most_p] <= pixel_positions[0, idx_down_neighbour]
and pixel_positions[1, idx_left_most_p] >= pixel_positions[1, idx_down_neighbour]
):
return (idx_left_most_p, idx_up_neighbour, idx_down_neighbour)
raise ValueError("Could not find trio of points matching conditions.")
idx_left_most, idx_up_neighbour, idx_down_neighbour = find_neighbours_trio()
# Get the angle required to align the left-most - down-most points line to the x axis
# That is the angle between the points line and the x-axis (not the reverse ! orientation is important)
#
# x left-most - x down-most
# .___________________________
# |
# | y left-most - y down-most
# |
# .
psi = np.arctan2(
pixel_positions[1, idx_left_most] - pixel_positions[1, idx_down_neighbour],
pixel_positions[0, idx_down_neighbour] - pixel_positions[0, idx_left_most],
)
# Apply rotation
pixel_positions = np.array([[np.cos(psi), -np.sin(psi)], [np.sin(psi), np.cos(psi)]]) @ pixel_positions
# Normalize the neighbours distance for both axis (make neighbours distance 1 on x and y)
pixel_positions[0, :] *= 1.0 / (pixel_positions[0, idx_down_neighbour] - pixel_positions[0, idx_left_most])
pixel_positions[1, :] *= 1.0 / (pixel_positions[1, idx_up_neighbour] - pixel_positions[1, idx_left_most])
# Apply linear deformation to the x coordinate based on y coordinate
pixel_positions[0, :] -= (pixel_positions[1, :] - pixel_positions[1, idx_left_most]) * (
pixel_positions[0, idx_up_neighbour] - pixel_positions[0, idx_left_most]
)
# Bin the pixel coordinates to find their index in 2D array
# After the normalization and deformation the pixels are placed at integer values on the grid
# To make sure the points are not near the bins edge, make the bins shifted by 0.5
# The values are then floored and converted to int by the astype
# Note: - what binning x would do: (x - (min(x) - shift)) / bin_size (with integer division)
# - bin size is 1.0, so no need to divide.
# - row indices are reversed, because the position in the matrix start in the top left, not
# bottom left as the y coordinates do.
row_indices = (pixel_positions[1] - (pixel_positions[1].min() - 0.5)).astype(np.int32)
row_indices = row_indices.max() - row_indices
col_indices = (pixel_positions[0] - (pixel_positions[0].min() - 0.5)).astype(np.int32)
return (row_indices, col_indices)
[docs]
def waveform1Dto2D(
maps1D_to_2D: Tuple[np.ndarray, np.ndarray], waveform: np.ndarray, waveform2D: np.ndarray = None
) -> np.ndarray:
"""Get the 2D array representation of waveform
If waveform2D is passed, the values of waveform2D will be updated with `waveform`'s.
Otherwise a new 2D array will be created.
Parameters
----------
maps1D_to_2D : tuple of np.ndarray
The tuple holding the indices to go from 1D to 2D representation
waveform : np.ndarray
Waveform values in 1D shape.
waveform2D : np.ndarray, optionnal
An already existing 2D representation of the waveform, to update with `waveform`'s values.
Returns
-------
np.ndarray
A 2D array representing the waveform. (The array is independent of waveform, not a view)
"""
# Initialize waveform2D it with 0s if needed.
if waveform2D is None:
waveform2D = np.zeros(
waveform.shape[:-1] + (maps1D_to_2D[0].max() + 1, maps1D_to_2D[1].max() + 1), dtype=waveform.dtype
)
waveform2D[..., *maps1D_to_2D] = waveform
return waveform2D
[docs]
def waveform2Dto1D(maps1D_to_2D: np.ndarray, waveform2D: np.ndarray) -> np.ndarray:
"""Get the 1D array representing the waveform from the 2D representation `waveform2D`.
Notes
-----
The result is a view in `waveform2D`! Modifying the view will modify waveform2D!
Parameters
----------
maps1D_to_2D : tuple of np.ndarray
The tuple holding the indices to go from 2D to 1D representation
waveform2D : np.ndarray
The 2D array holding the waveform values.
Returns
-------
np.ndarray
1D array holding the waveform values.
"""
return np.ascontiguousarray(waveform2D[..., *maps1D_to_2D])
[docs]
def neighbors_only_stencil() -> np.ndarray:
"""Return the stencil selecting __only the neighbors__ of a pixel in 2D array.
The pixel itself is not included: center element of the stencil is 0.
Notes
-----
The stencil values depend on the way the 1D to 2D conversion is performed. Depending on how the axis are deformed,
it could be transposed or reversed. The current implementation does
0 1 1 4
2 3 4 -> 0 3 6
5 6 2 5
so the stencil is 1 where the neighbors are and 0 otherwise (and 0 at the center to exclude current pixel).
"""
return np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]], dtype=bool)
[docs]
def leakage_pixel_masks(pixel_positions: np.ndarray) -> np.ndarray:
"""Returns a mask (boolean array) that multiplied with waveforms (1D) puts all pixels except the outer rings to 0.
Parameters
----------
pixel_positions : np.ndarray
2D array (shape: (2, N_pixels)) containing the pixels x and y coordinates. This should
be the actual pixels positions in the camera. It is important that the referential in
which the position are measured has equal metrics on both axis (1m in the real world
is 1 on x _and_ y axis)
Returns
-------
leakage_mask: np.ndarray
Mask array with values True at the indices of pixels on the outermost rings, False otherwise.
The shape is (2, nb_pixels). leakage_mask[0, :] selects the pixels of the outermost ring (leakage 1)
while leakage[1:] selects the pixels of the 1st and 2nd outermost ring (leakage 2).
"""
neighbors_indices = _get_neighbors_indices(pixel_positions)
# leakage 1 pixels are the one of the borders, so they have less than 6 neighbors.
# we use the fact that the neighbors_indices array has invalid indices for non-existing neighbors.
leakage_1_mask = (neighbors_indices >= pixel_positions.shape[1]).any(axis=1)
# leakage 2 pixels are the ones that are neighbors of leakage 1 pixels, but are not in leakage 1 themselves.
# find them by getting the indices of leakage 1 pixels with nonzero, then use isin to see if there are neighbors
# finally exclude the leakage 1 pixels
leakage_2_mask = np.isin(neighbors_indices, np.nonzero(leakage_1_mask)[0]).any(axis=1) & np.bitwise_not(
leakage_1_mask
)
return np.stack([leakage_1_mask, leakage_1_mask | leakage_2_mask], axis=0)