diff --git a/ctapipe/image/__init__.py b/ctapipe/image/__init__.py index 50c29065b71..a1eceb9554c 100644 --- a/ctapipe/image/__init__.py +++ b/ctapipe/image/__init__.py @@ -6,3 +6,4 @@ from .reducers import * from .muon import * from .geometry_converter import * +from .leakage import * diff --git a/ctapipe/image/leakage.py b/ctapipe/image/leakage.py new file mode 100644 index 00000000000..36ed895d9f6 --- /dev/null +++ b/ctapipe/image/leakage.py @@ -0,0 +1,50 @@ +""" +Leakage calculation +""" + +import numpy as np +from ..io.containers import LeakageContainer + + +__all__ = ['leakage'] + + +def leakage(geom, image, cleaning_mask): + ''' + Calculating the leakage-values for a given image. + Image must be cleaned for example with tailcuts_clean. + Leakage describes how strong a shower is on the edge of a telescope. + + Parameters + ---------- + geom: ctapipe.instrument.CameraGeometry + Camera geometry information + image: array + pixel values + cleaning_mask: array, dtype=bool + The pixel that survived cleaning, e.g. tailcuts_clean + + Returns + ------- + LeakageContainer + ''' + border1 = geom.get_border_pixel_mask(1) + border2 = geom.get_border_pixel_mask(2) + + mask1 = border1 & cleaning_mask + mask2 = border2 & cleaning_mask + + leakage_pixel1 = np.sum(mask1) + leakage_pixel2 = np.sum(mask2) + + leakage_intensity1 = np.sum(image[mask1]) + leakage_intensity2 = np.sum(image[mask2]) + + size = np.sum(image[cleaning_mask]) + + return LeakageContainer( + leakage1_pixel=leakage_pixel1 / geom.n_pixels, + leakage2_pixel=leakage_pixel2 / geom.n_pixels, + leakage1_intensity=leakage_intensity1 / size, + leakage2_intensity=leakage_intensity2 / size, + ) diff --git a/ctapipe/image/tests/test_leakage.py b/ctapipe/image/tests/test_leakage.py new file mode 100644 index 00000000000..59b3a3cb793 --- /dev/null +++ b/ctapipe/image/tests/test_leakage.py @@ -0,0 +1,21 @@ +from ctapipe.instrument.camera import CameraGeometry +import numpy as np + + +def test_leakage(): + from ctapipe.image.leakage import leakage + + geom = CameraGeometry.from_name('LSTCam') + + img = np.ones(geom.n_pixels) + mask = np.ones(len(geom), dtype=bool) + + l = leakage(geom, img, mask) + + ratio1 = np.sum(geom.get_border_pixel_mask(1)) / geom.n_pixels + ratio2 = np.sum(geom.get_border_pixel_mask(2)) / geom.n_pixels + + assert l.leakage1_intensity == ratio1 + assert l.leakage2_intensity == ratio2 + assert l.leakage1_pixel == ratio1 + assert l.leakage2_pixel == ratio2 diff --git a/ctapipe/instrument/camera.py b/ctapipe/instrument/camera.py index 53e700d42c3..9360147834d 100644 --- a/ctapipe/instrument/camera.py +++ b/ctapipe/instrument/camera.py @@ -15,6 +15,7 @@ from ctapipe.utils import get_table_dataset, find_all_matching_datasets from ctapipe.utils.linalg import rotation_matrix_2d + __all__ = ['CameraGeometry'] logger = logging.getLogger(__name__) @@ -90,6 +91,8 @@ def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type, pix_rotation="0d", cam_rotation="0d", neighbors=None, apply_derotation=True): + assert len(pix_x) == len(pix_y), 'pix_x and pix_y must have same length' + self.n_pixels = len(pix_x) self.cam_id = cam_id self.pix_id = pix_id self.pix_x = pix_x @@ -110,6 +113,9 @@ def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type, if len(pix_x.shape) == 1: self.rotate(cam_rotation) + # cache border pixel mask per instance + self.border_cache = {} + def __eq__(self, other): return ((self.cam_id == other.cam_id) and (self.pix_x == other.pix_x).all() @@ -119,6 +125,9 @@ def __eq__(self, other): and (self.pix_type == other.pix_type) ) + def __len__(self): + return self.n_pixels + def __getitem__(self, slice_): return CameraGeometry(cam_id=" ".join([self.cam_id, " sliced"]), pix_id=self.pix_id[slice_], @@ -450,6 +459,34 @@ def make_rectangular(cls, npix_x=40, npix_y=40, range_x=(-0.5, 0.5), neighbors=None, pix_type='rectangular') + def get_border_pixel_mask(self, width=1): + ''' + Get a mask for pixels at the border of the camera of arbitrary width + + Parameters + ---------- + width: int + The width of the border in pixels + + Returns + ------- + mask: array + A boolean mask, True if pixel is in the border of the specified width + ''' + if width in self.border_cache: + return self.border_cache[width] + + if width == 1: + n_neighbors = self.neighbor_matrix_sparse.sum(axis=1).A1 + max_neighbors = n_neighbors.max() + mask = n_neighbors < max_neighbors + else: + n = self.neighbor_matrix + mask = (n & self.get_border_pixel_mask(width - 1)).any(axis=1) + + self.border_cache[width] = mask + return mask + # ====================================================================== # utility functions: diff --git a/ctapipe/instrument/tests/test_camera.py b/ctapipe/instrument/tests/test_camera.py index e98f4c8c264..2bffe7ba551 100644 --- a/ctapipe/instrument/tests/test_camera.py +++ b/ctapipe/instrument/tests/test_camera.py @@ -150,3 +150,20 @@ def test_slicing(): assert sliced2.pix_id[0] == 5 assert sliced2.pix_id[1] == 7 assert len(sliced2.pix_x) == 5 + + +def test_border_pixels(): + from ctapipe.instrument.camera import CameraGeometry + + cam = CameraGeometry.from_name("LSTCam") + + assert np.sum(cam.get_border_pixel_mask(1)) == 168 + assert np.sum(cam.get_border_pixel_mask(2)) == 330 + + cam = CameraGeometry.from_name("ASTRICam") + assert np.sum(cam.get_border_pixel_mask(1)) == 212 + assert np.sum(cam.get_border_pixel_mask(2)) == 408 + + assert cam.get_border_pixel_mask(1)[0] + assert cam.get_border_pixel_mask(1)[2351] + assert not cam.get_border_pixel_mask(1)[521] diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py index 71361348844..a677bf4033b 100644 --- a/ctapipe/io/containers.py +++ b/ctapipe/io/containers.py @@ -39,7 +39,8 @@ 'DataContainer', 'TargetIODataContainer', 'SST1MDataContainer', - 'HillasParametersContainer'] + 'HillasParametersContainer', + 'LeakageContainer'] class SST1MCameraContainer(Container): @@ -669,3 +670,29 @@ class HillasParametersContainer(Container): skewness = Field(nan, 'measure of the asymmetry') kurtosis = Field(nan, 'measure of the tailedness') + + +class LeakageContainer(Container): + """ + Leakage + """ + leakage1_pixel = Field( + nan, + 'Percentage of pixels after cleaning' + ' that are in camera border of width=1' + ) + leakage2_pixel = Field( + nan, + 'Percentage of pixels after cleaning' + ' that are in camera border of width=2' + ) + leakage1_intensity = Field( + nan, + 'Percentage of photo-electrons after cleaning' + ' that are in the camera border of width=1' + ) + leakage2_intensity = Field( + nan, + 'Percentage of photo-electrons after cleaning' + ' that are in the camera border of width=2' + )