Add PixelResponseNonUniformity effect for fixed per-pixel gain variation (PRNU) (#886)

Author: oczoske

scopesim/effects/electronic/__init__.py

@@ -9,6 +9,7 @@
 - PoorMansHxRGReadoutNoise - simple readout noise for HAWAII detectors
 - BasicReadoutNoise - readout noise
 - ShotNoise - realisation of Poissonian photon noise
+- PixelResponseNonUniformity - per-pixel gain variation (PRNU)
 - DarkCurrent - add dark current
 - LinearityCurve - apply detector (non-)linearity and saturation
 - ReferencePixelBorder
@@ -23,7 +24,7 @@
 
 from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance
 from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise,
-                    ShotNoise, DarkCurrent)
+                    ShotNoise, DarkCurrent, PixelResponseNonUniformity)
 from .exposure import AutoExposure, ExposureIntegration, ExposureOutput
 from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage
 from .dmps import DetectorModePropertiesSetter

scopesim/effects/electronic/noise.py

@@ -134,6 +134,85 @@ def plot_hist(self, det, **kwargs):
         ax.hist(dtcr.data.flatten())
 
 
+class PixelResponseNonUniformity(Effect):
+    """Pixel Response Non-Uniformity (PRNU).
+
+    Models the fixed pattern of per-pixel gain variations across the detector
+    arising from manufacturing differences in quantum efficiency. Each pixel is
+    multiplied by a gain factor drawn from N(1, ``prnu_std``) keyed by detector ID.
+    The gain map is generated once per detector on first use and reused identically 
+    across all subsequent exposures.
+
+    Parameters
+    ----------
+    prnu_std : float or dict
+        Standard deviation of the per-pixel gain distribution. 
+
+    prnu_seed : int, fixed 
+
+    include:  "!DET.include_prnu" 
+
+    Example
+    --------
+
+      - name: prnu
+        description: Pixel response non-uniformity
+        class: PixelResponseNonUniformity
+        kwargs:
+          prnu_std: 0.001
+          prnu_seed: 42
+          include: "!DET.include_prnu" 
+
+    """
+
+    required_keys: ClassVar[set] = set()
+    z_order: ClassVar[tuple[int, ...]] = (805,)
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.meta.update(kwargs)
+        self._gain_maps = {}  # keyed by dtcr_id
+
+    def apply_to(self, obj, **kwargs):
+        if not isinstance(obj, Detector):
+            return obj
+
+        random_seed = from_currsys(self.meta.get("prnu_seed"), self.cmds)
+        id_key = real_colname("id", obj.meta)
+        dtcr_id = obj.meta[id_key] if id_key is not None else None
+
+        prnu_std_meta = from_currsys(self.meta["prnu_std"], self.cmds)
+        if isinstance(prnu_std_meta, dict):
+            prnu_std = float(from_currsys(prnu_std_meta[dtcr_id], self.cmds))
+        elif isinstance(prnu_std_meta, (int, float)):
+            prnu_std = float(prnu_std_meta)
+        else:
+            raise TypeError(
+                "<PixelResponseNonUniformity>.meta['prnu_std'] must be a float "
+                f"or a dict keyed by detector ID, got {type(prnu_std_meta)}"
+            )
+
+        shape = obj._hdu.data.shape
+        if dtcr_id not in self._gain_maps or self._gain_maps[dtcr_id].shape != shape:
+            self._gain_maps[dtcr_id] = np.random.default_rng(random_seed).normal(
+                loc=1.0, scale=prnu_std, size=shape)
+
+        obj._hdu.data = obj._hdu.data * self._gain_maps[dtcr_id]
+        return obj
+
+    def plot(self, det_id=None):
+        if not self._gain_maps:
+            raise RuntimeError("No gain map yet — run a simulation first.")
+        key = det_id if det_id in self._gain_maps else next(iter(self._gain_maps))
+        gain_map = self._gain_maps[key]
+        dev = np.max(np.abs(gain_map - 1.0))
+        fig, ax = figure_factory()
+        im = ax.imshow(gain_map, origin="lower", aspect="auto",
+                       vmin=1 - dev, vmax=1 + dev)
+        fig.colorbar(im, ax=ax, label="per-pixel gain")
+        return fig
+
+
 class ShotNoise(Effect):
     z_order: ClassVar[tuple[int, ...]] = (820,)
 

