better leakage parameter calculation (#783)

* Implement get_border_pixel_mask

* Add leakage image parameter

* Do not use lru_cache to fix pickling

* Fix size in leakage

ctapipe/image/__init__.py

@@ -6,3 +6,4 @@
 from .reducers import *
 from .muon import *
 from .geometry_converter import *
+from .leakage import *

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,
+    )

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

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:

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]

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'
+    )