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. diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index 6dc8ec54c62..e96b2bef48f 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -28,15 +28,16 @@ import numpy as np -from ctapipe.containers import MonitoringCameraContainer +from ctapipe.image.statistics import n_largest +from ..containers import MonitoringCameraContainer from ..core import TelescopeComponent from ..core.traits import ( BoolTelescopeParameter, FloatTelescopeParameter, IntTelescopeParameter, ) -from .morphology import brightest_island, number_of_islands +from .morphology import brightest_island, largest_island, number_of_islands def tailcuts_clean( @@ -59,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 @@ -112,6 +113,43 @@ def tailcuts_clean( ) +def bright_cleaning(image, threshold, fraction, n_pixels=3): + """ + Clean an image by removing pixels below a fraction of the mean charge + in the `n_pixels` brightest pixels. + + No pixels are removed instead if the mean charge of the brightest pixels + are below a certain threshold. + + Parameters + ---------- + image : np.ndarray + pixel charges + threshold : float + Minimum average charge in the `n_pixels` brightest pixels to apply + cleaning + fraction : float + 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 + ------- + A boolean mask of *clean* pixels. + + """ + mean_brightest_signal = np.mean(n_largest(n_pixels, image)) + + if mean_brightest_signal < threshold: + return np.ones(image.shape, dtype=bool) + + threshold_brightest = fraction * mean_brightest_signal + mask = image >= threshold_brightest + + return mask + + def mars_cleaning_1st_pass( geom, image, @@ -121,8 +159,9 @@ 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 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 @@ -131,21 +170,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 @@ -194,14 +233,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) @@ -216,21 +256,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 - Threshold to determine if a pixel survives cleaning steps. - These steps include checks of neighbor arrival time and value - time_limit: int or float + 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 | float arrival time limit for neighboring pixels Returns ------- - A boolean mask of *clean* pixels. """ pixels_to_keep = mask.copy() @@ -254,22 +293,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() @@ -310,21 +348,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 | np.ndarray 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 @@ -384,7 +422,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) @@ -395,22 +433,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 | np.ndarray threshold above which all pixels are retained - boundary_threshold: float or 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 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 @@ -459,6 +497,145 @@ def time_constrained_clean( return mask_core | mask_boundary +def nsb_image_cleaning( + geom, + image, + arrival_times, + picture_thresh_min=8, + boundary_thresh=4, + min_number_picture_neighbors=2, + 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, + pedestal_factor=2.5, + pedestal_std=None, +): + """ + Clean an image in 5 Steps: + + 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 - + `ctapipe.image.cleaning.apply_time_delta_cleaning` if time_limit is not None. + 4) Apply bright_cleaning - `ctapipe.image.cleaning.bright_cleaning` if + bright_cleaning_threshold is not None. + 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_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_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. + 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. + 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_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 * (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 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 + 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`. 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 + ------- + A boolean mask of *clean* pixels. + + """ + # 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_min, 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 np.count_nonzero(mask) == 0: + return mask + + # Step 3 + if time_limit is not None: + mask = apply_time_delta_cleaning( + geom, + mask, + arrival_times, + min_number_neighbors=time_num_neighbors, + time_limit=time_limit, + ) + + # Step 4 + if bright_cleaning_threshold is not None: + mask &= bright_cleaning( + image, + bright_cleaning_threshold, + bright_cleaning_fraction, + bright_cleaning_n_pixels, + ) + + # 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 @@ -479,14 +656,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. @@ -501,7 +678,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( @@ -544,6 +721,105 @@ def __call__( ) +class NSBImageCleaner(TailcutsImageCleaner): + """ + Clean images based on lstchains image cleaning technique [1]_. See + `ctapipe.image.nsb_image_cleaning` + + References + ---------- + .. [1] https://arxiv.org/pdf/2306.12960 + + """ + + 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.", + allow_none=True, + ).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) + + 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 * (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 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) + + 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 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__( + self, + tel_id: int, + image: np.ndarray, + arrival_times: np.ndarray = None, + *, + monitoring: MonitoringCameraContainer = None, + ) -> np.ndarray: + """ + Apply NSB image cleaning used by lstchain. See `ImageCleaner.__call__()` + """ + pedestal_std = None + if monitoring is not None: + pedestal_std = monitoring.tel[tel_id].pedestal.charge_std + + return nsb_image_cleaning( + self.subarray.tel[tel_id].camera.geometry, + image, + arrival_times, + 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], + 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], + 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`) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 2ba9281ee9c..ce772744c58 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,12 @@ def __call__( ).argmax(axis=-1) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort( + brightest = arg_n_largest( + n_pixels, waveforms.max( axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf - ) - )[..., -n_pixels:] + ), + ) # 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 ee8ddbd8aed..bd38f02a377 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -1,9 +1,18 @@ +from heapq import nlargest + import numpy as np -from numba import njit +from numba import float32, float64, guvectorize, int64, njit from ..containers import StatisticsContainer -__all__ = ["descriptive_statistics", "skewness", "kurtosis"] +__all__ = [ + "arg_n_largest", + "arg_n_largest_gu", + "n_largest", + "descriptive_statistics", + "skewness", + "kurtosis", +] @njit(cache=True) @@ -92,3 +101,41 @@ 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) + + +@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) + 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 diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py index 1e799a4d73e..159fca60692 100644 --- a/src/ctapipe/image/tests/test_cleaning.py +++ b/src/ctapipe/image/tests/test_cleaning.py @@ -457,3 +457,104 @@ 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, 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, n_pixels=3) + assert np.count_nonzero(~mask) == 0 + + +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) + + 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. + # These two islands does not overlap. + charge[neighbors_1] = 29 + charge[core_pixel_1] = 50 + charge[neighbors_2] = 29 + charge[core_pixel_2] = 50 + + args = { + "picture_thresh_min": 45, + "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, + "pedestal_factor": 2, + "pedestal_std": None, + } + + # return normal tailcuts cleaning result if every other step is None/False: + 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) + + # 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_std"] = pedestal_std + + 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 + + # 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.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 + # peak_time to 10 + peak_time[neighbors_1] = 10 + args["pedestal_std"] = None + args["time_limit"] = 3 + + 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 + # all boundaries should be deselected + args["time_limit"] = None + args["bright_cleaning_threshold"] = 40 + + 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 + # core_pixel_1 with its neighbors + charge[neighbors_2] = 0 + args["bright_cleaning_threshold"] = None + args["largest_island_only"] = True + + 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 1b8ef34274b..e0f34829a0e 100644 --- a/src/ctapipe/image/tests/test_image_cleaner_component.py +++ b/src/ctapipe/image/tests/test_image_cleaner_component.py @@ -16,6 +16,12 @@ def test_image_cleaner(method, prod5_mst_nectarcam, reference_location): "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, }, + "NSBImageCleaner": { + "boundary_threshold_pe": 5.0, + "picture_threshold_pe": 10.0, + "time_limit": None, + "bright_cleaning_threshold": None, + }, "MARSImageCleaner": { "boundary_threshold_pe": 5.0, "picture_threshold_pe": 10.0, diff --git a/src/ctapipe/image/tests/test_statistics.py b/src/ctapipe/image/tests/test_statistics.py index b9773edc9a7..006c9c0dbce 100644 --- a/src/ctapipe/image/tests/test_statistics.py +++ b/src/ctapipe/image/tests/test_statistics.py @@ -60,3 +60,25 @@ def test_return_type(): stats = descriptive_statistics(data, container_class=PeakTimeStatisticsContainer) assert isinstance(stats, PeakTimeStatisticsContainer) + + +def test_n_largest(): + from ctapipe.image.statistics import n_largest + + 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()