diff --git a/AGENTS.md b/AGENTS.md index 0b362da..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`. | @@ -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..fcf21f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,23 @@ 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 - Full-detector region-of-interest simulation through @@ -15,6 +32,61 @@ 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 + 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. +- `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. 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/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=" 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/docs/reference.md b/docs/reference.md index a57fa84..3b6cc28 100644 --- a/docs/reference.md +++ b/docs/reference.md @@ -86,6 +86,8 @@ ::: getframes.analysis.ptc +::: getframes.analysis.characterize + ## Datasets & scale ::: getframes.dataset diff --git a/examples/15_detector_characterization.py b/examples/15_detector_characterization.py new file mode 100644 index 0000000..10e0e03 --- /dev/null +++ b/examples/15_detector_characterization.py @@ -0,0 +1,259 @@ +# SPDX-License-Identifier: MIT +"""Characterise a detector from frame stacks, then rebuild it as a CameraConfig. + +Example 06 drives a *simulated* camera to build a photon transfer curve. This one +works the other way round, on stacks of frames that already exist --- which is what +you have when the frames came off real hardware. + +The workflow is: + + frames -> 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/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_cb1_0_5mp.toml b/src/getframes/presets/data/andor_cb1_0_5mp.toml index 9fe03d4..7aa1a5b 100644 --- a/src/getframes/presets/data/andor_cb1_0_5mp.toml +++ b/src/getframes/presets/data/andor_cb1_0_5mp.toml @@ -13,12 +13,27 @@ 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." +# 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. +# 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] body_dimensions_order = "length, width, height; confirm mounting/cable clearance in the Keck envelope" 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..1e64c29 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,32 @@ 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 +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 = "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/generic_scmos.toml b/src/getframes/presets/data/generic_scmos.toml index 183cf98..93b7578 100644 --- a/src/getframes/presets/data/generic_scmos.toml +++ b/src/getframes/presets/data/generic_scmos.toml @@ -15,7 +15,9 @@ 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 +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 80d525a..613250e 100644 --- a/src/getframes/presets/data/hamamatsu_orca_fusion.toml +++ b/src/getframes/presets/data/hamamatsu_orca_fusion.toml @@ -19,7 +19,9 @@ 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 +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 53b161c..7198d79 100644 --- a/src/getframes/presets/data/hamamatsu_orca_quest_2.toml +++ b/src/getframes/presets/data/hamamatsu_orca_quest_2.toml @@ -13,13 +13,27 @@ 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 +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. + +# 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/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/src/getframes/presets/data/tucsen_aries_6504_pro.toml b/src/getframes/presets/data/tucsen_aries_6504_pro.toml index 013ea00..7e4dd3a 100644 --- a/src/getframes/presets/data/tucsen_aries_6504_pro.toml +++ b/src/getframes/presets/data/tucsen_aries_6504_pro.toml @@ -18,11 +18,24 @@ 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 +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." +# 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. +# 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] body_dimensions_order = "width, height, depth; bare camera body only, excluding connectors and adapters" 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] diff --git a/tests/test_realism.py b/tests/test_realism.py index b86d081..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 +from getframes import Camera, CameraConfig, SensorType, available_presets, load_preset, noise def base_config(**overrides): @@ -84,12 +84,204 @@ 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_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 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 # ---------------------------------------------------------------------------