Repository: cta-observatory/ctapipe
Pull request #2037 (2022-07-22 12:48:43)

Fix pixel rotation for very large offsets in NominalFrame
The used vectors to estimate scale and translation were 1 * pixel.unit long, which e.g. for CHEC is 5 times the camera size.

To really get the local scale, I now use the pixel width instead of 1 * pixel.unit.


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/tool.py b/ctapipe/core/tool.py
index ab0fa1de..97e4e80b 100644
--- a/ctapipe/core/tool.py
+++ b/ctapipe/core/tool.py
@@ -6,10 +6,11 @@ import pathlib
 import re
 import textwrap
 from abc import abstractmethod
+from inspect import cleandoc
 from typing import Union
+
 import yaml
 from docutils.core import publish_parts
-from inspect import cleandoc
 
 try:
     import tomli as toml
@@ -18,7 +19,7 @@ try:
 except ImportError:
     HAS_TOML = False
 
-from traitlets import default, List
+from traitlets import List, default
 from traitlets.config import Application, Config, Configurable
 
 from .. import __version__ as version
@@ -312,7 +313,7 @@ class Tool(Application):
         after `Tool.start` when `Tool.run` is called."""
         self.log.info("Goodbye")
 
-    def run(self, argv=None):
+    def run(self, argv=None, raises=False):
         """Run the tool. This automatically calls `initialize()`,
         `start()` and `finish()`
 
@@ -360,6 +361,8 @@ class Tool(Application):
             self.log.exception(f"Caught unexpected exception: {err}")
             Provenance().finish_activity(activity_name=self.name, status="error")
             exit_status = 1  # any other error
+            if raises:
+                raise err
         finally:
             if not {"-h", "--help", "--help-all"}.intersection(self.argv):
                 self.write_provenance()
@@ -525,7 +528,7 @@ def export_tool_config_to_commented_yaml(tool_instance: Tool, classes=None):
     return "\n".join(lines)
 
 
-def run_tool(tool: Tool, argv=None, cwd=None):
+def run_tool(tool: Tool, argv=None, cwd=None, raises=False):
     """
     Utility run a certain tool in a python session without exitinig
 
@@ -539,7 +542,7 @@ def run_tool(tool: Tool, argv=None, cwd=None):
     try:
         # switch to cwd for running and back after
         os.chdir(cwd)
-        tool.run(argv or [])
+        tool.run(argv or [], raises=raises)
     except SystemExit as e:
         return e.code
     finally:
diff --git a/ctapipe/core/traits.py b/ctapipe/core/traits.py
index 793951bc..c82fac01 100644
--- a/ctapipe/core/traits.py
+++ b/ctapipe/core/traits.py
@@ -201,7 +201,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()
@@ -219,7 +219,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/instrument/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index afe69f0f..dac95b00 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -241,7 +241,10 @@ class CameraGeometry:
 
         # also transform the origin and unit vectors,
         # needed to account for translation, rotation / mirroring, scale
-        points = SkyCoord([0, 1, 0], [0, 0, 1], unit=cam.pix_x.unit, frame=cam.frame)
+        width = cam.pixel_width[0].to_value(cam.pix_x.unit)
+        points = SkyCoord(
+            [0, width, 0], [0, 0, width], unit=cam.pix_x.unit, frame=cam.frame
+        )
         points_trans = points.transform_to(frame)
 
         x_name, y_name = list(cam.frame.get_representation_component_names().keys())
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 6a95d8bf..98af3633 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -16,6 +16,7 @@ from ..containers import (
     ArrayEventContainer,
     EventIndexContainer,
     EventType,
+    PixelStatusContainer,
     PointingContainer,
     R0CameraContainer,
     R1CameraContainer,
@@ -597,6 +598,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,
@@ -620,6 +622,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/core/tests/test_tool.py b/ctapipe/core/tests/test_tool.py
index 8e402ecc..26df23c0 100644
--- a/ctapipe/core/tests/test_tool.py
+++ b/ctapipe/core/tests/test_tool.py
@@ -1,13 +1,14 @@
-import os
+import json
 import logging
+import os
 import tempfile
+from pathlib import Path
+
 import pytest
-import json
-from traitlets import Float, TraitError, Dict, Int
+from traitlets import Dict, Float, Int, TraitError
 from traitlets.config import Config
-from pathlib import Path
 
-from .. import Tool, Component
+from .. import Component, Tool
 from ..tool import export_tool_config_to_commented_yaml, run_tool
 
 
@@ -33,7 +34,7 @@ def test_tool_simple():
 
 
 def test_tool_version():
-    """ check that the tool gets an automatic version string"""
+    """check that the tool gets an automatic version string"""
 
     class MyTool(Tool):
         description = "test"
@@ -44,7 +45,7 @@ def test_tool_version():
 
 
 def test_provenance_dir():
-    """ check that the tool gets the provenance dir"""
+    """check that the tool gets the provenance dir"""
 
     class MyTool(Tool):
         description = "test"
@@ -57,7 +58,7 @@ def test_provenance_dir():
 
 
 def test_provenance_log_help(tmpdir):
-    """ check that the tool does not write a provenance log if only the help was run"""
+    """check that the tool does not write a provenance log if only the help was run"""
     from ctapipe.core.tool import run_tool
 
     class MyTool(Tool):
@@ -72,8 +73,9 @@ def test_provenance_log_help(tmpdir):
 
 
 def test_export_config_to_yaml():
-    """ test that we can export a Tool's config to YAML"""
+    """test that we can export a Tool's config to YAML"""
     import yaml
