Repository: cta-observatory/ctapipe
Pull request #2023 (2022-07-15 15:18:36)

Add abstract InvalidPixelHandler and NeighborAverage implementation, apply in CameraCalibrator and TwoPassWindowSum. Pass invalid pixels to extractors.
* Fill pixel_status in SimTelEventSource
* Ignore broken pixels in GlobalPeakWindowSum and NeighborPeakWindowSum extractors
* Interpolate images and peak_time of broken pixels in `CameraCalibrator` using a simple average of their neighbors

maxnoe: Before:

![before](https://user-images.githubusercontent.com/5488440/179255788-66854d79-51b6-49f4-afd4-fa54ec711f9c.png)

After:
![after](https://user-images.githubusercontent.com/5488440/179255841-4601b011-16df-4dfb-80d3-3ea3d27c11c1.png)

codecov[bot]: # [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
> Merging [#2023](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (f338b3e) into [master](https://codecov.io/gh/cta-observatory/ctapipe/commit/93527bf24a2107cf5c3613694f2ed01c0f657715?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (93527bf) will **increase** coverage by `0.03%`.
> The diff coverage is `98.10%`.

> :exclamation: Current head f338b3e differs from pull request most recent head 02ab60d. Consider uploading reports for the commit 02ab60d to get more accurate results

```diff
@@            Coverage Diff             @@
##           master    #2023      +/-   ##
==========================================
+ Coverage   92.33%   92.36%   +0.03%     
==========================================
  Files         193      195       +2     
  Lines       15708    15804      +96     
==========================================
+ Hits        14504    14598      +94     
- Misses       1204     1206       +2     
```


| [Impacted Files](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [ctapipe/image/extractor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9pbWFnZS9leHRyYWN0b3IucHk=) | `82.75% <85.00%> (-0.64%)` | :arrow_down: |
| [ctapipe/calib/camera/calibrator.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvY2FsaWJyYXRvci5weQ==) | `95.41% <100.00%> (+0.51%)` | :arrow_up: |
| [ctapipe/calib/camera/flatfield.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvZmxhdGZpZWxkLnB5) | `96.03% <100.00%> (+0.07%)` | :arrow_up: |
| [ctapipe/calib/camera/pedestals.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvcGVkZXN0YWxzLnB5) | `95.95% <100.00%> (+0.08%)` | :arrow_up: |
| [ctapipe/calib/camera/pixel\_interpolator.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvcGl4ZWxfaW50ZXJwb2xhdG9yLnB5) | `100.00% <100.00%> (ø)` | |
| [ctapipe/calib/camera/tests/test\_flatfield.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvdGVzdHMvdGVzdF9mbGF0ZmllbGQucHk=) | `92.85% <100.00%> (+0.17%)` | :arrow_up: |
| [ctapipe/calib/camera/tests/test\_pedestals.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvdGVzdHMvdGVzdF9wZWRlc3RhbHMucHk=) | `100.00% <100.00%> (ø)` | |
| [...pipe/calib/camera/tests/test\_pixel\_interpolator.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jYWxpYi9jYW1lcmEvdGVzdHMvdGVzdF9waXhlbF9pbnRlcnBvbGF0b3IucHk=) | `100.00% <100.00%> (ø)` | |
| [ctapipe/image/reducer.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9pbWFnZS9yZWR1Y2VyLnB5) | `92.30% <100.00%> (+0.15%)` | :arrow_up: |
| [ctapipe/image/tests/test\_extractor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9pbWFnZS90ZXN0cy90ZXN0X2V4dHJhY3Rvci5weQ==) | `100.00% <100.00%> (ø)` | |
| ... and [2 more](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Last update [93527bf...02ab60d](https://codecov.io/gh/cta-observatory/ctapipe/pull/2023?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).

maxnoe: > Those are normally computed after charge extraction, but before parameterization.

Those are computed from interleaved pedestals, and the information about current interleaved pedestal monitoring will be available in the current physics event exactly at this position (pedestal_failing_pixels), right?


kosack: > Ignore broken pixels in GlobalPeakWindowSum and NeighborPeakWindowSum extractors

Looks like also TwoPassWindowSum as well

kosack: > Those are computed from interleaved pedestals, and the information about current interleaved pedestal monitoring will be available in the current physics event exactly at this position (pedestal_failing_pixels), right?

Yes, I would imagine that was the intent there. Maybe the name is not great (it sounds more like pixels where a pedestal could not be computed, rather than where the pedestal variance is too high or too low).
maxnoe: > Looks like also TwoPassWindowSum as well

Not exactly. The `broken_pixels` argument had to be added to all extractors, because they share a common API. But only the two extractors mentioned make actual use of this information.

kosack: > Not exactly. The broken_pixels argument had to be added to all extractors, because they share a common API. But only the two extractors mentioned make actual use of this information.

Then shouldn't TwoPass also use it? In computing the time gradient in the first pass, broken pixels should be ignored.
maxnoe: > Then shouldn't TwoPass also use it? In computing the time gradient in the first pass, broken pixels should be ignored.

Ignored? Or should I add the interpolation?
maxnoe: > , especially when there are large regions of broken pixels like a full module or more where it could do more harm than good

Pixels where no valid neighbors are found are set to 0
maxnoe: > which to my mind means that this should be a separate step in the CameraCalibrator

It is, the interpolation happens in `CameraCalibrator`. The extractors combining multiple pixels however need to be able to ignore broken pixels.
maxnoe: > Perhaps BrokenPixelInpainter (since interpolation is only one way to do it). That way we can also include other types of broken pixel in the inpainting.


I started with a Component, then removed it again... The idea here was to start simple and adding another component that needs to have the API `__call__(event)` and does a loop over the telescope events was not something I found attractive.

It all should get better with the `ArrayEventContainer -> TelescopeEventContainer -> [DL0, DL1, DL2,MON]` structure, since then these components concerned with telescope calibration don't have a loop inside and can just take the single telescope event.
maxnoe: > I think it would only be necessary to interpolate the second pass. The first pass just determines the time gradient used for charge extraction - there just skipping broken pixels should be fine.

That's actually not how it works.

The first pass is just a simple integration for each pixel with some hard coded behavior to make it work exactly like MARS. Which makes it a bit ugly. But since no global or neighbor information is used here, nothing to worry about for the broken pixels.

The second pass starts with the parametrization of the first pass result. It applies a cleaning, thus using neighbor information.

So my strategy would be to add the interpolation in the second pass before the parametrization.

Like this, if only the first pass is used, the interpolation happens in `CameraCalibrator` for all pixels. If both passes are used, the cleaning and parametrization actually handle the broken pixels but the broken pixels are extracted again. Then the final result is interpolated in `CameraCalibrator`.

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3995f024..2b7695b0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -14,4 +14,4 @@ repos:
     rev: 4.0.1
     hooks:
       - id: flake8
-        args: [--max-line-length=88, --extend-ignore=E203]
+        args: [--max-line-length=88, "--extend-ignore=E203,E712"]
diff --git a/ctapipe/calib/camera/__init__.py b/ctapipe/calib/camera/__init__.py
index 861d256f..d576bafd 100644
--- a/ctapipe/calib/camera/__init__.py
+++ b/ctapipe/calib/camera/__init__.py
@@ -6,4 +6,7 @@ Camera calibration module.
 from .calibrator import CameraCalibrator
 from .gainselection import GainSelector
 
-__all__ = ["CameraCalibrator", "GainSelector"]
+__all__ = [
+    "CameraCalibrator",
+    "GainSelector",
+]
diff --git a/ctapipe/calib/camera/calibrator.py b/ctapipe/calib/camera/calibrator.py
index 505fe32d..0c722668 100644
--- a/ctapipe/calib/camera/calibrator.py
+++ b/ctapipe/calib/camera/calibrator.py
@@ -7,6 +7,8 @@ import warnings
 
 import astropy.units as u
 import numpy as np
+from numba import float32, float64, guvectorize, int64
+
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
@@ -15,12 +17,27 @@ from ctapipe.core.traits import (
     create_class_enum_trait,
 )
 from ctapipe.image.extractor import ImageExtractor
+from ctapipe.image.invalid_pixels import InvalidPixelHandler
 from ctapipe.image.reducer import DataVolumeReducer
-from numba import float32, float64, guvectorize, int64
 
 __all__ = ["CameraCalibrator"]
 
 
+def _get_invalid_pixels(n_pixels, pixel_status, selected_gain_channel):
+    broken_pixels = np.zeros(n_pixels, dtype=bool)
+    index = np.arange(n_pixels)
+    masks = (
+        pixel_status.hardware_failing_pixels,
+        pixel_status.pedestal_failing_pixels,
+        pixel_status.flatfield_failing_pixels,
+    )
+    for mask in masks:
+        if mask is not None:
+            broken_pixels |= mask[selected_gain_channel, index]
+
+    return broken_pixels
+
+
 class CameraCalibrator(TelescopeComponent):
     """
     Calibrator to handle the full camera calibration chain, in order to fill
@@ -48,6 +65,13 @@ class CameraCalibrator(TelescopeComponent):
         help="Name of the ImageExtractor subclass to be used.",
     ).tag(config=True)
 
+    invalid_pixel_handler_type = create_class_enum_trait(
+        InvalidPixelHandler,
+        default_value="NeighborAverage",
+        help="Name of the InvalidPixelHandler to use",
+        allow_none=True,
+    ).tag(config=True)
+
     apply_waveform_time_shift = BoolTelescopeParameter(
         default_value=False,
         help=(
@@ -123,6 +147,14 @@ class CameraCalibrator(TelescopeComponent):
         else:
             self.data_volume_reducer = data_volume_reducer
 
+        self.invalid_pixel_handler = None
+        if self.invalid_pixel_handler_type is not None:
+            self.invalid_pixel_handler = InvalidPixelHandler.from_name(
+                self.invalid_pixel_handler_type,
+                subarray=self.subarray,
+                parent=self,
+            )
+
     def _check_r1_empty(self, waveforms):
         if waveforms is None:
             if not self._r1_empty_warn:
@@ -162,19 +194,24 @@ class CameraCalibrator(TelescopeComponent):
         event.dl0.tel[telid].waveform = waveforms_copy
         event.dl0.tel[telid].selected_gain_channel = selected_gain_channel
 
-    def _calibrate_dl1(self, event, telid):
-        waveforms = event.dl0.tel[telid].waveform
-        selected_gain_channel = event.dl0.tel[telid].selected_gain_channel
-        dl1_calib = event.calibration.tel[telid].dl1
-
+    def _calibrate_dl1(self, event, tel_id):
+        waveforms = event.dl0.tel[tel_id].waveform
         if self._check_dl0_empty(waveforms):
             return
 
-        selected_gain_channel = event.r1.tel[telid].selected_gain_channel
-        time_shift = event.calibration.tel[telid].dl1.time_shift
-        readout = self.subarray.tel[telid].camera.readout
         n_pixels, n_samples = waveforms.shape
 
+        selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel
+        broken_pixels = _get_invalid_pixels(
+            n_pixels,
+            event.mon.tel[tel_id].pixel_status,
+            selected_gain_channel,
+        )
+
+        dl1_calib = event.calibration.tel[tel_id].dl1
+        time_shift = event.calibration.tel[tel_id].dl1.time_shift
+        readout = self.subarray.tel[tel_id].camera.readout
+
         # subtract any remaining pedestal before extraction
         if dl1_calib.pedestal_offset is not None:
             # this copies intentionally, we don't want to modify the dl0 data
@@ -197,7 +234,7 @@ class CameraCalibrator(TelescopeComponent):
 
             # shift waveforms if time_shift calibration is available
             if time_shift is not None:
-                if self.apply_waveform_time_shift.tel[telid]:
+                if self.apply_waveform_time_shift.tel[tel_id]:
                     sampling_rate = readout.sampling_rate.to_value(u.GHz)
                     time_shift_samples = time_shift * sampling_rate
                     waveforms, remaining_shift = shift_waveforms(
@@ -207,20 +244,32 @@ class CameraCalibrator(TelescopeComponent):
                 else:
                     remaining_shift = time_shift
 
-            extractor = self.image_extractors[self.image_extractor_type.tel[telid]]
+            extractor = self.image_extractors[self.image_extractor_type.tel[tel_id]]
             dl1 = extractor(
-                waveforms, telid=telid, selected_gain_channel=selected_gain_channel
+                waveforms,
+                telid=tel_id,
+                selected_gain_channel=selected_gain_channel,
+                broken_pixels=broken_pixels,
             )
 
             # correct non-integer remainder of the shift if given
-            if self.apply_peak_time_shift.tel[telid] and time_shift is not None:
+            if self.apply_peak_time_shift.tel[tel_id] and time_shift is not None:
                 dl1.peak_time -= remaining_shift
 
         # Calibrate extracted charge
         dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
 
+        # handle invalid pixels
+        if self.invalid_pixel_handler is not None:
+            dl1.image, dl1.peak_time = self.invalid_pixel_handler(
+                tel_id,
+                dl1.image,
+                dl1.peak_time,
+                broken_pixels,
+            )
+
         # store the results in the event structure
-        event.dl1.tel[telid] = dl1
+        event.dl1.tel[tel_id] = dl1
 
     def __call__(self, event):
         """
diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py
index 1e3263c6..0af81699 100644
--- a/ctapipe/calib/camera/flatfield.py
+++ b/ctapipe/calib/camera/flatfield.py
@@ -6,11 +6,14 @@ from abc import abstractmethod
 
 import numpy as np
 from astropy import units as u
+
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import Component
 from ctapipe.core.traits import Int, List, Unicode
 from ctapipe.image.extractor import ImageExtractor
 
+from .calibrator import _get_invalid_pixels
+
 __all__ = ["FlatFieldCalculator", "FlasherFlatFieldCalculator"]
 
 
@@ -181,10 +184,16 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
         waveforms = event.r1.tel[self.tel_id].waveform
         selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel
-
+        broken_pixels = _get_invalid_pixels(
+            n_pixels=waveforms.shape[-2],
+            pixel_status=event.mon.tel[self.tel_id].pixel_status,
+            selected_gain_channel=selected_gain_channel,
+        )
         # Extract charge and time
         if self.extractor:
-            return self.extractor(waveforms, self.tel_id, selected_gain_channel)
+            return self.extractor(
+                waveforms, self.tel_id, selected_gain_channel, broken_pixels
+            )
         else:
             return DL1CameraContainer(image=0, peak_pos=0, is_valid=False)
 
diff --git a/ctapipe/calib/camera/pedestals.py b/ctapipe/calib/camera/pedestals.py
index ed8c527d..0942d2c4 100644
--- a/ctapipe/calib/camera/pedestals.py
+++ b/ctapipe/calib/camera/pedestals.py
@@ -3,14 +3,16 @@ Factory for the estimation of the flat field coefficients
 """
 
 from abc import abstractmethod
+
 import numpy as np
 from astropy import units as u
+
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import Component
-
-
+from ctapipe.core.traits import Int, List, Unicode
 from ctapipe.image.extractor import ImageExtractor
-from ctapipe.core.traits import Int, Unicode, List
+
+from .calibrator import _get_invalid_pixels
 
 __all__ = ["calc_pedestals_from_traces", "PedestalCalculator", "PedestalIntegrator"]
 
@@ -210,10 +212,17 @@ class PedestalIntegrator(PedestalCalculator):
         """
         waveforms = event.r1.tel[self.tel_id].waveform
         selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel
+        broken_pixels = _get_invalid_pixels(
+            n_pixels=waveforms.shape[-2],
+            pixel_status=event.mon.tel[self.tel_id].pixel_status,
+            selected_gain_channel=selected_gain_channel,
+        )
 
         # Extract charge and time
         if self.extractor:
-            return self.extractor(waveforms, self.tel_id, selected_gain_channel)
+            return self.extractor(
+                waveforms, self.tel_id, selected_gain_channel, broken_pixels
+            )
         else:
             return DL1CameraContainer(image=0, peak_pos=0, is_valid=False)
 
diff --git a/ctapipe/core/traits.py b/ctapipe/core/traits.py
index c336e44c..083ca6d6 100644
--- a/ctapipe/core/traits.py
+++ b/ctapipe/core/traits.py
@@ -199,7 +199,7 @@ class Path(TraitType):
         return value
 
 
-def create_class_enum_trait(base_class, default_value, help=None):
+def create_class_enum_trait(base_class, default_value, help=None, allow_none=False):
     """create a configurable CaselessStrEnum traitlet from baseclass
 
     the enumeration should contain all names of non_abstract_children()
@@ -217,7 +217,10 @@ def create_class_enum_trait(base_class, default_value, help=None):
         raise ValueError(f"{default_value} is not in choices: {choices}")
 
     return CaselessStrEnum(
-        choices, default_value=default_value, allow_none=False, help=help
+        choices,
+        default_value=default_value,
+        help=help,
+        allow_none=allow_none,
     ).tag(config=True)
 
 
diff --git a/ctapipe/image/__init__.py b/ctapipe/image/__init__.py
index 5264cb0e..3f8f32fc 100644
--- a/ctapipe/image/__init__.py
+++ b/ctapipe/image/__init__.py
@@ -1,67 +1,66 @@
-from .hillas import (
-    hillas_parameters,
-    HillasParameterizationError,
-    camera_to_shower_coordinates,
-)
-from .timing import timing_parameters
-from .leakage import leakage_parameters
-from .concentration import concentration_parameters
-from .statistics import descriptive_statistics
-from .morphology import (
-    number_of_islands,
-    number_of_island_sizes,
-    morphology_parameters,
-    largest_island,
-    brightest_island,
-)
-
 from .cleaning import (
-    tailcuts_clean,
-    dilate,
-    mars_cleaning_1st_pass,
-    fact_image_cleaning,
-    apply_time_delta_cleaning,
     ImageCleaner,
     TailcutsImageCleaner,
+    apply_time_delta_cleaning,
+    dilate,
+    fact_image_cleaning,
+    mars_cleaning_1st_pass,
+    tailcuts_clean,
 )
-from .pixel_likelihood import (
-    neg_log_likelihood_approx,
-    neg_log_likelihood_numeric,
-    neg_log_likelihood,
-    mean_poisson_likelihood_gaussian,
-    mean_poisson_likelihood_full,
-    PixelLikelihoodError,
-    chi_squared,
-)
+from .concentration import concentration_parameters
 from .extractor import (
-    ImageExtractor,
-    FullWaveformSum,
+    BaselineSubtractedNeighborPeakWindowSum,
     FixedWindowSum,
+    FullWaveformSum,
     GlobalPeakWindowSum,
+    ImageExtractor,
     LocalPeakWindowSum,
-    SlidingWindowMaxSum,
     NeighborPeakWindowSum,
-    BaselineSubtractedNeighborPeakWindowSum,
+    SlidingWindowMaxSum,
     TwoPassWindowSum,
     extract_around_peak,
     extract_sliding_window,
+    integration_correction,
     neighbor_average_maximum,
     subtract_baseline,
-    integration_correction,
 )
-from .reducer import DataVolumeReducer, NullDataVolumeReducer, TailCutsDataVolumeReducer
+from .hillas import (
+    HillasParameterizationError,
+    camera_to_shower_coordinates,
+    hillas_parameters,
+)
+from .image_processor import ImageProcessor
+from .invalid_pixels import InvalidPixelHandler, NeighborAverage
+from .leakage import leakage_parameters
+from .modifications import ImageModifier
+from .morphology import (
+    brightest_island,
+    largest_island,
+    morphology_parameters,
+    number_of_island_sizes,
+    number_of_islands,
+)
 from .muon import (
     MuonIntensityFitter,
     MuonRingFitter,
+    intensity_ratio_inside_ring,
     kundu_chaudhuri_circle_fit,
     mean_squared_error,
-    intensity_ratio_inside_ring,
     ring_completeness,
     ring_containment,
 )
-from .modifications import ImageModifier
-from .image_processor import ImageProcessor
-
+from .pixel_likelihood import (
+    PixelLikelihoodError,
+    chi_squared,
+    mean_poisson_likelihood_full,
+    mean_poisson_likelihood_gaussian,
+    neg_log_likelihood,
+    neg_log_likelihood_approx,
+    neg_log_likelihood_numeric,
+)
+from .reducer import DataVolumeReducer, NullDataVolumeReducer, TailCutsDataVolumeReducer
+from .statistics import descriptive_statistics
+from .timing import timing_parameters
 
 __all__ = [
     "ImageModifier",
@@ -116,4 +115,6 @@ __all__ = [
     "DataVolumeReducer",
     "NullDataVolumeReducer",
     "TailCutsDataVolumeReducer",
+    "InvalidPixelHandler",
+    "NeighborAverage",
 ]
diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index a36c35b9..285b35dc 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -25,19 +25,23 @@ from functools import lru_cache
 from typing import Tuple
 
 import numpy as np
+from numba import float32, float64, guvectorize, int64, njit, prange
+from scipy.ndimage import convolve1d
+from traitlets import Bool, Int
+
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
     BoolTelescopeParameter,
     FloatTelescopeParameter,
     IntTelescopeParameter,
+    create_class_enum_trait,
 )
-from numba import float32, float64, guvectorize, int64, njit, prange
-from scipy.ndimage import convolve1d
-from traitlets import Bool, Int
 
-from . import brightest_island, number_of_islands, tailcuts_clean
+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 .timing import timing_parameters
 
 
@@ -195,7 +199,9 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time)
 
 
 @njit(cache=True)
-def neighbor_average_maximum(waveforms, neighbors_indices, neighbors_indptr, lwt):
+def neighbor_average_maximum(
+    waveforms, neighbors_indices, neighbors_indptr, local_weight, broken_pixels
+):
     """
     Obtain the average waveform built from the neighbors of each pixel
 
@@ -210,9 +216,12 @@ def neighbor_average_maximum(waveforms, neighbors_indices, neighbors_indptr, lwt
     neighbors_indptr : ndarray
         indptr of a scipy csr sparse matrix of neighbors, i.e.
         ``ctapipe.instrument.CameraGeometry.neighbor_matrix_sparse.indptr``.
-    lwt: int
+    local_weight : int
         Weight of the local pixel (0: peak from neighbors only,
         1: local pixel counts as much as any neighbor)
+    broken_pixels : ndarray
+        Mask of broken pixels. Broken pixels are ignored in the sum over the
+        neighbors.
 
     Returns
     -------
@@ -226,15 +235,18 @@ def neighbor_average_maximum(waveforms, neighbors_indices, neighbors_indptr, lwt
     indptr = neighbors_indptr
     indices = neighbors_indices
 
-    # initialize to waveforms weighted with lwt
+    # initialize to waveforms weighted with local_weight
     # so the value of the pixel itself is already taken into account
     peak_pos = np.empty(n_pixels, dtype=np.int64)
 
     for pixel in prange(n_pixels):
-        average = waveforms[pixel] * lwt
+        average = waveforms[pixel] * local_weight
         neighbors = indices[indptr[pixel] : indptr[pixel + 1]]
 
         for neighbor in neighbors:
+            if broken_pixels[neighbor]:
+                continue
+
             average += waveforms[neighbor]
 
         peak_pos[pixel] = np.argmax(average)
@@ -368,7 +380,9 @@ class ImageExtractor(TelescopeComponent):
         }
 
     @abstractmethod
-    def __call__(self, waveforms, telid, selected_gain_channel) -> DL1CameraContainer:
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         """
         Call the relevant functions to fully extract the charge and time
         for the particular extractor.
@@ -398,7 +412,9 @@ class FullWaveformSum(ImageExtractor):
     Extractor that sums the entire waveform.
     """
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         charge, peak_time = extract_around_peak(
             waveforms, 0, waveforms.shape[-1], 0, self.sampling_rate_ghz[telid]
         )
@@ -455,7 +471,9 @@ class FixedWindowSum(ImageExtractor):
             self.window_shift.tel[telid],
         )
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         charge, peak_time = extract_around_peak(
             waveforms,
             self.peak_index.tel[telid],
@@ -532,13 +550,17 @@ class GlobalPeakWindowSum(ImageExtractor):
             self.window_shift.tel[telid],
         )
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         if self.pixel_fraction.tel[telid] == 1.0:
             # average over pixels then argmax over samples
-            peak_index = waveforms.mean(axis=-2).argmax()
+            peak_index = waveforms[~broken_pixels].mean(axis=-2).argmax()
         else:
             n_pixels = int(self.pixel_fraction.tel[telid] * waveforms.shape[-2])
-            brightest = np.argsort(waveforms.max(axis=-1))[..., -n_pixels:]
+            brightest = np.argsort(waveforms.max(axis=-1))[~broken_pixels][
+                ..., -n_pixels:
+            ]
 
             # average over brightest pixels then argmax over samples
             peak_index = waveforms[brightest].mean(axis=-2).argmax()
@@ -604,7 +626,9 @@ class LocalPeakWindowSum(ImageExtractor):
             self.window_shift.tel[telid],
         )
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         peak_index = waveforms.argmax(axis=-1).astype(np.int64)
         charge, peak_time = extract_around_peak(
             waveforms,
@@ -684,7 +708,9 @@ class SlidingWindowMaxSum(ImageExtractor):
 
         return correction
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         charge, peak_time = extract_sliding_window(
             waveforms, self.window_width.tel[telid], self.sampling_rate_ghz[telid]
         )
@@ -709,7 +735,7 @@ class NeighborPeakWindowSum(ImageExtractor):
         "from the peak_index (peak_index - shift)",
     ).tag(config=True)
 
-    lwt = IntTelescopeParameter(
+    local_weight = IntTelescopeParameter(
         default_value=0,
         help="Weight of the local pixel (0: peak from neighbors only, "
         "1: local pixel counts as much as any neighbor)",
@@ -748,13 +774,16 @@ class NeighborPeakWindowSum(ImageExtractor):
             self.window_shift.tel[telid],
         )
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         neighbors = self.subarray.tel[telid].camera.geometry.neighbor_matrix_sparse
         peak_index = neighbor_average_maximum(
             waveforms,
             neighbors_indices=neighbors.indices,
             neighbors_indptr=neighbors.indptr,
-            lwt=self.lwt.tel[telid],
+            local_weight=self.local_weight.tel[telid],
+            broken_pixels=broken_pixels,
         )
         charge, peak_time = extract_around_peak(
             waveforms,
@@ -779,11 +808,15 @@ class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum):
     )
     baseline_end = Int(10, help="End sample for baseline estimation").tag(config=True)
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         baseline_corrected = subtract_baseline(
             waveforms, self.baseline_start, self.baseline_end
         )
-        return super().__call__(baseline_corrected, telid, selected_gain_channel)
+        return super().__call__(
+            baseline_corrected, telid, selected_gain_channel, broken_pixels
+        )
 
 
 class TwoPassWindowSum(ImageExtractor):
@@ -844,6 +877,23 @@ class TwoPassWindowSum(ImageExtractor):
         default_value=True, help="Apply the integration window correction"
     ).tag(config=True)
 
+    invalid_pixel_handler_type = create_class_enum_trait(
+        InvalidPixelHandler,
+        default_value="NeighborAverage",
+        help="Name of the InvalidPixelHandler to apply in the first pass.",
+        allow_none=True,
+    ).tag(config=True)
+
+    def __init__(self, subarray, **kwargs):
+        super().__init__(subarray=subarray, **kwargs)
+        self.invalid_pixel_handler = None
+        if self.invalid_pixel_handler_type is not None:
+            self.invalid_pixel_handler = InvalidPixelHandler.from_name(
+                self.invalid_pixel_handler_type,
+                self.subarray,
+                parent=self,
+            )
+
     @lru_cache(maxsize=4096)
     def _calculate_correction(self, telid, width, shift):
         """Obtain the correction for the integration window specified for each
@@ -959,6 +1009,7 @@ class TwoPassWindowSum(ImageExtractor):
         charge_1stpass_uncorrected,
         pulse_time_1stpass,
         correction,
+        broken_pixels,
     ) -> Tuple[np.ndarray, np.ndarray]:
         """
         Follow steps from 2 to 7.
@@ -997,16 +1048,23 @@ class TwoPassWindowSum(ImageExtractor):
             True=second-pass succeeded, False=second-pass failed, first pass used
         """
         # STEP 2
-
         # Apply correction to 1st pass charges
         charge_1stpass = charge_1stpass_uncorrected * correction[selected_gain_channel]
 
+        camera_geometry = self.subarray.tel[telid].camera.geometry
+        if self.invalid_pixel_handler is not None:
+            charge_1stpass, pulse_time_1stpass = self.invalid_pixel_handler(
+                telid,
+                charge_1stpass,
+                pulse_time_1stpass,
+                broken_pixels,
+            )
+
         # Set thresholds for core-pixels depending on telescope
         core_th = self.core_threshold.tel[telid]
         # Boundary thresholds will be half of core thresholds.
 
         # Preliminary image cleaning with simple two-level tail-cut
-        camera_geometry = self.subarray.tel[telid].camera.geometry
         mask_clean = tailcuts_clean(
             camera_geometry,
             charge_1stpass,
@@ -1201,7 +1259,9 @@ class TwoPassWindowSum(ImageExtractor):
 
         return charge_2ndpass, pulse_time_2ndpass, True
 
-    def __call__(self, waveforms, telid, selected_gain_channel):
+    def __call__(
+        self, waveforms, telid, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         charge1, pulse_time1, correction1 = self._apply_first_pass(waveforms, telid)
 
         # FIXME: properly make sure that output is 32Bit instead of downcasting here
@@ -1213,7 +1273,13 @@ class TwoPassWindowSum(ImageExtractor):
             )
 
         charge2, pulse_time2, is_valid = self._apply_second_pass(
-            waveforms, telid, selected_gain_channel, charge1, pulse_time1, correction1
+            waveforms,
+            telid,
+            selected_gain_channel,
+            charge1,
+            pulse_time1,
+            correction1,
+            broken_pixels,
         )
         # FIXME: properly make sure that output is 32Bit instead of downcasting here
         return DL1CameraContainer(
diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py
index 94450dda..fe5211f0 100644
--- a/ctapipe/image/image_processor.py
+++ b/ctapipe/image/image_processor.py
@@ -1,30 +1,28 @@
 """
 High level image processing  (ImageProcessor Component)
 """
-from ctapipe.coordinates import TelescopeFrame
 import numpy as np
 
+from ctapipe.coordinates import TelescopeFrame
+
 from ..containers import (
     ArrayEventContainer,
-    IntensityStatisticsContainer,
     ImageParametersContainer,
-    TimingParametersContainer,
+    IntensityStatisticsContainer,
     PeakTimeStatisticsContainer,
+    TimingParametersContainer,
 )
 from ..core import QualityQuery, TelescopeComponent
 from ..core.traits import Bool, BoolTelescopeParameter, List, create_class_enum_trait
 from ..instrument import SubarrayDescription
-from . import (
-    ImageCleaner,
-    ImageModifier,
-    concentration_parameters,
-    descriptive_statistics,
-    hillas_parameters,
-    leakage_parameters,
-    morphology_parameters,
-    timing_parameters,
-)
-
+from .cleaning import ImageCleaner
+from .concentration import concentration_parameters
+from .hillas import hillas_parameters
+from .leakage import leakage_parameters
+from .modifications import ImageModifier
+from .morphology import morphology_parameters
+from .statistics import descriptive_statistics
+from .timing import timing_parameters
 
 # avoid use of base containers for unparameterized images
 DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer()
diff --git a/ctapipe/image/invalid_pixels.py b/ctapipe/image/invalid_pixels.py
new file mode 100644
index 00000000..85dd2cce
--- /dev/null
+++ b/ctapipe/image/invalid_pixels.py
@@ -0,0 +1,106 @@
+"""
+Methods to interpolate broken pixels
+"""
+from abc import ABCMeta, abstractmethod
+from typing import Tuple
+
+import numpy as np
+
+from ..core import TelescopeComponent
+
+__all__ = [
+    "InvalidPixelHandler",
+    "NeighborAverage",
+]
+
+
+class InvalidPixelHandler(TelescopeComponent, metaclass=ABCMeta):
+    """
+    An abtract base class for algorithms treating invalid pixel data in images
+    """
+
+    @abstractmethod
+    def __call__(
+        self, tel_id, image, peak_time, pixel_mask
+    ) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Handle invalid (broken, high noise, ...) pixels.
+
+        E.g. replace by average over neighbors, see `NeighborAverage`).
+
+        Parameters
+        ----------
+        tel_id : int
+            telescope id
+        image : np.ndarray
+            Array of pixel image values
+        peak_time : np.ndarray
+            Array of pixel peak_time values
+        pixel_mask : np.ndarray
+            Boolean mask of the pixels to be interpolated
+
+        Returns
+        -------
+        image : np.ndarray
+            Image with interpolated values
+        peak_time : np.ndarray
+            peak_time with interpolated values
+        """
+
+
+class NeighborAverage(InvalidPixelHandler):
+    def __call__(self, tel_id, image, peak_time, pixel_mask):
+        """Interpolate pixels in dl1 images and peak_times
+
+        Pixels to be interpolated are replaced by the average value their neighbors.
+
+        Pixels where no valid neighbors are available are filled with zeros.
+
+        Parameters
+        ----------
+        tel_id : int
+            telescope id
+        image : np.ndarray
+            Array of pixel image values
+        peak_time : np.ndarray
+            Array of pixel peak_time values
+        pixel_mask : np.ndarray
+            Boolean mask of the pixels to be interpolated
+
+        Returns
+        -------
+        image : np.ndarray
+            Image with interpolated values
+        peak_time : np.ndarray
+            peak_time with interpolated values
+        """
+        geometry = self.subarray.tel[tel_id].camera.geometry
+
+        n_interpolated = np.count_nonzero(pixel_mask)
+        if n_interpolated == 0:
+            return image, peak_time
+
+        # exclude to-be-interpolated pixels from neighbors
+        neighbors = geometry.neighbor_matrix[pixel_mask] & ~pixel_mask
+
+        index, neighbor = np.nonzero(neighbors)
+        image_sum = np.zeros(n_interpolated, dtype=image.dtype)
+        count = np.zeros(n_interpolated, dtype=int)
+        peak_time_sum = np.zeros(n_interpolated, dtype=peak_time.dtype)
+
+        # calculate average of image and peak_time
+        np.add.at(count, index, 1)
+        np.add.at(image_sum, index, image[neighbor])
+        np.add.at(peak_time_sum, index, peak_time[neighbor])
+
+        valid = count > 0
+        np.divide(image_sum, count, out=image_sum, where=valid)
+        np.divide(peak_time_sum, count, out=peak_time_sum, where=valid)
+
+        peak_time = peak_time.copy()
+        peak_time[pixel_mask] = peak_time_sum
+
+        image = image.copy()
+        image[pixel_mask] = image_sum
+
+        return image, peak_time
diff --git a/ctapipe/image/reducer.py b/ctapipe/image/reducer.py
index 6386bb05..03b96020 100644
--- a/ctapipe/image/reducer.py
+++ b/ctapipe/image/reducer.py
@@ -2,18 +2,20 @@
 Algorithms for the data volume reduction.
 """
 from abc import abstractmethod
+
 import numpy as np
+
 from ctapipe.containers import DL1CameraContainer
-from ctapipe.image import TailcutsImageCleaner
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
-    IntTelescopeParameter,
     BoolTelescopeParameter,
+    IntTelescopeParameter,
     TelescopeParameter,
     create_class_enum_trait,
 )
-from ctapipe.image.extractor import ImageExtractor
+from ctapipe.image import TailcutsImageCleaner
 from ctapipe.image.cleaning import dilate
+from ctapipe.image.extractor import ImageExtractor
 
 __all__ = ["DataVolumeReducer", "NullDataVolumeReducer", "TailCutsDataVolumeReducer"]
 
@@ -181,8 +183,13 @@ class TailCutsDataVolumeReducer(DataVolumeReducer):
         camera_geom = self.subarray.tel[telid].camera.geometry
         # Pulse-integrate waveforms
         extractor = self.image_extractors[self.image_extractor_type.tel[telid]]
+        # do not treat broken pixels in data volume reduction
+        broken_pixels = np.zeros(camera_geom.n_pixels, dtype=bool)
         dl1: DL1CameraContainer = extractor(
-            waveforms, telid=telid, selected_gain_channel=selected_gain_channel
+            waveforms,
+            telid=telid,
+            selected_gain_channel=selected_gain_channel,
+            broken_pixels=broken_pixels,
         )
 
         # 1) Step: TailcutCleaning at first
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 196bda33..0def68eb 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -15,6 +15,7 @@ from ..containers import (
     ArrayEventContainer,
     EventIndexContainer,
     EventType,
+    PixelStatusContainer,
     PointingContainer,
     R0CameraContainer,
     R1CameraContainer,
@@ -505,6 +506,7 @@ class SimTelEventSource(EventSource):
                 mon = data.mon.tel[tel_id]
                 mon.calibration.dc_to_pe = dc_to_pe
                 mon.calibration.pedestal_per_sample = pedestal
+                mon.pixel_status = self._get_pixels_status(tel_id)
 
                 r1_waveform, selected_gain_channel = apply_simtel_r1_calibration(
                     adc_samples,
@@ -528,6 +530,21 @@ class SimTelEventSource(EventSource):
 
             yield data
 
+    def _get_pixels_status(self, tel_id):
+        tel = self.file_.telescope_descriptions[tel_id]
+        n_pixels = tel["camera_organization"]["n_pixels"]
+        n_gains = tel["camera_organization"]["n_gains"]
+
+        disabled_ids = tel["disabled_pixels"]["HV_disabled"]
+        disabled_pixels = np.zeros((n_gains, n_pixels), dtype=bool)
+        disabled_pixels[:, disabled_ids] = True
+
+        return PixelStatusContainer(
+            hardware_failing_pixels=disabled_pixels,
+            pedestal_failing_pixels=disabled_pixels.copy(),
+            flatfield_failing_pixels=disabled_pixels.copy(),
+        )
+
     @staticmethod
     def _fill_event_pointing(tracking_position):
         azimuth_raw = tracking_position["azimuth_raw"]
diff --git a/docs/ctapipe_api/image/index.rst b/docs/ctapipe_api/image/index.rst
index 662231fb..0e02ec1b 100644
--- a/docs/ctapipe_api/image/index.rst
+++ b/docs/ctapipe_api/image/index.rst
@@ -28,6 +28,7 @@ This module contains the following sub-modules, but the most important functions
   extractor
   reducers
   modifications
+  invalid_pixels
 
 
 
diff --git a/docs/ctapipe_api/image/invalid_pixels.rst b/docs/ctapipe_api/image/invalid_pixels.rst
new file mode 100644
index 00000000..17e53173
--- /dev/null
+++ b/docs/ctapipe_api/image/invalid_pixels.rst
@@ -0,0 +1,12 @@
+.. _invalid_pixels:
+
+=======================
+Invalid Pixels Handling
+=======================
+
+
+
+API Reference
+-------------
+
+.. automodapi:: ctapipe.image.invalid_pixels
diff --git a/setup.cfg b/setup.cfg
index 0b1d6c94..2f26e090 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -39,7 +39,7 @@ exclude=
 max-line-length=90
 #ignore=W291,E303,W391,F403,F401,W503,W1202
 select = C,E,F,W,B,B950
-ignore = E501,W503,E203,W201,E731
+ignore = E501,W503,E203,W201,E712,E731
 per-file-ignores =
     */__init__.py: F401, F403
 


Test patch:
diff --git a/ctapipe/calib/camera/tests/test_calibrator.py b/ctapipe/calib/camera/tests/test_calibrator.py
index fae1b051..ba9346f9 100644
--- a/ctapipe/calib/camera/tests/test_calibrator.py
+++ b/ctapipe/calib/camera/tests/test_calibrator.py
@@ -1,22 +1,23 @@
 """
 Tests for CameraCalibrator and related functions
 """
+from copy import deepcopy
+
 import astropy.units as u
 import numpy as np
 import pytest
 from scipy.stats import norm
-from traitlets.config.configurable import Config
+from traitlets.config import Config
 
 from ctapipe.calib.camera.calibrator import CameraCalibrator
+from ctapipe.containers import ArrayEventContainer
 from ctapipe.image.extractor import (
-    NeighborPeakWindowSum,
-    LocalPeakWindowSum,
     FullWaveformSum,
     GlobalPeakWindowSum,
+    LocalPeakWindowSum,
+    NeighborPeakWindowSum,
 )
 from ctapipe.image.reducer import NullDataVolumeReducer, TailCutsDataVolumeReducer
-from copy import deepcopy
-from ctapipe.containers import ArrayEventContainer
 
 
 def test_camera_calibrator(example_event, example_subarray):
@@ -263,3 +264,48 @@ def test_shift_waveforms():
     assert shifted_waveforms[2, 12] == 1
     assert shifted_waveforms[3, 7] == 1
     assert shifted_waveforms[4, 14] == 1
+
+
+def test_invalid_pixels(example_event, example_subarray):
+
+    # switching off the corrections makes it easier to test for
+    # the exact value of 1.0
+    config = Config(
+        {
+            "CameraCalibrator": {
+                "image_extractor_type": "LocalPeakWindowSum",
+                "apply_peak_time_shift": False,
+                "LocalPeakWindowSum": {
+                    "apply_integration_correction": False,
+                },
+            }
+        }
+    )
+    # going to modify this
+    event = deepcopy(example_event)
+    tel_id = list(event.r0.tel)[0]
+    camera = example_subarray.tel[tel_id].camera
+    sampling_rate = camera.readout.sampling_rate.to_value(u.GHz)
+
+    event.mon.tel[tel_id].pixel_status.flatfield_failing_pixels[:, 0] = True
+    event.r1.tel[tel_id].waveform.fill(0.0)
+    event.r1.tel[tel_id].waveform[1:, 20] = 1.0
+    event.r1.tel[tel_id].waveform[0, 10] = 9999
+
+    calibrator = CameraCalibrator(
+        subarray=example_subarray,
+        config=config,
+    )
+    calibrator(event)
+    assert np.all(event.dl1.tel[tel_id].image == 1.0)
+    assert np.all(event.dl1.tel[tel_id].peak_time == 20.0 / sampling_rate)
+
+    # test we can set the invalid pixel handler to None
+    config.CameraCalibrator.invalid_pixel_handler_type = None
+    calibrator = CameraCalibrator(
+        subarray=example_subarray,
+        config=config,
+    )
+    calibrator(event)
+    assert event.dl1.tel[tel_id].image[0] == 9999
+    assert event.dl1.tel[tel_id].peak_time[0] == 10.0 / sampling_rate
diff --git a/ctapipe/calib/camera/tests/test_flatfield.py b/ctapipe/calib/camera/tests/test_flatfield.py
index f9387c25..70cbed30 100644
--- a/ctapipe/calib/camera/tests/test_flatfield.py
+++ b/ctapipe/calib/camera/tests/test_flatfield.py
@@ -3,7 +3,7 @@ from copy import deepcopy
 import astropy.units as u
 import numpy as np
 from astropy.time import Time
-from traitlets.config.loader import Config
+from traitlets.config import Config
 
 from ctapipe.calib.camera.flatfield import FlasherFlatFieldCalculator
 from ctapipe.containers import ArrayEventContainer
@@ -52,6 +52,7 @@ def test_flasherflatfieldcalculator(prod5_sst):
         (n_gain, n_pixels), dtype=bool
     )
     data.r1.tel[tel_id].waveform = np.zeros((n_gain, n_pixels, 40))
+    data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8)
 
     # flat-field signal put == delta function of height ff_level at sample 20
     data.r1.tel[tel_id].waveform[:, :, 20] = ff_level
diff --git a/ctapipe/calib/camera/tests/test_pedestals.py b/ctapipe/calib/camera/tests/test_pedestals.py
index 5cec5ec0..32deb714 100644
--- a/ctapipe/calib/camera/tests/test_pedestals.py
+++ b/ctapipe/calib/camera/tests/test_pedestals.py
@@ -15,7 +15,7 @@ from ctapipe.containers import ArrayEventContainer
 from ctapipe.instrument import SubarrayDescription
 
 
-def test_pedestal_calculator(prod5_sst):
+def test_pedestal_integrator(prod5_sst):
     """test of PedestalIntegrator"""
 
     tel_id = 0
@@ -48,6 +48,7 @@ def test_pedestal_calculator(prod5_sst):
         (n_gain, n_pixels), dtype=bool
     )
     data.r1.tel[tel_id].waveform = np.full((2, n_pixels, 40), ped_level)
+    data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8)
 
     while ped_calculator.num_events_seen < n_events:
         if ped_calculator.calculate_pedestals(data):
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index b93071b2..460db2bf 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -191,11 +191,13 @@ def test_extract_around_peak_charge_expected(toymodel):
 def test_neighbor_average_peakpos(toymodel):
     waveforms, subarray, telid, _, _, _ = toymodel
     neighbors = subarray.tel[telid].camera.geometry.neighbor_matrix_sparse
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
     peak_pos = neighbor_average_maximum(
         waveforms,
         neighbors_indices=neighbors.indices,
         neighbors_indptr=neighbors.indptr,
-        lwt=0,
+        local_weight=0,
+        broken_pixels=broken_pixels,
     )
 
     pixel = 0
@@ -204,17 +206,18 @@ def test_neighbor_average_peakpos(toymodel):
     expected_peak_pos = np.argmax(expected_average, axis=-1)
     assert (peak_pos[pixel] == expected_peak_pos).all()
 
-    lwt = 4
+    local_weight = 4
     peak_pos = neighbor_average_maximum(
         waveforms,
         neighbors_indices=neighbors.indices,
         neighbors_indptr=neighbors.indptr,
-        lwt=lwt,
+        local_weight=local_weight,
+        broken_pixels=broken_pixels,
     )
 
     pixel = 1
     _, nei_pixel = np.where(neighbors[pixel].A)
-    nei_pixel = np.concatenate([nei_pixel, [pixel] * lwt])
+    nei_pixel = np.concatenate([nei_pixel, [pixel] * local_weight])
     expected_average = waveforms[nei_pixel].sum(0) / len(nei_pixel)
     expected_peak_pos = np.argmax(expected_average, axis=-1)
     assert (peak_pos[pixel] == expected_peak_pos).all()
@@ -292,7 +295,8 @@ def test_integration_correction_outofbounds(subarray):
 def test_extractors(Extractor, toymodel):
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
     extractor = Extractor(subarray=subarray)
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
     assert_allclose(dl1.image, true_charge, rtol=0.1)
     assert_allclose(dl1.peak_time, true_time, rtol=0.1)
     assert dl1.is_valid == True
@@ -306,7 +310,8 @@ def test_integration_correction_off(Extractor, toymodel):
 
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
     extractor = Extractor(subarray=subarray, apply_integration_correction=False)
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
 
     assert dl1.is_valid == True
 
@@ -320,7 +325,8 @@ def test_integration_correction_off(Extractor, toymodel):
 def test_fixed_window_sum(toymodel):
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
     extractor = FixedWindowSum(subarray=subarray, peak_index=47)
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
     assert_allclose(dl1.image, true_charge, rtol=0.1)
     assert_allclose(dl1.peak_time, true_time, rtol=0.1)
     assert dl1.is_valid == True
@@ -329,18 +335,20 @@ def test_fixed_window_sum(toymodel):
 def test_sliding_window_max_sum(toymodel):
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
     extractor = SlidingWindowMaxSum(subarray=subarray)
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
     print(true_charge, dl1.image, true_time, dl1.peak_time)
     assert_allclose(dl1.image, true_charge, rtol=0.1)
     assert_allclose(dl1.peak_time, true_time, rtol=0.1)
     assert dl1.is_valid == True
 
 
-def test_neighbor_peak_window_sum_lwt(toymodel):
+def test_neighbor_peak_window_sum_local_weight(toymodel):
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
-    extractor = NeighborPeakWindowSum(subarray=subarray, lwt=4)
-    assert extractor.lwt.tel[telid] == 4
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+    extractor = NeighborPeakWindowSum(subarray=subarray, local_weight=4)
+    assert extractor.local_weight.tel[telid] == 4
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
     assert_allclose(dl1.image, true_charge, rtol=0.1)
     assert_allclose(dl1.peak_time, true_time, rtol=0.1)
     assert dl1.is_valid == True
@@ -417,7 +425,8 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST):
 
     # Test only the 1st pass
     extractor.disable_second_pass = True
-    dl1_pass1 = extractor(waveforms, 1, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1_pass1 = extractor(waveforms, 1, selected_gain_channel, broken_pixels)
     assert_allclose(
         dl1_pass1.image[signal_pixels & integration_window_inside],
         true_charge[signal_pixels & integration_window_inside],
@@ -431,7 +440,7 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST):
 
     # Test also the 2nd pass
     extractor.disable_second_pass = False
-    dl1_pass2 = extractor(waveforms, 1, selected_gain_channel)
+    dl1_pass2 = extractor(waveforms, 1, selected_gain_channel, broken_pixels)
 
     # Check that we have gained signal charge by using the 2nd pass
     # This also checks that the test image has triggered the 2nd pass
@@ -462,7 +471,9 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST):
 def test_waveform_extractor_factory(toymodel):
     waveforms, subarray, telid, selected_gain_channel, true_charge, true_time = toymodel
     extractor = ImageExtractor.from_name("LocalPeakWindowSum", subarray=subarray)
-    dl1 = extractor(waveforms, telid, selected_gain_channel)
+
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, telid, selected_gain_channel, broken_pixels)
     assert_allclose(dl1.image, true_charge, rtol=0.1)
     assert_allclose(dl1.peak_time, true_time, rtol=0.1)
 
@@ -521,7 +532,8 @@ def test_dtype(Extractor, subarray):
 
     waveforms = np.ones((n_pixels, 50), dtype="float64")
     extractor = Extractor(subarray=subarray)
-    dl1 = extractor(waveforms, tel_id, selected_gain_channel)
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
     assert dl1.image.dtype == np.float32
     assert dl1.peak_time.dtype == np.float32
 
@@ -555,9 +567,8 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
         apply_integration_correction=False,
     )
 
-    dl1 = extractor(
-        waveforms=waveforms, telid=tel_id, selected_gain_channel=selected_gain_channel
-    )
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
 
     assert np.allclose(dl1.image[bright_pixels], 18)
     assert np.allclose(dl1.image[~bright_pixels], 0)
diff --git a/ctapipe/image/tests/test_invalid_pixels.py b/ctapipe/image/tests/test_invalid_pixels.py
new file mode 100644
index 00000000..280c84f8
--- /dev/null
+++ b/ctapipe/image/tests/test_invalid_pixels.py
@@ -0,0 +1,61 @@
+import numpy as np
+
+from ctapipe.io.eventsource import EventSource
+
+
+def test_neighbor_average(prod5_gamma_simtel_path):
+    from ctapipe.image.invalid_pixels import NeighborAverage
+
+    with EventSource(prod5_gamma_simtel_path) as source:
+        subarray = source.subarray
+        geometry = subarray.tel[1].camera.geometry
+
+    neighbor_average = NeighborAverage(subarray)
+
+    broken_pixels = np.zeros(geometry.n_pixels, dtype=bool)
+
+    # one border pixel, one non-border pixel
+    border = geometry.get_border_pixel_mask(1)
+    broken_pixels[np.nonzero(border)[0][0]] = True
+    broken_pixels[np.nonzero(~border)[0][0]] = True
+
+    image = np.ones(geometry.n_pixels)
+    image[broken_pixels] = 9999
+    peak_time = np.full(geometry.n_pixels, 20.0)
+    peak_time[broken_pixels] = -1
+
+    interpolated_image, interpolated_peaktime = neighbor_average(
+        tel_id=1, image=image, peak_time=peak_time, pixel_mask=broken_pixels
+    )
+
+    assert np.allclose(interpolated_image[broken_pixels], 1.0)
+    assert np.allclose(interpolated_peaktime[broken_pixels], 20.0)
+
+
+def test_negative_image(prod5_gamma_simtel_path):
+    """Test with negative charges"""
+    from ctapipe.image.invalid_pixels import NeighborAverage
+
+    with EventSource(prod5_gamma_simtel_path) as source:
+        subarray = source.subarray
+        geometry = subarray.tel[1].camera.geometry
+
+    neighbor_average = NeighborAverage(subarray)
+
+    broken_pixels = np.zeros(geometry.n_pixels, dtype=bool)
+
+    # one border pixel, one non-border pixel
+    border = geometry.get_border_pixel_mask(1)
+    broken_pixels[np.nonzero(border)[0][0]] = True
+    broken_pixels[np.nonzero(~border)[0][0]] = True
+
+    image = np.full(geometry.n_pixels, -10.0)
+    image[broken_pixels] = 9999
+    peak_time = np.full(geometry.n_pixels, 20.0)
+    peak_time[broken_pixels] = -1
+
+    interpolated_image, interpolated_peaktime = neighbor_average(
+        tel_id=1, image=image, peak_time=peak_time, pixel_mask=broken_pixels
+    )
+    assert np.allclose(interpolated_image[broken_pixels], -10)
+    assert np.allclose(interpolated_peaktime[broken_pixels], 20.0)
diff --git a/ctapipe/image/tests/test_sliding_window_correction.py b/ctapipe/image/tests/test_sliding_window_correction.py
index a93fec28..c46acf5c 100644
--- a/ctapipe/image/tests/test_sliding_window_correction.py
+++ b/ctapipe/image/tests/test_sliding_window_correction.py
@@ -52,7 +52,10 @@ def test_sw_pulse_lst(prod5_lst):
         "SlidingWindowMaxSum", subarray=subarray, config=config
     )
 
-    dl1: DL1CameraContainer = extractor(waveform, telid, selected_gain_channel)
+    broken_pixels = np.zeros(n_pixels, dtype=bool)
+    dl1: DL1CameraContainer = extractor(
+        waveform, telid, selected_gain_channel, broken_pixels
+    )
     print(dl1.image / charge_true)
     assert_allclose(dl1.image, charge_true, rtol=0.02)
     assert dl1.is_valid