Cleaning for biggest cluster

ctapipe/image/cleaning.py

@@ -218,7 +218,7 @@ def number_of_islands(geom, mask):
         Total number of clusters
     island_labels: ndarray
         Contains cluster membership of each pixel.
-        Dimesion equals input mask.
+        Dimension equals input geometry.
         Entries range from 0 (not in the pixel mask) to num_islands.
     """
     # compress sparse neighbor matrix
@@ -365,6 +365,30 @@ def fact_image_cleaning(
     return pixels_to_keep
 
 
+def largest_island(islands_labels):
+    """Find the biggest island and filter it from the image.
+
+    This function takes a list of islands in an image and isolates the largest one
+    for later parametrization.
+
+    Parameters
+    ----------
+    islands_labels : array
+        Flattened array containing a list of labelled islands from a cleaned image.
+        Second value returned by the function 'number_of_islands'.
+
+    Returns
+    -------
+    islands_labels : array
+        A boolean mask created from the input labels and filtered for the largest island.
+        If no islands survived the cleaning the array is all False.
+
+    """
+    if np.count_nonzero(islands_labels) == 0:
+        return np.zeros(islands_labels.shape, dtype="bool")
+    return islands_labels == np.argmax(np.bincount(islands_labels[islands_labels > 0]))
+
+
 class ImageCleaner(TelescopeComponent):
     """
     Abstract class for all configurable Image Cleaning algorithms.   Use
@@ -403,15 +427,15 @@ class TailcutsImageCleaner(ImageCleaner):
     """
 
     picture_threshold_pe = FloatTelescopeParameter(
-        help="top-level threshold in photoelectrons", default_value=10.0,
+        help="top-level threshold in photoelectrons", default_value=10.0
     ).tag(config=True)
 
     boundary_threshold_pe = FloatTelescopeParameter(
-        help="second-level threshold in photoelectrons", default_value=5.0,
+        help="second-level threshold in photoelectrons", default_value=5.0
     ).tag(config=True)
 
     min_picture_neighbors = IntTelescopeParameter(
-        help="Minimum number of neighbors above threshold to consider", default_value=2,
+        help="Minimum number of neighbors above threshold to consider", default_value=2
     ).tag(config=True)
 
     def __call__(

ctapipe/image/tests/test_cleaning.py

@@ -235,6 +235,54 @@ def test_number_of_islands():
     assert island_mask.dtype == np.int32
 
 
+def test_largest_island():
+    """Test selection of largest island in imagea with given cleaning masks."""
+
+    # Create a simple rectangular camera made of 17 pixels
+    camera = CameraGeometry.make_rectangular(17, 1)
+
+    # Assume a simple image (flattened array) made of 0, 1 or 2 photoelectrons
+    # [2, 1, 1, 1, 1, 2, 2, 1, 0, 2, 1, 1, 1, 0, 2, 2, 2]
+    # Assume a virtual tailcut cleaning that requires:
+    # - picture_threshold = 2 photoelectrons,
+    # - boundary_threshold = 1 photoelectron,
+    # - min_number_picture_neighbors = 0
+    # There will be 4 islands left after cleaning:
+    clean_mask = np.zeros(camera.n_pixels).astype("bool")  # initialization
+    clean_mask[[0, 1]] = 1
+    clean_mask[[4, 5, 6, 7]] = 2  # this is the biggest
+    clean_mask[[9, 10]] = 3
+    clean_mask[[14, 15, 16]] = 4
+    # Label islands (their number is not important here)
+    _, islands_labels = cleaning.number_of_islands(camera, clean_mask)
+    # Create the true mask which takes into account only the biggest island
+    # Pixels with no signal are labelled with a 0
+    true_mask_largest = np.zeros(camera.n_pixels).astype("bool")
+    true_mask_largest[[4, 5, 6, 7]] = 1
+    # Apply the function to test
+    mask_largest = cleaning.largest_island(islands_labels)
+
+    # Now the degenerate case of only one island surviving
+    # Same process as before
+    clean_mask_one = np.zeros(camera.n_pixels).astype("bool")
+    clean_mask_one[[0, 1]] = 1
+    _, islands_labels_one = cleaning.number_of_islands(camera, clean_mask_one)
+    true_mask_largest_one = np.zeros(camera.n_pixels).astype("bool")
+    true_mask_largest_one[[0, 1]] = 1
+    mask_largest_one = cleaning.largest_island(islands_labels_one)
+
+    # Last the case of no islands surviving
+    clean_mask_0 = np.zeros(camera.n_pixels).astype("bool")
+    _, islands_labels_0 = cleaning.number_of_islands(camera, clean_mask_0)
+    true_mask_largest_0 = np.zeros(camera.n_pixels).astype("bool")
+    mask_largest_0 = cleaning.largest_island(islands_labels_0)
+
+    # Check if the function recovers the ground truth in all of the three cases
+    assert (mask_largest_one == true_mask_largest_one).all()
+    assert (mask_largest_0 == true_mask_largest_0).all()
+    assert_allclose(mask_largest, true_mask_largest)
+
+
 def test_fact_image_cleaning():
     # use LST pixel geometry
     geom = CameraGeometry.from_name("LSTCam")
@@ -279,11 +327,7 @@ def test_apply_time_delta_cleaning():
 
     # Test unchanged
     td_mask = cleaning.apply_time_delta_cleaning(
-        geom,
-        mask,
-        pulse_time,
-        min_number_neighbors=1,
-        time_limit=5,
+        geom, mask, pulse_time, min_number_neighbors=1, time_limit=5
     )
     test_mask = mask.copy()
     assert (test_mask == td_mask).all()
@@ -292,23 +336,15 @@ def test_apply_time_delta_cleaning():
     noise_neighbour = neighbours[0]
     pulse_time[noise_neighbour] += 10
     td_mask = cleaning.apply_time_delta_cleaning(
-        geom,
-        mask,
-        pulse_time,
-        min_number_neighbors=1,
-        time_limit=5,
+        geom, mask, pulse_time, min_number_neighbors=1, time_limit=5
     )
     test_mask = mask.copy()
     test_mask[noise_neighbour] = 0
     assert (test_mask == td_mask).all()
 
     # Test min_number_neighbours
     td_mask = cleaning.apply_time_delta_cleaning(
-        geom,
-        mask,
-        pulse_time,
-        min_number_neighbors=4,
-        time_limit=5,
+        geom, mask, pulse_time, min_number_neighbors=4, time_limit=5
     )
     test_mask = mask.copy()
     test_mask[neighbours] = 0