Coverage for src/pyhiperta/hillas.py: 100%
25 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"""Hillas (and some additional) parameters computation from integrated waveforms."""
7import numpy as np
9from pyhiperta.utils.nanify import nanify
12def hillas(
13 waveform: np.ndarray, position: np.ndarray, nan_threshold: float = 1.0e-4
14) -> tuple[np.ndarray, np.ndarray]:
15 """Compute hillas parameters (ellipse position, orientation, skewness, kurtosis) from shower images.
17 The Hillas parametrization consists in estimating the ellipse formed by a gamma-ray shower over the
18 background noise in the telescope integrated and cleaned images (R1 data level).
20 The ellipse parameters are found by performing a principal component analysis (PCA) of the images pixels
21 weighted by their charge. The ellipse center is defined as the image center of gravity ; the ellipse
22 axis and angle with respect to the x-axis are given by the PCA eigenvectors ; the ellipse length and
23 width are defined as twice the standard deviation along the ellipse axis.
24 In addition, the skewness and kurtosis of the image are computed along the major axis of the ellipse.
26 This method works with a single image, of with a batch of images. The implementation first computes
27 the 1st and 2nd order image moments to build the covariance matrix (PCA is computed via the
28 covariance matrix). The ellipse parameters are deduced from the covariance matrix and finally the
29 skewness and kurtosis are computed as 3rd and 4rd along the ellipse major axis.
31 This function also returns the longitudinal coordinate of the pixels along the ellipsis greater axis
32 since it is usefull for timing parameters computation.
34 Notes
35 -----
36 A gamma-ray shower is radially symmetric in the direction of arrival, therefore the distribution
37 along the small ellipse axis should be gaussian and the skewness and kurtosis are not computed
38 along this axis.
40 An ellipse with given orientation and skewness is equivalent to another ellipse with
41 orientation + pi and opposite skewness (psi + pi, -skewness). This method consistently returns
42 the orientation in the ]-pi/2, pi/2[ range and corresponding skewness.
44 The skewness and kurtosis are computed using the 3rd and 4th moments, which are biased estimators.
46 Hillas, A. (1985). Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei.
47 In 19th International Cosmic Ray Conference (ICRC19), Volume 3 (pp. 445).
49 Parameters
50 ----------
51 waveform : np.ndarray
52 1D shower image(s). If a single image is provided the shape must be (N_pixels,). If a
53 batch of images is provided, the shape should be (N_batch, N_pixels).
54 position : np.ndarray
55 Image pixel position in (x, y) coordinates. The shape must be (2, N_pixels): position[0,:]
56 is the x coordinate, and position[1, :] the y coordinate, of the image pixels.
57 nan_threshold : float
58 Threshold to set some computed values to NaN instead of leaving them with problematic values.
59 Images with intensity smaller than the threshold will have all output set to NaN.
60 Images with lengths smaller than the threshold will have length, skewness and kurtosis set
61 to NaN.
63 Returns
64 -------
65 hillas_parameters : np.ndarray
66 Image(s) hillas parameters. This is a 1D array with shape (10,) if a single image was
67 provided, or a 2D array with shape (N_batch, 10) if a batch of images was provided. The
68 parameters are ordered like so:
69 hillas_parameters[..., 0]: intensity
70 hillas_parameters[..., 1]: center of gravity x coordinate
71 hillas_parameters[..., 2]: center of gravity y coordinate
72 hillas_parameters[..., 3]: center of gravity radial coordinate radius (r)
73 hillas_parameters[..., 4]: center of gravity radial coordinate angle (phi)
74 hillas_parameters[..., 5]: ellipse length
75 hillas_parameters[..., 6]: ellipse width
76 hillas_parameters[..., 7]: ellipse angle with respect to x-axis (psi) in ]-pi/2, pi/2[
77 hillas_parameters[..., 8]: ellipse skewness
78 hillas_parameters[..., 9]: ellipse kurtosis
79 longitudinal_coordinates : np.ndarray
80 Coordinates of the pixels along the long axis of the ellipsis.
81 """
82 # some computation could have division by 0 or sqrt(negative_value)
83 # we don't want warnings and actually want the NaNs in the result, so
84 # simply disable these warnings for this block
85 with np.errstate(divide="ignore", invalid="ignore"):
86 # intensity is simply the sum of all pixels (for each image in batch)
87 intensity = waveform.sum(axis=-1)
88 # If intensity is too small, set it to NaN (all results will be meaningless)
89 intensity = nanify(intensity, intensity < nan_threshold)
91 # The PCA can be computed with the singular value decomposition (SVD) or covariance matrix. Here we
92 # use the covariance matrix.
93 # We want to compute the variance of x and y, as well as the covariance(x,y) along the distribution
94 # given by the image charge (pixel values). The distribution will be normalized by the total charge.
95 # To do so we will compute the 1st and 2nd order moments of x and y.
96 # (Notation: E[x] is expected value of x, w is pixel value)
98 # Compute all moments: E[w*x], E[w*y], E[w*x**2], E[w*y**2], E[w*x*y] and normalize
99 # (stack all positions powers: x, y, x**2, y**2, x*y) and broadcast multiply with waveform
100 moments = (
101 waveform[..., np.newaxis, :] # shape (N_batch, 1, N_pixels) to broadcast for all moments
102 * np.stack([position[0], position[1], position[0] ** 2, position[1] ** 2, position.prod(axis=0)])
103 ).sum(axis=-1) / intensity[..., np.newaxis]
105 # E[x]**2 and E[y]**2 (will be re-used for center of gravity radial coordinates computation)
106 E_x_y_square = moments[..., 0:2] ** 2
107 # E[w*x**2] - E[w*x]**2 and E[w*y**2] - E[w*y]**2 (shape: (N_batch, 2))
108 variance = moments[..., 2:4] - E_x_y_square
109 # E[w*x*y] - E[w*x]*E[w*y] (shape: (N_batch,))
110 covariance = moments[..., 4] - moments[..., 0:2].prod(axis=-1)
112 # The covariance matrix would be: (Var(x) Cov(x,y))
113 # (Cov(x,y) Var(y) )
114 # Now we compute the eigenvalues lambda as the roots of the characteristic polynomial
115 # (lambda - Var(x)) * (lambda - Var(y)) - Cov(x,y)**2 = 0
116 varx_plus_vary = variance.sum(axis=-1)
117 char_pol_delta = np.sqrt(varx_plus_vary**2 + 4.0 * (covariance**2 - variance.prod(axis=-1)))
118 eigenValueHigh = (varx_plus_vary + char_pol_delta) / 2.0
119 eigenValueLow = (varx_plus_vary - char_pol_delta) / 2.0
121 # The length and width are twice the square roots of the eigenvalues (*2 is performed later)
122 # (<=> std along the major and minor axis)
123 # Note: an ellipse parametrized by A and B is defined by
124 # x^T (A B/2) x = 1
125 # (B/2 A )
126 # the semi-axis are given by the square root of the ellipse matrix eigenvalues.
127 length = np.sqrt(eigenValueHigh)
128 width = np.sqrt(eigenValueLow)
129 # If the length is too small, set it to NaN to avoid crazy values later on.
130 length = nanify(length, length < nan_threshold)
132 # The ellipse axis are along the eigenvectors of the covariance matrix, and the ellipse
133 # orientation angle (psi) is defined as the angle between the major axis and the x-axis.
134 # However, we can define the axis with 2 possible directions: for a given eigenvector e,
135 # -e is also an eigenvector with the same direction, but reversed sense.
136 # The convention here is to chose the eigenvector so that psi is in ]-pi/2, pi/2[
137 # Note: eigenvectors e are defined by lamba e = (Var(x) Cov(x, y)) e
138 # (Cov(x, y) Var(y) )
139 # To solve this system e[0] and e[1] must respect a fixed ratio, that is most simply
140 # expressed as e = (Cov(x, y), (lambda - Var(x))) or e = (lambda - Var(y), Cov(x,y))
141 # We are free to chose the direction of e, so we chose the one giving psi in ]-pi/2,pi/2[
142 psi = np.arctan2(covariance, eigenValueHigh - variance[..., 1])
144 # Finally, we compute the 3rd and 4th order moments along the ellipse major axis
145 # First we compute the pixels coordinates in the PCA coordinates: (cos(psi), sin(psi)) (p - cog)
146 # (-sin(psi), cos(psi))
147 # but we only care about the major axis, so only the 1st coordinate is computed
148 # longitudinal shape (N_batch, N_pixels) (the sum is along coordinate axis to do p * cos + p * sin)
149 longitudinal = (
150 (position - moments[..., 0:2, np.newaxis]) # shape (N_batch, 2, N_pixels)
151 * np.stack([np.cos(psi), np.sin(psi)], axis=-1)[..., np.newaxis] # shape (N_batch, 2, 1)
152 ).sum(axis=-2)
154 # we want to compute longitudinal **3 and **4
155 # to do this with broadcasting, we build an array containing [3, 4] that broadcasts
156 # to (2, *longitudinal.shape). See np.power()
157 exponents = np.array([3, 4], dtype=np.float32)[:, *([np.newaxis] * longitudinal.ndim)]
158 # Compute 3rd and 4th moments E[w * l**3] and E[w * l**4], shape: (2, N_batch)
159 longitudinal_moments = ((longitudinal**exponents) * waveform).sum(axis=-1)
160 # Normalize by ellipse length
161 skewness, kurtosis = longitudinal_moments / ((length ** exponents.squeeze(-1)) * intensity)
163 # Radial coordinates of the center of gravity
164 r = np.sqrt(E_x_y_square.sum(axis=-1))
165 phi = np.arctan2(moments[..., 1], moments[..., 0])
167 # Stack everything in columns so it is ready for scikit-learn
168 return np.stack(
169 [
170 intensity,
171 moments[..., 0],
172 moments[..., 1],
173 r,
174 phi,
175 length * 2.0,
176 width * 2.0,
177 psi,
178 skewness,
179 kurtosis,
180 ],
181 axis=-1,
182 ), longitudinal