From c2667882d829db3ea68300362204c9365c21e0f9 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 4 Apr 2024 12:26:36 +0200 Subject: [PATCH 01/21] Initial commit --- src/ctapipe/image/cleaning.py | 227 +++++++++++++++++++++++++++++++++- 1 file changed, 226 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 6dc8ec54c62..95aa1998354 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -36,7 +36,7 @@ FloatTelescopeParameter, IntTelescopeParameter, ) -from .morphology import brightest_island, number_of_islands +from .morphology import brightest_island, largest_island, number_of_islands def tailcuts_clean( @@ -112,6 +112,41 @@ def tailcuts_clean( ) +def bright_cleaning(image, threshold=267, fraction=0.03): + """ + Clean an image by removing pixels below + fraction * (mean charge in the 3 brightest pixels). Select all pixels + instead if the mean charge of the brightest pixels are below a + certain threshold. + + Parameters + ---------- + image: array + pixel values + threshold: `float` + Minimum average charge in the 3 brightest pixels to apply + cleaning + fraction: `float` + Pixels below fraction * (average charge in the 3 brightest pixels) + will be removed from the cleaned image + + Returns + ------- + A boolean mask of *clean* pixels. + + """ + max_3_value_index = np.argsort(image)[-3:] + mean_3_max_signal = np.mean(image[max_3_value_index]) + + if mean_3_max_signal < threshold: + return np.ones(image.shape, dtype=bool) + + threshold_brightest = fraction * mean_3_max_signal + mask = image >= threshold_brightest + + return mask + + def mars_cleaning_1st_pass( geom, image, @@ -459,6 +494,122 @@ def time_constrained_clean( return mask_core | mask_boundary +def lst_image_cleaning( + geom, + image, + arrival_times, + picture_thresh=8, + boundary_thresh=4, + min_number_picture_neighbors=2, + keep_isolated_pixels=False, + time_limit=2, + apply_bright_cleaning=True, + fraction_bright_cleaning=0.03, + threshold_bright_cleaning=267, + largest_island_only=False, + pedestal_factor=2.5, + pedestal_std=None, +): + """ + Clean an image in 5 Steps: + + 1) Get picture threshold for `tailcuts_clean` in step 2) from interleaved + pedestal events if `pedestal_factor` and `pedestal_std` is not None. + 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. + 3) Apply time_delta_cleaning algorithm - + `ctapipe.image.cleaning.apply_time_delta_cleaning` if `time_limit` is not None. + 4) Apply bright_cleaning - `ctapipe.image.cleaning.bright_cleaning` if + `apply_bright_cleaning` is set to true. + 5) Get only largest island - `ctapipe.image.morphology.largest_island` if + `largest_island_only` is set to true. + + Parameters + ---------- + geom: `ctapipe.instrument.CameraGeometry` + Camera geometry information + image: `np.ndarray` + pixel charges + arrival_times: `np.ndarray` + Pixel timing information + picture_thresh: `float` or `np.ndarray` + threshold above which all pixels are retained. Used for `tailcuts_clean`. + boundary_thresh: `float` or `np.ndarray` + threshold above which pixels are retained if they have a neighbor + already above the picture_thresh. Used for `tailcuts_clean`. + min_number_picture_neighbors: `int` + A picture pixel survives cleaning only if it has at least this number + of picture neighbors. This has no effect in case keep_isolated_pixels is True. + Used for `tailcuts_clean`. + keep_isolated_pixels: `bool` + If True, pixels above the picture threshold will be included always, + if not they are only included if a neighbor is in the picture or + boundary. Used for `tailcuts_clean`. + time_limit: `float` + Time limit for the `time_delta_cleaning`. Set to None if no + `time_delta_cleaning` should be applied. + apply_bright_cleaning: `bool` + Set to true if `bright_cleaning` should be applied. + fraction_bright_cleaning: `float` + `fraction` parameter for `bright_cleaning`. Pixels below + fraction * (average charge in the 3 brightest pixels) will be removed from + the cleaned image. + threshold_bright_cleaning: `float` + `threshold` parameter for `bright_cleaning`. Minimum average charge + in the 3 brightest pixels to apply the cleaning. + largest_island_only: `bool` + Set to true to get only largest island. + pedestal_factor: `float` + Factor for interleaved pedestal cleaning. It is multiplied by the + pedestal standard deviation for each pixel to calculate pixelwise picture + threshold parameters for `tailcuts_clean` considering the current background. + Set to None if no pedestal cleaning should be applied. + pedestal_std: `np.ndarray` + Pedestal standard deviation for each pixel. See + `ctapipe.containers.PedestalContainer` + + Returns + ------- + A boolean mask of *clean* pixels. + + """ + # Step 1 + if pedestal_factor is not None and pedestal_std is not None: + pedestal_threshold = pedestal_std * pedestal_factor + picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) + # Step 2 + mask = tailcuts_clean( + geom, + image, + picture_thresh=picture_thresh, + boundary_thresh=boundary_thresh, + min_number_picture_neighbors=min_number_picture_neighbors, + keep_isolated_pixels=keep_isolated_pixels, + ) + # Check that at least one pixel survives tailcuts_clean + if ~mask.any(): + return mask + # Step 3 + if time_limit is not None: + mask = apply_time_delta_cleaning( + geom, + mask, + arrival_times, + min_number_neighbors=1, + time_limit=time_limit, + ) + # Step 4 + if apply_bright_cleaning: + mask = bright_cleaning( + image, mask, threshold_bright_cleaning, fraction_bright_cleaning + ) + # Step 5 + if largest_island_only: + _, island_labels = number_of_islands(geom, mask) + mask = largest_island(island_labels) + + return mask + + class ImageCleaner(TelescopeComponent): """ Abstract class for all configurable Image Cleaning algorithms. Use @@ -544,6 +695,80 @@ def __call__( ) +class LSTImageCleaner(TailcutsImageCleaner): + """ + Clean images using lstchains image cleaning technique. See + `ctapipe.image.lst_image_cleaning` + """ + + time_limit = FloatTelescopeParameter( + default_value=2, + help="Time limit for the `time_delta_cleaning`. Set to None if no" + "`time_delta_cleaning` should be applied", + ).tag(config=True) + + apply_bright_cleaning = BoolTelescopeParameter( + default_value=True, help="Set to true if `bright_cleaning` should be applied" + ).tag(config=True) + + fraction_bright_cleaning = FloatTelescopeParameter( + default_value=0.03, + help="`fraction` parameter for `bright_cleaning`. Pixels below " + "fraction * (average charge in the 3 brightest pixels) will be removed from " + "the cleaned image", + ).tag(config=True) + + threshold_bright_cleaning = FloatTelescopeParameter( + default_value=267, + help="`threshold` parameter for `bright_cleaning`. Minimum average charge " + "in the 3 brightest pixels to apply the cleaning", + ).tag(config=True) + + largest_island_only = BoolTelescopeParameter( + default_value=False, help="Set to true to get only largest island" + ).tag(config=True) + + pedestal_factor = FloatTelescopeParameter( + default_value=2.5, + help="Factor for interleaved pedestal cleaning. It is multiplied by the " + "pedestal standard deviation for each pixel to calculate pixelwise picture " + "threshold parameters for `tailcuts_clean` considering the current background. " + "Set to None if no pedestal cleaning should be applied.", + ).tag(config=True) + + def __call__( + self, + tel_id: int, + image: np.ndarray, + arrival_times: np.ndarray = None, + *, + monitoring: MonitoringCameraContainer = None, + ) -> np.ndarray: + """ + Apply LST image cleaning. See `ImageCleaner.__call__()` + """ + pedestal_std = None + if monitoring is not None: + pedestal_std = monitoring.tel[tel_id].pedestal.charge_std + + return lst_image_cleaning( + self.subarray.tel[tel_id].camera.geometry, + image, + arrival_times, + picture_thresh=self.picture_threshold_pe.tel[tel_id], + boundary_thresh=self.boundary_threshold_pe.tel[tel_id], + min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id], + keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], + time_limit=self.time_limit.tel[tel_id], + apply_bright_cleaning=self.apply_bright_cleaning.tel[tel_id], + fraction_bright_cleaning=self.fraction_bright_cleaning.tel[tel_id], + threshold_bright_cleaning=self.threshold_bright_cleaning.tel[tel_id], + largest_island_only=self.largest_island_only.tel[tel_id], + pedestal_factor=self.pedestal_factor.tel[tel_id], + pedestal_std=pedestal_std, + ) + + class MARSImageCleaner(TailcutsImageCleaner): """ 1st-pass MARS-like Image cleaner (See `ctapipe.image.mars_cleaning_1st_pass`) From 2db9ab7458ff3bcc9c1009cc153b7c6e8797ae34 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 12 Apr 2024 11:59:39 +0200 Subject: [PATCH 02/21] Fix bright cleaning --- src/ctapipe/image/cleaning.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 95aa1998354..7590134a272 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -115,7 +115,7 @@ def tailcuts_clean( def bright_cleaning(image, threshold=267, fraction=0.03): """ Clean an image by removing pixels below - fraction * (mean charge in the 3 brightest pixels). Select all pixels + fraction * (mean charge in the 3 brightest pixels). Select no pixels instead if the mean charge of the brightest pixels are below a certain threshold. @@ -139,7 +139,7 @@ def bright_cleaning(image, threshold=267, fraction=0.03): mean_3_max_signal = np.mean(image[max_3_value_index]) if mean_3_max_signal < threshold: - return np.ones(image.shape, dtype=bool) + return np.zeros(image.shape, dtype=bool) threshold_brightest = fraction * mean_3_max_signal mask = image >= threshold_brightest @@ -576,6 +576,7 @@ def lst_image_cleaning( if pedestal_factor is not None and pedestal_std is not None: pedestal_threshold = pedestal_std * pedestal_factor picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) + # Step 2 mask = tailcuts_clean( geom, @@ -586,8 +587,9 @@ def lst_image_cleaning( keep_isolated_pixels=keep_isolated_pixels, ) # Check that at least one pixel survives tailcuts_clean - if ~mask.any(): + if np.count_nonzero(mask) == 0: return mask + # Step 3 if time_limit is not None: mask = apply_time_delta_cleaning( @@ -597,11 +599,13 @@ def lst_image_cleaning( min_number_neighbors=1, time_limit=time_limit, ) + # Step 4 if apply_bright_cleaning: - mask = bright_cleaning( - image, mask, threshold_bright_cleaning, fraction_bright_cleaning + mask &= bright_cleaning( + image, threshold_bright_cleaning, fraction_bright_cleaning ) + # Step 5 if largest_island_only: _, island_labels = number_of_islands(geom, mask) From 4ff15bf9f38d063a75cdd1cb95b4bd2c4bf72294 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 12 Apr 2024 16:01:10 +0200 Subject: [PATCH 03/21] Adding unit tests --- src/ctapipe/image/cleaning.py | 71 +++++++----- src/ctapipe/image/tests/test_cleaning.py | 102 ++++++++++++++++++ .../tests/test_image_cleaner_component.py | 7 ++ 3 files changed, 151 insertions(+), 29 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 7590134a272..6ceeb4eb795 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -258,8 +258,8 @@ def apply_time_delta_cleaning( arrival_times: array pixel timing information min_number_neighbors: int - Threshold to determine if a pixel survives cleaning steps. - These steps include checks of neighbor arrival time and value + a selected pixel needs at least this number of (already selected) neighbors + that arrived within a given time_limit to itself to survive the cleaning. time_limit: int or float arrival time limit for neighboring pixels @@ -503,9 +503,10 @@ def lst_image_cleaning( min_number_picture_neighbors=2, keep_isolated_pixels=False, time_limit=2, + time_num_neighbors=1, apply_bright_cleaning=True, - fraction_bright_cleaning=0.03, - threshold_bright_cleaning=267, + bright_cleaning_fraction=0.03, + bright_cleaning_threshold=267, largest_island_only=False, pedestal_factor=2.5, pedestal_std=None, @@ -514,10 +515,10 @@ def lst_image_cleaning( Clean an image in 5 Steps: 1) Get picture threshold for `tailcuts_clean` in step 2) from interleaved - pedestal events if `pedestal_factor` and `pedestal_std` is not None. + pedestal events if `pedestal_factor` and `pedestal_std` is not 0. 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. 3) Apply time_delta_cleaning algorithm - - `ctapipe.image.cleaning.apply_time_delta_cleaning` if `time_limit` is not None. + `ctapipe.image.cleaning.apply_time_delta_cleaning` if `time_limit` is not 0. 4) Apply bright_cleaning - `ctapipe.image.cleaning.bright_cleaning` if `apply_bright_cleaning` is set to true. 5) Get only largest island - `ctapipe.image.morphology.largest_island` if @@ -545,15 +546,19 @@ def lst_image_cleaning( if not they are only included if a neighbor is in the picture or boundary. Used for `tailcuts_clean`. time_limit: `float` - Time limit for the `time_delta_cleaning`. Set to None if no + Time limit for the `time_delta_cleaning`. Set to 0 if no `time_delta_cleaning` should be applied. + time_num_neighbors: int + Used for `time_delta_cleaning`. + A selected pixel needs at least this number of (already selected) neighbors + that arrived within a given time_limit to itself to survive this cleaning. apply_bright_cleaning: `bool` Set to true if `bright_cleaning` should be applied. - fraction_bright_cleaning: `float` + bright_cleaning_fraction: `float` `fraction` parameter for `bright_cleaning`. Pixels below fraction * (average charge in the 3 brightest pixels) will be removed from the cleaned image. - threshold_bright_cleaning: `float` + bright_cleaning_threshold: `float` `threshold` parameter for `bright_cleaning`. Minimum average charge in the 3 brightest pixels to apply the cleaning. largest_island_only: `bool` @@ -562,7 +567,7 @@ def lst_image_cleaning( Factor for interleaved pedestal cleaning. It is multiplied by the pedestal standard deviation for each pixel to calculate pixelwise picture threshold parameters for `tailcuts_clean` considering the current background. - Set to None if no pedestal cleaning should be applied. + Set to 0 if no pedestal cleaning should be applied. pedestal_std: `np.ndarray` Pedestal standard deviation for each pixel. See `ctapipe.containers.PedestalContainer` @@ -573,7 +578,7 @@ def lst_image_cleaning( """ # Step 1 - if pedestal_factor is not None and pedestal_std is not None: + if pedestal_factor != 0 and pedestal_std is not None: pedestal_threshold = pedestal_std * pedestal_factor picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) @@ -591,19 +596,19 @@ def lst_image_cleaning( return mask # Step 3 - if time_limit is not None: + if time_limit != 0: mask = apply_time_delta_cleaning( geom, mask, arrival_times, - min_number_neighbors=1, + min_number_neighbors=time_num_neighbors, time_limit=time_limit, ) # Step 4 if apply_bright_cleaning: mask &= bright_cleaning( - image, threshold_bright_cleaning, fraction_bright_cleaning + image, bright_cleaning_threshold, bright_cleaning_fraction ) # Step 5 @@ -707,25 +712,32 @@ class LSTImageCleaner(TailcutsImageCleaner): time_limit = FloatTelescopeParameter( default_value=2, - help="Time limit for the `time_delta_cleaning`. Set to None if no" - "`time_delta_cleaning` should be applied", + help="Time limit for the `time_delta_cleaning`. Set to 0 if no" + " `time_delta_cleaning` should be applied", + ).tag(config=True) + + time_num_neighbors = IntTelescopeParameter( + default_value=1, + help="Used for `time_delta_cleaning`." + " A selected pixel needs at least this number of (already selected) neighbors" + " that arrived within a given time_limit to itself to survive this cleaning.", ).tag(config=True) apply_bright_cleaning = BoolTelescopeParameter( default_value=True, help="Set to true if `bright_cleaning` should be applied" ).tag(config=True) - fraction_bright_cleaning = FloatTelescopeParameter( + bright_cleaning_fraction = FloatTelescopeParameter( default_value=0.03, - help="`fraction` parameter for `bright_cleaning`. Pixels below " - "fraction * (average charge in the 3 brightest pixels) will be removed from " - "the cleaned image", + help="`fraction` parameter for `bright_cleaning`. Pixels below" + " fraction * (average charge in the 3 brightest pixels) will be removed from" + " the cleaned image", ).tag(config=True) - threshold_bright_cleaning = FloatTelescopeParameter( + bright_cleaning_threshold = FloatTelescopeParameter( default_value=267, - help="`threshold` parameter for `bright_cleaning`. Minimum average charge " - "in the 3 brightest pixels to apply the cleaning", + help="`threshold` parameter for `bright_cleaning`. Minimum average charge" + " in the 3 brightest pixels to apply the cleaning", ).tag(config=True) largest_island_only = BoolTelescopeParameter( @@ -734,10 +746,10 @@ class LSTImageCleaner(TailcutsImageCleaner): pedestal_factor = FloatTelescopeParameter( default_value=2.5, - help="Factor for interleaved pedestal cleaning. It is multiplied by the " - "pedestal standard deviation for each pixel to calculate pixelwise picture " - "threshold parameters for `tailcuts_clean` considering the current background. " - "Set to None if no pedestal cleaning should be applied.", + help="Factor for interleaved pedestal cleaning. It is multiplied by the" + " pedestal standard deviation for each pixel to calculate pixelwise picture" + " threshold parameters for `tailcuts_clean` considering the current background." + " Set to 0 if no pedestal cleaning should be applied.", ).tag(config=True) def __call__( @@ -764,9 +776,10 @@ def __call__( min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id], keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], time_limit=self.time_limit.tel[tel_id], + time_num_neighbors=self.time_num_neighbors.tel[tel_id], apply_bright_cleaning=self.apply_bright_cleaning.tel[tel_id], - fraction_bright_cleaning=self.fraction_bright_cleaning.tel[tel_id], - threshold_bright_cleaning=self.threshold_bright_cleaning.tel[tel_id], + bright_cleaning_fraction=self.bright_cleaning_fraction.tel[tel_id], + bright_cleaning_threshold=self.bright_cleaning_threshold.tel[tel_id], largest_island_only=self.largest_island_only.tel[tel_id], pedestal_factor=self.pedestal_factor.tel[tel_id], pedestal_std=pedestal_std, diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 1e799a4d73e..b7c158ebc97 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -457,3 +457,105 @@ def test_time_constrained_clean(prod5_lst): ) test_mask[noise_boundary] = 0 assert (test_mask == mask_reco).all() + + +def test_bright_cleaning(): + n_pixels = 1855 + fraction = 0.5 + image = np.zeros(n_pixels) + # set 3 bright pixels to 100 so mean of them is 100 as well + image[:3] = 100 + # 10 pixels above fraction * mean + image[10:20] = 60 + # 15 pixels below fraction * mean + image[50:65] = 30 + + threshold = 90 + mask = cleaning.bright_cleaning(image, threshold, fraction) + assert np.count_nonzero(mask) == 3 + 10 + # test that it doesn't select any pixels if mean of the 3 brightest pixels + # is below threshold + threshold = 110 + mask = cleaning.bright_cleaning(image, threshold, fraction) + assert np.count_nonzero(mask) == 0 + + +def test_lst_image_cleaning(prod5_lst): + geom = prod5_lst.camera.geometry + charge = np.zeros(geom.n_pixels, dtype=np.float32) + peak_time = np.zeros(geom.n_pixels, dtype=np.float32) + + core_pixel_1 = 100 + neighbors_1 = geom.neighbors[core_pixel_1] + core_pixel_2 = 1100 + neighbors_2 = geom.neighbors[core_pixel_2] + + # Set core pixel to 50 and first row of neighbors to 29. + charge[neighbors_1] = 29 + charge[core_pixel_1] = 50 + charge[neighbors_2] = 29 + charge[core_pixel_2] = 50 + + args = { + "picture_thresh": 45, + "boundary_thresh": 20, + "keep_isolated_pixels": True, + "time_limit": 0, + "apply_bright_cleaning": False, + "bright_cleaning_fraction": 0.5, + "bright_cleaning_threshold": 40, + "largest_island_only": False, + "pedestal_factor": 0, + "pedestal_std": None, + } + + # return normal tailcuts cleaning result if every other step is None/False: + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + # 2 * (1 core and 6 boundary_pixels) should be selected + assert np.count_nonzero(mask) == 2 * (1 + 6) + + # Test that tailcuts core threshold is adapted correctly with the pedestal + # charge std. + pedestal_std = np.ones(geom.n_pixels) + pedestal_std[core_pixel_1] = 30 + args["pedestal_factor"] = 2 + args["pedestal_std"] = pedestal_std + + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + # only core_pixel_2 with its boundaries should be selected since + # pedestal_std[core_pixel_1] * pedestal_factor > charge[core_pixel_1] + assert np.count_nonzero(mask) == 1 + 6 + + # if no pixel survives tailcuts cleaning it should not select any pixel + pedestal_std[core_pixel_2] = 30 + args["pedestal_std"] = pedestal_std + + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + assert np.count_nonzero(mask) == 0 + + # deselect core_pixel_1 with time_delta_cleaning by setting all its neighbors + # peak_time to 10 + peak_time[neighbors_1] = 10 + args["pedestal_factor"] = 0 + args["time_limit"] = 3 + + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + assert np.count_nonzero(mask) == 1 + 6 + 6 + + # 3 brightest pixels are [50, 50, 29], so mean is 43. With a fraction of 0.9 + # all boundaries should be deselected + args["time_limit"] = 0 + args["apply_bright_cleaning"] = True + args["bright_cleaning_fraction"] = 0.9 + + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + assert np.count_nonzero(mask) == 1 + 1 + + # Set neighbors of core_pixel_2 to 0 so largest_island_only should select only + # core_pixel_1 with its neighbors + charge[neighbors_2] = 0 + args["apply_bright_cleaning"] = False + args["largest_island_only"] = True + + mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + assert np.count_nonzero(mask) == 1 + 6 diff --git a/src/ctapipe/image/tests/test_image_cleaner_component.py b/src/ctapipe/image/tests/test_image_cleaner_component.py index 1b8ef34274b..e927fe93d51 100644 --- a/src/ctapipe/image/tests/test_image_cleaner_component.py +++ b/src/ctapipe/image/tests/test_image_cleaner_component.py @@ -16,6 +16,13 @@ def test_image_cleaner(method, prod5_mst_nectarcam, reference_location): "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, }, + "LSTImageCleaner": { + "boundary_threshold_pe": 5.0, + "picture_threshold_pe": 10.0, + "time_limit": 0, + "apply_bright_cleaning": False, + "pedestal_factor": 0, + }, "MARSImageCleaner": { "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, From 47dda603e0df934cad3a54b41f97076ce457a33e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 12 Apr 2024 17:19:47 +0200 Subject: [PATCH 04/21] Add apply_... traits for pedestal cleaning and time_delta_cleaning --- src/ctapipe/image/cleaning.py | 37 +++++++++++++------ src/ctapipe/image/tests/test_cleaning.py | 17 +++++---- .../tests/test_image_cleaner_component.py | 4 +- 3 files changed, 37 insertions(+), 21 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 6ceeb4eb795..f8ad36602a2 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -502,12 +502,14 @@ def lst_image_cleaning( boundary_thresh=4, min_number_picture_neighbors=2, keep_isolated_pixels=False, + apply_time_cleaning=True, time_limit=2, time_num_neighbors=1, apply_bright_cleaning=True, bright_cleaning_fraction=0.03, bright_cleaning_threshold=267, largest_island_only=False, + apply_pedestal_cleaning=True, pedestal_factor=2.5, pedestal_std=None, ): @@ -515,10 +517,12 @@ def lst_image_cleaning( Clean an image in 5 Steps: 1) Get picture threshold for `tailcuts_clean` in step 2) from interleaved - pedestal events if `pedestal_factor` and `pedestal_std` is not 0. + pedestal events if `apply_pedestal_cleaning` is set to true + and `pedestal_std` is not None. 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. 3) Apply time_delta_cleaning algorithm - - `ctapipe.image.cleaning.apply_time_delta_cleaning` if `time_limit` is not 0. + `ctapipe.image.cleaning.apply_time_delta_cleaning` if `apply_time_cleaning` is + set to true. 4) Apply bright_cleaning - `ctapipe.image.cleaning.bright_cleaning` if `apply_bright_cleaning` is set to true. 5) Get only largest island - `ctapipe.image.morphology.largest_island` if @@ -545,9 +549,10 @@ def lst_image_cleaning( If True, pixels above the picture threshold will be included always, if not they are only included if a neighbor is in the picture or boundary. Used for `tailcuts_clean`. + apply_time_cleaning: `bool` + Set to true if `time_delta_cleaning` should be applied. time_limit: `float` - Time limit for the `time_delta_cleaning`. Set to 0 if no - `time_delta_cleaning` should be applied. + Time limit for the `time_delta_cleaning`. time_num_neighbors: int Used for `time_delta_cleaning`. A selected pixel needs at least this number of (already selected) neighbors @@ -563,11 +568,12 @@ def lst_image_cleaning( in the 3 brightest pixels to apply the cleaning. largest_island_only: `bool` Set to true to get only largest island. + apply_pedestal_cleaning: `bool` + Set to true if pedestal cleaning should be applied. pedestal_factor: `float` Factor for interleaved pedestal cleaning. It is multiplied by the pedestal standard deviation for each pixel to calculate pixelwise picture threshold parameters for `tailcuts_clean` considering the current background. - Set to 0 if no pedestal cleaning should be applied. pedestal_std: `np.ndarray` Pedestal standard deviation for each pixel. See `ctapipe.containers.PedestalContainer` @@ -578,7 +584,7 @@ def lst_image_cleaning( """ # Step 1 - if pedestal_factor != 0 and pedestal_std is not None: + if apply_pedestal_cleaning and pedestal_std is not None: pedestal_threshold = pedestal_std * pedestal_factor picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) @@ -596,7 +602,7 @@ def lst_image_cleaning( return mask # Step 3 - if time_limit != 0: + if apply_time_cleaning: mask = apply_time_delta_cleaning( geom, mask, @@ -710,10 +716,14 @@ class LSTImageCleaner(TailcutsImageCleaner): `ctapipe.image.lst_image_cleaning` """ + apply_time_cleaning = BoolTelescopeParameter( + default_value=True, + help="Set to true if `time_delta_cleaning` should be applied", + ).tag(config=True) + time_limit = FloatTelescopeParameter( default_value=2, - help="Time limit for the `time_delta_cleaning`. Set to 0 if no" - " `time_delta_cleaning` should be applied", + help="Time limit for the `time_delta_cleaning`.", ).tag(config=True) time_num_neighbors = IntTelescopeParameter( @@ -744,12 +754,15 @@ class LSTImageCleaner(TailcutsImageCleaner): default_value=False, help="Set to true to get only largest island" ).tag(config=True) + apply_pedestal_cleaning = BoolTelescopeParameter( + default_value=True, help="Set to true if pedestal cleaning should be applied" + ).tag(config=True) + pedestal_factor = FloatTelescopeParameter( default_value=2.5, help="Factor for interleaved pedestal cleaning. It is multiplied by the" " pedestal standard deviation for each pixel to calculate pixelwise picture" - " threshold parameters for `tailcuts_clean` considering the current background." - " Set to 0 if no pedestal cleaning should be applied.", + " threshold parameters for `tailcuts_clean` considering the current background.", ).tag(config=True) def __call__( @@ -775,12 +788,14 @@ def __call__( boundary_thresh=self.boundary_threshold_pe.tel[tel_id], min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id], keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], + apply_time_cleaning=self.apply_time_cleaning.tel[tel_id], time_limit=self.time_limit.tel[tel_id], time_num_neighbors=self.time_num_neighbors.tel[tel_id], apply_bright_cleaning=self.apply_bright_cleaning.tel[tel_id], bright_cleaning_fraction=self.bright_cleaning_fraction.tel[tel_id], bright_cleaning_threshold=self.bright_cleaning_threshold.tel[tel_id], largest_island_only=self.largest_island_only.tel[tel_id], + apply_pedestal_cleaning=self.apply_pedestal_cleaning.tel[tel_id], pedestal_factor=self.pedestal_factor.tel[tel_id], pedestal_std=pedestal_std, ) diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index b7c158ebc97..06bed795951 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -500,12 +500,14 @@ def test_lst_image_cleaning(prod5_lst): "picture_thresh": 45, "boundary_thresh": 20, "keep_isolated_pixels": True, - "time_limit": 0, + "apply_time_cleaning": False, + "time_limit": 3, "apply_bright_cleaning": False, - "bright_cleaning_fraction": 0.5, + "bright_cleaning_fraction": 0.9, "bright_cleaning_threshold": 40, "largest_island_only": False, - "pedestal_factor": 0, + "apply_pedestal_cleaning": False, + "pedestal_factor": 2, "pedestal_std": None, } @@ -518,7 +520,7 @@ def test_lst_image_cleaning(prod5_lst): # charge std. pedestal_std = np.ones(geom.n_pixels) pedestal_std[core_pixel_1] = 30 - args["pedestal_factor"] = 2 + args["apply_pedestal_cleaning"] = True args["pedestal_std"] = pedestal_std mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) @@ -536,17 +538,16 @@ def test_lst_image_cleaning(prod5_lst): # deselect core_pixel_1 with time_delta_cleaning by setting all its neighbors # peak_time to 10 peak_time[neighbors_1] = 10 - args["pedestal_factor"] = 0 - args["time_limit"] = 3 + args["apply_pedestal_cleaning"] = False + args["apply_time_cleaning"] = True mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 6 + 6 # 3 brightest pixels are [50, 50, 29], so mean is 43. With a fraction of 0.9 # all boundaries should be deselected - args["time_limit"] = 0 + args["apply_time_cleaning"] = False args["apply_bright_cleaning"] = True - args["bright_cleaning_fraction"] = 0.9 mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 1 diff --git a/src/ctapipe/image/tests/test_image_cleaner_component.py b/src/ctapipe/image/tests/test_image_cleaner_component.py index e927fe93d51..cb3bb532c87 100644 --- a/src/ctapipe/image/tests/test_image_cleaner_component.py +++ b/src/ctapipe/image/tests/test_image_cleaner_component.py @@ -19,9 +19,9 @@ def test_image_cleaner(method, prod5_mst_nectarcam, reference_location): "LSTImageCleaner": { "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, - "time_limit": 0, + "apply_time_cleaning": False, "apply_bright_cleaning": False, - "pedestal_factor": 0, + "apply_pedestal_cleaning": False, }, "MARSImageCleaner": { "boundary_threshold_pe": 5.0, From 81d873f54d6703010915c3724111729164bcf62b Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 09:51:44 +0200 Subject: [PATCH 05/21] Add some comment --- src/ctapipe/image/tests/test_cleaning.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 06bed795951..55c6cb0a336 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -491,6 +491,7 @@ def test_lst_image_cleaning(prod5_lst): neighbors_2 = geom.neighbors[core_pixel_2] # Set core pixel to 50 and first row of neighbors to 29. + # These two islands does not overlap. charge[neighbors_1] = 29 charge[core_pixel_1] = 50 charge[neighbors_2] = 29 From 4e57616dcbc249a1f38f0916012a4b067f380b0b Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 11:00:03 +0200 Subject: [PATCH 06/21] Add changelog --- docs/changes/2541.feature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2541.feature.rst diff --git a/docs/changes/2541.feature.rst b/docs/changes/2541.feature.rst new file mode 100644 index 00000000000..30a27dbea21 --- /dev/null +++ b/docs/changes/2541.feature.rst @@ -0,0 +1 @@ +Add lstchains image cleaning procedure including its pedestal cleaning method. From 146c353b1e676e5d70f75bd6c8a02e880d3e34d1 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 11:08:56 +0200 Subject: [PATCH 07/21] Delete default values from bright cleaning --- src/ctapipe/image/cleaning.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index f8ad36602a2..f2eb6bfc12e 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -112,7 +112,7 @@ def tailcuts_clean( ) -def bright_cleaning(image, threshold=267, fraction=0.03): +def bright_cleaning(image, threshold, fraction): """ Clean an image by removing pixels below fraction * (mean charge in the 3 brightest pixels). Select no pixels From 4df3d94ab5ad5008dfa27dadf3e760071db0016c Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 12:13:03 +0200 Subject: [PATCH 08/21] Adapt bright cleaning docstring --- src/ctapipe/image/cleaning.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index f2eb6bfc12e..d2dbd38dc71 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -114,10 +114,11 @@ def tailcuts_clean( def bright_cleaning(image, threshold, fraction): """ - Clean an image by removing pixels below - fraction * (mean charge in the 3 brightest pixels). Select no pixels - instead if the mean charge of the brightest pixels are below a - certain threshold. + Clean an image by removing pixels below a fraction of the mean charge + in the 3 brightest pixels. + + Select no pixels instead if the mean charge of the brightest pixels + are below a certain threshold. Parameters ---------- @@ -156,7 +157,7 @@ def mars_cleaning_1st_pass( min_number_picture_neighbors=0, ): """ - Clean an image by selection pixels that pass a three-threshold tail-cuts + Clean an image by select pixels that pass a three-threshold tail-cuts procedure. All thresholds are defined with respect to the pedestal dispersion. All pixels that have a signal higher than the core (picture) From 9dca99b7115b6c5c16ee7d963de0f231e97fb4bf Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 14:06:49 +0200 Subject: [PATCH 09/21] Delete apply-traits - instead allow time_limit and bright_cleaning_threshold to be None --- src/ctapipe/image/cleaning.py | 61 ++++++------------- src/ctapipe/image/tests/test_cleaning.py | 18 +++--- .../tests/test_image_cleaner_component.py | 5 +- 3 files changed, 29 insertions(+), 55 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index d2dbd38dc71..78afed826c4 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -28,8 +28,7 @@ import numpy as np -from ctapipe.containers import MonitoringCameraContainer - +from ..containers import MonitoringCameraContainer from ..core import TelescopeComponent from ..core.traits import ( BoolTelescopeParameter, @@ -503,29 +502,24 @@ def lst_image_cleaning( boundary_thresh=4, min_number_picture_neighbors=2, keep_isolated_pixels=False, - apply_time_cleaning=True, time_limit=2, time_num_neighbors=1, - apply_bright_cleaning=True, bright_cleaning_fraction=0.03, bright_cleaning_threshold=267, largest_island_only=False, - apply_pedestal_cleaning=True, pedestal_factor=2.5, pedestal_std=None, ): """ Clean an image in 5 Steps: - 1) Get picture threshold for `tailcuts_clean` in step 2) from interleaved - pedestal events if `apply_pedestal_cleaning` is set to true - and `pedestal_std` is not None. + 1) Get picture thresholds for `tailcuts_clean` in step 2) from interleaved + pedestal events `pedestal_std` is not None. 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. 3) Apply time_delta_cleaning algorithm - - `ctapipe.image.cleaning.apply_time_delta_cleaning` if `apply_time_cleaning` is - set to true. + `ctapipe.image.cleaning.apply_time_delta_cleaning` if time_limit is not None. 4) Apply bright_cleaning - `ctapipe.image.cleaning.bright_cleaning` if - `apply_bright_cleaning` is set to true. + bright_cleaning_threshold is not None. 5) Get only largest island - `ctapipe.image.morphology.largest_island` if `largest_island_only` is set to true. @@ -550,27 +544,24 @@ def lst_image_cleaning( If True, pixels above the picture threshold will be included always, if not they are only included if a neighbor is in the picture or boundary. Used for `tailcuts_clean`. - apply_time_cleaning: `bool` - Set to true if `time_delta_cleaning` should be applied. time_limit: `float` - Time limit for the `time_delta_cleaning`. + Time limit for the `time_delta_cleaning`. Set to None if no `time_delta_cleaning` + should be applied. time_num_neighbors: int Used for `time_delta_cleaning`. A selected pixel needs at least this number of (already selected) neighbors that arrived within a given time_limit to itself to survive this cleaning. - apply_bright_cleaning: `bool` - Set to true if `bright_cleaning` should be applied. bright_cleaning_fraction: `float` `fraction` parameter for `bright_cleaning`. Pixels below fraction * (average charge in the 3 brightest pixels) will be removed from - the cleaned image. + the cleaned image. Set `bright_cleaning_threshold` to None if no + `bright_cleaning` should be applied. bright_cleaning_threshold: `float` `threshold` parameter for `bright_cleaning`. Minimum average charge - in the 3 brightest pixels to apply the cleaning. + in the 3 brightest pixels to apply the cleaning. Set to None if no + `bright_cleaning` should be applied. largest_island_only: `bool` Set to true to get only largest island. - apply_pedestal_cleaning: `bool` - Set to true if pedestal cleaning should be applied. pedestal_factor: `float` Factor for interleaved pedestal cleaning. It is multiplied by the pedestal standard deviation for each pixel to calculate pixelwise picture @@ -585,7 +576,7 @@ def lst_image_cleaning( """ # Step 1 - if apply_pedestal_cleaning and pedestal_std is not None: + if pedestal_std is not None: pedestal_threshold = pedestal_std * pedestal_factor picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) @@ -603,7 +594,7 @@ def lst_image_cleaning( return mask # Step 3 - if apply_time_cleaning: + if time_limit is not None: mask = apply_time_delta_cleaning( geom, mask, @@ -613,7 +604,7 @@ def lst_image_cleaning( ) # Step 4 - if apply_bright_cleaning: + if bright_cleaning_threshold is not None: mask &= bright_cleaning( image, bright_cleaning_threshold, bright_cleaning_fraction ) @@ -717,14 +708,11 @@ class LSTImageCleaner(TailcutsImageCleaner): `ctapipe.image.lst_image_cleaning` """ - apply_time_cleaning = BoolTelescopeParameter( - default_value=True, - help="Set to true if `time_delta_cleaning` should be applied", - ).tag(config=True) - time_limit = FloatTelescopeParameter( default_value=2, - help="Time limit for the `time_delta_cleaning`.", + help="Time limit for the `time_delta_cleaning`. Set to None if no" + " `time_delta_cleaning` should be applied.", + allow_none=True, ).tag(config=True) time_num_neighbors = IntTelescopeParameter( @@ -734,10 +722,6 @@ class LSTImageCleaner(TailcutsImageCleaner): " that arrived within a given time_limit to itself to survive this cleaning.", ).tag(config=True) - apply_bright_cleaning = BoolTelescopeParameter( - default_value=True, help="Set to true if `bright_cleaning` should be applied" - ).tag(config=True) - bright_cleaning_fraction = FloatTelescopeParameter( default_value=0.03, help="`fraction` parameter for `bright_cleaning`. Pixels below" @@ -748,17 +732,15 @@ class LSTImageCleaner(TailcutsImageCleaner): bright_cleaning_threshold = FloatTelescopeParameter( default_value=267, help="`threshold` parameter for `bright_cleaning`. Minimum average charge" - " in the 3 brightest pixels to apply the cleaning", + " in the 3 brightest pixels to apply the cleaning. Set to None if no" + " `bright_cleaning` should be applied.", + allow_none=True, ).tag(config=True) largest_island_only = BoolTelescopeParameter( default_value=False, help="Set to true to get only largest island" ).tag(config=True) - apply_pedestal_cleaning = BoolTelescopeParameter( - default_value=True, help="Set to true if pedestal cleaning should be applied" - ).tag(config=True) - pedestal_factor = FloatTelescopeParameter( default_value=2.5, help="Factor for interleaved pedestal cleaning. It is multiplied by the" @@ -789,14 +771,11 @@ def __call__( boundary_thresh=self.boundary_threshold_pe.tel[tel_id], min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id], keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], - apply_time_cleaning=self.apply_time_cleaning.tel[tel_id], time_limit=self.time_limit.tel[tel_id], time_num_neighbors=self.time_num_neighbors.tel[tel_id], - apply_bright_cleaning=self.apply_bright_cleaning.tel[tel_id], bright_cleaning_fraction=self.bright_cleaning_fraction.tel[tel_id], bright_cleaning_threshold=self.bright_cleaning_threshold.tel[tel_id], largest_island_only=self.largest_island_only.tel[tel_id], - apply_pedestal_cleaning=self.apply_pedestal_cleaning.tel[tel_id], pedestal_factor=self.pedestal_factor.tel[tel_id], pedestal_std=pedestal_std, ) diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 55c6cb0a336..56ab663cc40 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -501,13 +501,10 @@ def test_lst_image_cleaning(prod5_lst): "picture_thresh": 45, "boundary_thresh": 20, "keep_isolated_pixels": True, - "apply_time_cleaning": False, - "time_limit": 3, - "apply_bright_cleaning": False, + "time_limit": None, "bright_cleaning_fraction": 0.9, - "bright_cleaning_threshold": 40, + "bright_cleaning_threshold": None, "largest_island_only": False, - "apply_pedestal_cleaning": False, "pedestal_factor": 2, "pedestal_std": None, } @@ -521,7 +518,6 @@ def test_lst_image_cleaning(prod5_lst): # charge std. pedestal_std = np.ones(geom.n_pixels) pedestal_std[core_pixel_1] = 30 - args["apply_pedestal_cleaning"] = True args["pedestal_std"] = pedestal_std mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) @@ -539,16 +535,16 @@ def test_lst_image_cleaning(prod5_lst): # deselect core_pixel_1 with time_delta_cleaning by setting all its neighbors # peak_time to 10 peak_time[neighbors_1] = 10 - args["apply_pedestal_cleaning"] = False - args["apply_time_cleaning"] = True + args["pedestal_std"] = None + args["time_limit"] = 3 mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 6 + 6 # 3 brightest pixels are [50, 50, 29], so mean is 43. With a fraction of 0.9 # all boundaries should be deselected - args["apply_time_cleaning"] = False - args["apply_bright_cleaning"] = True + args["time_limit"] = None + args["bright_cleaning_threshold"] = 40 mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 1 @@ -556,7 +552,7 @@ def test_lst_image_cleaning(prod5_lst): # Set neighbors of core_pixel_2 to 0 so largest_island_only should select only # core_pixel_1 with its neighbors charge[neighbors_2] = 0 - args["apply_bright_cleaning"] = False + args["bright_cleaning_threshold"] = None args["largest_island_only"] = True mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) diff --git a/src/ctapipe/image/tests/test_image_cleaner_component.py b/src/ctapipe/image/tests/test_image_cleaner_component.py index cb3bb532c87..d2d443c4bbf 100644 --- a/src/ctapipe/image/tests/test_image_cleaner_component.py +++ b/src/ctapipe/image/tests/test_image_cleaner_component.py @@ -19,9 +19,8 @@ def test_image_cleaner(method, prod5_mst_nectarcam, reference_location): "LSTImageCleaner": { "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, - "apply_time_cleaning": False, - "apply_bright_cleaning": False, - "apply_pedestal_cleaning": False, + "time_limit": None, + "bright_cleaning_threshold": None, }, "MARSImageCleaner": { "boundary_threshold_pe": 5.0, From 3006fe00dd5950610d30ce3955b31435396020c0 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 26 Apr 2024 14:28:10 +0200 Subject: [PATCH 10/21] Substitute argsort with n_largest --- src/ctapipe/image/cleaning.py | 5 +++-- src/ctapipe/image/statistics.py | 10 +++++++++- src/ctapipe/image/tests/test_statistics.py | 10 ++++++++++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 78afed826c4..4bdf0cb456b 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -28,6 +28,8 @@ import numpy as np +from ctapipe.image.statistics import n_largest + from ..containers import MonitoringCameraContainer from ..core import TelescopeComponent from ..core.traits import ( @@ -135,8 +137,7 @@ def bright_cleaning(image, threshold, fraction): A boolean mask of *clean* pixels. """ - max_3_value_index = np.argsort(image)[-3:] - mean_3_max_signal = np.mean(image[max_3_value_index]) + mean_3_max_signal = np.mean(n_largest(3, image)) if mean_3_max_signal < threshold: return np.zeros(image.shape, dtype=bool) diff --git a/src/ctapipe/image/statistics.py b/src/ctapipe/image/statistics.py index ee8ddbd8aed..7ad6c7b2317 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -1,9 +1,11 @@ +from heapq import nlargest + import numpy as np from numba import njit from ..containers import StatisticsContainer -__all__ = ["descriptive_statistics", "skewness", "kurtosis"] +__all__ = ["n_largest", "descriptive_statistics", "skewness", "kurtosis"] @njit(cache=True) @@ -92,3 +94,9 @@ def descriptive_statistics( skewness=skewness(values, mean=mean, std=std), kurtosis=kurtosis(values, mean=mean, std=std), ) + + +@njit +def n_largest(n, array): + """return the n largest values of an array""" + return nlargest(n, array) diff --git a/src/ctapipe/image/tests/test_statistics.py b/src/ctapipe/image/tests/test_statistics.py index b9773edc9a7..600dbf6c2e2 100644 --- a/src/ctapipe/image/tests/test_statistics.py +++ b/src/ctapipe/image/tests/test_statistics.py @@ -60,3 +60,13 @@ def test_return_type(): stats = descriptive_statistics(data, container_class=PeakTimeStatisticsContainer) assert isinstance(stats, PeakTimeStatisticsContainer) + + +def test_nlargest(): + from ctapipe.image.statistics import n_largest + + image = np.random.rand(1855) + image[-3:] = 10 + + largest_3 = n_largest(3, image) + assert largest_3 == [10, 10, 10] From a40553dfba132dc7467a1be24bd5de0a4f80b770 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 26 Apr 2024 14:43:54 +0200 Subject: [PATCH 11/21] Set some advanced options to None by default --- src/ctapipe/image/cleaning.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 4bdf0cb456b..f09b00edf89 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -503,10 +503,10 @@ def lst_image_cleaning( boundary_thresh=4, min_number_picture_neighbors=2, keep_isolated_pixels=False, - time_limit=2, + time_limit=None, time_num_neighbors=1, bright_cleaning_fraction=0.03, - bright_cleaning_threshold=267, + bright_cleaning_threshold=None, largest_island_only=False, pedestal_factor=2.5, pedestal_std=None, From 794b2011e1f58c310e23663e4be4ccbbfb1735d6 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 26 Apr 2024 15:02:01 +0200 Subject: [PATCH 12/21] Rename and add reference --- src/ctapipe/image/cleaning.py | 15 ++++++++++----- src/ctapipe/image/tests/test_cleaning.py | 14 +++++++------- .../image/tests/test_image_cleaner_component.py | 2 +- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index f09b00edf89..0f56a38c4e2 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -495,7 +495,7 @@ def time_constrained_clean( return mask_core | mask_boundary -def lst_image_cleaning( +def nsb_image_cleaning( geom, image, arrival_times, @@ -660,7 +660,7 @@ def __call__( class TailcutsImageCleaner(ImageCleaner): """ Clean images using the standard picture/boundary technique. See - `ctapipe.image.tailcuts_clean` + `ctapipe.image.tailcuts_clean`. """ picture_threshold_pe = FloatTelescopeParameter( @@ -703,10 +703,15 @@ def __call__( ) -class LSTImageCleaner(TailcutsImageCleaner): +class NSBImageCleaner(TailcutsImageCleaner): """ - Clean images using lstchains image cleaning technique. See + Clean images based on lstchains image cleaning technique [1]_. See `ctapipe.image.lst_image_cleaning` + + References + ---------- + .. [1] https://arxiv.org/pdf/2306.12960 + """ time_limit = FloatTelescopeParameter( @@ -764,7 +769,7 @@ def __call__( if monitoring is not None: pedestal_std = monitoring.tel[tel_id].pedestal.charge_std - return lst_image_cleaning( + return nsb_image_cleaning( self.subarray.tel[tel_id].camera.geometry, image, arrival_times, diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 56ab663cc40..92e4b0b2019 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -480,7 +480,7 @@ def test_bright_cleaning(): assert np.count_nonzero(mask) == 0 -def test_lst_image_cleaning(prod5_lst): +def test_nsb_image_cleaning(prod5_lst): geom = prod5_lst.camera.geometry charge = np.zeros(geom.n_pixels, dtype=np.float32) peak_time = np.zeros(geom.n_pixels, dtype=np.float32) @@ -510,7 +510,7 @@ def test_lst_image_cleaning(prod5_lst): } # return normal tailcuts cleaning result if every other step is None/False: - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) # 2 * (1 core and 6 boundary_pixels) should be selected assert np.count_nonzero(mask) == 2 * (1 + 6) @@ -520,7 +520,7 @@ def test_lst_image_cleaning(prod5_lst): pedestal_std[core_pixel_1] = 30 args["pedestal_std"] = pedestal_std - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) # only core_pixel_2 with its boundaries should be selected since # pedestal_std[core_pixel_1] * pedestal_factor > charge[core_pixel_1] assert np.count_nonzero(mask) == 1 + 6 @@ -529,7 +529,7 @@ def test_lst_image_cleaning(prod5_lst): pedestal_std[core_pixel_2] = 30 args["pedestal_std"] = pedestal_std - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 0 # deselect core_pixel_1 with time_delta_cleaning by setting all its neighbors @@ -538,7 +538,7 @@ def test_lst_image_cleaning(prod5_lst): args["pedestal_std"] = None args["time_limit"] = 3 - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 6 + 6 # 3 brightest pixels are [50, 50, 29], so mean is 43. With a fraction of 0.9 @@ -546,7 +546,7 @@ def test_lst_image_cleaning(prod5_lst): args["time_limit"] = None args["bright_cleaning_threshold"] = 40 - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 1 # Set neighbors of core_pixel_2 to 0 so largest_island_only should select only @@ -555,5 +555,5 @@ def test_lst_image_cleaning(prod5_lst): args["bright_cleaning_threshold"] = None args["largest_island_only"] = True - mask = cleaning.lst_image_cleaning(geom, charge, peak_time, **args) + mask = cleaning.nsb_image_cleaning(geom, charge, peak_time, **args) assert np.count_nonzero(mask) == 1 + 6 diff --git a/src/ctapipe/image/tests/test_image_cleaner_component.py b/src/ctapipe/image/tests/test_image_cleaner_component.py index d2d443c4bbf..e0f34829a0e 100644 --- a/src/ctapipe/image/tests/test_image_cleaner_component.py +++ b/src/ctapipe/image/tests/test_image_cleaner_component.py @@ -16,7 +16,7 @@ def test_image_cleaner(method, prod5_mst_nectarcam, reference_location): "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, }, - "LSTImageCleaner": { + "NSBImageCleaner": { "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, "time_limit": None, From 6e1e60877dab07e81969d4b4839e13aa7d832387 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 26 Apr 2024 16:58:55 +0200 Subject: [PATCH 13/21] Add arg_n_largest function - substitute it with the argsort term in the GlobalPeakWindowSum Co-authored-by: Maximilian Linhoff --- src/ctapipe/image/extractor.py | 13 ++++++++----- src/ctapipe/image/statistics.py | 13 +++++++++++++ src/ctapipe/image/tests/test_statistics.py | 16 ++++++++++++++-- 3 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 2ba9281ee9c..8dbd9dc11a7 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -47,6 +47,7 @@ from .hillas import camera_to_shower_coordinates, hillas_parameters from .invalid_pixels import InvalidPixelHandler from .morphology import brightest_island, number_of_islands +from .statistics import arg_n_largest from .timing import timing_parameters @@ -585,11 +586,13 @@ def __call__( ).argmax(axis=-1) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort( - waveforms.max( - axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf - ) - )[..., -n_pixels:] + n_channels = waveforms.shape[0] + brightest = np.empty((n_channels, n_pixels), dtype=np.int64) + waveforms_max = waveforms.max( + axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf + ) + for channel in range(n_channels): + brightest[channel] = arg_n_largest(n_pixels, waveforms_max[channel]) # average over brightest pixels then argmax over samples peak_index = ( diff --git a/src/ctapipe/image/statistics.py b/src/ctapipe/image/statistics.py index 7ad6c7b2317..d3b9091087c 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -100,3 +100,16 @@ def descriptive_statistics( def n_largest(n, array): """return the n largest values of an array""" return nlargest(n, array) + + +@njit +def arg_n_largest(n, array): + """return the indices of the n largest values of an array""" + values = [] + for i in range(len(array)): + values.append((array[i], i)) + largest = n_largest(n, values) + idx = np.empty(n, dtype=np.int64) + for i in range(n): + idx[i] = largest[i][1] + return idx diff --git a/src/ctapipe/image/tests/test_statistics.py b/src/ctapipe/image/tests/test_statistics.py index 600dbf6c2e2..006c9c0dbce 100644 --- a/src/ctapipe/image/tests/test_statistics.py +++ b/src/ctapipe/image/tests/test_statistics.py @@ -62,11 +62,23 @@ def test_return_type(): assert isinstance(stats, PeakTimeStatisticsContainer) -def test_nlargest(): +def test_n_largest(): from ctapipe.image.statistics import n_largest - image = np.random.rand(1855) + rng = np.random.default_rng(0) + image = rng.random(1855) image[-3:] = 10 largest_3 = n_largest(3, image) assert largest_3 == [10, 10, 10] + + +def test_arg_n_largest(): + from ctapipe.image.statistics import arg_n_largest + + rng = np.random.default_rng(0) + image = rng.random(1855) + image[-3:] = 10 + + largest_3 = arg_n_largest(3, image) + assert (largest_3 == [1854, 1853, 1852]).all() From 4bfb532230ad9ea704bd728be1f83b9f6b057841 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 2 May 2024 13:21:45 +0200 Subject: [PATCH 14/21] Change to guvectorize - Adding a wrapper --- src/ctapipe/image/extractor.py | 11 +++++------ src/ctapipe/image/statistics.py | 29 ++++++++++++++++++++++++----- 2 files changed, 29 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 8dbd9dc11a7..ce772744c58 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -586,13 +586,12 @@ def __call__( ).argmax(axis=-1) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - n_channels = waveforms.shape[0] - brightest = np.empty((n_channels, n_pixels), dtype=np.int64) - waveforms_max = waveforms.max( - axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf + brightest = arg_n_largest( + n_pixels, + waveforms.max( + axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf + ), ) - for channel in range(n_channels): - brightest[channel] = arg_n_largest(n_pixels, waveforms_max[channel]) # average over brightest pixels then argmax over samples peak_index = ( diff --git a/src/ctapipe/image/statistics.py b/src/ctapipe/image/statistics.py index d3b9091087c..37b5166be18 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -1,7 +1,7 @@ from heapq import nlargest import numpy as np -from numba import njit +from numba import float32, float64, guvectorize, int64, njit from ..containers import StatisticsContainer @@ -102,14 +102,33 @@ def n_largest(n, array): return nlargest(n, array) -@njit -def arg_n_largest(n, array): - """return the indices of the n largest values of an array""" +@guvectorize( + [ + (int64[:], float64[:], int64[:]), + (int64[:], float32[:], int64[:]), + ], + "(n),(p)->(n)", + nopython=True, + cache=True, +) +def arg_n_largest_gu(dummy, array, idx): + """ + Returns the indices of the len(dummy) largest values of an array. + + Used by the `arg_n_largest` wrapper to return the indices of the n largest + values of an array. + """ values = [] for i in range(len(array)): values.append((array[i], i)) + n = len(dummy) largest = n_largest(n, values) - idx = np.empty(n, dtype=np.int64) for i in range(n): idx[i] = largest[i][1] + + +def arg_n_largest(n, array): + """Return the indices of the n largest values of an array""" + dummy = np.zeros(n, dtype=np.int64) + idx = arg_n_largest_gu(dummy, array) return idx From 83d3ca20f504f5400e9c5939767a0114a753d39c Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 6 May 2024 16:16:15 +0200 Subject: [PATCH 15/21] Minor fix --- src/ctapipe/image/cleaning.py | 3 ++- src/ctapipe/image/statistics.py | 9 ++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 0f56a38c4e2..084c237029b 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -157,8 +157,9 @@ def mars_cleaning_1st_pass( min_number_picture_neighbors=0, ): """ - Clean an image by select pixels that pass a three-threshold tail-cuts + Clean an image by selecting pixels that pass a three-threshold tail-cuts procedure. + All thresholds are defined with respect to the pedestal dispersion. All pixels that have a signal higher than the core (picture) threshold will be retained, along with all those above the boundary diff --git a/src/ctapipe/image/statistics.py b/src/ctapipe/image/statistics.py index 37b5166be18..bd38f02a377 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -5,7 +5,14 @@ from ..containers import StatisticsContainer -__all__ = ["n_largest", "descriptive_statistics", "skewness", "kurtosis"] +__all__ = [ + "arg_n_largest", + "arg_n_largest_gu", + "n_largest", + "descriptive_statistics", + "skewness", + "kurtosis", +] @njit(cache=True) From 5c2e00c5a712bd7fd570a82db37f43529a4d0b87 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 3 Jun 2024 14:03:17 +0200 Subject: [PATCH 16/21] Update docstrings --- src/ctapipe/image/cleaning.py | 137 +++++++++++++++++----------------- 1 file changed, 68 insertions(+), 69 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 084c237029b..7f399b7e244 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -60,20 +60,20 @@ def tailcuts_clean( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - image: array - pixel values - picture_thresh: float or array + image : np.ndarray + pixel charges + picture_thresh : float | np.ndarray threshold above which all pixels are retained - boundary_thresh: float or array + boundary_thresh : float | np.ndarray threshold above which pixels are retained if they have a neighbor already above the picture_thresh - keep_isolated_pixels: bool + keep_isolated_pixels : bool If True, pixels above the picture threshold will be included always, if not they are only included if a neighbor is in the picture or boundary - min_number_picture_neighbors: int + min_number_picture_neighbors : int A picture pixel survives cleaning only if it has at least this number of picture neighbors. This has no effect in case keep_isolated_pixels is True @@ -123,12 +123,12 @@ def bright_cleaning(image, threshold, fraction): Parameters ---------- - image: array - pixel values - threshold: `float` + image : np.ndarray + pixel charges + threshold : float Minimum average charge in the 3 brightest pixels to apply cleaning - fraction: `float` + fraction : float Pixels below fraction * (average charge in the 3 brightest pixels) will be removed from the cleaned image @@ -168,21 +168,21 @@ def mars_cleaning_1st_pass( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - image: array - pixel values - picture_thresh: float + image : np.ndarray + pixel charges + picture_thresh : float threshold above which all pixels are retained - boundary_thresh: float + boundary_thresh : float threshold above which pixels are retained if they have a neighbor already above the picture_thresh; it is then reapplied iteratively to the neighbor of the neighbor - keep_isolated_pixels: bool + keep_isolated_pixels : bool If True, pixels above the picture threshold will be included always, if not they are only included if a neighbor is in the picture or boundary - min_number_picture_neighbors: int + min_number_picture_neighbors : int A picture pixel survives cleaning only if it has at least this number of picture neighbors. This has no effect in case keep_isolated_pixels is True @@ -231,14 +231,15 @@ def dilate(geom, mask): """ Add one row of neighbors to the True values of a pixel mask and return the new mask. + This can be used to include extra rows of pixels in a mask that was pre-computed, e.g. via `tailcuts_clean`. Parameters ---------- - geom: `~ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - mask: ndarray + mask : np.ndarray input mask (array of booleans) to be dilated """ return mask | geom.neighbor_matrix_sparse.dot(mask) @@ -253,21 +254,20 @@ def apply_time_delta_cleaning( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - mask: array, boolean + mask : np.ndarray boolean mask of *clean* pixels before time_delta_cleaning - arrival_times: array + arrival_times : np.ndarray pixel timing information - min_number_neighbors: int + min_number_neighbors : int a selected pixel needs at least this number of (already selected) neighbors that arrived within a given time_limit to itself to survive the cleaning. - time_limit: int or float + time_limit : int | float arrival time limit for neighboring pixels Returns ------- - A boolean mask of *clean* pixels. """ pixels_to_keep = mask.copy() @@ -291,22 +291,21 @@ def apply_time_average_cleaning( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - mask: array, boolean + mask : np.ndarray boolean mask of *clean* pixels before time_delta_cleaning - image: array - pixel values - arrival_times: array + image : np.ndarray + pixel charges + arrival_times : np.ndarray pixel timing information - picture_thresh: float + picture_thresh : float threshold above which time limit is extended twice its value - time_limit: int or float + time_limit : int | float arrival time limit w.r.t. the average time of the core pixels Returns ------- - A boolean mask of *clean* pixels. """ mask = mask.copy() @@ -347,21 +346,21 @@ def fact_image_cleaning( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - image: array - pixel values - arrival_times: array + image : np.ndarray + pixel charges + arrival_times : np.ndarray pixel timing information - picture_threshold: float or array + picture_threshold : float | np.ndarray threshold above which all pixels are retained - boundary_threshold: float or array + boundary_threshold : float | array threshold above which pixels are retained if they have a neighbor already above the picture_thresh - min_number_neighbors: int + min_number_neighbors : int Threshold to determine if a pixel survives cleaning steps. These steps include checks of neighbor arrival time and value - time_limit: int or float + time_limit : int | float arrival time limit for neighboring pixels Returns @@ -421,7 +420,7 @@ def time_constrained_clean( min_number_picture_neighbors=1, ): """ - time constrained cleaning by MAGIC + Time constrained cleaning by MAGIC Cleaning contains the following steps: - Find core pixels (containing more photons than a picture threshold) @@ -432,22 +431,22 @@ def time_constrained_clean( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - image: array - pixel values - arrival_times: array + image : np.ndarray + pixel charges + arrival_times : np.ndarray pixel timing information - picture_threshold: float or array + picture_threshold : float | array threshold above which all pixels are retained - boundary_threshold: float or array + boundary_threshold : float | array threshold above which pixels are retained if they have a neighbor already above the picture_thresh - time_limit_core: int or float + time_limit_core : int | float arrival time limit of core pixels w.r.t the average time - time_limit_boundary: int or float + time_limit_boundary : int | float arrival time limit of boundary pixels w.r.t neighboring core pixels - min_number_neighbors: int + min_number_neighbors : int Threshold to determine if a pixel survives cleaning steps. These steps include checks of neighbor arrival time and value @@ -527,48 +526,48 @@ def nsb_image_cleaning( Parameters ---------- - geom: `ctapipe.instrument.CameraGeometry` + geom : ctapipe.instrument.CameraGeometry Camera geometry information - image: `np.ndarray` + image : np.ndarray pixel charges - arrival_times: `np.ndarray` + arrival_times : np.ndarray Pixel timing information - picture_thresh: `float` or `np.ndarray` + picture_thresh : float | np.ndarray threshold above which all pixels are retained. Used for `tailcuts_clean`. - boundary_thresh: `float` or `np.ndarray` + boundary_thresh : float | np.ndarray threshold above which pixels are retained if they have a neighbor already above the picture_thresh. Used for `tailcuts_clean`. - min_number_picture_neighbors: `int` + min_number_picture_neighbors : int A picture pixel survives cleaning only if it has at least this number of picture neighbors. This has no effect in case keep_isolated_pixels is True. Used for `tailcuts_clean`. - keep_isolated_pixels: `bool` + keep_isolated_pixels : bool If True, pixels above the picture threshold will be included always, if not they are only included if a neighbor is in the picture or boundary. Used for `tailcuts_clean`. - time_limit: `float` + time_limit : float Time limit for the `time_delta_cleaning`. Set to None if no `time_delta_cleaning` should be applied. - time_num_neighbors: int + time_num_neighbors : int Used for `time_delta_cleaning`. A selected pixel needs at least this number of (already selected) neighbors that arrived within a given time_limit to itself to survive this cleaning. - bright_cleaning_fraction: `float` + bright_cleaning_fraction : float `fraction` parameter for `bright_cleaning`. Pixels below fraction * (average charge in the 3 brightest pixels) will be removed from the cleaned image. Set `bright_cleaning_threshold` to None if no `bright_cleaning` should be applied. - bright_cleaning_threshold: `float` + bright_cleaning_threshold : float `threshold` parameter for `bright_cleaning`. Minimum average charge in the 3 brightest pixels to apply the cleaning. Set to None if no `bright_cleaning` should be applied. - largest_island_only: `bool` + largest_island_only : bool Set to true to get only largest island. - pedestal_factor: `float` + pedestal_factor : float Factor for interleaved pedestal cleaning. It is multiplied by the pedestal standard deviation for each pixel to calculate pixelwise picture threshold parameters for `tailcuts_clean` considering the current background. - pedestal_std: `np.ndarray` + pedestal_std : np.ndarray | None Pedestal standard deviation for each pixel. See `ctapipe.containers.PedestalContainer` @@ -639,14 +638,14 @@ def __call__( Parameters ---------- - tel_id: int + tel_id : int which telescope id in the subarray is being used (determines which cut is used) image : np.ndarray image pixel data corresponding to the camera geometry - arrival_times: np.ndarray + arrival_times : np.ndarray image of arrival time (not used in this method) - monitoring: `ctapipe.containers.MonitoringCameraContainer` + monitoring : ctapipe.containers.MonitoringCameraContainer MonitoringCameraContainer to make use of additional parameters from monitoring data e.g. pedestal std. From dd2dcafdc4721a55dbb4020fd13e114199b23375 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 3 Jun 2024 14:45:00 +0200 Subject: [PATCH 17/21] Rename function parameter --- src/ctapipe/image/cleaning.py | 39 ++++++++++++++---------- src/ctapipe/image/tests/test_cleaning.py | 2 +- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 7f399b7e244..e9ae632bf29 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -499,7 +499,7 @@ def nsb_image_cleaning( geom, image, arrival_times, - picture_thresh=8, + picture_thresh_min=8, boundary_thresh=4, min_number_picture_neighbors=2, keep_isolated_pixels=False, @@ -515,7 +515,7 @@ def nsb_image_cleaning( Clean an image in 5 Steps: 1) Get picture thresholds for `tailcuts_clean` in step 2) from interleaved - pedestal events `pedestal_std` is not None. + pedestal events if `pedestal_std` is not None. 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. 3) Apply time_delta_cleaning algorithm - `ctapipe.image.cleaning.apply_time_delta_cleaning` if time_limit is not None. @@ -529,14 +529,16 @@ def nsb_image_cleaning( geom : ctapipe.instrument.CameraGeometry Camera geometry information image : np.ndarray - pixel charges + Pixel charges arrival_times : np.ndarray Pixel timing information - picture_thresh : float | np.ndarray - threshold above which all pixels are retained. Used for `tailcuts_clean`. + picture_thresh_min : float | np.ndarray + Defines the minimum value used for the picture threshold for `tailcuts_clean`. + The threshold used will be at least this value, or higher if `pedestal_std` + and `pedestal_factor` are set. boundary_thresh : float | np.ndarray - threshold above which pixels are retained if they have a neighbor - already above the picture_thresh. Used for `tailcuts_clean`. + Threshold above which pixels are retained if they have a neighbor + already above the picture_thresh_min. Used for `tailcuts_clean`. min_number_picture_neighbors : int A picture pixel survives cleaning only if it has at least this number of picture neighbors. This has no effect in case keep_isolated_pixels is True. @@ -553,12 +555,12 @@ def nsb_image_cleaning( A selected pixel needs at least this number of (already selected) neighbors that arrived within a given time_limit to itself to survive this cleaning. bright_cleaning_fraction : float - `fraction` parameter for `bright_cleaning`. Pixels below + Fraction parameter for `bright_cleaning`. Pixels below fraction * (average charge in the 3 brightest pixels) will be removed from the cleaned image. Set `bright_cleaning_threshold` to None if no `bright_cleaning` should be applied. bright_cleaning_threshold : float - `threshold` parameter for `bright_cleaning`. Minimum average charge + Threshold parameter for `bright_cleaning`. Minimum average charge in the 3 brightest pixels to apply the cleaning. Set to None if no `bright_cleaning` should be applied. largest_island_only : bool @@ -567,9 +569,13 @@ def nsb_image_cleaning( Factor for interleaved pedestal cleaning. It is multiplied by the pedestal standard deviation for each pixel to calculate pixelwise picture threshold parameters for `tailcuts_clean` considering the current background. + Has no effect if `pedestal_std` is set to None. pedestal_std : np.ndarray | None Pedestal standard deviation for each pixel. See - `ctapipe.containers.PedestalContainer` + `ctapipe.containers.PedestalContainer`. Used to calculate pixelwise picture + threshold parameters for `tailcuts_clean` by multiplying it with `pedestal_factor` + and clip (limit) the product with picture_thresh_min. If set to None, only + `picture_thresh_min` is used to set the picture threshold for `tailcuts_clean`. Returns ------- @@ -577,9 +583,10 @@ def nsb_image_cleaning( """ # Step 1 + picture_thresh = picture_thresh_min if pedestal_std is not None: pedestal_threshold = pedestal_std * pedestal_factor - picture_thresh = np.clip(pedestal_threshold, picture_thresh, None) + picture_thresh = np.clip(pedestal_threshold, picture_thresh_min, None) # Step 2 mask = tailcuts_clean( @@ -706,7 +713,7 @@ def __call__( class NSBImageCleaner(TailcutsImageCleaner): """ Clean images based on lstchains image cleaning technique [1]_. See - `ctapipe.image.lst_image_cleaning` + `ctapipe.image.nsb_image_cleaning` References ---------- @@ -725,19 +732,19 @@ class NSBImageCleaner(TailcutsImageCleaner): default_value=1, help="Used for `time_delta_cleaning`." " A selected pixel needs at least this number of (already selected) neighbors" - " that arrived within a given time_limit to itself to survive this cleaning.", + " that arrived within a given `time_limit` to itself to survive this cleaning.", ).tag(config=True) bright_cleaning_fraction = FloatTelescopeParameter( default_value=0.03, - help="`fraction` parameter for `bright_cleaning`. Pixels below" + help="Fraction parameter for `bright_cleaning`. Pixels below" " fraction * (average charge in the 3 brightest pixels) will be removed from" " the cleaned image", ).tag(config=True) bright_cleaning_threshold = FloatTelescopeParameter( default_value=267, - help="`threshold` parameter for `bright_cleaning`. Minimum average charge" + help="Threshold parameter for `bright_cleaning`. Minimum average charge" " in the 3 brightest pixels to apply the cleaning. Set to None if no" " `bright_cleaning` should be applied.", allow_none=True, @@ -773,7 +780,7 @@ def __call__( self.subarray.tel[tel_id].camera.geometry, image, arrival_times, - picture_thresh=self.picture_threshold_pe.tel[tel_id], + picture_thresh_min=self.picture_threshold_pe.tel[tel_id], boundary_thresh=self.boundary_threshold_pe.tel[tel_id], min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id], keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 92e4b0b2019..e773e4c0b37 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -498,7 +498,7 @@ def test_nsb_image_cleaning(prod5_lst): charge[core_pixel_2] = 50 args = { - "picture_thresh": 45, + "picture_thresh_min": 45, "boundary_thresh": 20, "keep_isolated_pixels": True, "time_limit": None, From 98616079bac09e5f0ab228cf497851526e76d152 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 3 Jun 2024 15:28:56 +0200 Subject: [PATCH 18/21] Fix bright cleaning --- src/ctapipe/image/cleaning.py | 14 +++++++------- src/ctapipe/image/tests/test_cleaning.py | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index e9ae632bf29..8ba0d158e4c 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -113,12 +113,12 @@ def tailcuts_clean( ) -def bright_cleaning(image, threshold, fraction): +def bright_cleaning(image, threshold, fraction, n_pixels=3): """ Clean an image by removing pixels below a fraction of the mean charge - in the 3 brightest pixels. + in the `n_pixels` brightest pixels. - Select no pixels instead if the mean charge of the brightest pixels + No pixels are removed instead if the mean charge of the brightest pixels are below a certain threshold. Parameters @@ -137,12 +137,12 @@ def bright_cleaning(image, threshold, fraction): A boolean mask of *clean* pixels. """ - mean_3_max_signal = np.mean(n_largest(3, image)) + mean_brightest_signal = np.mean(n_largest(n_pixels, image)) - if mean_3_max_signal < threshold: - return np.zeros(image.shape, dtype=bool) + if mean_brightest_signal < threshold: + return np.ones(image.shape, dtype=bool) - threshold_brightest = fraction * mean_3_max_signal + threshold_brightest = fraction * mean_brightest_signal mask = image >= threshold_brightest return mask diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index e773e4c0b37..899b81b6208 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -471,13 +471,13 @@ def test_bright_cleaning(): image[50:65] = 30 threshold = 90 - mask = cleaning.bright_cleaning(image, threshold, fraction) + mask = cleaning.bright_cleaning(image, threshold, fraction, n_pixels=3) assert np.count_nonzero(mask) == 3 + 10 # test that it doesn't select any pixels if mean of the 3 brightest pixels # is below threshold threshold = 110 - mask = cleaning.bright_cleaning(image, threshold, fraction) - assert np.count_nonzero(mask) == 0 + mask = cleaning.bright_cleaning(image, threshold, fraction, n_pixels=3) + assert np.count_nonzero(~mask) == 0 def test_nsb_image_cleaning(prod5_lst): From e466d4e648d208dffcff17ab2d3dc16d810ffecd Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 5 Jun 2024 10:50:06 +0200 Subject: [PATCH 19/21] Add bright_cleaning_n_pixels traitlet - adapt docstrings --- src/ctapipe/image/cleaning.py | 63 +++++++++++++++++------- src/ctapipe/image/tests/test_cleaning.py | 1 + 2 files changed, 45 insertions(+), 19 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 8ba0d158e4c..e5e340cccaa 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -126,11 +126,13 @@ def bright_cleaning(image, threshold, fraction, n_pixels=3): image : np.ndarray pixel charges threshold : float - Minimum average charge in the 3 brightest pixels to apply + Minimum average charge in the `n_pixels` brightest pixels to apply cleaning fraction : float - Pixels below fraction * (average charge in the 3 brightest pixels) + Pixels below fraction * (average charge in the `n_pixels` brightest pixels) will be removed from the cleaned image + n_pixels : int + Consider this number of the brightest pixels to calculate the mean charge Returns ------- @@ -505,6 +507,7 @@ def nsb_image_cleaning( keep_isolated_pixels=False, time_limit=None, time_num_neighbors=1, + bright_cleaning_n_pixels=3, bright_cleaning_fraction=0.03, bright_cleaning_threshold=None, largest_island_only=False, @@ -514,7 +517,7 @@ def nsb_image_cleaning( """ Clean an image in 5 Steps: - 1) Get picture thresholds for `tailcuts_clean` in step 2) from interleaved + 1) Get pixelwise picture thresholds for `tailcuts_clean` in step 2) from interleaved pedestal events if `pedestal_std` is not None. 2) Apply tailcuts image cleaning algorithm - `ctapipe.image.cleaning.tailcuts_clean`. 3) Apply time_delta_cleaning algorithm - @@ -554,15 +557,20 @@ def nsb_image_cleaning( Used for `time_delta_cleaning`. A selected pixel needs at least this number of (already selected) neighbors that arrived within a given time_limit to itself to survive this cleaning. + bright_cleaning_n_pixels : int + Consider this number of the brightest pixels for calculating the mean charge. + Pixels below fraction * (mean charge in the `bright_cleaning_n_pixels` + brightest pixels) will be removed from the cleaned image. Set + `bright_cleaning_threshold` to None if no `bright_cleaning` should be applied. bright_cleaning_fraction : float Fraction parameter for `bright_cleaning`. Pixels below - fraction * (average charge in the 3 brightest pixels) will be removed from - the cleaned image. Set `bright_cleaning_threshold` to None if no - `bright_cleaning` should be applied. + fraction * (mean charge in the `bright_cleaning_n_pixels` brightest pixels) + will be removed from the cleaned image. Set `bright_cleaning_threshold` to None + if no `bright_cleaning` should be applied. bright_cleaning_threshold : float - Threshold parameter for `bright_cleaning`. Minimum average charge - in the 3 brightest pixels to apply the cleaning. Set to None if no - `bright_cleaning` should be applied. + Threshold parameter for `bright_cleaning`. Minimum mean charge + in the `bright_cleaning_n_pixels` brightest pixels to apply the cleaning. + Set to None if no `bright_cleaning` should be applied. largest_island_only : bool Set to true to get only largest island. pedestal_factor : float @@ -574,7 +582,7 @@ def nsb_image_cleaning( Pedestal standard deviation for each pixel. See `ctapipe.containers.PedestalContainer`. Used to calculate pixelwise picture threshold parameters for `tailcuts_clean` by multiplying it with `pedestal_factor` - and clip (limit) the product with picture_thresh_min. If set to None, only + and clip (limit) the product with `picture_thresh_min`. If set to None, only `picture_thresh_min` is used to set the picture threshold for `tailcuts_clean`. Returns @@ -614,7 +622,10 @@ def nsb_image_cleaning( # Step 4 if bright_cleaning_threshold is not None: mask &= bright_cleaning( - image, bright_cleaning_threshold, bright_cleaning_fraction + image, + bright_cleaning_threshold, + bright_cleaning_fraction, + bright_cleaning_n_pixels, ) # Step 5 @@ -735,18 +746,28 @@ class NSBImageCleaner(TailcutsImageCleaner): " that arrived within a given `time_limit` to itself to survive this cleaning.", ).tag(config=True) + bright_cleaning_n_pixels = IntTelescopeParameter( + default_value=3, + help="Consider this number of the brightest pixels for calculating the" + " mean charge. Pixels below fraction * (mean charge in the" + " `bright_cleaning_n_pixels` brightest pixels) will be removed from the" + " cleaned image. Set `bright_cleaning_threshold` to None if no" + " `bright_cleaning` should be applied.", + ).tag(config=True) + bright_cleaning_fraction = FloatTelescopeParameter( default_value=0.03, help="Fraction parameter for `bright_cleaning`. Pixels below" - " fraction * (average charge in the 3 brightest pixels) will be removed from" - " the cleaned image", + " fraction * (mean charge in the `bright_cleaning_n_pixels` brightest pixels)" + " will be removed from the cleaned image. Set `bright_cleaning_threshold` to" + " None if no `bright_cleaning` should be applied.", ).tag(config=True) bright_cleaning_threshold = FloatTelescopeParameter( default_value=267, - help="Threshold parameter for `bright_cleaning`. Minimum average charge" - " in the 3 brightest pixels to apply the cleaning. Set to None if no" - " `bright_cleaning` should be applied.", + help="Threshold parameter for `bright_cleaning`. Minimum mean charge" + " in the `bright_cleaning_n_pixels` brightest pixels to apply the cleaning." + " Set to None if no `bright_cleaning` should be applied.", allow_none=True, ).tag(config=True) @@ -757,8 +778,12 @@ class NSBImageCleaner(TailcutsImageCleaner): pedestal_factor = FloatTelescopeParameter( default_value=2.5, help="Factor for interleaved pedestal cleaning. It is multiplied by the" - " pedestal standard deviation for each pixel to calculate pixelwise picture" - " threshold parameters for `tailcuts_clean` considering the current background.", + " pedestal standard deviation for each pixel to calculate pixelwise upper limit" + " picture thresholds and clip them with `picture_thresh_pe` of" + " `TailcutsImageCleaner` for `tailcuts_clean` considering the current background." + " If no pedestal standard deviation is given, interleaved pedestal cleaning is" + " not applied and `picture_thresh_pe` of `TailcutsImageCleaner` is used alone" + " instead.", ).tag(config=True) def __call__( @@ -770,7 +795,7 @@ def __call__( monitoring: MonitoringCameraContainer = None, ) -> np.ndarray: """ - Apply LST image cleaning. See `ImageCleaner.__call__()` + Apply NSB image cleaning used by lstchain. See `ImageCleaner.__call__()` """ pedestal_std = None if monitoring is not None: diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 899b81b6208..159fca60692 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -502,6 +502,7 @@ def test_nsb_image_cleaning(prod5_lst): "boundary_thresh": 20, "keep_isolated_pixels": True, "time_limit": None, + "bright_cleaning_n_pixels": 3, "bright_cleaning_fraction": 0.9, "bright_cleaning_threshold": None, "largest_island_only": False, From 5f64042b6e52c65df2ebd26e7ec793b71ba46864 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 5 Jun 2024 11:10:41 +0200 Subject: [PATCH 20/21] Adapt docstrings --- src/ctapipe/image/cleaning.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index e5e340cccaa..a2be122d36d 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -356,7 +356,7 @@ def fact_image_cleaning( pixel timing information picture_threshold : float | np.ndarray threshold above which all pixels are retained - boundary_threshold : float | array + boundary_threshold : float | np.ndarray threshold above which pixels are retained if they have a neighbor already above the picture_thresh min_number_neighbors : int @@ -439,9 +439,9 @@ def time_constrained_clean( pixel charges arrival_times : np.ndarray pixel timing information - picture_threshold : float | array + picture_threshold : float | np.ndarray threshold above which all pixels are retained - boundary_threshold : float | array + boundary_threshold : float | np.ndarray threshold above which pixels are retained if they have a neighbor already above the picture_thresh time_limit_core : int | float From e58699eecd6427db1786cb3d12e696b8abff322e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 5 Jun 2024 14:08:48 +0200 Subject: [PATCH 21/21] Add missing parameter --- src/ctapipe/image/cleaning.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index a2be122d36d..e96b2bef48f 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -811,6 +811,7 @@ def __call__( keep_isolated_pixels=self.keep_isolated_pixels.tel[tel_id], time_limit=self.time_limit.tel[tel_id], time_num_neighbors=self.time_num_neighbors.tel[tel_id], + bright_cleaning_n_pixels=self.bright_cleaning_n_pixels.tel[tel_id], bright_cleaning_fraction=self.bright_cleaning_fraction.tel[tel_id], bright_cleaning_threshold=self.bright_cleaning_threshold.tel[tel_id], largest_island_only=self.largest_island_only.tel[tel_id],