From 8428ee0b0149a5e81deb7a62e099fcd787133a18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 6 Sep 2017 17:53:57 +0200 Subject: [PATCH 1/4] Implement hillas features usen eigh --- ctapipe/image/hillas.py | 78 +++++++++++++++++++++++++++++- ctapipe/image/tests/test_hillas.py | 8 ++- 2 files changed, 83 insertions(+), 3 deletions(-) diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index c333b0ef57c..73dd7fc55d2 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -14,6 +14,7 @@ from ctapipe.instrument import CameraGeometry from ..io.containers import HillasParametersContainer + __all__ = [ 'MomentParameters', 'hillas_parameters', @@ -21,6 +22,7 @@ '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,79 @@ 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 + + cov = np.cov(delta_x, delta_y, aweights=image) + 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 diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index ced39260f85..e02143707ee 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -2,7 +2,7 @@ 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 @@ -79,6 +79,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: @@ -132,9 +133,12 @@ def test_hillas_failure(): with pytest.raises(HillasParameterizationError): hillas_parameters_4(geom, blank_image) + with pytest.raises(HillasParameterizationError): + hillas_parameters_5(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)) From 07bffe6c3392ac01afb0863f8a1696b5c983d4f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Mon, 18 Jun 2018 15:06:17 +0200 Subject: [PATCH 2/4] Use ddof=0 to make it comparable with the other methods --- ctapipe/image/hillas.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 73dd7fc55d2..10ac365bf1e 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -704,7 +704,7 @@ def hillas_parameters_5(geom: CameraGeometry, image): delta_x = pix_x - cog_x delta_y = pix_y - cog_y - cov = np.cov(delta_x, delta_y, aweights=image) + 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 From e289b11f74f3551bc712a5a99f580882abcf7da9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Mon, 18 Jun 2018 15:23:50 +0200 Subject: [PATCH 3/4] Add comment about ddof --- ctapipe/image/hillas.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 10ac365bf1e..ed2b20d4dda 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -704,6 +704,9 @@ def hillas_parameters_5(geom: CameraGeometry, image): 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) From 9c358459840ceb913144ff08d98c2b6d78266343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Mon, 18 Jun 2018 16:06:53 +0200 Subject: [PATCH 4/4] Fix tests --- ctapipe/image/tests/test_hillas.py | 31 +++++++++++++----------------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index e02143707ee..e8429ec3b23 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -10,6 +10,14 @@ 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'): @@ -90,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) @@ -112,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) @@ -121,21 +127,9 @@ 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) - - with pytest.raises(HillasParameterizationError): - hillas_parameters_5(geom, blank_image) - + for method in methods: + with pytest.raises(HillasParameterizationError): + method(geom, blank_image) def test_hillas_api_change(): @@ -145,5 +139,6 @@ def test_hillas_api_change(): 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)