Implement hillas features usen eigh

ctapipe/image/hillas.py

@@ -14,13 +14,15 @@
 from ctapipe.instrument import CameraGeometry
 from ..io.containers import HillasParametersContainer
 
+
 __all__ = [
     'MomentParameters',
     'hillas_parameters',
     'hillas_parameters_1',
     'hillas_parameters_2',
     'hillas_parameters_3',
     'hillas_parameters_4',
+    'hillas_parameters_5',
     'HillasParameterizationError',
 ]
 
@@ -518,7 +520,7 @@ def hillas_parameters_4(geom: CameraGeometry, image, container=False):
     unit = geom.pix_x.unit
 
     # MP: Actually, I don't know why we need to strip the units... shouldn'
-    # the calculations all work with them?'''
+    # the calculations all work with them?
 
     pix_x = Quantity(np.asanyarray(geom.pix_x, dtype=np.float64)).value
     pix_y = Quantity(np.asanyarray(geom.pix_y, dtype=np.float64)).value
@@ -658,5 +660,82 @@ def hillas_parameters_4(geom: CameraGeometry, image, container=False):
                                 skewness=skewness, kurtosis=kurtosis)
 
 
+def hillas_parameters_5(geom: CameraGeometry, image):
+    """
+    Compute Hillas parameters for a given shower image.
+
+    Implementation uses a PCA analogous to the implementation in
+    src/main/java/fact/features/HillasParameters.java
+    from
+    https://github.com/fact-project/fact-tools
+
+    Parameters
+    ----------
+    geom: ctapipe.instrument.CameraGeometry
+        Camera geometry
+    image : array_like
+        Pixel values
+
+    Returns
+    -------
+    HillasParametersContainer:
+        container of hillas parametesr
+    """
+    unit = geom.pix_x.unit
+    pix_x = Quantity(np.asanyarray(geom.pix_x, dtype=np.float64)).value
+    pix_y = Quantity(np.asanyarray(geom.pix_y, dtype=np.float64)).value
+    image = np.asanyarray(image, dtype=np.float64)
+    assert pix_x.shape == pix_y.shape == image.shape, 'Image and pixel shape do not match'
+
+    size = np.sum(image)
+
+    if size == 0.0:
+        raise HillasParameterizationError('size=0, cannot calculate HillasParameters')
+
+    # calculate the cog as the mean of the coordinates weighted with the image
+    cog_x = np.average(pix_x, weights=image)
+    cog_y = np.average(pix_y, weights=image)
+
+    # polar coordinates of the cog
+    cog_r = np.linalg.norm([cog_x, cog_y])
+    cog_phi = np.arctan2(cog_y, cog_x)
+
+    # do the PCA for the hillas parameters
+    delta_x = pix_x - cog_x
+    delta_y = pix_y - cog_y
+
+    # The ddof=0 makes this comparable to the other methods,
+    # but ddof=1 should be more correct, mostly affects small showers
+    # on a percent level
+    cov = np.cov(delta_x, delta_y, aweights=image, ddof=0)
+    eig_vals, eig_vecs = np.linalg.eigh(cov)
+
+    # width and length are eigen values of the PCA
+    width, length = np.sqrt(eig_vals)
+
+    # psi is the angle of the eigenvector to length to the x-axis
+    psi = np.arctan2(eig_vecs[1, 1], eig_vecs[0, 1])
+
+    # calculate higher order moments along shower axes
+    longitudinal = delta_x * np.cos(psi) + delta_y * np.sin(psi)
+
+    m3_long = np.average(longitudinal**3, weights=image)
+    skewness_long = m3_long / length**3
+
+    m4_long = np.average(longitudinal**4, weights=image)
+    kurtosis_long = m4_long / length**4
+
+    return HillasParametersContainer(
+        x=cog_x * unit, y=cog_y * unit,
+        r=cog_r * unit, phi=Angle(cog_phi * u.rad),
+        intensity=size,
+        length=length * unit,
+        width=width * unit,
+        psi=Angle(psi * u.rad),
+        skewness=skewness_long,
+        kurtosis=kurtosis_long,
+    )
+
+
 # use the 4 version by default.
 hillas_parameters = hillas_parameters_4

ctapipe/image/tests/test_hillas.py

@@ -2,14 +2,22 @@
 from ctapipe.image import tailcuts_clean, toymodel
 from ctapipe.image.hillas import (hillas_parameters_1, hillas_parameters_2,
                                   hillas_parameters_3, hillas_parameters_4,
-                                  HillasParameterizationError)
+                                  hillas_parameters_5, HillasParameterizationError)
 from ctapipe.io.containers import HillasParametersContainer
 from astropy import units as u
 from numpy import isclose, zeros_like, arange
 from numpy.random import seed
 from numpy.ma import masked_array
 import pytest
 
+methods = (
+    hillas_parameters_1,
+    hillas_parameters_2,
+    hillas_parameters_3,
+    hillas_parameters_4,
+    hillas_parameters_5
+)
+
 
 def create_sample_image(psi='-30d'):
 
@@ -79,6 +87,7 @@ def test_hillas():
         results['v2'] = hillas_parameters_2(geom, image)
         results['v3'] = hillas_parameters_3(geom, image)
         results['v4'] = hillas_parameters_4(geom, image)
+        results['v5'] = hillas_parameters_5(geom, image)
         # compare each method's output
         for aa in results:
             for bb in results:
@@ -89,7 +98,6 @@ def test_hillas():
                     compare_result(results[aa].r, results[bb].r)
                     compare_result(results[aa].phi.deg, results[bb].phi.deg)
                     compare_result(results[aa].psi.deg, results[bb].psi.deg)
-                    compare_result(results[aa].miss, results[bb].miss)
                     compare_result(results[aa].skewness, results[bb].skewness)
                     # compare_result(results[aa].kurtosis, results[bb].kurtosis)
 
@@ -111,7 +119,6 @@ def test_hillas_masked():
     compare_result(results.r, results_ma.r)
     compare_result(results.phi.deg, results_ma.phi.deg)
     compare_result(results.psi.deg, results_ma.psi.deg)
-    compare_result(results.miss, results_ma.miss)
     compare_result(results.skewness, results_ma.skewness)
     # compare_result(results.kurtosis, results_ma.kurtosis)
 
@@ -120,26 +127,18 @@ def test_hillas_failure():
     geom, image = create_sample_image_zeros(psi='0d')
     blank_image = zeros_like(image)
 
-    with pytest.raises(HillasParameterizationError):
-        hillas_parameters_1(geom, blank_image)
-
-    with pytest.raises(HillasParameterizationError):
-        hillas_parameters_2(geom, blank_image)
-
-    with pytest.raises(HillasParameterizationError):
-        hillas_parameters_3(geom, blank_image)
-
-    with pytest.raises(HillasParameterizationError):
-        hillas_parameters_4(geom, blank_image)
+    for method in methods:
+        with pytest.raises(HillasParameterizationError):
+            method(geom, blank_image)
 
 
 def test_hillas_api_change():
-    import numpy as np
     with pytest.raises(ValueError):
         hillas_parameters_4(arange(10), arange(10), arange(10))
 
 
 def test_hillas_container():
     geom, image = create_sample_image_zeros(psi='0d')
+
     params = hillas_parameters_4(geom, image, container=True)
-    assert type(params) is HillasParametersContainer
+    assert isinstance(params, HillasParametersContainer)