scopesim/tests/tests_effects/test_prnu.py

@@ -0,0 +1,73 @@
+from matplotlib import pyplot as plt
+import numpy as np
+import pytest
+
+from scopesim.detector import Detector
+from scopesim.effects.electronic import PixelResponseNonUniformity
+from scopesim.optics.image_plane_utils import header_from_list_of_xy
+
+PLOTS = False
+
+def make_detector(value=1000, size=10):
+    hdr = header_from_list_of_xy([-size/2, size/2], [-size/2, size/2], 1, "D")
+    dtcr = Detector(hdr)
+    dtcr._hdu.data[:] = value
+    return dtcr
+
+
+class TestApplyTo:
+    def test_output_std_matches_prnu_std(self):
+        """Relative std of output should be close to prnu_std."""
+        prnu_std = 0.05
+        dtcr = make_detector(value=1000, size=100)
+        PixelResponseNonUniformity(prnu_std=prnu_std, prnu_seed=42).apply_to(dtcr)
+        rel_std = dtcr._hdu.data.std() / dtcr._hdu.data.mean()
+        assert abs(rel_std - prnu_std) < 0.01
+
+    def test_gain_map_is_reused(self):
+        """Same map should be applied on repeated calls (fixed pattern)."""
+        dtcr1 = make_detector()
+        dtcr2 = make_detector()
+        prnu = PixelResponseNonUniformity(prnu_std=0.01, prnu_seed=42)
+        prnu.apply_to(dtcr1)
+        prnu.apply_to(dtcr2)
+        np.testing.assert_array_equal(dtcr1._hdu.data, dtcr2._hdu.data)
+
+    def test_dict_prnu_std(self):
+        """Dict mode: different amplitude per detector ID."""
+        hdr = header_from_list_of_xy([-5, 5], [-5, 5], 1, "D")
+        dtcr = Detector(hdr)
+        dtcr.meta["id"] = "H2RG"
+        dtcr._hdu.data[:] = 1000
+        prnu = PixelResponseNonUniformity(
+            prnu_std={"H2RG": 0.005, "GeoSnap": 0.020}, prnu_seed=42)
+        prnu.apply_to(dtcr)
+        assert dtcr._hdu.data.std() > 0
+
+    def test_multiplicative_zero_signal(self):
+        """Zero signal should remain zero — confirms multiplicative (not additive)."""
+        dtcr = make_detector(value=0)
+        PixelResponseNonUniformity(prnu_std=0.01, prnu_seed=42).apply_to(dtcr)
+        assert dtcr._hdu.data.sum() == 0
+
+    def test_non_detector_passthrough(self):
+        """Non-Detector objects should be returned unchanged."""
+        prnu = PixelResponseNonUniformity(prnu_std=0.01)
+        result = prnu.apply_to("not a detector")
+        assert result == "not a detector"
+
+    def test_invalid_prnu_std_raises(self):
+        prnu = PixelResponseNonUniformity(prnu_std="not a number nor a dict")
+        with pytest.raises(TypeError):
+            prnu.apply_to(make_detector())
+
+    def test_plot_raises_before_simulation(self):
+        prnu = PixelResponseNonUniformity(prnu_std=0.01)
+        with pytest.raises(RuntimeError):
+            prnu.plot()
+
+    def test_plot_returns_figure(self):
+        prnu = PixelResponseNonUniformity(prnu_std=0.01, prnu_seed=42)
+        prnu.apply_to(make_detector())
+        if PLOTS:
+            assert isinstance(prnu.plot(), plt.Figure)