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