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.

Functions:

Name	Description
leakage_pixel_masks	Returns a mask (boolean array) that multiplied with waveforms (1D) puts all pixels except the outer rings to 0.
neighbors_only_stencil	Return the stencil selecting only the neighbors of a pixel in 2D array.
waveform1Dto2D	Get the 2D array representation of waveform
waveform2Dto1D	Get the 1D array representing the waveform from the 2D representation waveform2D.
waveform_1D_to_2D_maps	Create the 2D waveform array with "squared geometry" and 1D to 2D indices maps.

leakage_pixel_masks

leakage_pixel_masks(pixel_positions)

Returns a mask (boolean array) that multiplied with waveforms (1D) puts all pixels except the outer rings to 0.

Parameters:

Name	Type	Description	Default
pixel_positions	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)	required

Returns:

Name	Type	Description
leakage_mask	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).

Source code in src/pyhiperta/waveform_indexing.py

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)

neighbors_only_stencil

neighbors_only_stencil()

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).

Source code in src/pyhiperta/waveform_indexing.py

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)

waveform1Dto2D

waveform1Dto2D(maps1D_to_2D, waveform, waveform2D=None)

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:

Name	Type	Description	Default
maps1D_to_2D	tuple of np.ndarray	The tuple holding the indices to go from 1D to 2D representation	required
waveform	ndarray	Waveform values in 1D shape.	required
waveform2D	(ndarray, optionnal)	An already existing 2D representation of the waveform, to update with waveform's values.	None

Returns:

Type	Description
ndarray	A 2D array representing the waveform. (The array is independent of waveform, not a view)

Source code in src/pyhiperta/waveform_indexing.py

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

waveform2Dto1D

waveform2Dto1D(maps1D_to_2D, waveform2D)

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:

Name	Type	Description	Default
maps1D_to_2D	tuple of np.ndarray	The tuple holding the indices to go from 2D to 1D representation	required
waveform2D	ndarray	The 2D array holding the waveform values.	required

Returns:

Type	Description
ndarray	1D array holding the waveform values.

Source code in src/pyhiperta/waveform_indexing.py

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])

waveform_1D_to_2D_maps

waveform_1D_to_2D_maps(pixel_positions)

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:

Name	Type	Description	Default
pixel_positions	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)	required

Returns:

Type	Description
Tuple[ndarray, ndarray]	Row and column indices (shape (N_pixels,) each) to use as 1D to 2D indexing maps.

Source code in src/pyhiperta/waveform_indexing.py

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)