"""Implements the "Hillas" (and some additional) parameters computation from integrated waveforms."""
from typing import Tuple
import numpy as np
from pyhiperta.utils.nanify import nanify
[docs]
def hillas(waveform: np.ndarray, position: np.ndarray, nan_threshold: float = 1.0e-4) -> Tuple[np.ndarray, np.ndarray]:
"""Compute hillas parameters (ellipse position, orientation, skewness, kurtosis) from shower images.
The Hillas parametrization consists in estimating the ellipse formed by a gamma-ray shower over the
background noise in the telescope integrated and cleaned images (R1 data level).
The ellipse parameters are found by performing a principal component analysis (PCA) of the images pixels
weighted by their charge. The ellipse center is defined as the image center of gravity ; the ellipse
axis and angle with respect to the x-axis are given by the PCA eigenvectors ; the ellipse length and
width are defined as twice the standard deviation along the ellipse axis.
In addition, the skewness and kurtosis of the image are computed along the major axis of the ellipse.
This method works with a single image, of with a batch of images. The implementation first computes
the 1st and 2nd order image moments to build the covariance matrix (PCA is computed via the
covariance matrix). The ellipse parameters are deduced from the covariance matrix and finally the
skewness and kurtosis are computed as 3rd and 4rd along the ellipse major axis.
This function also returns the longitudinal coordinate of the pixels along the ellipsis greater axis
since it is usefull for timing parameters computation.
Notes
-----
A gamma-ray shower is radially symmetric in the direction of arrival, therefore the distribution
along the small ellipse axis should be gaussian and the skewness and kurtosis are not computed
along this axis.
An ellipse with given orientation and skewness is equivalent to another ellipse with
orientation + pi and opposite skewness (psi + pi, -skewness). This method consistently returns
the orientation in the ]-pi/2, pi/2[ range and corresponding skewness.
The skewness and kurtosis are computed using the 3rd and 4th moments, which are biased estimators.
Hillas, A. (1985). Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei.
In 19th International Cosmic Ray Conference (ICRC19), Volume 3 (pp. 445).
Parameters
----------
waveform : np.ndarray
1D shower image(s). If a single image is provided the shape must be (N_pixels,). If a
batch of images is provided, the shape should be (N_batch, N_pixels).
position : np.ndarray
Image pixel position in (x, y) coordinates. The shape must be (2, N_pixels): position[0,:]
is the x coordinate, and position[1, :] the y coordinate, of the image pixels.
nan_threshold : float
Threshold to set some computed values to NaN instead of leaving them with problematic values.
Images with intensity smaller than the threshold will have all output set to NaN.
Images with lengths smaller than the threshold will have length, skewness and kurtosis set
to NaN.
Returns
-------
hillas_parameters : np.ndarray
Image(s) hillas parameters. This is a 1D array with shape (10,) if a single image was
provided, or a 2D array with shape (N_batch, 10) if a batch of images was provided. The
parameters are ordered like so:
hillas_parameters[..., 0]: intensity
hillas_parameters[..., 1]: center of gravity x coordinate
hillas_parameters[..., 2]: center of gravity y coordinate
hillas_parameters[..., 3]: center of gravity radial coordinate radius (r)
hillas_parameters[..., 4]: center of gravity radial coordinate angle (phi)
hillas_parameters[..., 5]: ellipse length
hillas_parameters[..., 6]: ellipse width
hillas_parameters[..., 7]: ellipse angle with respect to x-axis (psi) in ]-pi/2, pi/2[
hillas_parameters[..., 8]: ellipse skewness
hillas_parameters[..., 9]: ellipse kurtosis
longitudinal_coordinates : np.ndarray
Coordinates of the pixels along the long axis of the ellipsis.
"""
# some computation could have division by 0 or sqrt(negative_value)
# we don't want warnings and actually want the NaNs in the result, so
# simply disable these warnings for this block
with np.errstate(divide="ignore", invalid="ignore"):
# intensity is simply the sum of all pixels (for each image in batch)
intensity = waveform.sum(axis=-1)
# If intensity is too small, set it to NaN (all results will be meaningless)
intensity = nanify(intensity, intensity < nan_threshold)
# The PCA can be computed with the singular value decomposition (SVD) or covariance matrix. Here we
# use the covariance matrix.
# We want to compute the variance of x and y, as well as the covariance(x,y) along the distribution
# given by the image charge (pixel values). The distribution will be normalized by the total charge.
# To do so we will compute the 1st and 2nd order moments of x and y.
# (Notation: E[x] is expected value of x, w is pixel value)
# Compute all moments: E[w*x], E[w*y], E[w*x**2], E[w*y**2], E[w*x*y] and normalize
# (stack all positions powers: x, y, x**2, y**2, x*y) and broadcast multiply with waveform
moments = (
waveform[..., np.newaxis, :] # shape (N_batch, 1, N_pixels) to broadcast for all moments
* np.stack([position[0], position[1], position[0] ** 2, position[1] ** 2, position.prod(axis=0)])
).sum(axis=-1) / intensity[..., np.newaxis]
# E[x]**2 and E[y]**2 (will be re-used for center of gravity radial coordinates computation)
E_x_y_square = moments[..., 0:2] ** 2
# E[w*x**2] - E[w*x]**2 and E[w*y**2] - E[w*y]**2 (shape: (N_batch, 2))
variance = moments[..., 2:4] - E_x_y_square
# E[w*x*y] - E[w*x]*E[w*y] (shape: (N_batch,))
covariance = moments[..., 4] - moments[..., 0:2].prod(axis=-1)
# The covariance matrix would be: (Var(x) Cov(x,y))
# (Cov(x,y) Var(y) )
# Now we compute the eigenvalues lambda as the roots of the characteristic polynomial
# (lambda - Var(x)) * (lambda - Var(y)) - Cov(x,y)**2 = 0
varx_plus_vary = variance.sum(axis=-1)
char_pol_delta = np.sqrt(varx_plus_vary**2 + 4.0 * (covariance**2 - variance.prod(axis=-1)))
eigenValueHigh = (varx_plus_vary + char_pol_delta) / 2.0
eigenValueLow = (varx_plus_vary - char_pol_delta) / 2.0
# The length and width are twice the square roots of the eigenvalues (*2 is performed later)
# (<=> std along the major and minor axis)
# Note: an ellipse parametrized by A and B is defined by
# x^T (A B/2) x = 1
# (B/2 A )
# the semi-axis are given by the square root of the ellipse matrix eigenvalues.
length = np.sqrt(eigenValueHigh)
width = np.sqrt(eigenValueLow)
# If the length is too small, set it to NaN to avoid crazy values later on.
length = nanify(length, length < nan_threshold)
# The ellipse axis are along the eigenvectors of the covariance matrix, and the ellipse
# orientation angle (psi) is defined as the angle between the major axis and the x-axis.
# However, we can define the axis with 2 possible directions: for a given eigenvector e,
# -e is also an eigenvector with the same direction, but reversed sense.
# The convention here is to chose the eigenvector so that psi is in ]-pi/2, pi/2[
# Note: eigenvectors e are defined by lamba e = (Var(x) Cov(x, y)) e
# (Cov(x, y) Var(y) )
# To solve this system e[0] and e[1] must respect a fixed ratio, that is most simply
# expressed as e = (Cov(x, y), (lambda - Var(x))) or e = (lambda - Var(y), Cov(x,y))
# We are free to chose the direction of e, so we chose the one giving psi in ]-pi/2,pi/2[
psi = np.arctan2(covariance, eigenValueHigh - variance[..., 1])
# Finally, we compute the 3rd and 4th order moments along the ellipse major axis
# First we compute the pixels coordinates in the PCA coordinates: (cos(psi), sin(psi)) (p - cog)
# (-sin(psi), cos(psi))
# but we only care about the major axis, so only the 1st coordinate is computed
# longitudinal shape (N_batch, N_pixels) (the sum is along coordinate axis to do p * cos + p * sin)
longitudinal = (
(position - moments[..., 0:2, np.newaxis]) # shape (N_batch, 2, N_pixels)
* np.stack([np.cos(psi), np.sin(psi)], axis=-1)[..., np.newaxis] # shape (N_batch, 2, 1)
).sum(axis=-2)
# we want to compute longitudinal **3 and **4
# to do this with broadcasting, we build an array containing [3, 4] that broadcasts
# to (2, *longitudinal.shape). See np.power()
exponents = np.array([3, 4], dtype=np.float32)[:, *([np.newaxis] * longitudinal.ndim)]
# Compute 3rd and 4th moments E[w * l**3] and E[w * l**4], shape: (2, N_batch)
longitudinal_moments = ((longitudinal**exponents) * waveform).sum(axis=-1)
# Normalize by ellipse length
skewness, kurtosis = longitudinal_moments / ((length ** exponents.squeeze(-1)) * intensity)
# Radial coordinates of the center of gravity
r = np.sqrt(E_x_y_square.sum(axis=-1))
phi = np.arctan2(moments[..., 1], moments[..., 0])
# Stack everything in columns so it is ready for scikit-learn
return np.stack(
[
intensity,
moments[..., 0],
moments[..., 1],
r,
phi,
length * 2.0,
width * 2.0,
psi,
skewness,
kurtosis,
],
axis=-1,
), longitudinal