Repository: cta-observatory/ctapipe
Pull request #2541 (2024-04-23 08:57:12)

Implement lstchains image cleaning
Changing the `ImageCleaner`s API in #2511 now allows to make use of monitoring information in the cleaning process.

This enables implementing lstchains image cleaning procedure from its [r0_to_dl1 script](https://github.com/cta-observatory/cta-lstchain/blob/main/lstchain/reco/r0_to_dl1.py#L219) and [dl1ab script](https://github.com/cta-observatory/cta-lstchain/blob/main/lstchain/scripts/lstchain_dl1ab.py) including the 'pedestal cleaning' method.



ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 98.90% Coverage (92.70% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.70% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 98.90% Coverage (92.70% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.70% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 98.00% Coverage (94.10% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 91.90% Coverage (94.10% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 91.90% Coverage (94.10% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 92.70% Coverage (94.10% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 92.70% Coverage (94.20% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2541)



Patch:
diff --git a/docs/changes/2541.feature.rst b/docs/changes/2541.feature.rst
new file mode 100644
index 00000000..30a27dbe
--- /dev/null
+++ b/docs/changes/2541.feature.rst
@@ -0,0 +1,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 6dc8ec54..e96b2bef 100644
--- a/src/ctapipe/image/cleaning.py
+++ b/src/ctapipe/image/cleaning.py
@@ -28,15 +28,16 @@ from abc import abstractmethod
 
 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 @@ class ImageCleaner(TelescopeComponent):
 
         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 @@ class ImageCleaner(TelescopeComponent):
 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 @@ class TailcutsImageCleaner(ImageCleaner):
         )
 
 
+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 2ba9281e..ce772744 100644
--- a/src/ctapipe/image/extractor.py
+++ b/src/ctapipe/image/extractor.py
@@ -47,6 +47,7 @@ from .cleaning import tailcuts_clean
 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 @@ class GlobalPeakWindowSum(ImageExtractor):
             ).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 ee8ddbd8..bd38f02a 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


Test patch:
diff --git a/src/ctapipe/image/tests/test_cleaning.py b/src/ctapipe/image/tests/test_cleaning.py
index 1e799a4d..159fca60 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 1b8ef342..e0f34829 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 b9773edc..006c9c0d 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()