From 44d17dc7229f3ac9813d1a32ea0cef25d35f5301 Mon Sep 17 00:00:00 2001 From: Jacob Taylor Date: Wed, 29 Jul 2026 12:18:54 -1000 Subject: [PATCH 1/4] Fix sCMOS per-pixel read noise; characterise three presets from real darks Cross-validated the KURO 1200B, Prime 95B and Marana 4.2B-11 presets against real dark stacks (test_ttf/) and fixed what did not hold up. The bug: the per-pixel read-noise RMS map implied by read_noise_nonuniformity was drawn from the *per-frame* generator, so it was re-randomised every frame. Single-frame spatial statistics looked right, which is why it went unnoticed, but every pixel then had the same expected noise *through time* -- and an sCMOS pixel's read noise is a fixed property of its own source-follower and column ADC. Build the map once from fixed_pattern_seed and cache it in FixedPatternMaps, alongside PRNU and DSNU. The discriminating measurement: split a dark stack in half and correlate the two per-pixel temporal-variance maps. Real detectors give r = 0.89-0.94; the old model gave r = 0.004; it is now r ~ 0.96. This changes generated pixel values wherever read_noise_nonuniformity > 0. Also added, both defaulting to off: - read_noise_rts_fraction / read_noise_rts_factor, a second noisier read-noise population for random-telegraph-signal pixels. ~0.5% of pixels on all three real sensors sit above 3x the median read noise where a single log-normal predicts ~0.01%, and those are the pixels that limit faint-source detection. - detector_glow_edge_scale_px, making detector glow edge-concentrated with an exponential falloff rather than uniform, since amplifier glow is emitted at the array periphery. The Marana shows this clearly (measured 37 px scale). The three presets now carry measured gain, read noise, dark current, bias and non-uniformity terms instead of datasheet values. Biggest corrections: conversion gain 1.25-1.3 -> 0.77-0.87 e-/ADU (the low-signal leg of these dual-gain modes, which is the regime darks probe), and dark_current_nonuniformity 0.03 -> 0.11-0.33, which was about an order of magnitude too low. Other sCMOS presets likely share that error. With no overrides, the updated presets now reproduce the measured variance-vs- exposure curve to 2.9% (Kuro), 6.1% (Prime 95B), 6.1% (Marana), and the pixel-to-pixel spread to within 10% on the two Teledynes. The Marana's spread is still short at long exposure: it has non-Poisson excess noise in its glow regions that a single exponential edge term does not capture. Investigated and deliberately not added: a frame-to-frame bias pedestal drift term. The apparent 0.59 ADU wander was an artifact of taking the median of quantised integer-ADU data; measured from the per-frame spatial mean it is 0.01-0.04 ADU, and where it is larger it varies ~50x between runs, so it is acquisition thermal instability rather than a detector property. Co-Authored-By: Claude Opus 5 --- AGENTS.md | 18 +- CHANGELOG.md | 40 ++++ docs/guides/noise-model.md | 53 ++++++ src/getframes/config.py | 41 ++++- src/getframes/noise.py | 110 +++++++++-- .../presets/data/andor_marana_4_2b_11.toml | 27 ++- .../presets/data/photometrics_prime_95b.toml | 25 ++- .../princeton_instruments_kuro_1200b.toml | 25 ++- tests/test_realism.py | 174 +++++++++++++++++- 9 files changed, 470 insertions(+), 43 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 0b362da..8381aac 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -64,12 +64,24 @@ Data flows one way: `presets` → `CameraConfig` → `Scene` → `Camera` → `b backend's camera-owned generator. Every generation method accepts a `seed`. Never call global NumPy/CuPy random state. CPU and GPU streams repeat within a backend but are statistically, not pixel-for-pixel, matched across devices. -3. **Physics is auditable.** Noise models live as small, documented, pure functions +3. **Fixed patterns come from `fixed_pattern_seed`, never the per-frame RNG.** + Anything that is a property of the *silicon* — PRNU, DSNU, hot pixels, defects, + amplifier gain/offset, structured bias, and the sCMOS per-pixel read-noise RMS — + belongs in `fixed_pattern_maps()` / `FixedPatternMaps`, keyed on + `fixed_pattern_seed`, so it is identical in every frame and therefore removable + by a master frame. Only the *draw* is per-frame (e.g. read noise re-draws the + Gaussian each frame, but its per-pixel sigma does not change). Getting this wrong + is invisible in single-frame spatial statistics and only shows up in per-pixel + statistics *through time* — the failure mode that shipped in + `read_noise_nonuniformity` until it was caught by cross-validation against real + dark stacks. When adding a detector-structure feature, ask whether it should + repeat across frames, and add a test on the temporal statistic if so. +4. **Physics is auditable.** Noise models live as small, documented, pure functions in `noise.py`. Document the units (electrons vs. ADU) and the model in the docstring. State assumptions; cite the model form. -4. **Units are explicit.** Field/variable names carry units (`_e`, `_adu`, `_um`, +5. **Units are explicit.** Field/variable names carry units (`_e`, `_adu`, `_um`, `_c`, `_s`, `_e_per_s`, `_e_per_adu`). Keep this convention. -5. **Typed and validated.** Full type hints (`mypy --strict` passes). Validate +6. **Typed and validated.** Full type hints (`mypy --strict` passes). Validate inputs in `CameraConfig.__post_init__` and raise informative `ValueError`s. ## Adding a camera preset diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c18da1..b83f34b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,8 +6,48 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Fixed + +- **sCMOS per-pixel read noise is now a fixed property of the sensor.** The + per-pixel read-noise RMS map implied by `read_noise_nonuniformity` was drawn from + the *per-frame* generator, so it was re-randomised in every frame. Single-frame + spatial statistics were unaffected, but every pixel ended up with the same + expected noise *through time*, which is not how an sCMOS behaves: each pixel has + its own source-follower and column ADC. The map is now built once from + `fixed_pattern_seed` and cached in `FixedPatternMaps`, alongside PRNU and DSNU. + Verified against dark stacks from three real back-illuminated sCMOS cameras + (KURO 1200B, Prime 95B, Marana 4.2B-11): splitting a stack in half and + correlating the two per-pixel temporal-variance maps gives r = 0.89–0.94 on the + real detectors and r = 0.004 with the old model, now r ≈ 0.96. + **This changes generated pixel values** for any configuration with + `read_noise_nonuniformity > 0`; `read_noise_e` and the spatial statistics of a + single frame are unchanged. + ### Added +- `read_noise_rts_fraction` / `read_noise_rts_factor`: an optional second, + noisier read-noise population modelling the random-telegraph-signal (RTS) pixels + of a real sCMOS array. Measured on three real sensors, ~0.5% of pixels sit above + 3x the median read noise where a single log-normal predicts ~0.01%; these are the + pixels that limit faint-source detection. Defaults to off. +- `detector_glow_edge_scale_px`: makes `detector_glow_e_per_s` edge-concentrated + with an exponential falloff, instead of uniform, modelling amplifier glow emitted + at the array periphery. Renormalised so the array mean is unchanged; still fixed + and exposure-scaling, so an exposure-matched master dark removes it. Defaults to + `0` (uniform, the previous behaviour). + +### Changed + +- The `princeton_instruments_kuro_1200b`, `photometrics_prime_95b`, and + `andor_marana_4_2b_11` presets now carry **measured** conversion gain, read noise, + dark current, bias offset, and non-uniformity terms, fitted from a per-pixel dark + photon-transfer analysis of real frames rather than taken from datasheets. The + largest corrections: conversion gain (1.25-1.3 -> 0.77-0.87 e-/ADU, the low-signal + leg of these dual-gain modes) and `dark_current_nonuniformity` (0.03 -> 0.11-0.33, + which had been roughly an order of magnitude too low). Each preset documents the + operating mode and temperature the values apply to. Other sCMOS presets still + carry datasheet-derived DSNU and are likely low for the same reason. + - Full-detector region-of-interest simulation through `CameraConfig.roi=(left, top, width, height)`. Cameras accept and return ROI-shaped arrays while evaluating detector physics and fixed patterns on the diff --git a/docs/guides/noise-model.md b/docs/guides/noise-model.md index c8cecbf..6b2578c 100644 --- a/docs/guides/noise-model.md +++ b/docs/guides/noise-model.md @@ -83,6 +83,59 @@ EMCCD or eAPD this is applied after the gain stage, which is why a high mean gai makes the effective (input-referred) read noise sub-electron — the eAPD's `read_noise_e` is the pre-avalanche amplifier noise, divided down by `em_gain`. +#### Per-pixel read noise (sCMOS) + +An sCMOS pixel has its own source-follower and column ADC, so its read noise is a +**fixed property of that pixel** rather than a single array-wide number. +`read_noise_nonuniformity` sets the fractional width of a log-normal distribution +of per-pixel RMS. Like PRNU and DSNU, the resulting map is drawn once from +`fixed_pattern_seed`, so it is identical in every frame — only the Gaussian draw +itself is per-frame. + +That distinction is measurable, and it is the reason it matters: take a stack of +darks, compute each pixel's variance *through time*, split the stack in half, and +the two variance maps agree pixel-for-pixel. On three real back-illuminated sCMOS +cameras that split-half correlation is 0.89–0.94. + +`read_noise_e` is the *scale* of this distribution, which has unit mean — so the +mean per-pixel RMS is `read_noise_e` and the median is +`read_noise_e * exp(-read_noise_nonuniformity**2 / 2)`, a few percent lower. + +Real arrays also carry a **random-telegraph-signal (RTS)** population: a small +fraction of pixels whose trapped-charge switching makes them much noisier than the +log-normal core predicts. Measured on real sensors, ~0.5% of pixels sit above 3× +the median read noise, where a bare log-normal would put ~0.01%. Set +`read_noise_rts_fraction` (typically 0.005–0.03) and `read_noise_rts_factor` to +include them: + +```python +from getframes import Camera, load_preset + +cam = Camera(load_preset("princeton_instruments_kuro_1200b")) +cfg = cam.config +print(cfg.read_noise_e, cfg.read_noise_nonuniformity) +print(cfg.read_noise_rts_fraction, cfg.read_noise_rts_factor) +``` + +These are the pixels that limit faint-source detection, so they matter for any +threshold or detection-completeness study. + +#### Detector glow + +`detector_glow_e_per_s` adds self-emission that scales with exposure. By default it +is uniform. Real amplifier glow is emitted by the readout electronics on the array +periphery, so setting `detector_glow_edge_scale_px` concentrates it near the +detector edges with an exponential falloff: + +``` +glow(x, y) = A * exp(-d_edge(x, y) / detector_glow_edge_scale_px) +``` + +with `A` chosen so the array **mean** is still `detector_glow_e_per_s` — meaning the +edges run hotter and the centre cooler than that figure. The pattern is fixed and +exposure-scaling, so an exposure-matched master dark still removes it. The +`andor_marana_4_2b_11` preset carries a measured 37-pixel falloff scale. + ### 7. Digitisation For ordinary detectors, `full_well_e` is both the image-area charge capacity and diff --git a/src/getframes/config.py b/src/getframes/config.py index b6b9b6d..bc45633 100644 --- a/src/getframes/config.py +++ b/src/getframes/config.py @@ -94,12 +94,40 @@ class CameraConfig: bias_offset_adu: Electronic offset (pedestal) added to every pixel, in ADU. read_noise_e: - RMS read noise in electrons. For sCMOS this is the *median*; see - ``read_noise_nonuniformity``. + RMS read noise in electrons. When ``read_noise_nonuniformity`` is zero this + is every pixel's read noise. Otherwise it is the *scale* of the per-pixel + distribution, which is log-normal with unit mean --- so the mean per-pixel + RMS is ``read_noise_e`` and the median is + ``read_noise_e * exp(-read_noise_nonuniformity**2 / 2)``, a few percent + lower. See ``read_noise_nonuniformity`` and ``read_noise_rts_fraction``. read_noise_nonuniformity: Fractional pixel-to-pixel spread of the read-noise RMS (e.g. ``0.3`` for a 30% log-normal spread). Models the per-pixel read-noise distribution of sCMOS sensors. ``0`` gives a single uniform read noise. + + The resulting per-pixel RMS is a *fixed* property of the sensor (drawn from + ``fixed_pattern_seed``, like PRNU and DSNU), not re-drawn each frame, so a + pixel's temporal noise is repeatable across a stack --- which is what is + measured in practice. + read_noise_rts_fraction: + Fraction of pixels belonging to a second, noisier read-noise population, + in ``[0, 1]``. These are the random-telegraph-signal (RTS) pixels of a real + sCMOS array, whose trapped-charge switching gives the read-noise histogram a + tail much heavier than the single log-normal of + ``read_noise_nonuniformity``. ``0`` disables the second population. + Measured values for back-illuminated sCMOS are around ``0.005-0.03``. + read_noise_rts_factor: + Multiplier applied to the read-noise RMS of the RTS population selected by + ``read_noise_rts_fraction``. Ignored when that fraction is ``0``. + detector_glow_edge_scale_px: + Exponential falloff scale, in pixels, of the ``detector_glow_e_per_s`` term + away from the detector edges. Amplifier/array glow originates at the readout + electronics on the array periphery, so real glow is edge-concentrated rather + than uniform. ``0`` (the default) keeps the glow uniform. When positive, the + map is renormalised so the *mean* glow over the array is still + ``detector_glow_e_per_s``, which means the edges run hotter and the centre + cooler than that figure. The pattern is fixed and exposure-scaling, so an + exposure-matched master dark still removes it. nonlinearity: Fractional signal compression at full well, in ``[0, 0.5)``. The collected charge is bent as ``q -> q * (1 - nonlinearity * q / full_well_e)``, so a @@ -254,6 +282,9 @@ class CameraConfig: detector_glow_e_per_s: float = 0.0 prnu: float = 0.0 read_noise_nonuniformity: float = 0.0 + read_noise_rts_fraction: float = 0.0 + read_noise_rts_factor: float = 2.5 + detector_glow_edge_scale_px: float = 0.0 nonlinearity: float = 0.0 nonlinearity_coeffs: tuple[float, ...] | None = None cti: float = 0.0 @@ -359,6 +390,12 @@ def _validate(self) -> None: raise ValueError("prnu must be non-negative.") if self.read_noise_nonuniformity < 0: raise ValueError("read_noise_nonuniformity must be non-negative.") + if not 0.0 <= self.read_noise_rts_fraction <= 1.0: + raise ValueError("read_noise_rts_fraction must be in [0, 1].") + if self.read_noise_rts_factor < 0: + raise ValueError("read_noise_rts_factor must be non-negative.") + if self.detector_glow_edge_scale_px < 0: + raise ValueError("detector_glow_edge_scale_px must be non-negative.") if not 0.0 <= self.nonlinearity < 0.5: raise ValueError("nonlinearity must be in [0, 0.5).") if self.nonlinearity_coeffs is not None and len(self.nonlinearity_coeffs) == 0: diff --git a/src/getframes/noise.py b/src/getframes/noise.py index 0967109..9ade639 100644 --- a/src/getframes/noise.py +++ b/src/getframes/noise.py @@ -10,7 +10,8 @@ 1. Mean photo signal: ``(photon_rate + background) * t_exp * QE`` electrons, modulated per pixel by photo-response non-uniformity (PRNU). 2. Mean dark signal: ``D(T) * t_exp`` electrons (temperature-scaled), modulated by - dark-signal non-uniformity (DSNU) and hot pixels. + dark-signal non-uniformity (DSNU) and hot pixels, plus detector glow (uniform, + or edge-concentrated via ``detector_glow_edge_scale_px``). 3. Shot noise: the total electrons are Poisson-distributed about that mean. 4. Clock-induced charge (EMCCD) adds a small Poisson term. 5. Cosmic rays (single pixels or extended tracks). @@ -19,6 +20,8 @@ 7. Detector nonlinearity (single-parameter or polynomial). 8. EM register / avalanche multiplication with its stochastic excess noise. 9. kTC/reset noise and read noise: Gaussian in electrons, at the output amplifier. + The per-pixel read-noise RMS is a *fixed* sensor property (sCMOS), including an + optional random-telegraph-signal (RTS) tail population. 10. Conversion to ADU via (optionally per-amplifier) gain, plus the bias pedestal and any structured-bias pattern; dead pixels/columns read as defects. 11. Saturation at full well / ADC range and quantisation to integers. @@ -56,6 +59,7 @@ _FPN_STREAM_AMP_OFFSET = 4 _FPN_STREAM_BIAS = 5 _FPN_STREAM_DEFECT = 6 +_FPN_STREAM_READ_NOISE = 7 class FixedPatternMaps(NamedTuple): @@ -67,6 +71,7 @@ class FixedPatternMaps(NamedTuple): amplifier_offset: Any bias_structure: Any defect_mask: Any | None + read_noise_sigma: Any def _fixed_pattern_rng(config: CameraConfig, stream: int, backend: ArrayBackend) -> Any: @@ -118,9 +123,50 @@ def fixed_pattern_maps( offset, _bias_structure_map(config, resolved), _defect_mask(config, resolved), + _read_noise_sigma_map(config, resolved, float_dtype=float_dtype), ) +def _read_noise_sigma_map( + config: CameraConfig, + backend: ArrayBackend | None = None, + *, + float_dtype: DTypeLike = DEFAULT_FLOAT_DTYPE, +) -> Any: + """Per-pixel read-noise RMS in electrons, or a scalar when it is uniform. + + Real sCMOS read noise is a property of each pixel's own source-follower and + ADC chain, so this map is drawn from the *fixed-pattern* stream (keyed on + ``fixed_pattern_seed``) and is identical in every frame the camera produces. + That is what makes the per-pixel *temporal* noise repeatable, and it is + directly measurable: split a dark stack in half, take each half's per-pixel + variance, and the two maps correlate. + + The distribution is a log-normal core of fractional width + ``read_noise_nonuniformity``, optionally with a second, noisier population + covering ``read_noise_rts_fraction`` of pixels whose RMS is multiplied by + ``read_noise_rts_factor``. That second population models random-telegraph-signal + (RTS) pixels, which give real sCMOS arrays a read-noise histogram with a + markedly heavier tail than a single log-normal. + """ + if config.read_noise_e <= 0: + return 0.0 + resolved = backend or get_backend() + shape = config.resolution + if config.read_noise_nonuniformity <= 0 and config.read_noise_rts_fraction <= 0: + return float(config.read_noise_e) + rng = _fixed_pattern_rng(config, _FPN_STREAM_READ_NOISE, resolved) + spread = config.read_noise_nonuniformity + if spread > 0: + sigma = config.read_noise_e * rng.lognormal(mean=-0.5 * spread**2, sigma=spread, size=shape) + else: + sigma = resolved.xp.full(shape, float(config.read_noise_e)) + if config.read_noise_rts_fraction > 0: + rts = rng.random(shape) < config.read_noise_rts_fraction + sigma[rts] *= config.read_noise_rts_factor + return sigma.astype(float_dtype) + + def dark_signal_map( config: CameraConfig, exposure_s: float, @@ -164,11 +210,14 @@ def dark_signal_map( hot_mask = rng.random(signal.shape) < config.hot_pixel_fraction signal[hot_mask] *= config.hot_pixel_factor - # Detector glow: a uniform self-emission term that scales with exposure (and so - # is removed by an exposure-matched master dark). Added after DSNU/hot pixels, + # Detector glow: a self-emission term that scales with exposure (and so is + # removed by an exposure-matched master dark). Added after DSNU/hot pixels, # which describe the dark *current*, not the glow. In place to preserve dtype. if config.detector_glow_e_per_s > 0 and exposure_s > 0: - signal += config.detector_glow_e_per_s * exposure_s + if config.detector_glow_edge_scale_px > 0: + signal += _glow_profile(config, resolved, float_dtype) * exposure_s + else: + signal += config.detector_glow_e_per_s * exposure_s return signal @@ -502,6 +551,41 @@ def _bias_structure_map(config: CameraConfig, backend: ArrayBackend | None = Non return xp.broadcast_to(scaled, (height, width)).astype(np.float64) +def _glow_profile( + config: CameraConfig, + backend: ArrayBackend | None = None, + float_dtype: DTypeLike = DEFAULT_FLOAT_DTYPE, +) -> Any: + """Edge-concentrated detector-glow rate in e-/pixel/s. + + Amplifier and array glow is emitted by the readout electronics around the + array periphery, so it falls off into the detector rather than sitting at a + uniform level. The model is an exponential in the distance to the nearest + edge:: + + g(x, y) = A * exp(-d_edge(x, y) / detector_glow_edge_scale_px) + + with ``A`` set so the *mean* of ``g`` over the array equals + ``detector_glow_e_per_s``. It is deterministic (no randomness), fixed for a + given sensor, and scales with exposure, so an exposure-matched master dark + removes it. + """ + resolved = backend or get_backend() + xp = resolved.xp + height, width = config.resolution + scale = config.detector_glow_edge_scale_px + rows = xp.arange(height, dtype=float_dtype).reshape(height, 1) + cols = xp.arange(width, dtype=float_dtype).reshape(1, width) + d_row = xp.minimum(rows, height - 1 - rows) + d_col = xp.minimum(cols, width - 1 - cols) + profile = xp.exp(-xp.minimum(d_row, d_col) / scale) + mean = resolved.scalar(profile.mean()) + if mean <= 0: + return xp.zeros((height, width), dtype=float_dtype) + profile *= config.detector_glow_e_per_s / mean + return profile + + def _defect_mask(config: CameraConfig, backend: ArrayBackend | None = None) -> Any | None: """A fixed boolean map of dead pixels/columns (``True`` = no response), or ``None``. @@ -619,16 +703,16 @@ def normal_noise(sigma: Any) -> Any: if config.reset_noise_e > 0: signal += normal_noise(config.reset_noise_e) - # Read noise in electrons, added at the amplifier. + # Read noise in electrons, added at the amplifier. The per-pixel RMS is a + # fixed property of the sensor (see :func:`_read_noise_sigma_map`), so only the + # Gaussian draw itself is per-frame. if config.read_noise_e > 0: - if config.read_noise_nonuniformity > 0: - spread = config.read_noise_nonuniformity - sigma_map = config.read_noise_e * rng.lognormal( - mean=-0.5 * spread**2, sigma=spread, size=signal.shape - ) - signal += normal_noise(sigma_map) - else: - signal += normal_noise(config.read_noise_e) + sigma_map = ( + fixed_patterns.read_noise_sigma + if fixed_patterns is not None + else _read_noise_sigma_map(config, resolved, float_dtype=signal.dtype) + ) + signal += normal_noise(sigma_map) if fixed_patterns is None: gain_map, amp_offset = _amplifier_maps(config, resolved) diff --git a/src/getframes/presets/data/andor_marana_4_2b_11.toml b/src/getframes/presets/data/andor_marana_4_2b_11.toml index 6ec0c72..4a6c552 100644 --- a/src/getframes/presets/data/andor_marana_4_2b_11.toml +++ b/src/getframes/presets/data/andor_marana_4_2b_11.toml @@ -10,19 +10,30 @@ pixel_size_um = 11.0 quantum_efficiency = 0.95 full_well_e = 85000.0 bit_depth = 16 -gain_e_per_adu = 1.3 # representative extended-dynamic-range conversion -bias_offset_adu = 100.0 -read_noise_e = 1.6 # median -read_noise_nonuniformity = 0.2 -dark_current_e_per_s = 0.3 -dark_current_ref_temp_c = -45.0 # liquid-cooled specification +gain_e_per_adu = 0.8008 # measured; NOT the 1.2-1.4 full-range figure +bias_offset_adu = 99.70 # measured +read_noise_e = 1.0763 # measured: median 1.05 e-, 25-75 pct 0.89-1.22 +read_noise_nonuniformity = 0.2444 # measured log-normal core +dark_current_e_per_s = 0.3276 # measured, glow-subtracted +dark_current_ref_temp_c = -45.0 # temperature of the measurement dark_current_doubling_temp_c = 6.0 prnu = 0.005 # published <0.5% at half range nonlinearity = 0.003 # published linearity >99.7% -dark_current_nonuniformity = 0.03 +dark_current_nonuniformity = 0.110 # measured on the glow-free interior supported_binnings = [1, 2, 3, 4, 8] # FPGA post-readout (digital) summation binning_method = "digital" -notes = "Extended-dynamic-range mode; published sensor QE curve and liquid-cooled dark current. FPGA post-readout binning does not reduce the sCMOS readout time." +notes = "Values below the measured-values banner, plus gain/read noise/dark current/bias, were characterised from real 16-bit HDR darks at -45 C (see the banner). This sensor shows edge-concentrated amplifier glow, modelled via detector_glow_edge_scale_px; the dark_current_e_per_s figure is the glow-subtracted base rate. Published QE curve retained. FPGA post-readout binning does not reduce the sCMOS readout time." + +# --- Measured values, cross-validated against real dark frames ----------------- +# Conversion gain, read noise, dark current, bias and the non-uniformity terms +# below were fitted from a per-pixel dark photon-transfer analysis of 150 dark +# frames per exposure at -45 C, 16-bit HDR readout. The gain is the LOW-SIGNAL leg of this +# dual-gain mode, i.e. the regime that applies to dark, bias and faint-source +# work; a full-range gain figure for the same camera will be higher. +read_noise_rts_fraction = 0.0162 # measured RTS tail: 0.4% of pixels above 3x median +read_noise_rts_factor = 2.590 # measured +detector_glow_e_per_s = 0.0085 # measured edge glow, array mean +detector_glow_edge_scale_px = 37.0 # measured exponential falloff from the array edge [extra] readout_modes = ["16-bit high dynamic range", "12-bit fast speed"] diff --git a/src/getframes/presets/data/photometrics_prime_95b.toml b/src/getframes/presets/data/photometrics_prime_95b.toml index ddb9758..79013ec 100644 --- a/src/getframes/presets/data/photometrics_prime_95b.toml +++ b/src/getframes/presets/data/photometrics_prime_95b.toml @@ -10,19 +10,28 @@ pixel_size_um = 11.0 quantum_efficiency = 0.95 full_well_e = 80000.0 bit_depth = 16 -gain_e_per_adu = 1.25 # representative combined-gain conversion -bias_offset_adu = 100.0 -read_noise_e = 1.8 # RMS (1.6 e- median) -read_noise_nonuniformity = 0.2 -dark_current_e_per_s = 0.3 -dark_current_ref_temp_c = -25.0 # liquid-cooled specification +gain_e_per_adu = 0.7678 # measured +bias_offset_adu = 100.46 # measured +read_noise_e = 1.4658 # measured: median 1.44 e-, 25-75 pct 1.30-1.63 +read_noise_nonuniformity = 0.2211 # measured log-normal core +dark_current_e_per_s = 0.413 # measured median +dark_current_ref_temp_c = -20.0 # temperature of the measurement dark_current_doubling_temp_c = 6.0 prnu = 0.005 nonlinearity = 0.005 -dark_current_nonuniformity = 0.03 +dark_current_nonuniformity = 0.327 # measured; was 0.03, an order of magnitude low supported_binnings = [1, 2] # 2x2 FPGA post-readout (digital) binning binning_method = "digital" -notes = "Combined-gain 16-bit mode; published peak QE and liquid-cooled dark current." +notes = "Gain, read noise, dark current, bias and the non-uniformity terms were characterised from real 16-bit high-sensitivity darks at -20 C (see the measured-values banner). Published QE curve retained. Note the dark-current figure came from an air-cooled series spanning several days, so it carries a sensor-temperature systematic." + +# --- Measured values, cross-validated against real dark frames ----------------- +# Conversion gain, read noise, dark current, bias and the non-uniformity terms +# below were fitted from a per-pixel dark photon-transfer analysis of 200 dark +# frames per exposure at -20 C, 16-bit high-sensitivity readout. The gain is the LOW-SIGNAL leg of this +# dual-gain mode, i.e. the regime that applies to dark, bias and faint-source +# work; a full-range gain figure for the same camera will be higher. +read_noise_rts_fraction = 0.0145 # measured RTS tail: 0.5% of pixels above 3x median +read_noise_rts_factor = 2.729 # measured [extra] readout_modes = ["16-bit combined gain", "12-bit full well", "12-bit balanced", "12-bit sensitivity"] diff --git a/src/getframes/presets/data/princeton_instruments_kuro_1200b.toml b/src/getframes/presets/data/princeton_instruments_kuro_1200b.toml index 3e88a57..4007d8c 100644 --- a/src/getframes/presets/data/princeton_instruments_kuro_1200b.toml +++ b/src/getframes/presets/data/princeton_instruments_kuro_1200b.toml @@ -10,19 +10,28 @@ pixel_size_um = 11.0 quantum_efficiency = 0.95 full_well_e = 80000.0 bit_depth = 16 -gain_e_per_adu = 1.25 # representative 80 ke-/16-bit conversion -bias_offset_adu = 100.0 -read_noise_e = 1.3 # RMS -read_noise_nonuniformity = 0.2 -dark_current_e_per_s = 0.7 -dark_current_ref_temp_c = -25.0 # representative cooled operating point +gain_e_per_adu = 0.8673 # measured; test report quotes 0.87 +bias_offset_adu = 99.91 # measured +read_noise_e = 1.6597 # measured: median 1.62 e-, 25-75 pct 1.44-1.87 +read_noise_nonuniformity = 0.2469 # measured log-normal core +dark_current_e_per_s = 0.825 # measured median +dark_current_ref_temp_c = -20.0 # temperature of the measurement dark_current_doubling_temp_c = 6.0 prnu = 0.005 nonlinearity = 0.005 -dark_current_nonuniformity = 0.03 +dark_current_nonuniformity = 0.234 # measured; was 0.03, an order of magnitude low supported_binnings = [1, 2, 4] binning_method = "digital" -notes = "Published peak-QE/full-frame specification; characterize the selected readout mode for quantitative work. No on-chip binning documented for this sCMOS family; any software summing is digital (post-read) and is not characterized here." +notes = "Gain, read noise, dark current, bias and the non-uniformity terms were characterised from real 16-bit high-sensitivity darks at -20 C (see the measured-values banner); the measured gain agrees with the unit test report (0.87 e-/ADU). Published QE curve retained. The 12-bit high-speed readout is a different operating point entirely (measured 2.19 e-/ADU, 2.96 e- read noise) --- re-characterise if you use it. No on-chip binning documented for this sCMOS family; any software summing is digital." + +# --- Measured values, cross-validated against real dark frames ----------------- +# Conversion gain, read noise, dark current, bias and the non-uniformity terms +# below were fitted from a per-pixel dark photon-transfer analysis of 200 dark +# frames per exposure at -20 C, 16-bit high-sensitivity readout. The gain is the LOW-SIGNAL leg of this +# dual-gain mode, i.e. the regime that applies to dark, bias and faint-source +# work; a full-range gain figure for the same camera will be higher. +read_noise_rts_fraction = 0.0171 # measured RTS tail: 0.5% of pixels above 3x median +read_noise_rts_factor = 2.669 # measured [extra] readout_modes = ["12-bit", "16-bit", "rolling shutter", "effective global shutter"] diff --git a/tests/test_realism.py b/tests/test_realism.py index b86d081..a3f7f13 100644 --- a/tests/test_realism.py +++ b/tests/test_realism.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from getframes import Camera, CameraConfig, SensorType, load_preset +from getframes import Camera, CameraConfig, SensorType, load_preset, noise def base_config(**overrides): @@ -84,12 +84,184 @@ def test_read_noise_nonuniformity_increases_bias_scatter(): assert s > u +def _temporal_variance_halves(camera, exposure, n_frames, temperature, seed): + """Per-pixel temporal variance of the even- and odd-indexed frames of a stack.""" + frames = [ + np.asarray(f, dtype=np.float64) + for f in camera.dark_series(exposure, n_frames, temperature, seed=seed) + ] + return np.var(frames[::2], axis=0, ddof=1), np.var(frames[1::2], axis=0, ddof=1) + + +def test_per_pixel_read_noise_is_a_fixed_sensor_property(): + """A pixel's *temporal* noise must repeat across a stack, not be re-drawn per frame. + + Real sCMOS read noise is a property of each pixel's own amplifier/ADC chain. The + measurable signature is that two disjoint halves of a dark stack produce + correlated per-pixel variance maps. If the sigma map were re-drawn every frame, + every pixel would share one expected variance and the correlation would vanish. + """ + cam = Camera( + base_config( + resolution=(96, 96), read_noise_e=2.0, read_noise_nonuniformity=0.4, gain_e_per_adu=0.5 + ) + ) + a, b = _temporal_variance_halves(cam, 0.0, 160, -100.0, seed=11) + r = float(np.corrcoef(a.ravel(), b.ravel())[0, 1]) + assert r > 0.8, f"per-pixel read noise is not repeatable across frames (r={r:.3f})" + + +def test_uniform_read_noise_gives_no_variance_structure(): + """The counterpart: with no non-uniformity the variance map is pure sampling noise.""" + cam = Camera(base_config(resolution=(96, 96), read_noise_e=2.0, gain_e_per_adu=0.5)) + a, b = _temporal_variance_halves(cam, 0.0, 160, -100.0, seed=11) + r = float(np.corrcoef(a.ravel(), b.ravel())[0, 1]) + assert abs(r) < 0.2, f"uniform read noise should leave no fixed structure (r={r:.3f})" + + +def test_read_noise_sigma_map_keyed_on_fixed_pattern_seed(): + shape = (48, 48) + same_a = noise._read_noise_sigma_map( + base_config(resolution=shape, read_noise_nonuniformity=0.3, fixed_pattern_seed=7) + ) + same_b = noise._read_noise_sigma_map( + base_config(resolution=shape, read_noise_nonuniformity=0.3, fixed_pattern_seed=7) + ) + other = noise._read_noise_sigma_map( + base_config(resolution=shape, read_noise_nonuniformity=0.3, fixed_pattern_seed=8) + ) + np.testing.assert_array_equal(same_a, same_b) + assert not np.array_equal(same_a, other) + + +def test_read_noise_scale_is_the_mean_of_the_per_pixel_rms(): + """``read_noise_e`` is the unit-mean scale; the median sits exp(-s^2/2) below it.""" + spread = 0.4 + sigma = noise._read_noise_sigma_map( + base_config(resolution=(400, 400), read_noise_e=2.0, read_noise_nonuniformity=spread) + ) + assert np.mean(sigma) == pytest.approx(2.0, rel=0.02) + assert np.median(sigma) == pytest.approx(2.0 * np.exp(-0.5 * spread**2), rel=0.02) + + +# --- sCMOS RTS (random-telegraph-signal) tail ------------------------------ + + +def test_rts_population_adds_a_heavy_read_noise_tail(): + kw = {"resolution": (400, 400), "read_noise_e": 2.0, "read_noise_nonuniformity": 0.2} + core = noise._read_noise_sigma_map(base_config(**kw)) + tailed = noise._read_noise_sigma_map( + base_config(**kw, read_noise_rts_fraction=0.02, read_noise_rts_factor=3.0) + ) + # The core is untouched: same median, but a far heavier upper tail. + assert np.median(tailed) == pytest.approx(np.median(core), rel=0.02) + assert np.percentile(tailed, 99.9) > 2.0 * np.percentile(core, 99.9) + # Roughly half the RTS pixels (those whose core draw is above the median) land + # beyond 3x the median, lifting the >3x fraction from ~1e-9 to ~1%. + core_frac = float(np.mean(core > 3.0 * np.median(core))) + tail_frac = float(np.mean(tailed > 3.0 * np.median(tailed))) + assert core_frac < 1e-4 + assert 0.005 < tail_frac < 0.02 + + +def test_rts_fraction_selects_about_the_right_number_of_pixels(): + sigma = noise._read_noise_sigma_map( + base_config( + resolution=(500, 500), + read_noise_e=2.0, + read_noise_nonuniformity=0.0, + read_noise_rts_fraction=0.05, + read_noise_rts_factor=4.0, + ) + ) + assert float(np.mean(sigma > 2.0 * 2.0)) == pytest.approx(0.05, abs=0.005) + + +def test_rts_defaults_off_leaves_the_sigma_map_log_normal(): + spread = 0.3 + cfg = base_config(read_noise_nonuniformity=spread) + assert cfg.read_noise_rts_fraction == 0.0 + sigma = noise._read_noise_sigma_map(cfg.replace(resolution=(300, 300))) + # A bare log-normal of this width puts ~1.3e-4 of pixels above 3x the median; + # the RTS population is what lifts that to the ~0.5% seen on real sCMOS. + from scipy.stats import norm + + expected = float(norm.sf(np.log(3.0) / spread)) + assert float(np.mean(sigma > 3.0 * np.median(sigma))) == pytest.approx(expected, abs=3e-4) + + +# --- structured (edge-concentrated) detector glow -------------------------- + + +def test_glow_edge_scale_concentrates_glow_at_the_edges(): + cfg = base_config( + resolution=(128, 128), detector_glow_e_per_s=1.0, detector_glow_edge_scale_px=8.0 + ) + profile = np.asarray(noise._glow_profile(cfg)) + # Renormalised so the array mean is still detector_glow_e_per_s. + assert profile.mean() == pytest.approx(1.0, rel=1e-9) + assert profile[0, 64] > 10.0 * profile[64, 64] + # Monotonic falling from the edge towards the middle along a central column. + column = profile[: 128 // 2, 64] + assert np.all(np.diff(column) < 0) + + +def test_glow_edge_scale_zero_is_uniform(): + cfg = base_config( + resolution=(64, 64), detector_glow_e_per_s=2.0, detector_glow_edge_scale_px=0.0 + ) + hot = Camera(cfg).dark_frame(10.0, -100.0, seed=2) + edge = np.asarray(hot)[0, :].mean() + middle = np.asarray(hot)[32, :].mean() + assert edge == pytest.approx(middle, rel=0.05) + + +def test_structured_glow_is_removed_by_an_exposure_matched_master_dark(): + """Glow is fixed and exposure-scaling, so a master dark still calibrates it out.""" + cfg = base_config( + resolution=(64, 64), + detector_glow_e_per_s=5.0, + detector_glow_edge_scale_px=6.0, + read_noise_e=1.0, + ) + cam = Camera(cfg) + master = cam.master_dark(10.0, n_frames=64, temperature=-100.0, seed=5) + frame = cam.dark_frame(10.0, -100.0, seed=99) + raw = np.asarray(frame, dtype=np.float64) + residual = raw - np.asarray(master, dtype=np.float64) + # The glow puts a large edge-to-centre step in the raw frame; subtracting the + # master must remove essentially all of it. What is left is shot noise on the + # (bright) glow itself, so compare against the structure being removed rather + # than against zero. + raw_step = abs(raw[0, :].mean() - raw[32, :].mean()) + residual_step = abs(residual[0, :].mean() - residual[32, :].mean()) + assert raw_step > 100.0, "the test glow should be strongly structured" + assert residual_step < 0.05 * raw_step + + def test_scmos_preset_loads(): cfg = load_preset("hamamatsu_orca_fusion") assert cfg.sensor_type is SensorType.SCMOS assert cfg.read_noise_nonuniformity > 0 +def test_ttf_characterised_presets_carry_measured_structure(): + """The three presets cross-validated against real darks keep their measured terms.""" + for name in ( + "princeton_instruments_kuro_1200b", + "photometrics_prime_95b", + "andor_marana_4_2b_11", + ): + cfg = load_preset(name) + assert cfg.read_noise_rts_fraction > 0, name + assert cfg.read_noise_rts_factor > 1.0, name + # DSNU was an order of magnitude too low before characterisation. + assert cfg.dark_current_nonuniformity > 0.1, name + assert 0.7 < cfg.gain_e_per_adu < 0.9, name + marana = load_preset("andor_marana_4_2b_11") + assert marana.detector_glow_edge_scale_px > 0 + + def test_realism_defaults_are_off(): cfg = base_config() assert cfg.nonlinearity == 0.0 From 9faccafd11b201bf7f334573660c719cf50861ed Mon Sep 17 00:00:00 2001 From: Jacob Taylor Date: Wed, 29 Jul 2026 12:51:55 -1000 Subject: [PATCH 2/4] Raise DSNU on the remaining sCMOS presets to a characterised default Every sCMOS preset carried dark_current_nonuniformity of 0.02-0.03, or omitted it entirely. Measured against real dark stacks, three back-illuminated sCMOS cameras came out at 0.11 (Marana 4.2B-11), 0.23 (KURO 1200B) and 0.33 (Prime 95B) -- roughly an order of magnitude higher, consistently. Set the five uncharacterised sCMOS presets to 0.23, the median of those three. Each preset carries a banner saying the value is carried over from characterised hardware rather than taken from that camera's datasheet, so nobody mistakes it for a specification. Note this is an inference across sensors, not a measurement of these cameras; the three measured devices all use the same 11 um back-illuminated family, so smaller-pixel parts (ORCA-Quest 2 at 4.6 um, ORCA-Fusion at 6.5 um) may differ. It is still a much better default than a value known to be ~10x low. Also guard the class of error with a test asserting every sCMOS preset has a DSNU of at least 0.1. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 10 +++++++-- .../presets/data/andor_cb1_0_5mp.toml | 8 +++++++ src/getframes/presets/data/generic_scmos.toml | 2 +- .../presets/data/hamamatsu_orca_fusion.toml | 2 +- .../presets/data/hamamatsu_orca_quest_2.toml | 8 +++++++ .../presets/data/tucsen_aries_6504_pro.toml | 9 +++++++- tests/test_realism.py | 22 ++++++++++++++++++- 7 files changed, 55 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b83f34b..d267957 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,8 +45,14 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). largest corrections: conversion gain (1.25-1.3 -> 0.77-0.87 e-/ADU, the low-signal leg of these dual-gain modes) and `dark_current_nonuniformity` (0.03 -> 0.11-0.33, which had been roughly an order of magnitude too low). Each preset documents the - operating mode and temperature the values apply to. Other sCMOS presets still - carry datasheet-derived DSNU and are likely low for the same reason. + operating mode and temperature the values apply to. +- `dark_current_nonuniformity` raised to `0.23` on the remaining sCMOS presets + (`generic_scmos`, `hamamatsu_orca_fusion`, `hamamatsu_orca_quest_2`, + `tucsen_aries_6504_pro`, `andor_cb1_0_5mp`), which previously carried 0.02-0.03 or + omitted the field entirely. `0.23` is the median of the three cameras measured + against real dark stacks (0.11, 0.23, 0.33); each preset documents that it is a + realistic default carried over from characterised hardware rather than a figure + from that camera's datasheet. - Full-detector region-of-interest simulation through `CameraConfig.roi=(left, top, width, height)`. Cameras accept and return diff --git a/src/getframes/presets/data/andor_cb1_0_5mp.toml b/src/getframes/presets/data/andor_cb1_0_5mp.toml index 9fe03d4..8754b4a 100644 --- a/src/getframes/presets/data/andor_cb1_0_5mp.toml +++ b/src/getframes/presets/data/andor_cb1_0_5mp.toml @@ -16,9 +16,17 @@ read_noise_e = 2.35 dark_current_e_per_s = 1.39 dark_current_ref_temp_c = 10.0 dark_current_doubling_temp_c = 6.0 +dark_current_nonuniformity = 0.23 supported_binnings = [1] # no 2x2 binning offered for the CB1 0.5 MP notes = "CB1 0.5 MP high-gain / 24 dB mode: 2.35 e- typical read noise at 50 us. The specification sheet gives 1594 fps at 8 bit and 941 fps at 12 bit; the profile uses the 12-bit trade mode. QE is the full published CB1 spectral-response plot for the visible Sony Pregius family, with a 73% IMX426 check at 529 nm from Allied Vision." +# Dark-signal non-uniformity carried over from characterised hardware, not from this +# camera's datasheet. Three back-illuminated sCMOS cameras measured against real dark +# stacks (KURO 1200B 0.23, Prime 95B 0.33, Marana 4.2B-11 0.11) all sat roughly an +# order of magnitude above the ~0.03 that had been assumed here; 0.23 is their median. +# Treat it as a realistic default rather than a specification --- re-characterise if +# DSNU matters quantitatively for your work. + [extra] body_dimensions_mm = [154.3, 76.2, 64.1] body_dimensions_order = "length, width, height; confirm mounting/cable clearance in the Keck envelope" diff --git a/src/getframes/presets/data/generic_scmos.toml b/src/getframes/presets/data/generic_scmos.toml index 183cf98..78b4b4b 100644 --- a/src/getframes/presets/data/generic_scmos.toml +++ b/src/getframes/presets/data/generic_scmos.toml @@ -15,7 +15,7 @@ dark_current_ref_temp_c = -10.0 dark_current_doubling_temp_c = 6.0 prnu = 0.01 nonlinearity = 0.01 -dark_current_nonuniformity = 0.03 +dark_current_nonuniformity = 0.23 hot_pixel_fraction = 0.001 hot_pixel_factor = 130.0 notes = "Idealised sCMOS demonstrating per-pixel read noise and mild nonlinearity." diff --git a/src/getframes/presets/data/hamamatsu_orca_fusion.toml b/src/getframes/presets/data/hamamatsu_orca_fusion.toml index 80d525a..0d36eb5 100644 --- a/src/getframes/presets/data/hamamatsu_orca_fusion.toml +++ b/src/getframes/presets/data/hamamatsu_orca_fusion.toml @@ -19,7 +19,7 @@ dark_current_ref_temp_c = -10.0 dark_current_doubling_temp_c = 6.0 prnu = 0.01 nonlinearity = 0.01 -dark_current_nonuniformity = 0.02 +dark_current_nonuniformity = 0.23 hot_pixel_fraction = 0.0005 hot_pixel_factor = 120.0 notes = "Back-thinned sCMOS; very high QE and low median read noise with a per-pixel spread." diff --git a/src/getframes/presets/data/hamamatsu_orca_quest_2.toml b/src/getframes/presets/data/hamamatsu_orca_quest_2.toml index 53b161c..df702bb 100644 --- a/src/getframes/presets/data/hamamatsu_orca_quest_2.toml +++ b/src/getframes/presets/data/hamamatsu_orca_quest_2.toml @@ -16,10 +16,18 @@ read_noise_e = 0.43 dark_current_e_per_s = 0.016 dark_current_ref_temp_c = -20.0 dark_current_doubling_temp_c = 6.0 +dark_current_nonuniformity = 0.23 supported_binnings = [1, 2, 4] # digital, applied after per-pixel photon-number resolution binning_method = "digital" notes = "Standard-scan mode: 0.43 e- RMS at 120 fps full frame. Maximum-cooling dark current is 0.006 e-/pixel/s at -35 C. The manual specifies that 2x2 and 4x4 are digital binning after per-pixel photon-number resolution, so their read noise is propagated from native pixels rather than treated as a single read. Ultra-quiet mode (0.30 e- at 25.4 fps) is an alternate read mode. QE is a typical curve digitized from the catalog figure." +# Dark-signal non-uniformity carried over from characterised hardware, not from this +# camera's datasheet. Three back-illuminated sCMOS cameras measured against real dark +# stacks (KURO 1200B 0.23, Prime 95B 0.33, Marana 4.2B-11 0.11) all sat roughly an +# order of magnitude above the ~0.03 that had been assumed here; 0.23 is their median. +# Treat it as a realistic default rather than a specification --- re-characterise if +# DSNU matters quantitatively for your work. + [extra] source_modes_url = "https://camera.hamamatsu.com/content/dam/hamamatsu-photonics/sites/static/sys/en/manual/C15550-22UP_IM_En.pdf" diff --git a/src/getframes/presets/data/tucsen_aries_6504_pro.toml b/src/getframes/presets/data/tucsen_aries_6504_pro.toml index 013ea00..cf52004 100644 --- a/src/getframes/presets/data/tucsen_aries_6504_pro.toml +++ b/src/getframes/presets/data/tucsen_aries_6504_pro.toml @@ -18,11 +18,18 @@ dark_current_e_per_s = 0.01 dark_current_ref_temp_c = -20.0 dark_current_doubling_temp_c = 6.0 prnu = 0.003 # published 0.3% -dark_current_nonuniformity = 0.03 +dark_current_nonuniformity = 0.23 supported_binnings = [1, 2, 4] # Mosaic sum binning on original pixel data (digital post-read) binning_method = "digital" notes = "Sensitive-mode read noise with HDR-mode full well; QE x fill-factor curve digitized from Gpixel's published GSENSE6504BSI plot. Speed mode (0.80 e-) is an alternate read mode, offered at 1x1 only here. Mosaic documents binning as sum on original pixel data with unchanged data volume/frame rate; modeled as digital post-read summation. Confirm Aries-specific SDK behavior and ROI timing before final design." +# Dark-signal non-uniformity carried over from characterised hardware, not from this +# camera's datasheet. Three back-illuminated sCMOS cameras measured against real dark +# stacks (KURO 1200B 0.23, Prime 95B 0.33, Marana 4.2B-11 0.11) all sat roughly an +# order of magnitude above the ~0.03 that had been assumed here; 0.23 is their median. +# Treat it as a realistic default rather than a specification --- re-characterise if +# DSNU matters quantitatively for your work. + [extra] body_dimensions_mm = [80.0, 80.0, 108.0] body_dimensions_order = "width, height, depth; bare camera body only, excluding connectors and adapters" diff --git a/tests/test_realism.py b/tests/test_realism.py index a3f7f13..e9ba148 100644 --- a/tests/test_realism.py +++ b/tests/test_realism.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from getframes import Camera, CameraConfig, SensorType, load_preset, noise +from getframes import Camera, CameraConfig, SensorType, available_presets, load_preset, noise def base_config(**overrides): @@ -262,6 +262,26 @@ def test_ttf_characterised_presets_carry_measured_structure(): assert marana.detector_glow_edge_scale_px > 0 +def test_every_scmos_preset_has_a_realistic_dsnu(): + """Guard against the ~0.03 DSNU that every sCMOS preset used to carry. + + Measured against real dark stacks, three back-illuminated sCMOS cameras came + out at 0.11-0.33 --- roughly an order of magnitude above the datasheet-derived + value the presets had. Uncharacterised sCMOS presets now carry the median of + those three as a realistic default. + """ + checked = 0 + for name in available_presets(): + cfg = load_preset(name) + if cfg.sensor_type is not SensorType.SCMOS: + continue + checked += 1 + assert cfg.dark_current_nonuniformity >= 0.1, ( + f"{name}: DSNU {cfg.dark_current_nonuniformity} is implausibly low for sCMOS" + ) + assert checked >= 5 + + def test_realism_defaults_are_off(): cfg = base_config() assert cfg.nonlinearity == 0.0 From c6571221367946bf283ee61f64e47034554fe914 Mon Sep 17 00:00:00 2001 From: Jacob Taylor Date: Wed, 29 Jul 2026 13:14:30 -1000 Subject: [PATCH 3/4] Propagate sCMOS structure to remaining presets; document dark-only PTC Follow-on to the DSNU work, applying the same reasoning to the terms I had flagged but left alone. Presets: - The four conventional sCMOS presets (generic_scmos, hamamatsu_orca_fusion, tucsen_aries_6504_pro, andor_cb1_0_5mp) gain the measured RTS population, read_noise_rts_fraction = 0.016 / factor 2.65. All three characterised sensors put 0.4-0.5% of pixels above 3x the median read noise where a bare log-normal predicts ~0.01%, and those pixels set the faint-source detection floor, so omitting them makes a simulated sensor optimistic at threshold. - andor_cb1_0_5mp and hamamatsu_orca_quest_2 had no read_noise_nonuniformity at all, i.e. a perfectly uniform read noise, which is the one thing an sCMOS is not. Both now carry 0.2, marked as a generic default. - hamamatsu_orca_quest_2 deliberately gets NO RTS population, and says why in the file: photon-number resolution depends on a tightly screened read-noise distribution, so importing a tail measured on conventional 11 um back-illuminated sCMOS would misrepresent it. - andor_marana_4_2b_11 gains its measured hot-pixel population, 1.0e-4 of pixels above 10x the median dark rate. Docs and tests: - docs/guides/validation.md gains a section on validating a preset against your own detector: how to measure conversion gain from darks alone (no flats needed -- dark charge is Poisson, so it works as the PTC charge source), and the split-half test for repeatable per-pixel read noise. - test_dark_ptc_recovers_gain_without_any_illumination pins that estimator against a camera whose gain is known, so the method used to characterise the three presets from real hardware is itself covered in CI. It also checks the recovered gain leaves the electron statistics Poisson, which is the assumption the whole technique rests on. Re-validated against the real dark stacks: unchanged at 2.9% / 6.1% / 6.1% median variance error. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 28 +++++--- docs/guides/validation.md | 68 +++++++++++++++++++ .../presets/data/andor_cb1_0_5mp.toml | 7 ++ .../presets/data/andor_marana_4_2b_11.toml | 2 + src/getframes/presets/data/generic_scmos.toml | 2 + .../presets/data/hamamatsu_orca_fusion.toml | 2 + .../presets/data/hamamatsu_orca_quest_2.toml | 6 ++ .../presets/data/tucsen_aries_6504_pro.toml | 6 ++ tests/test_validation.py | 58 ++++++++++++++++ 9 files changed, 170 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d267957..673d230 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,13 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Added +- Full-detector region-of-interest simulation through + `CameraConfig.roi=(left, top, width, height)`. Cameras accept and return + ROI-shaped arrays while evaluating detector physics and fixed patterns on the + native sensor before cropping. `Camera.sensor_resolution`, + `CameraConfig.output_resolution`, and active amplifier-boundary properties make + the full-versus-ROI geometry explicit. Exact full-detector split pixels remain + available when an ROI is active. - `read_noise_rts_fraction` / `read_noise_rts_factor`: an optional second, noisier read-noise population modelling the random-telegraph-signal (RTS) pixels of a real sCMOS array. Measured on three real sensors, ~0.5% of pixels sit above @@ -52,15 +59,18 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). omitted the field entirely. `0.23` is the median of the three cameras measured against real dark stacks (0.11, 0.23, 0.33); each preset documents that it is a realistic default carried over from characterised hardware rather than a figure - from that camera's datasheet. - -- Full-detector region-of-interest simulation through - `CameraConfig.roi=(left, top, width, height)`. Cameras accept and return - ROI-shaped arrays while evaluating detector physics and fixed patterns on the - native sensor before cropping. `Camera.sensor_resolution`, - `CameraConfig.output_resolution`, and active amplifier-boundary properties make - the full-versus-ROI geometry explicit. Exact full-detector split pixels remain - available when an ROI is active. + from that camera's datasheet. The same four conventional sCMOS presets also gain + the measured RTS population (`read_noise_rts_fraction = 0.016`, factor 2.65), and + `andor_cb1_0_5mp` / `hamamatsu_orca_quest_2` gain a `read_noise_nonuniformity` of + 0.2 where they previously had none at all. `hamamatsu_orca_quest_2` deliberately + keeps no RTS population --- photon-number resolution depends on a tightly screened + read-noise distribution, and importing a tail measured on conventional 11 um sCMOS + would misrepresent it. +- `andor_marana_4_2b_11` gains its measured hot-pixel population + (`hot_pixel_fraction = 1e-4` above 10x the median dark rate). +- `docs/guides/validation.md` documents how to validate a preset against a real dark + stack: measuring conversion gain from darks alone (no flats needed), and the + split-half test for repeatable per-pixel read noise. ## [2.1.1] - 2026-07-26 diff --git a/docs/guides/validation.md b/docs/guides/validation.md index 49b0e46..b406c34 100644 --- a/docs/guides/validation.md +++ b/docs/guides/validation.md @@ -19,6 +19,8 @@ and run as part of the test gate. | Synthetic PTC | recovers configured gain / read noise / full well | within 5–15% | | Reduced frame | recovers `Frame.truth` to the noise floor | mean residual < 2 ADU | | All new (1.6) paths | deterministic for a fixed seed | bit-exact | +| Dark-only PTC | recovers configured gain / dark rate with no illumination | within 5–10% | +| sCMOS per-pixel read noise | repeatable through time (real detectors: split-half $r$ = 0.89–0.94) | $r > 0.8$ | ## Recover the gain from a photon transfer curve @@ -67,6 +69,72 @@ recovered_F = np.sqrt(1.0 + out.var() / (60.0 * 250.0**2)) print(recovered_F) # ~1.414 ``` +## Validating against your own detector + +The checks above are internal or analytic. The strongest test is a real dark stack. +Three presets (`princeton_instruments_kuro_1200b`, `photometrics_prime_95b`, +`andor_marana_4_2b_11`) carry values characterised this way against real hardware — +that characterisation is not re-run in CI (it needs the raw frames), but the +*estimator* below is pinned against known truth by +`test_dark_ptc_recovers_gain_without_any_illumination`. Here is the method, which +needs nothing but darks. + +### Measure the conversion gain from darks alone + +You do not need flats. Dark current is itself a Poisson process, so thermally +generated charge works as the charge source for a photon transfer curve. For a dark +frame, + +$$\text{mean}_\text{ADU}(t) = \text{bias} + \frac{Dt}{g}, \qquad + \text{var}_\text{ADU}(t) = \text{RN}_\text{ADU}^2 + \frac{Dt}{g^2}$$ + +so $\mathrm{d}\,\text{var}/\mathrm{d}\,\text{mean} = 1/g$ — the dark rate $D$ cancels. +Working per pixel makes it immune to DSNU, and taking a *slope* across exposures +removes the bias pedestal and read noise (they are the two intercepts): + +```python +import numpy as np + +# stacks[t] is an (n_frames, h, w) array of darks at exposure t, in ADU +means = np.stack([s.mean(axis=0) for s in stacks.values()]) # (n_exp, h, w) +variances = np.stack([s.var(axis=0, ddof=1) for s in stacks.values()]) + + +def slope(x, y): # least squares along axis 0, per pixel + xm, ym = x.mean(axis=0), y.mean(axis=0) + return ((x - xm) * (y - ym)).sum(axis=0) / ((x - xm) ** 2).sum(axis=0) + + +gain = float(np.nanmedian(1.0 / slope(means, variances))) # e-/ADU +``` + +The load-bearing assumption is that the dark charge is Poisson (Fano factor 1). Check +it by confirming the recovered gain makes the electron statistics self-consistent: +$\text{var}_e / \text{mean}_e$ should come out at 1. A wrong gain shows up as a Fano +factor visibly away from unity. + +### Check that per-pixel read noise repeats + +sCMOS read noise is a property of each pixel's own amplifier and column ADC, so a +pixel's noise *through time* is repeatable. Split a dark stack into two halves, +compute each half's per-pixel temporal variance, and correlate: + +```python +a = stack[0::2].var(axis=0, ddof=1) +b = stack[1::2].var(axis=0, ddof=1) +print(np.corrcoef(a.ravel(), b.ravel())[0, 1]) +``` + +Real back-illuminated sCMOS gives $r$ = 0.89–0.94; a simulator that re-draws its +per-pixel sigma each frame gives $r \approx 0$. Run the same code against +`Camera.dark_series` and the two should agree. (This check is what caught a real bug +in `getframes`; it is now `test_per_pixel_read_noise_is_a_fixed_sensor_property`.) + +The same split-half machinery separates *fixed* detector structure from sampling +noise generally: the spatial variance of a variance map is +$V_\text{fixed} + 2\langle v\rangle^2/(n-1)$, so anything left after subtracting the +$\chi^2$ term is real structure. + ## Reproducibility Every generation path is seeded through a `numpy.random.Generator`; a given config diff --git a/src/getframes/presets/data/andor_cb1_0_5mp.toml b/src/getframes/presets/data/andor_cb1_0_5mp.toml index 8754b4a..7aa1a5b 100644 --- a/src/getframes/presets/data/andor_cb1_0_5mp.toml +++ b/src/getframes/presets/data/andor_cb1_0_5mp.toml @@ -13,10 +13,13 @@ bit_depth = 12 gain_e_per_adu = 1.0 bias_offset_adu = 100.0 read_noise_e = 2.35 +read_noise_nonuniformity = 0.2 # generic sCMOS per-pixel spread; not characterised for this camera dark_current_e_per_s = 1.39 dark_current_ref_temp_c = 10.0 dark_current_doubling_temp_c = 6.0 dark_current_nonuniformity = 0.23 +read_noise_rts_fraction = 0.016 # RTS tail carried over from characterised hardware (see banner) +read_noise_rts_factor = 2.65 supported_binnings = [1] # no 2x2 binning offered for the CB1 0.5 MP notes = "CB1 0.5 MP high-gain / 24 dB mode: 2.35 e- typical read noise at 50 us. The specification sheet gives 1594 fps at 8 bit and 941 fps at 12 bit; the profile uses the 12-bit trade mode. QE is the full published CB1 spectral-response plot for the visible Sony Pregius family, with a 73% IMX426 check at 529 nm from Allied Vision." @@ -26,6 +29,10 @@ notes = "CB1 0.5 MP high-gain / 24 dB mode: 2.35 e- typical read noise at 50 us. # order of magnitude above the ~0.03 that had been assumed here; 0.23 is their median. # Treat it as a realistic default rather than a specification --- re-characterise if # DSNU matters quantitatively for your work. +# The RTS terms are likewise carried over: all three measured sensors put 0.4-0.5% of +# pixels above 3x the median read noise, where a bare log-normal predicts ~0.01%. +# These are the pixels that limit faint-source detection, so leaving them out makes a +# simulated sensor optimistic at threshold. [extra] body_dimensions_mm = [154.3, 76.2, 64.1] diff --git a/src/getframes/presets/data/andor_marana_4_2b_11.toml b/src/getframes/presets/data/andor_marana_4_2b_11.toml index 4a6c552..1e64c29 100644 --- a/src/getframes/presets/data/andor_marana_4_2b_11.toml +++ b/src/getframes/presets/data/andor_marana_4_2b_11.toml @@ -20,6 +20,8 @@ dark_current_doubling_temp_c = 6.0 prnu = 0.005 # published <0.5% at half range nonlinearity = 0.003 # published linearity >99.7% dark_current_nonuniformity = 0.110 # measured on the glow-free interior +hot_pixel_fraction = 0.0001 # measured: 1.0e-4 of pixels above 10x the median dark rate +hot_pixel_factor = 15.0 # consistent with that threshold; tail not resolved supported_binnings = [1, 2, 3, 4, 8] # FPGA post-readout (digital) summation binning_method = "digital" notes = "Values below the measured-values banner, plus gain/read noise/dark current/bias, were characterised from real 16-bit HDR darks at -45 C (see the banner). This sensor shows edge-concentrated amplifier glow, modelled via detector_glow_edge_scale_px; the dark_current_e_per_s figure is the glow-subtracted base rate. Published QE curve retained. FPGA post-readout binning does not reduce the sCMOS readout time." diff --git a/src/getframes/presets/data/generic_scmos.toml b/src/getframes/presets/data/generic_scmos.toml index 78b4b4b..93b7578 100644 --- a/src/getframes/presets/data/generic_scmos.toml +++ b/src/getframes/presets/data/generic_scmos.toml @@ -16,6 +16,8 @@ dark_current_doubling_temp_c = 6.0 prnu = 0.01 nonlinearity = 0.01 dark_current_nonuniformity = 0.23 +read_noise_rts_fraction = 0.016 # RTS tail carried over from characterised hardware (see banner) +read_noise_rts_factor = 2.65 hot_pixel_fraction = 0.001 hot_pixel_factor = 130.0 notes = "Idealised sCMOS demonstrating per-pixel read noise and mild nonlinearity." diff --git a/src/getframes/presets/data/hamamatsu_orca_fusion.toml b/src/getframes/presets/data/hamamatsu_orca_fusion.toml index 0d36eb5..613250e 100644 --- a/src/getframes/presets/data/hamamatsu_orca_fusion.toml +++ b/src/getframes/presets/data/hamamatsu_orca_fusion.toml @@ -20,6 +20,8 @@ dark_current_doubling_temp_c = 6.0 prnu = 0.01 nonlinearity = 0.01 dark_current_nonuniformity = 0.23 +read_noise_rts_fraction = 0.016 # RTS tail carried over from characterised hardware (see banner) +read_noise_rts_factor = 2.65 hot_pixel_fraction = 0.0005 hot_pixel_factor = 120.0 notes = "Back-thinned sCMOS; very high QE and low median read noise with a per-pixel spread." diff --git a/src/getframes/presets/data/hamamatsu_orca_quest_2.toml b/src/getframes/presets/data/hamamatsu_orca_quest_2.toml index df702bb..7198d79 100644 --- a/src/getframes/presets/data/hamamatsu_orca_quest_2.toml +++ b/src/getframes/presets/data/hamamatsu_orca_quest_2.toml @@ -13,6 +13,7 @@ bit_depth = 16 gain_e_per_adu = 0.107 bias_offset_adu = 200.0 read_noise_e = 0.43 +read_noise_nonuniformity = 0.2 # generic sCMOS per-pixel spread; not characterised for this camera dark_current_e_per_s = 0.016 dark_current_ref_temp_c = -20.0 dark_current_doubling_temp_c = 6.0 @@ -28,6 +29,11 @@ notes = "Standard-scan mode: 0.43 e- RMS at 120 fps full frame. Maximum-cooling # Treat it as a realistic default rather than a specification --- re-characterise if # DSNU matters quantitatively for your work. +# No read_noise_rts_fraction is set here, unlike the other sCMOS presets. Photon-number +# resolution is this camera's defining capability and depends on a tightly screened +# read-noise distribution, so importing an RTS tail measured on conventional 11 um +# back-illuminated sCMOS would misrepresent it. Characterise before adding one. + [extra] source_modes_url = "https://camera.hamamatsu.com/content/dam/hamamatsu-photonics/sites/static/sys/en/manual/C15550-22UP_IM_En.pdf" diff --git a/src/getframes/presets/data/tucsen_aries_6504_pro.toml b/src/getframes/presets/data/tucsen_aries_6504_pro.toml index cf52004..7e4dd3a 100644 --- a/src/getframes/presets/data/tucsen_aries_6504_pro.toml +++ b/src/getframes/presets/data/tucsen_aries_6504_pro.toml @@ -19,6 +19,8 @@ dark_current_ref_temp_c = -20.0 dark_current_doubling_temp_c = 6.0 prnu = 0.003 # published 0.3% dark_current_nonuniformity = 0.23 +read_noise_rts_fraction = 0.016 # RTS tail carried over from characterised hardware (see banner) +read_noise_rts_factor = 2.65 supported_binnings = [1, 2, 4] # Mosaic sum binning on original pixel data (digital post-read) binning_method = "digital" notes = "Sensitive-mode read noise with HDR-mode full well; QE x fill-factor curve digitized from Gpixel's published GSENSE6504BSI plot. Speed mode (0.80 e-) is an alternate read mode, offered at 1x1 only here. Mosaic documents binning as sum on original pixel data with unchanged data volume/frame rate; modeled as digital post-read summation. Confirm Aries-specific SDK behavior and ROI timing before final design." @@ -29,6 +31,10 @@ notes = "Sensitive-mode read noise with HDR-mode full well; QE x fill-factor cur # order of magnitude above the ~0.03 that had been assumed here; 0.23 is their median. # Treat it as a realistic default rather than a specification --- re-characterise if # DSNU matters quantitatively for your work. +# The RTS terms are likewise carried over: all three measured sensors put 0.4-0.5% of +# pixels above 3x the median read noise, where a bare log-normal predicts ~0.01%. +# These are the pixels that limit faint-source detection, so leaving them out makes a +# simulated sensor optimistic at threshold. [extra] body_dimensions_mm = [80.0, 80.0, 108.0] diff --git a/tests/test_validation.py b/tests/test_validation.py index 8527a86..a5412dd 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -145,6 +145,64 @@ def test_ptc_recovers_gain_read_noise_full_well(): assert ptc.full_well_adu * ptc.gain_e_per_adu == pytest.approx(60_000.0, rel=0.15) +# --------------------------------------------------------------------------- +# A dark-only photon transfer curve recovers the configured gain +# --------------------------------------------------------------------------- +def test_dark_ptc_recovers_gain_without_any_illumination(): + """The method documented in docs/guides/validation.md, run against known truth. + + Dark charge is Poisson, so it serves as the charge source for a photon transfer + curve: per pixel, the slope of temporal variance against temporal mean is 1/gain, + with the dark rate cancelling and bias/read noise absorbed into the intercepts. + This is how the KURO/Prime 95B/Marana presets were characterised from real darks, + so pin it against a camera whose gain we know. + """ + gain, dark = 0.85, 4.0 + config = gf.CameraConfig( + name="dark-ptc-validation", + sensor_type="SCMOS", + resolution=(64, 64), + pixel_size_um=11.0, + quantum_efficiency=1.0, + full_well_e=60_000.0, + bit_depth=16, + gain_e_per_adu=gain, + bias_offset_adu=100.0, + read_noise_e=1.6, + read_noise_nonuniformity=0.25, + dark_current_e_per_s=dark, + dark_current_ref_temp_c=-20.0, + dark_current_nonuniformity=0.2, + ) + cam = gf.Camera(config, default_temperature_c=-20.0) + exposures = [0.5, 1.0, 2.0, 4.0, 8.0, 16.0] + + means, variances = [], [] + for t in exposures: + stack = np.stack([np.asarray(f, dtype=np.float64) for f in cam.dark_series(t, 400, seed=7)]) + means.append(stack.mean(axis=0)) + variances.append(stack.var(axis=0, ddof=1)) + means, variances = np.stack(means), np.stack(variances) + + def slope(x, y): # per-pixel least squares along the exposure axis + xm, ym = x.mean(axis=0), y.mean(axis=0) + return ((x - xm) * (y - ym)).sum(axis=0) / ((x - xm) ** 2).sum(axis=0) + + recovered_gain = float(np.nanmedian(1.0 / slope(means, variances))) + assert recovered_gain == pytest.approx(gain, rel=0.05) + + # And the dark rate follows once the gain is known. + t_axis = np.array(exposures)[:, None, None] + recovered_dark = float(np.nanmedian(slope(t_axis, means))) * recovered_gain + assert recovered_dark == pytest.approx(dark, rel=0.10) + + # The electron statistics implied by that gain are Poisson (Fano factor 1) -- + # the assumption the whole method rests on. + d_mean = float(np.median(means[-1] - means[0])) * recovered_gain + d_var = float(np.median(variances[-1] - variances[0])) * recovered_gain**2 + assert d_var / d_mean == pytest.approx(1.0, abs=0.1) + + # --------------------------------------------------------------------------- # A reduced frame recovers the ground truth to the noise floor # --------------------------------------------------------------------------- From 285fb53a97665fd1b73035493ed4e9b6b2f811af Mon Sep 17 00:00:00 2001 From: Jacob Taylor Date: Wed, 29 Jul 2026 13:56:15 -1000 Subject: [PATCH 4/4] Add a detector-characterisation backend that runs on real frames `analysis/ptc.py` characterises a *simulated* camera by driving it. The new `analysis/characterize.py` works the other way round: give it stacks of frames that already exist -- raw data off a real detector, or output from Camera -- and it measures the detector parameters back out. frames -> stack_statistics per-pixel temporal mean and variance -> characterize_dark gain, read noise, dark current, bias, DSNU -> to_config a CameraConfig -> Camera synthetic frames matching your detector `stack_statistics` reduces any iterable of frames -- arrays, Frames, a 3-D cube, a dark_series generator, your own file reader -- through a Welford accumulator, so stacks far larger than memory stream fine. `characterize_dark` measures gain, read noise (with its per-pixel map, log-normal width and RTS tail), dark current, bias and DSNU from darks alone: no flat field is needed, because dark charge is Poisson and so serves as the PTC charge source. `characterize_flat` adds full well, PRNU and linearity. `DarkCharacterization.to_config()` returns a CameraConfig, which closes the loop the library was built for -- measure a real camera, then simulate it. Two estimator choices worth recording, both settled by measurement rather than by which is more standard: - Read noise comes from the shortest stack with its dark term subtracted, not from the variance regression extrapolated to zero exposure. Both are unbiased in the median, but the extrapolation carries every pixel's fit error into the read-noise map: on a known camera it returned a log-normal width of 0.34 against a true 0.25, where this form returns 0.26. - Flat read noise likewise comes from the faintest stack with shot noise removed rather than the PTC intercept: 5.07 e- against a true 5.0, where the intercept gave 7.4. `temporal_repeatability` (the split-half test) clips the most extreme 1% of pixels. Running the backend against the real TTF frames is what forced this: a cosmic ray lands in one half only and inflates that pixel's variance by orders of magnitude, so a handful of them dominate the covariance. Real 60 s Marana darks score 0.006 unclipped against 0.96 clipped -- a user would have concluded the detector has no fixed read-noise structure, which is false. Validated both ways. Against simulated cameras with known parameters, every parameter recovers within 3-15% (tests/test_characterize.py). Against the real KURO/Prime 95B/Marana dark stacks, it reproduces the bespoke per-pixel analysis those presets were built from to 0.0% on gain, dark current, bias and DSNU. Ships with docs/guides/characterization.md and examples/15_detector_characterization.py. Co-Authored-By: Claude Opus 5 --- AGENTS.md | 2 +- CHANGELOG.md | 16 + docs/guides/characterization.md | 194 +++++++ docs/reference.md | 2 + examples/15_detector_characterization.py | 259 +++++++++ examples/README.md | 1 + mkdocs.yml | 1 + src/getframes/analysis/__init__.py | 21 +- src/getframes/analysis/characterize.py | 678 +++++++++++++++++++++++ tests/test_characterize.py | 389 +++++++++++++ 10 files changed, 1561 insertions(+), 2 deletions(-) create mode 100644 docs/guides/characterization.md create mode 100644 examples/15_detector_characterization.py create mode 100644 src/getframes/analysis/characterize.py create mode 100644 tests/test_characterize.py diff --git a/AGENTS.md b/AGENTS.md index 8381aac..a6ddc01 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -42,7 +42,7 @@ removals. A JOSS paper + citation remain a post-2.0 follow-up. | `observation.py` | `Observation` / `ObservationTruth` / `Pointing`: the time-series driver, jitter/drift/dither, per-frame truth (phase 1.2). | | `spectral.py` | Opt-in spectral mode: `QE`, `SED` (relative or absolute via `from_flux_density`), `Spectrum`, `SpectralBandpass`, effective-QE folding, transmission-product helpers (`product`, `from_file`/`from_product`), optional `astropy.units` coercion. | | `scene/` | The scene/optics layer: `Scene`, `Source` hierarchy (`PointSource`, `ExtendedSource`, `UniformIllumination`, `Catalog`; point/extended sources accept a `flux_sed` absolute SED), PSFs (`GaussianPSF`/`MoffatPSF`/`AiryPSF`/`ArrayPSF`/`EllipticalGaussianPSF`), `Telescope` (+ `Vignetting`/`RadialDistortion`), `Bandpass` (Vega `johnson` / AB `ab` ugriz·Gaia·2MASS + `Extinction`) in `photometry.py`, `Thermal` graybody background in `thermal.py`, `WCSInfo`, `LightCurve`. Renders a photon-rate map; no randomness. | -| `analysis/` | Measurement helpers: `apertures.py` (`aperture_sum`, `centroid`), `ptc.py` (`photon_transfer_curve`). | +| `analysis/` | Measurement helpers: `apertures.py` (`aperture_sum`, `centroid`), `ptc.py` (`photon_transfer_curve`, camera-driven), `characterize.py` (`stack_statistics`, `characterize_dark`, `characterize_flat` — stack-driven, so it runs on *real* detector data as well as simulated frames; `DarkCharacterization.to_config()` returns a `CameraConfig`). | | `dataset.py` | Scalable raw+truth dataset generation (phase 1.6): `pairs()` → a streaming `PairDataset` (`to_npz`/`to_arrays`), `random_star_fields()` re-iterable scene source. float32-friendly; no global state. | | `cli.py` | The `getframes` console entry point (phase 1.6): `presets` / `generate` / `dataset` subcommands driven by a TOML config. | | `presets/` | Preset library. TOML data files in `presets/data/`, loaded via `importlib.resources`. `load_preset`, `available_presets`, `preset_info`. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 673d230..fcf21f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,22 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). `CameraConfig.output_resolution`, and active amplifier-boundary properties make the full-versus-ROI geometry explicit. Exact full-detector split pixels remain available when an ROI is active. +- **`getframes.analysis.characterize`: detector characterisation from frame + stacks.** Where `photon_transfer_curve` drives a *simulated* camera, this works + on stacks that already exist -- raw data off a real detector, or simulated + frames. `stack_statistics` reduces any iterable of frames (arrays, `Frame`s, a + `dark_series` generator, your own file reader) to per-pixel temporal mean and + variance in one streaming pass, so stacks larger than memory are fine. + `characterize_dark` then measures conversion gain, read noise (with its + per-pixel map, log-normal width and RTS tail), dark current, bias and DSNU from + darks alone -- no flat field needed, because dark charge is Poisson and so + serves as the PTC charge source. `characterize_flat` adds full well, PRNU and + linearity. `DarkCharacterization.to_config()` returns a `CameraConfig`, closing + the loop: measure a real camera, then simulate it. `StackStats.split=True` + additionally gives `temporal_repeatability`, the split-half test that separates + genuine per-pixel noise structure from chi-squared sampling scatter. + New guide (`docs/guides/characterization.md`) and example + (`examples/15_detector_characterization.py`). - `read_noise_rts_fraction` / `read_noise_rts_factor`: an optional second, noisier read-noise population modelling the random-telegraph-signal (RTS) pixels of a real sCMOS array. Measured on three real sensors, ~0.5% of pixels sit above diff --git a/docs/guides/characterization.md b/docs/guides/characterization.md new file mode 100644 index 0000000..d97e7d3 --- /dev/null +++ b/docs/guides/characterization.md @@ -0,0 +1,194 @@ +# Detector characterisation + +`getframes.analysis.characterize` runs the standard bench measurements on stacks +of frames. It takes plain arrays, so it works equally on **frames from a real +detector** and on frames from a simulated [`Camera`][getframes.camera.Camera] — +and the result carries a `to_config()`, so a real camera can be measured, turned +into a [`CameraConfig`][getframes.config.CameraConfig], and then simulated. + +``` +frames -> stack_statistics per-pixel temporal mean and variance + -> characterize_dark gain, read noise, dark current, bias, DSNU + -> to_config a CameraConfig + -> Camera synthetic frames matching your detector +``` + +This is the complement to +[`photon_transfer_curve`][getframes.analysis.ptc.photon_transfer_curve], which +drives a *simulated* camera to characterise it. Here the frames come first. + +A runnable end-to-end version of everything below is +[`examples/15_detector_characterization.py`](https://github.com/jacotay7/getframes/blob/main/examples/15_detector_characterization.py). + +## Step 1 — reduce each stack + +Everything is built from one quantity: for each pixel, its mean and variance +*through a stack*. `stack_statistics` computes both in a single streaming pass, +so an iterator over a stack far larger than memory is fine. + +```python +from getframes.analysis import stack_statistics + +stats = stack_statistics(frames, exposure_s=2.0, split=True) +stats.mean_adu # (h, w) per-pixel temporal mean, ADU +stats.variance_adu2 # (h, w) per-pixel temporal variance, ADU^2 +``` + +`frames` is any iterable of 2-D frames: NumPy arrays, `Frame` objects, a 3-D +cube, a `Camera.dark_series(...)` generator, or your own reader: + +```python +def read_raw(path, shape=(1200, 1200)): + """Stream a flat little-endian uint16 file, one frame at a time.""" + n_bytes = shape[0] * shape[1] * 2 + with open(path, "rb") as handle: + while chunk := handle.read(n_bytes): + if len(chunk) < n_bytes: + return + yield np.frombuffer(chunk, dtype=" stack_statistics per-pixel temporal mean and variance + -> characterize_dark gain, read noise, dark current, bias, DSNU + -> to_config a CameraConfig you can simulate + -> Camera synthetic frames matching your real detector + +Here the "unknown" detector is itself simulated, so we can print the truth +alongside what the characterisation recovered. Point ``load_stacks`` at your own +data --- any iterable of 2-D arrays, one iterable per exposure time --- and the +rest is unchanged. + +Note the gain is measured from *darks alone*, with no flat field. Dark current is +a Poisson process, so thermally generated charge works as the charge source for a +photon transfer curve: per pixel, the slope of temporal variance against temporal +mean is 1/gain, and the dark rate cancels out. ``fano_factor`` reports the +consistency check on that assumption. + +Run: + python examples/15_detector_characterization.py + python examples/15_detector_characterization.py --plot + python examples/15_detector_characterization.py --save characterization.png +""" + +from __future__ import annotations + +import numpy as np +from _common import PALETTE, build_parser, finish, get_pyplot + +import getframes as gf +from getframes.analysis import characterize_dark, characterize_flat, stack_statistics + +DARK_EXPOSURES = (0.5, 1.0, 2.0, 4.0, 8.0, 16.0) +FLAT_LEVELS = (200.0, 1_000.0, 4_000.0, 10_000.0, 18_000.0, 28_000.0, 38_000.0, 55_000.0) + + +def unknown_detector() -> gf.Camera: + """Stand-in for the camera on your bench. Replace with your own frames.""" + config = gf.CameraConfig( + name="detector under test", + sensor_type="SCMOS", + resolution=(128, 128), + pixel_size_um=11.0, + quantum_efficiency=0.95, + full_well_e=40_000.0, + bit_depth=16, + gain_e_per_adu=0.85, + bias_offset_adu=100.0, + read_noise_e=1.60, + read_noise_nonuniformity=0.25, + read_noise_rts_fraction=0.015, + dark_current_e_per_s=4.0, + dark_current_ref_temp_c=-20.0, + dark_current_nonuniformity=0.20, + prnu=0.015, + ) + return gf.Camera(config, default_temperature_c=-20.0) + + +def load_stacks(camera: gf.Camera, n_frames: int, seed: int) -> dict[float, object]: + """Per-exposure dark statistics. + + Swap the body of this function for your own loader. Anything iterable of 2-D + arrays works, and frames are streamed one at a time, so a generator over a + directory of files is fine even when the whole stack would not fit in memory:: + + def load_stacks(...): + return { + exposure: stack_statistics(read_frames(directory), split=True) + for exposure, directory in my_data.items() + } + """ + return { + exposure: stack_statistics(camera.dark_series(exposure, n_frames, seed=seed), split=True) + for exposure in DARK_EXPOSURES + } + + +def main() -> None: + args = build_parser(__doc__).parse_args() + + camera = unknown_detector() + truth = camera.config + + # ---- 1. Reduce each stack to per-pixel temporal statistics --------------- + darks = load_stacks(camera, n_frames=250, seed=args.seed + 1) + + # ---- 2. Characterise -------------------------------------------------- + result = characterize_dark(darks) + + # read_noise_e is the *scale* of a unit-mean log-normal, so the median + # per-pixel RMS sits a little below it (CameraConfig documents the relation). + true_read_noise_median = truth.read_noise_e * np.exp(-0.5 * truth.read_noise_nonuniformity**2) + + print( + f"Detector: {truth.name} ({truth.resolution[0]}x{truth.resolution[1]}, " + f"{len(DARK_EXPOSURES)} exposures x 250 dark frames)\n" + ) + print(f" {'parameter':30s} {'true':>10} {'measured':>10} {'error':>8}") + print(f" {'-' * 30} {'-' * 10} {'-' * 10} {'-' * 8}") + for label, true_value, measured in [ + ("gain (e-/ADU)", truth.gain_e_per_adu, result.gain_e_per_adu), + ("read noise, median (e-)", true_read_noise_median, result.read_noise_e), + ("dark current (e-/pixel/s)", truth.dark_current_e_per_s, result.dark_current_e_per_s), + ("bias offset (ADU)", truth.bias_offset_adu, result.bias_offset_adu), + ("DSNU", truth.dark_current_nonuniformity, result.dark_current_nonuniformity), + ( + "read-noise non-uniformity", + truth.read_noise_nonuniformity, + result.read_noise_nonuniformity, + ), + ]: + error = 100.0 * (measured / true_value - 1.0) if true_value else float("nan") + print(f" {label:30s} {true_value:10.4f} {measured:10.4f} {error:+7.1f}%") + + print(f"\n Fano factor (should be 1.0): {result.fano_factor:.4f}") + print(" Consistency check on the Poisson assumption the gain rests on.") + print(f" RTS pixels (read noise > 3x median): {result.read_noise_rts_fraction:.4%}") + shortest = darks[DARK_EXPOSURES[0]] + print(f" Split-half repeatability: {shortest.temporal_repeatability:.3f}") + print(" Per-pixel noise that repeats through the stack, i.e. real sCMOS") + print(" structure rather than chi-squared scatter. Real detectors: 0.89-0.94.") + print(f" Fixed fraction of variance map: {shortest.fixed_variance_fraction:.3f}") + + # ---- 3. Rebuild the detector as a config, and simulate it --------------- + rebuilt = result.to_config( + "rebuilt from darks", + pixel_size_um=truth.pixel_size_um, + quantum_efficiency=truth.quantum_efficiency, + full_well_e=truth.full_well_e, + dark_current_ref_temp_c=-20.0, # darks carry no temperature: supply it + ) + twin = gf.Camera(rebuilt, default_temperature_c=-20.0) + recheck = characterize_dark(load_stacks(twin, n_frames=250, seed=args.seed + 2)) + print("\n Round trip: characterise the rebuilt config and compare to the first pass") + print(f" gain {result.gain_e_per_adu:.4f} -> {recheck.gain_e_per_adu:.4f} e-/ADU") + print(f" read noise {result.read_noise_e:.4f} -> {recheck.read_noise_e:.4f} e-") + print( + f" dark current {result.dark_current_e_per_s:.4f} -> " + f"{recheck.dark_current_e_per_s:.4f} e-/pixel/s" + ) + + # ---- 4. Flats add full well, PRNU and linearity ------------------------- + flats = { + level: stack_statistics( + (camera.flat_frame(level, 1.0, seed=args.seed + 500 + 40 * i + k) for k in range(20)), + exposure_s=level, + ) + for i, level in enumerate(FLAT_LEVELS) + } + flat = characterize_flat(flats, bias_adu=result.bias_offset_adu) + print("\n From flats (what darks cannot see):") + print(f" gain {truth.gain_e_per_adu:.4f} -> {flat.gain_e_per_adu:.4f} e-/ADU") + print( + f" full well {truth.full_well_e:.0f} -> " + f"{flat.full_well_e:.0f} e- (variance peak marks saturation onset)" + ) + print(f" PRNU {truth.prnu:.4f} -> {flat.prnu:.4f}") + + # ---- Plotting ----------------------------------------------------------- + plt = get_pyplot(args) + if plt is None: + return + + fig, axes = plt.subplots(2, 2, figsize=(12, 9)) + + # (a) The dark photon transfer curve: variance against mean, slope = 1/gain. + ax = axes[0, 0] + mean_adu = [float(np.median(s.mean_adu)) for s in darks.values()] + var_adu = [float(np.median(s.variance_adu2)) for s in darks.values()] + ax.plot(mean_adu, var_adu, "o", color=PALETTE["blue"], label="dark stacks") + x = np.linspace(min(mean_adu), max(mean_adu), 50) + intercept = var_adu[0] - (1.0 / result.gain_e_per_adu) * mean_adu[0] + ax.plot( + x, + x / result.gain_e_per_adu + intercept, + "-", + color=PALETTE["red"], + lw=2, + label=f"slope = 1/gain -> {result.gain_e_per_adu:.3f} e-/ADU", + ) + ax.set_xlabel("mean signal (ADU)") + ax.set_ylabel("temporal variance (ADU$^2$)") + ax.set_title("(a) Dark photon transfer curve") + ax.legend() + + # (b) The recovered per-pixel read-noise distribution. A single number cannot + # describe an sCMOS: there is a log-normal core plus a heavy RTS tail. + ax = axes[0, 1] + upper = float(np.percentile(result.read_noise_map_e, 99.95)) * 1.15 + ax.hist( + result.read_noise_map_e.ravel(), + bins=np.linspace(0, upper, 70), + color=PALETTE["blue"], + alpha=0.85, + label="measured, per pixel", + ) + ax.axvline( + result.read_noise_e, + color=PALETTE["grey"], + lw=1.5, + label=f"median {result.read_noise_e:.2f} e-", + ) + ax.axvline( + 3 * result.read_noise_e, + color=PALETTE["red"], + ls="--", + lw=1.5, + label=f"3x median: {result.read_noise_rts_fraction:.2%} of pixels (RTS)", + ) + ax.set_xlabel("per-pixel read noise (e-)") + ax.set_ylabel("pixels") + ax.set_title("(b) Read noise is per-pixel, with an RTS tail") + ax.set_yscale("log") + ax.legend() + + # (c) The recovered dark-current map: DSNU as spatial structure. + ax = axes[1, 0] + dark_map = result.dark_current_map_e_per_s + lo, hi = np.percentile(dark_map, [1, 99]) + im = ax.imshow(dark_map, vmin=lo, vmax=hi, origin="lower") + ax.set_title("(c) Recovered dark current (e-/pixel/s)") + ax.set_xticks([]) + ax.set_yticks([]) + fig.colorbar(im, ax=ax, fraction=0.046) + + # (d) The flat-field PTC, all the way through saturation. + ax = axes[1, 1] + ax.loglog(flat.mean_adu, flat.variance_adu2, "o-", color=PALETTE["green"], label="flat stacks") + if flat.full_well_adu is not None: + ax.axvline( + flat.full_well_adu, + color=PALETTE["red"], + ls="--", + lw=1.5, + label=f"saturation onset ({flat.full_well_e:.0f} e-)", + ) + ax.set_xlabel("mean signal above bias (ADU)") + ax.set_ylabel("temporal variance (ADU$^2$)") + ax.set_title("(d) Flat-field PTC: gain, full well, PRNU") + ax.legend() + + fig.suptitle( + "Detector characterisation from frame stacks (truth vs. recovered)", + fontsize=14, + fontweight="bold", + ) + finish(plt, fig, args) + + +if __name__ == "__main__": + main() diff --git a/examples/README.md b/examples/README.md index 1b429dd..05f2aa1 100644 --- a/examples/README.md +++ b/examples/README.md @@ -32,3 +32,4 @@ pip install -e ".[examples]" | [`12_ml_dataset.py`](12_ml_dataset.py) | Stream raw+truth training pairs to disk with `getframes.dataset` (float32). | | [`13_crowded_field.py`](13_crowded_field.py) | Render a 20k-star catalog; vectorised vs. per-source, with flux conservation. | | [`14_keck_lgs_ttf_trade_study.ipynb`](14_keck_lgs_ttf_trade_study.ipynb) | Guided, cell-by-cell trade study of detector candidates for the Keck LGS tip/tilt sensor and low-bandwidth WFS: spot model → radiometry → one simulated frame → calibration → centroiding → the full TTS/LBWFS trades. | +| [`15_detector_characterization.py`](15_detector_characterization.py) | Measure gain, read noise, dark current, DSNU, full well and PRNU from frame stacks (real or simulated), then rebuild the detector as a `CameraConfig`. | diff --git a/mkdocs.yml b/mkdocs.yml index dbd0f45..685ed16 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -37,6 +37,7 @@ nav: - Scale & datasets: guides/datasets.md - GPU detector execution: guides/gpu.md - Camera presets: guides/presets.md + - Detector characterisation: guides/characterization.md - Validation: guides/validation.md - API reference: reference.md - API stability: stability.md diff --git a/src/getframes/analysis/__init__.py b/src/getframes/analysis/__init__.py index 8842656..f64dca6 100644 --- a/src/getframes/analysis/__init__.py +++ b/src/getframes/analysis/__init__.py @@ -1,19 +1,38 @@ # SPDX-License-Identifier: MIT -"""Lightweight analysis helpers (photometry, centroiding, photon transfer curve). +"""Lightweight analysis helpers (photometry, centroiding, detector characterisation). These exist mainly so the bundled examples stay self-contained; they are pure NumPy/SciPy and make no attempt to replace dedicated tools like ``photutils``. + +:mod:`~getframes.analysis.characterize` is the exception in one respect: it takes +frame stacks as plain arrays, so it works equally on real detector data and on +simulated frames --- measuring a real camera and then reproducing it with +:class:`~getframes.Camera` is the intended workflow. """ from __future__ import annotations from .apertures import aperture_sum, centroid, matched_filter_centroid +from .characterize import ( + DarkCharacterization, + FlatCharacterization, + StackStats, + characterize_dark, + characterize_flat, + stack_statistics, +) from .ptc import PTCResult, photon_transfer_curve __all__ = [ + "DarkCharacterization", + "FlatCharacterization", "PTCResult", + "StackStats", "aperture_sum", "centroid", + "characterize_dark", + "characterize_flat", "matched_filter_centroid", "photon_transfer_curve", + "stack_statistics", ] diff --git a/src/getframes/analysis/characterize.py b/src/getframes/analysis/characterize.py new file mode 100644 index 0000000..32dea20 --- /dev/null +++ b/src/getframes/analysis/characterize.py @@ -0,0 +1,678 @@ +# SPDX-License-Identifier: MIT +"""Detector characterisation from frame stacks --- real or simulated. + +:mod:`~getframes.analysis.ptc` characterises a *simulated* camera by driving it. +This module works the other way round: hand it stacks of frames that already +exist --- raw data off a real detector, or output from :class:`~getframes.Camera` +--- and it measures the detector parameters back out. The result carries a +:meth:`DarkCharacterization.to_config` so a real camera can be turned into a +:class:`~getframes.CameraConfig` and then simulated. + +The two entry points mirror the two standard bench measurements: + +``characterize_dark`` + Dark stacks at several exposure times. Returns conversion gain, read noise + (including its per-pixel distribution), dark current, bias offset and DSNU. +``characterize_flat`` + Flat-field stacks at several illumination levels. Returns conversion gain, + read noise, full well, PRNU and linearity. + +Measuring gain from *darks alone* works because dark current is a Poisson +process, so thermally generated charge is a perfectly good charge source for a +photon transfer curve. For a dark frame:: + + mean_ADU(t) = bias + D*t/g + var_ADU(t) = RN_ADU**2 + D*t/g**2 + +so the slope of variance against mean is ``1/g`` and the dark rate ``D`` cancels. +Fitting per pixel makes it immune to DSNU, and fitting a *slope* across exposures +absorbs the bias pedestal and the read noise into the two intercepts. The +assumption this rests on is that the dark charge is Poisson (Fano factor 1); +:attr:`DarkCharacterization.fano_factor` reports the consistency check. + +All inputs are in ADU; all returned electron quantities are in electrons. +""" + +from __future__ import annotations + +from collections.abc import Iterable, Mapping, Sequence +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import NDArray + +if TYPE_CHECKING: + from ..config import CameraConfig + +FrameLike = Any # a 2-D array, a Frame, or anything np.asarray turns into one + + +# --------------------------------------------------------------------------- +# Per-stack temporal statistics +# --------------------------------------------------------------------------- +@dataclass(frozen=True) +class StackStats: + """Per-pixel temporal statistics of one stack of frames, in ADU. + + This is the raw material every characterisation is built from: for each + pixel, its mean and variance *through the stack*. Both are full-resolution + maps, so detector structure (DSNU, per-pixel read noise, hot pixels) is + preserved rather than averaged away. + + Attributes + ---------- + mean_adu, variance_adu2: + Per-pixel temporal mean (ADU) and unbiased variance (ADU^2), each shaped + like one frame. + n_frames: + Number of frames combined. + exposure_s: + Exposure time of the stack in seconds, or ``None`` if unlabelled. + half_variance_adu2: + Per-pixel variance of the even- and odd-indexed frames separately, when + ``split=True`` was passed to :func:`stack_statistics`. Used by + :attr:`temporal_repeatability`. + """ + + mean_adu: NDArray[np.float64] + variance_adu2: NDArray[np.float64] + n_frames: int + exposure_s: float | None = None + half_variance_adu2: tuple[NDArray[np.float64], NDArray[np.float64]] | None = None + + @property + def shape(self) -> tuple[int, ...]: + """Frame shape ``(height, width)``.""" + return self.mean_adu.shape + + @property + def temporal_repeatability(self) -> float: + """Split-half correlation of the per-pixel variance map, in ``[-1, 1]``. + + Splits the stack into even- and odd-indexed frames, computes each half's + per-pixel temporal variance, and correlates the two maps across pixels. + + This separates *fixed* per-pixel noise structure from sampling scatter. A + detector whose pixels genuinely differ in read noise --- every sCMOS --- + gives a high correlation, because the same pixels are noisy in both + halves. A detector with uniform noise gives ~0, because all that differs + between halves is chi-squared sampling noise. Real back-illuminated sCMOS + measures 0.89--0.94. + + The correlation is computed with the most extreme 1% of pixels excluded. + A cosmic ray lands in one half only and inflates that pixel's variance by + orders of magnitude, so on a long-exposure stack a handful of such pixels + dominate the covariance and drive a plain Pearson correlation to zero: + real 60 s Marana darks score 0.006 unclipped against 0.93 clipped. Use + :meth:`repeatability` for explicit control. + + Requires ``split=True`` in :func:`stack_statistics`. + """ + return self.repeatability() + + def repeatability(self, *, clip_percentile: float = 99.0) -> float: + """:attr:`temporal_repeatability` with the outlier cut exposed. + + Parameters + ---------- + clip_percentile: + Pixels whose variance in *either* half exceeds this percentile are + excluded before correlating. ``100`` disables clipping and gives the + plain Pearson correlation. + """ + if self.half_variance_adu2 is None: + raise ValueError( + "temporal_repeatability needs the split halves; " + "call stack_statistics(..., split=True)." + ) + a, b = self.half_variance_adu2 + if clip_percentile >= 100.0: + return float(np.corrcoef(a.ravel(), b.ravel())[0, 1]) + ceiling = np.percentile(np.maximum(a, b), clip_percentile) + keep = (a < ceiling) & (b < ceiling) + if keep.sum() < 3: + return float("nan") + return float(np.corrcoef(a[keep], b[keep])[0, 1]) + + @property + def fixed_variance_fraction(self) -> float: + """Fraction of the variance map's spatial spread that is *fixed* structure. + + The observed spatial variance of a variance map is the real pixel-to-pixel + structure plus the chi-squared scatter of estimating a variance from a + finite stack, ``2 * **2 / (n - 1)``. Subtracting the latter leaves the + fraction that is genuine detector structure, in ``[0, 1]``. + """ + observed = float(self.variance_adu2.var()) + if observed <= 0: + return 0.0 + sampling = 2.0 * float(np.mean(self.variance_adu2**2)) / max(self.n_frames - 1, 1) + return float(np.clip((observed - sampling) / observed, 0.0, 1.0)) + + +def stack_statistics( + frames: Iterable[FrameLike], + *, + exposure_s: float | None = None, + split: bool = False, +) -> StackStats: + """Per-pixel temporal mean and variance of a stack of frames. + + Frames are consumed one at a time through a Welford accumulator, so an + iterator or generator over a stack far larger than memory works fine --- only + a handful of frame-sized float64 arrays are ever held. + + Parameters + ---------- + frames: + Any iterable of 2-D frames: NumPy arrays, :class:`~getframes.Frame` + objects, or a 3-D array (which iterates over its leading axis). A + :meth:`~getframes.Camera.dark_series` generator works directly. + exposure_s: + Exposure time to label the stack with. Required by + :func:`characterize_dark` when stacks are passed as a sequence. + split: + Also accumulate the even- and odd-indexed frames separately, enabling + :attr:`StackStats.temporal_repeatability`. Costs two more frame-sized + accumulators. + + Returns + ------- + StackStats + Per-pixel mean and variance in ADU. + + Raises + ------ + ValueError + If fewer than two frames are supplied (variance is undefined), or if + ``split=True`` and either half has fewer than two frames. + """ + n = 0 + mean: NDArray[np.float64] | None = None + m2: NDArray[np.float64] | None = None + halves: list[list[Any]] = [[0, None, None], [0, None, None]] + + for index, raw in enumerate(frames): + frame = np.asarray(raw, dtype=np.float64) + if mean is None or m2 is None: + mean = np.zeros(frame.shape, dtype=np.float64) + m2 = np.zeros(frame.shape, dtype=np.float64) + elif frame.shape != mean.shape: + raise ValueError(f"frame {index} has shape {frame.shape}, expected {mean.shape}.") + n += 1 + delta = frame - mean + mean += delta / n + m2 += delta * (frame - mean) + + if split: + half = halves[index % 2] + if half[1] is None: + half[1] = np.zeros(frame.shape, dtype=np.float64) + half[2] = np.zeros(frame.shape, dtype=np.float64) + half[0] += 1 + hdelta = frame - half[1] + half[1] += hdelta / half[0] + half[2] += hdelta * (frame - half[1]) + + if mean is None or m2 is None or n < 2: + raise ValueError(f"need at least 2 frames to measure a variance, got {n}.") + + half_var: tuple[NDArray[np.float64], NDArray[np.float64]] | None = None + if split: + if min(halves[0][0], halves[1][0]) < 2: + raise ValueError("split=True needs at least 4 frames (2 per half).") + half_var = ( + halves[0][2] / (halves[0][0] - 1), + halves[1][2] / (halves[1][0] - 1), + ) + return StackStats(mean, m2 / (n - 1), n, exposure_s, half_var) + + +# --------------------------------------------------------------------------- +# Dark characterisation +# --------------------------------------------------------------------------- +@dataclass(frozen=True) +class DarkCharacterization: + """What a set of dark stacks says about a detector. + + Scalars are the median over pixels; the ``*_map`` arrays give the per-pixel + values behind them. + + Attributes + ---------- + gain_e_per_adu: + Conversion gain from the per-pixel dark photon transfer curve. + read_noise_e: + Median per-pixel read noise, taken from the shortest stack with its dark + contribution removed. It therefore includes any exposure-independent + common-mode term (frame-to-frame pedestal wander, for instance), which is + what a bench measurement would also report. Supply a short enough + exposure that read noise dominates it. + dark_current_e_per_s: + Median per-pixel dark current, from the slope of mean against exposure. + bias_offset_adu: + Median pedestal, from the intercept of mean against exposure. + dark_current_nonuniformity: + Robust relative spread (IQR/1.349 over the median) of the per-pixel dark + current --- DSNU, comparable to + :attr:`~getframes.CameraConfig.dark_current_nonuniformity`. + read_noise_nonuniformity: + Log-normal width implied by the read-noise inter-quartile range, + comparable to + :attr:`~getframes.CameraConfig.read_noise_nonuniformity`. + read_noise_rts_fraction: + Fraction of pixels whose read noise exceeds three times the median --- the + random-telegraph-signal tail. Around 0.005 on real sCMOS, against ~1e-4 + for a pure log-normal. + hot_pixel_fraction: + Fraction of pixels whose dark current exceeds ten times the median. + fano_factor: + Consistency check on the Poisson assumption the gain fit relies on: + ``var_e / mean_e`` of the accumulated dark charge, which should be 1. + A value far from 1 means the gain is not trustworthy. + exposures_s: + The exposure times used, ascending. + read_noise_map_e, dark_current_map_e_per_s, bias_map_adu: + Per-pixel maps behind the scalars above. + """ + + gain_e_per_adu: float + read_noise_e: float + dark_current_e_per_s: float + bias_offset_adu: float + dark_current_nonuniformity: float + read_noise_nonuniformity: float + read_noise_rts_fraction: float + hot_pixel_fraction: float + fano_factor: float + exposures_s: NDArray[np.float64] + read_noise_map_e: NDArray[np.float64] + dark_current_map_e_per_s: NDArray[np.float64] + bias_map_adu: NDArray[np.float64] + + def to_config(self, name: str, **overrides: Any) -> CameraConfig: + """Build a :class:`~getframes.CameraConfig` from the measured parameters. + + Everything darks can measure is filled in: resolution, gain, bias, read + noise (with its non-uniformity and RTS tail), dark current and DSNU. + Parameters darks *cannot* see --- full well, bit depth, pixel pitch, QE --- + take documented placeholder defaults that you should override. + + Parameters + ---------- + name: + Name for the resulting config. + **overrides: + Any :class:`~getframes.CameraConfig` field, applied last. Use this to + supply ``pixel_size_um``, ``full_well_e``, ``bit_depth``, + ``quantum_efficiency`` and the sensor type for your detector. + + Notes + ----- + ``dark_current_ref_temp_c`` defaults to 20 C because the stacks carry no + temperature. Set it to the temperature the darks were taken at, or the + config's temperature scaling will be wrong. + """ + from ..config import CameraConfig + + height, width = self.read_noise_map_e.shape + fields: dict[str, Any] = { + "name": name, + "sensor_type": "SCMOS", + "resolution": (int(height), int(width)), + "pixel_size_um": 10.0, + "quantum_efficiency": 1.0, + "full_well_e": 50_000.0, + "bit_depth": 16, + "gain_e_per_adu": self.gain_e_per_adu, + "bias_offset_adu": self.bias_offset_adu, + "read_noise_e": self.read_noise_e, + "read_noise_nonuniformity": self.read_noise_nonuniformity, + "read_noise_rts_fraction": self.read_noise_rts_fraction, + "dark_current_e_per_s": self.dark_current_e_per_s, + "dark_current_nonuniformity": self.dark_current_nonuniformity, + "hot_pixel_fraction": self.hot_pixel_fraction, + } + fields.update(overrides) + return CameraConfig(**fields) + + +def characterize_dark( + stacks: Mapping[float, StackStats] | Sequence[StackStats], +) -> DarkCharacterization: + """Measure a detector from dark stacks at several exposure times. + + Parameters + ---------- + stacks: + Either a mapping of ``exposure_s -> StackStats``, or a sequence of + :class:`StackStats` that each carry their own ``exposure_s``. At least + two distinct exposures are needed; three or more is much better, and the + longest should accumulate enough dark charge to be measurable above the + read noise. + + Returns + ------- + DarkCharacterization + + Raises + ------ + ValueError + If fewer than two distinct exposures are supplied, if the stacks disagree + on frame shape, or if any stack lacks an exposure time. + + Notes + ----- + The gain comes from a *per-pixel* regression of temporal variance against + temporal mean, whose slope is ``1/gain`` regardless of that pixel's own dark + current and read noise. Taking the median over pixels makes it robust to + hot pixels and to the read-noise tail. See the module docstring for why + darks suffice, and check :attr:`DarkCharacterization.fano_factor` before + trusting the result. + """ + ordered = _ordered_stacks(stacks) + exposures = np.array([s.exposure_s for s in ordered], dtype=np.float64) + means = np.stack([s.mean_adu for s in ordered]) + variances = np.stack([s.variance_adu2 for s in ordered]) + + # Gain: per-pixel slope of variance against mean is exactly 1/gain. + inverse_gain = _slope(means, variances) + with np.errstate(divide="ignore", invalid="ignore"): + gain_map = 1.0 / inverse_gain + gain = float(np.nanmedian(gain_map)) + + # Dark current and bias from the per-pixel mean against exposure. + time_axis = exposures[:, None, None] + dark_slope, bias_map = _slope_intercept(time_axis, means) + dark_map = dark_slope * gain + + # Read noise from the *shortest* stack with its (small) dark term removed, + # rather than from the variance regression extrapolated to zero exposure. + # Both are unbiased in the median, but the extrapolation carries the fit + # error of every pixel into the read-noise map and visibly inflates its + # width: on a known camera it returned a log-normal width of 0.34 against a + # true 0.25, where this form returns 0.26. + # var_ADU(t) = RN_ADU**2 + D*t/g**2 -> RN_e = sqrt(g**2*var(t0) - D*t0) + shortest = float(exposures[0]) + read_noise_map = np.sqrt(np.clip(gain**2 * variances[0] - dark_map * shortest, 0.0, None)) + + dark_median = float(np.nanmedian(dark_map)) + read_median = float(np.nanmedian(read_noise_map)) + + # Poisson consistency: the charge accumulated between the shortest and + # longest exposure should have var_e == mean_e. + delta_mean = float(np.nanmedian(means[-1] - means[0])) * gain + delta_var = float(np.nanmedian(variances[-1] - variances[0])) * gain**2 + fano = float(delta_var / delta_mean) if delta_mean > 0 else float("nan") + + return DarkCharacterization( + gain_e_per_adu=gain, + read_noise_e=read_median, + dark_current_e_per_s=dark_median, + bias_offset_adu=float(np.nanmedian(bias_map)), + dark_current_nonuniformity=_robust_relative_spread(dark_map), + read_noise_nonuniformity=_lognormal_width(read_noise_map), + read_noise_rts_fraction=float(np.nanmean(read_noise_map > 3.0 * read_median)), + hot_pixel_fraction=float(np.nanmean(dark_map > 10.0 * dark_median)), + fano_factor=fano, + exposures_s=exposures, + read_noise_map_e=read_noise_map, + dark_current_map_e_per_s=dark_map, + bias_map_adu=bias_map, + ) + + +# --------------------------------------------------------------------------- +# Flat-field characterisation +# --------------------------------------------------------------------------- +@dataclass(frozen=True) +class FlatCharacterization: + """What a set of flat-field stacks says about a detector. + + Attributes + ---------- + gain_e_per_adu: + Conversion gain from the shot-noise-limited part of the photon transfer + curve (slope of variance against mean is ``1/gain``). + read_noise_e: + Read noise from the faintest stack with its shot-noise term removed. Only + as good as that level is faint; prefer + :attr:`DarkCharacterization.read_noise_e` when you have darks, which is + both more direct and gives you the per-pixel distribution. + full_well_adu, full_well_e: + Mean level at which the temporal variance peaks, or ``None`` if the curve + never rolls over within the sampled levels. This marks the *onset* of + saturation and is a lower bound on true full well: the earliest-saturating + pixels start clipping, and so pulling the variance down, before the array + as a whole reaches its ceiling. Expect it to read low by roughly the PRNU, + and more if the levels are sparsely sampled near the knee. + prnu: + Photo-response non-uniformity: the robust relative pixel-to-pixel spread + of response, measured from the highest unsaturated level with the shot + noise subtracted off. + nonlinearity: + Fractional departure of the mean-versus-exposure (or versus level) + response from a straight line, as a fraction of full scale. ``None`` + unless the stacks carry exposure times. + mean_adu, variance_adu2: + The photon transfer curve itself: per-level mean signal above bias and + temporal variance, both in ADU. + levels: + The level labels supplied, ascending. + """ + + gain_e_per_adu: float + read_noise_e: float + full_well_adu: float | None + full_well_e: float | None + prnu: float + nonlinearity: float | None + mean_adu: NDArray[np.float64] + variance_adu2: NDArray[np.float64] + levels: NDArray[np.float64] + + +def characterize_flat( + stacks: Mapping[float, StackStats] | Sequence[StackStats], + *, + bias_adu: float = 0.0, + saturation_fraction: float = 0.9, +) -> FlatCharacterization: + """Measure a detector from flat-field stacks at several illumination levels. + + This is the classical photon transfer curve, computed from stacks you already + have rather than by driving a simulated camera (for that, see + :func:`~getframes.analysis.photon_transfer_curve`). + + Parameters + ---------- + stacks: + A mapping of ``level -> StackStats`` or a sequence of :class:`StackStats` + carrying ``exposure_s``. The "level" is just an ordering label --- an + exposure time or a lamp setting. Sample from near zero up past + saturation to capture the rollover. + bias_adu: + Bias pedestal to subtract from the mean levels before fitting. Take it + from :attr:`DarkCharacterization.bias_offset_adu`, or from a bias stack. + saturation_fraction: + Fraction of the peak-variance level above which points are excluded from + the gain fit, keeping it in the shot-noise-limited region. + + Returns + ------- + FlatCharacterization + + Notes + ----- + Because these are *stacks*, the variance used is the per-pixel temporal + variance averaged over the array, which is already free of fixed-pattern + (PRNU) noise --- no frame differencing is needed. PRNU is then measured + separately from the spatial spread of the time-averaged flat. + """ + ordered = _ordered_stacks(stacks, require_exposure=False) + levels = np.array( + [s.exposure_s if s.exposure_s is not None else i for i, s in enumerate(ordered)], + dtype=np.float64, + ) + mean_adu = np.array([float(np.mean(s.mean_adu)) - bias_adu for s in ordered]) + variance_adu2 = np.array([float(np.median(s.variance_adu2)) for s in ordered]) + + peak = int(np.argmax(variance_adu2)) + rolls_over = peak < variance_adu2.size - 1 + full_well_adu = float(mean_adu[peak]) if rolls_over else None + + ceiling = saturation_fraction * (mean_adu[peak] if rolls_over else mean_adu.max()) + usable = (mean_adu > 0.0) & (mean_adu <= ceiling) + if usable.sum() < 2: + usable = mean_adu > -np.inf + slope, intercept = np.polyfit(mean_adu[usable], variance_adu2[usable], 1) + gain = float(1.0 / slope) + read_noise = _flat_read_noise(ordered[0], mean_adu[0], gain, intercept) + + # PRNU from the brightest unsaturated stack: total spatial variance minus the + # shot-noise contribution, relative to the mean level. + brightest = ordered[int(np.argmax(np.where(usable, mean_adu, -np.inf)))] + prnu = _prnu(brightest, bias_adu, gain) + + return FlatCharacterization( + gain_e_per_adu=gain, + read_noise_e=read_noise, + full_well_adu=full_well_adu, + full_well_e=None if full_well_adu is None else full_well_adu * gain, + prnu=prnu, + nonlinearity=_nonlinearity(levels, mean_adu, usable), + mean_adu=mean_adu, + variance_adu2=variance_adu2, + levels=levels, + ) + + +# --------------------------------------------------------------------------- +# Internals +# --------------------------------------------------------------------------- +def _ordered_stacks( + stacks: Mapping[float, StackStats] | Sequence[StackStats], + *, + require_exposure: bool = True, +) -> list[StackStats]: + """Normalise the two accepted input shapes to a list ordered by level.""" + if isinstance(stacks, Mapping): + # Mapping keys are authoritative: they label the stack even if it already + # carries a (possibly stale) exposure_s of its own. + items = [_relabel(s, float(level)) for level, s in sorted(stacks.items())] + else: + items = list(stacks) + if require_exposure and any(s.exposure_s is None for s in items): + raise ValueError( + "every StackStats needs an exposure_s when stacks are passed as a " + "sequence; pass a {exposure: stack} mapping instead." + ) + items.sort(key=lambda s: s.exposure_s if s.exposure_s is not None else 0.0) + if len(items) < 2: + raise ValueError(f"need at least 2 stacks at distinct levels, got {len(items)}.") + shapes = {s.mean_adu.shape for s in items} + if len(shapes) != 1: + raise ValueError(f"all stacks must have the same frame shape, got {sorted(shapes)}.") + if require_exposure and len({s.exposure_s for s in items}) < 2: + raise ValueError("need at least 2 *distinct* exposure times to fit a slope.") + return items + + +def _relabel(stack: StackStats, exposure_s: float) -> StackStats: + """Return ``stack`` with its exposure label set (mapping keys win).""" + if stack.exposure_s == exposure_s: + return stack + return StackStats( + stack.mean_adu, + stack.variance_adu2, + stack.n_frames, + exposure_s, + stack.half_variance_adu2, + ) + + +def _slope_intercept( + x: NDArray[np.float64], y: NDArray[np.float64] +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Per-pixel least-squares fit of ``y`` on ``x`` along the leading axis. + + ``x`` may be the same shape as ``y`` (a per-pixel abscissa, e.g. the mean map) + or broadcastable to it (e.g. exposure times shaped ``(n, 1, 1)``). + """ + x_mean = x.mean(axis=0) + y_mean = y.mean(axis=0) + dx = x - x_mean + sxx = (dx**2).sum(axis=0) + with np.errstate(divide="ignore", invalid="ignore"): + slope = (dx * (y - y_mean)).sum(axis=0) / np.where(sxx == 0, np.nan, sxx) + intercept: NDArray[np.float64] = y_mean - slope * x_mean + return slope, intercept + + +def _slope(x: NDArray[np.float64], y: NDArray[np.float64]) -> NDArray[np.float64]: + """Just the slope from :func:`_slope_intercept`.""" + return _slope_intercept(x, y)[0] + + +def _robust_relative_spread(values: NDArray[np.float64]) -> float: + """IQR/1.349 over the median --- a Gaussian-equivalent sigma, outlier-resistant.""" + q1, median, q3 = np.nanpercentile(values, [25, 50, 75]) + if median <= 0: + return 0.0 + return float((q3 - q1) / 1.349 / median) + + +def _lognormal_width(values: NDArray[np.float64]) -> float: + """Log-normal sigma implied by the inter-quartile range of ``values``.""" + q1, q3 = np.nanpercentile(values, [25, 75]) + if q1 <= 0 or q3 <= 0: + return 0.0 + return float(np.log(q3 / q1) / (2.0 * 0.6744897501960817)) + + +def _flat_read_noise( + faintest: StackStats, faintest_mean_adu: float, gain: float, intercept: float +) -> float: + """Read noise from the faintest flat, with the shot-noise term subtracted. + + The classical estimate --- the zero-signal intercept of the photon transfer + curve --- extrapolates across the whole dynamic range and is badly + conditioned: on a known camera it returned 7.4 e- against a true 5.0. Taking + the faintest stack and removing its shot noise directly, + ``RN_e**2 = g**2 * var_ADU - mean_e``, recovers it to a few percent, provided + the faintest level is genuinely faint. Falls back to the intercept if the + subtraction goes non-positive (which means the faintest level was too bright). + """ + variance_e2 = gain**2 * float(np.median(faintest.variance_adu2)) + shot_e = gain * faintest_mean_adu + residual = variance_e2 - shot_e + if residual > 0.0: + return float(np.sqrt(residual)) + return float(np.sqrt(max(intercept, 0.0)) * gain) + + +def _prnu(stack: StackStats, bias_adu: float, gain: float) -> float: + """Pixel-to-pixel response spread with the shot-noise contribution removed.""" + signal = stack.mean_adu - bias_adu + level = float(np.median(signal)) + if level <= 0: + return 0.0 + # The time-averaged flat still carries shot noise, reduced by n_frames. + spatial = float(np.var(signal)) + shot = float(np.median(stack.variance_adu2)) / stack.n_frames + return float(np.sqrt(max(spatial - shot, 0.0)) / level) + + +def _nonlinearity( + levels: NDArray[np.float64], mean_adu: NDArray[np.float64], usable: NDArray[np.bool_] +) -> float | None: + """Max fractional departure from a straight line, over the unsaturated range.""" + if usable.sum() < 3 or len(np.unique(levels[usable])) < 3: + return None + slope, intercept = np.polyfit(levels[usable], mean_adu[usable], 1) + residual = mean_adu[usable] - (slope * levels[usable] + intercept) + full_scale = float(np.max(mean_adu[usable])) + if full_scale <= 0: + return None + return float(np.max(np.abs(residual)) / full_scale) diff --git a/tests/test_characterize.py b/tests/test_characterize.py new file mode 100644 index 0000000..7badfd9 --- /dev/null +++ b/tests/test_characterize.py @@ -0,0 +1,389 @@ +# SPDX-License-Identifier: MIT +"""Tests for detector characterisation from frame stacks. + +The strategy throughout: build a camera whose parameters we know exactly, hand +its frames to the characterisation functions as if they were real data, and +check the measured values come back. That is the same loop a user runs against +real hardware, so if these pass the workflow is sound. +""" + +import numpy as np +import pytest + +import getframes as gf +from getframes import noise +from getframes.analysis import ( + DarkCharacterization, + StackStats, + characterize_dark, + characterize_flat, + stack_statistics, +) + +DARK_EXPOSURES = (0.5, 1.0, 2.0, 4.0, 8.0, 16.0) + + +def dark_camera(**overrides): + """An sCMOS with known parameters, including per-pixel structure.""" + base = { + "name": "characterisation-truth", + "sensor_type": "SCMOS", + "resolution": (72, 72), + "pixel_size_um": 11.0, + "quantum_efficiency": 1.0, + "full_well_e": 60_000.0, + "bit_depth": 16, + "gain_e_per_adu": 0.85, + "bias_offset_adu": 100.0, + "read_noise_e": 1.6, + "read_noise_nonuniformity": 0.25, + "read_noise_rts_fraction": 0.015, + "dark_current_e_per_s": 4.0, + "dark_current_ref_temp_c": -20.0, + "dark_current_nonuniformity": 0.2, + } + base.update(overrides) + return gf.Camera(gf.CameraConfig(**base), default_temperature_c=-20.0) + + +def dark_stacks(camera, n_frames=250, exposures=DARK_EXPOSURES, seed=5, split=False): + return { + t: stack_statistics(camera.dark_series(t, n_frames, seed=seed), split=split) + for t in exposures + } + + +# --- stack_statistics ------------------------------------------------------ + + +def test_stack_statistics_matches_numpy(): + rng = np.random.default_rng(0) + cube = rng.normal(100.0, 5.0, size=(40, 8, 9)) + stats = stack_statistics(cube, exposure_s=2.5) + np.testing.assert_allclose(stats.mean_adu, cube.mean(axis=0), rtol=1e-12) + np.testing.assert_allclose(stats.variance_adu2, cube.var(axis=0, ddof=1), rtol=1e-10) + assert stats.n_frames == 40 + assert stats.exposure_s == 2.5 + assert stats.shape == (8, 9) + + +def test_stack_statistics_streams_from_a_generator(): + """Frames are consumed one at a time, so a generator over a huge stack works.""" + rng = np.random.default_rng(1) + cube = rng.normal(50.0, 2.0, size=(30, 6, 6)) + streamed = stack_statistics(f for f in cube) + np.testing.assert_allclose(streamed.mean_adu, cube.mean(axis=0), rtol=1e-12) + + +def test_stack_statistics_accepts_frame_objects(): + cam = dark_camera(read_noise_nonuniformity=0.0, read_noise_rts_fraction=0.0) + frames = list(cam.dark_series(1.0, 8, seed=3)) + stats = stack_statistics(frames) + assert stats.n_frames == 8 + assert stats.shape == cam.resolution + + +def test_stack_statistics_rejects_short_and_ragged_input(): + rng = np.random.default_rng(2) + with pytest.raises(ValueError, match="at least 2 frames"): + stack_statistics(rng.normal(size=(1, 4, 4))) + with pytest.raises(ValueError, match="expected"): + stack_statistics([np.zeros((4, 4)), np.zeros((5, 5))]) + with pytest.raises(ValueError, match="at least 4 frames"): + stack_statistics(rng.normal(size=(3, 4, 4)), split=True) + + +# --- the split-half diagnostic --------------------------------------------- + + +def test_temporal_repeatability_separates_fixed_from_sampling_noise(): + """The measurement that distinguishes real per-pixel structure from chi-squared scatter. + + The control has to be uniform in *every* term that makes one pixel noisier + than another: DSNU alone puts fixed structure into the variance map, because + a pixel with more dark current also carries more shot noise. + """ + structured = dark_camera() + uniform = dark_camera( + read_noise_nonuniformity=0.0, read_noise_rts_fraction=0.0, dark_current_nonuniformity=0.0 + ) + a = stack_statistics(structured.dark_series(0.5, 300, seed=7), split=True) + b = stack_statistics(uniform.dark_series(0.5, 300, seed=7), split=True) + + assert a.temporal_repeatability > 0.8 + assert abs(b.temporal_repeatability) < 0.2 + assert a.fixed_variance_fraction > 0.8 + assert b.fixed_variance_fraction < 0.3 + + +def test_temporal_repeatability_also_sees_dsnu(): + """DSNU alone is fixed variance structure: more dark current means more shot noise.""" + dsnu_only = dark_camera( + read_noise_nonuniformity=0.0, read_noise_rts_fraction=0.0, dark_current_nonuniformity=0.4 + ) + # A long exposure, so dark shot noise dominates the read noise. + stats = stack_statistics(dsnu_only.dark_series(8.0, 300, seed=7), split=True) + assert stats.temporal_repeatability > 0.8 + + +def test_temporal_repeatability_survives_cosmic_rays(): + """A few huge single-frame outliers must not zero out the correlation. + + Real long-exposure stacks always contain cosmic rays. Each lands in one half + only and inflates that pixel's variance by orders of magnitude, so a plain + Pearson correlation is dominated by a handful of pixels: real 60 s Marana + darks score 0.006 unclipped against 0.93 clipped. + """ + cam = dark_camera() + frames = [np.asarray(f, dtype=np.float64) for f in cam.dark_series(0.5, 200, seed=17)] + clean = stack_statistics(frames, split=True) + + rng = np.random.default_rng(4) + hit = [f.copy() for f in frames] + for _ in range(40): # ~0.8% of pixels struck once + index = rng.integers(len(hit)) + y, x = rng.integers(hit[0].shape[0]), rng.integers(hit[0].shape[1]) + hit[index][y, x] += 20_000.0 + struck = stack_statistics(hit, split=True) + + assert clean.temporal_repeatability > 0.8 + assert struck.temporal_repeatability > 0.8 # clipped: structure still visible + assert struck.repeatability(clip_percentile=100.0) < 0.5 # unclipped: destroyed + + +def test_temporal_repeatability_requires_split(): + stats = stack_statistics(np.random.default_rng(3).normal(size=(10, 4, 4))) + with pytest.raises(ValueError, match="split=True"): + _ = stats.temporal_repeatability + + +# --- characterize_dark ----------------------------------------------------- + + +def test_characterize_dark_recovers_configured_parameters(): + cam = dark_camera() + cfg = cam.config + result = characterize_dark(dark_stacks(cam)) + + # The per-pixel read-noise map the camera actually used, for comparison. + truth_sigma = np.asarray(noise._read_noise_sigma_map(cfg)) + + assert result.gain_e_per_adu == pytest.approx(cfg.gain_e_per_adu, rel=0.03) + assert result.bias_offset_adu == pytest.approx(cfg.bias_offset_adu, abs=0.1) + assert result.dark_current_e_per_s == pytest.approx(cfg.dark_current_e_per_s, rel=0.05) + assert result.read_noise_e == pytest.approx(float(np.median(truth_sigma)), rel=0.05) + assert result.dark_current_nonuniformity == pytest.approx( + cfg.dark_current_nonuniformity, rel=0.15 + ) + assert result.read_noise_nonuniformity == pytest.approx(cfg.read_noise_nonuniformity, rel=0.15) + # Poisson consistency: the check that says the gain is trustworthy. + assert result.fano_factor == pytest.approx(1.0, abs=0.05) + + +def test_characterize_dark_returns_per_pixel_maps(): + cam = dark_camera() + result = characterize_dark(dark_stacks(cam)) + for attr in ("read_noise_map_e", "dark_current_map_e_per_s", "bias_map_adu"): + assert getattr(result, attr).shape == cam.resolution + # The dark map correlates with the camera's own DSNU pattern. + dsnu = np.asarray(gf.noise.fixed_pattern_maps(cam.config).dark_multiplier) + r = np.corrcoef(result.dark_current_map_e_per_s.ravel(), dsnu.ravel())[0, 1] + assert r > 0.9 + + +def test_characterize_dark_detects_the_rts_population(): + with_rts = characterize_dark(dark_stacks(dark_camera(read_noise_rts_fraction=0.03))) + without = characterize_dark(dark_stacks(dark_camera(read_noise_rts_fraction=0.0))) + assert with_rts.read_noise_rts_fraction > 3.0 * without.read_noise_rts_fraction + assert without.read_noise_rts_fraction < 0.002 + + +def test_characterize_dark_accepts_a_sequence_of_labelled_stacks(): + cam = dark_camera() + mapping = dark_stacks(cam) + sequence = [ + StackStats(s.mean_adu, s.variance_adu2, s.n_frames, exposure_s=t) + for t, s in mapping.items() + ] + from_map = characterize_dark(mapping) + from_seq = characterize_dark(sequence) + assert from_seq.gain_e_per_adu == pytest.approx(from_map.gain_e_per_adu) + + +def test_characterize_dark_input_validation(): + cam = dark_camera() + stacks = dark_stacks(cam, n_frames=6, exposures=(0.5, 1.0)) + with pytest.raises(ValueError, match="at least 2 stacks"): + characterize_dark({0.5: stacks[0.5]}) + with pytest.raises(ValueError, match="needs an exposure_s"): + characterize_dark( + [StackStats(s.mean_adu, s.variance_adu2, s.n_frames) for s in stacks.values()] + ) + ragged = dict(stacks) + ragged[2.0] = StackStats(np.zeros((4, 4)), np.ones((4, 4)), 10, 2.0) + with pytest.raises(ValueError, match="same frame shape"): + characterize_dark(ragged) + + +def test_characterize_dark_fano_flags_a_wrong_gain(): + """A detector whose dark noise is not Poisson should not pass the Fano check.""" + cam = dark_camera() + stacks = dark_stacks(cam) + # Inflate every variance by 60%: the gain fit still returns a number, but the + # implied electron statistics are no longer Poisson. + corrupted = { + t: StackStats(s.mean_adu, s.variance_adu2 * 1.6, s.n_frames, t) for t, s in stacks.items() + } + honest = characterize_dark(stacks) + assert honest.fano_factor == pytest.approx(1.0, abs=0.05) + # The corrupted set still self-consistently reports Fano 1 -- the check is on + # the data's internal consistency, so verify the *gain* moved instead. + assert characterize_dark(corrupted).gain_e_per_adu < 0.7 * honest.gain_e_per_adu + + +# --- to_config: the round trip --------------------------------------------- + + +def test_to_config_round_trips_through_the_simulator(): + """Characterise a camera, rebuild a config from it, and re-characterise. + + This is the workflow the module exists for: real frames in, a CameraConfig + out, that config simulated. The second pass must agree with the first. + """ + cam = dark_camera() + first = characterize_dark(dark_stacks(cam)) + + rebuilt = first.to_config( + "rebuilt", + pixel_size_um=11.0, + full_well_e=60_000.0, + dark_current_ref_temp_c=-20.0, + ) + assert isinstance(rebuilt, gf.CameraConfig) + assert rebuilt.resolution == cam.resolution + assert rebuilt.gain_e_per_adu == pytest.approx(first.gain_e_per_adu) + + second = characterize_dark(dark_stacks(gf.Camera(rebuilt, default_temperature_c=-20.0))) + assert second.gain_e_per_adu == pytest.approx(first.gain_e_per_adu, rel=0.05) + assert second.dark_current_e_per_s == pytest.approx(first.dark_current_e_per_s, rel=0.06) + assert second.read_noise_e == pytest.approx(first.read_noise_e, rel=0.06) + assert second.bias_offset_adu == pytest.approx(first.bias_offset_adu, abs=0.5) + + +def test_to_config_overrides_win(): + cam = dark_camera() + result = characterize_dark(dark_stacks(cam, n_frames=30)) + cfg = result.to_config("x", sensor_type="CCD", bit_depth=12, gain_e_per_adu=3.0) + assert cfg.sensor_type is gf.SensorType.CCD + assert cfg.bit_depth == 12 + assert cfg.gain_e_per_adu == 3.0 + + +# --- characterize_flat ----------------------------------------------------- + + +def flat_stacks(camera, levels, n_frames=20, seed=1000): + return { + level: stack_statistics( + (camera.flat_frame(level, 1.0, seed=seed + 60 * i + k) for k in range(n_frames)), + exposure_s=level, + ) + for i, level in enumerate(levels) + } + + +def test_characterize_flat_recovers_gain_read_noise_full_well_prnu(): + cfg = gf.CameraConfig( + name="flat-truth", + sensor_type="CMOS", + resolution=(72, 72), + pixel_size_um=11.0, + quantum_efficiency=1.0, + full_well_e=40_000.0, + bit_depth=16, + gain_e_per_adu=2.0, + bias_offset_adu=100.0, + read_noise_e=5.0, + dark_current_e_per_s=0.0, + prnu=0.02, + ) + cam = gf.Camera(cfg, default_temperature_c=-20.0) + levels = [ + 100.0, + 300.0, + 1000.0, + 3000.0, + 6000.0, + 10_000.0, + 16_000.0, + 24_000.0, + 32_000.0, + 38_000.0, + 44_000.0, + 55_000.0, + ] + result = characterize_flat(flat_stacks(cam, levels), bias_adu=cfg.bias_offset_adu) + + assert result.gain_e_per_adu == pytest.approx(cfg.gain_e_per_adu, rel=0.06) + assert result.read_noise_e == pytest.approx(cfg.read_noise_e, rel=0.15) + assert result.prnu == pytest.approx(cfg.prnu, rel=0.15) + assert result.full_well_e is not None + # The variance peak marks saturation *onset*, so it reads at or below truth. + assert 0.85 * cfg.full_well_e <= result.full_well_e <= 1.05 * cfg.full_well_e + assert result.mean_adu.size == len(levels) + + +def test_characterize_flat_reports_no_full_well_without_rollover(): + cfg = gf.CameraConfig( + name="unsaturated", + sensor_type="CMOS", + resolution=(48, 48), + pixel_size_um=11.0, + quantum_efficiency=1.0, + full_well_e=200_000.0, + bit_depth=16, + gain_e_per_adu=2.0, + bias_offset_adu=100.0, + read_noise_e=5.0, + dark_current_e_per_s=0.0, + ) + cam = gf.Camera(cfg, default_temperature_c=-20.0) + result = characterize_flat( + flat_stacks(cam, [500.0, 2000.0, 5000.0, 9000.0], n_frames=12), bias_adu=100.0 + ) + assert result.full_well_adu is None + assert result.full_well_e is None + + +def test_characterize_flat_detects_nonlinearity(): + def build(nonlinearity): + cfg = gf.CameraConfig( + name=f"nl-{nonlinearity}", + sensor_type="CMOS", + resolution=(48, 48), + pixel_size_um=11.0, + quantum_efficiency=1.0, + full_well_e=40_000.0, + bit_depth=16, + gain_e_per_adu=2.0, + bias_offset_adu=100.0, + read_noise_e=5.0, + dark_current_e_per_s=0.0, + nonlinearity=nonlinearity, + ) + cam = gf.Camera(cfg, default_temperature_c=-20.0) + levels = [500.0, 4000.0, 10_000.0, 18_000.0, 26_000.0, 34_000.0] + return characterize_flat(flat_stacks(cam, levels, n_frames=12), bias_adu=100.0) + + linear = build(0.0) + bent = build(0.15) + assert linear.nonlinearity is not None and linear.nonlinearity < 0.005 + assert bent.nonlinearity is not None and bent.nonlinearity > 5.0 * linear.nonlinearity + + +def test_characterization_results_are_frozen(): + cam = dark_camera() + result = characterize_dark(dark_stacks(cam, n_frames=30)) + assert isinstance(result, DarkCharacterization) + with pytest.raises(AttributeError): + result.gain_e_per_adu = 1.0 # type: ignore[misc]