diff --git a/scopesim/effects/electronic/__init__.py b/scopesim/effects/electronic/__init__.py index 9ae7c4ff..5d10e0d4 100644 --- a/scopesim/effects/electronic/__init__.py +++ b/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 diff --git a/scopesim/effects/electronic/noise.py b/scopesim/effects/electronic/noise.py index d2ce522a..51db1e17 100644 --- a/scopesim/effects/electronic/noise.py +++ b/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( + ".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,) diff --git a/scopesim/tests/tests_effects/test_prnu.py b/scopesim/tests/tests_effects/test_prnu.py new file mode 100644 index 00000000..2c7bbcbd --- /dev/null +++ b/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)