+
     from ctapipe.tools.process import ProcessorTool
 
     tool = ProcessorTool()
@@ -87,7 +89,7 @@ def test_export_config_to_yaml():
 
 
 def test_tool_html_rep(tmp_path):
-    """ check that the HTML rep for Jupyter notebooks works"""
+    """check that the HTML rep for Jupyter notebooks works"""
 
     class MyTool(Tool):
         description = "test"
@@ -100,7 +102,7 @@ def test_tool_html_rep(tmp_path):
         val = Float(1.0, help="val").tag(config=True)
 
     class MyTool2(Tool):
-        """ A docstring description"""
+        """A docstring description"""
 
         userparam = Float(5.0, help="parameter").tag(config=True)
 
@@ -121,7 +123,7 @@ def test_tool_html_rep(tmp_path):
 
 
 def test_tool_current_config():
-    """ Check that we can get the full instance configuration """
+    """Check that we can get the full instance configuration"""
 
     class MyTool(Tool):
         description = "test"
@@ -137,7 +139,7 @@ def test_tool_current_config():
 
 
 def test_tool_current_config_subcomponents():
-    """ Check that we can get the full instance configuration """
+    """Check that we can get the full instance configuration"""
     from ctapipe.core.component import Component
 
     class SubComponent(Component):
@@ -172,7 +174,7 @@ def test_tool_current_config_subcomponents():
 
 
 def test_tool_exit_code():
-    """ Check that we can get the full instance configuration """
+    """Check that we can get the full instance configuration"""
 
     class MyTool(Tool):
 
@@ -340,3 +342,26 @@ def test_invalid_traits(tmp_path, caplog):
         json.dump({"MyTool": {"foo": 5}}, f)
 
     assert run_tool(MyTool(), [f"--config={config}"]) == 2
+
+
+def test_tool_raises():
+    class ToolGood(Tool):
+        name = "Good"
+        description = "This runs without an exception."
+
+        def start(self):
+            print("All good.")
+
+    class ToolBad(Tool):
+        name = "Bad"
+        description = "This tool raises an exception."
+
+        def start(self):
+            raise ValueError("1 does not equal 0.")
+
+    assert run_tool(ToolGood(), raises=True) == 0
+
+    assert run_tool(ToolBad(), raises=False) == 1
+
+    with pytest.raises(ValueError):
+        run_tool(ToolBad(), raises=True)
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