2D detector performance characterization

2D CCD and sCMOS detectors are used to measure scattered X-rays in RSoXS experiments. They are the primary detectors of the RSoXS station. This page details the performance characterization of these detectors to ensure reliable data.



Background

The detector converts incident photons into electrons that are then converted into a digital signal. Detctor performance characterization involves quantifying various metrics throughout this process. The metrics describing the photon-to-electron conversion are related to the quantum perforamance of the detector. The metrics describing the electron-to-digitial signal conversion are related to the electronic performance of the detector.

A combination of the photon transfer curve (PTC) and X-ray transfer curve is the ideal approach to characterize the electronic and quantum performance of the detector, respectively. The PTC alone falls short because low-energy visible light photons are absorbed near the surface of the sensor and are not representative of the deeper-penetrating high-energy X-ray photons. As a result, the PTC is insufficient for measuring charge collection efficiency, charge transfer efficiency, and energy resolution.


The PTC is a graph of the noise (\(\sigma\)) versus signal (\(S_{DN}\), where DN is digital units, as opposed to electrons or photons) of the detector plotted across the full operational exposure time range of the detector. A log-log representation of this plot will show the following regimes:

  • At very low exposure times (low signal), the curve flatlines to a slope of ~0. The fixed electronic read noise dominates the noise.
  • At low-moderate exposure times and increasing signal, the curve becomes linear with a slope of 1/2. The random Poisson-statistical arrival of photons is the dominant noise source.
  • At moderate-high exposure times, the curve will bend up and have a slope of 1. The spatial fixed pattern noise related to pixel-to-pixel manufacturing non-uniformities will become larger than the random shot noise. This regime will only be present in cases where the image data is not temporally averaged.
  • At very high exposure times and high signal, the variance will either drastically skew upward or go to zero. The pixels have reached full-well capacity and are saturated.


References:

  • James R. Janesick, “Scientific Charge-Coupled Devices” (textbook)
  • https://www.teledynevisionsolutions.com/learn/learning-center/imaging-fundamentals/camera-test-protocol/
  • https://camera.hamamatsu.com/content/dam/hamamatsu-photonics/sites/static/sys/en/documents/FOM2017_Poster.pdf
  • https://www.researching.cn/ArticlePdf/m00092/2016/9/3/JIOHS-2016-03-1630008.pdf
  • https://www.adimec.com/how-to-measure-the-photon-transfer-curve-for-ccd-or-cmos-cameras/
  • https://www.adimec.com/how-to-measure-the-photon-transfer-curve-for-ccd-or-cmos-cameras-2/
  • https://pubmed.ncbi.nlm.nih.gov/36256233/
  • https://blog.kasson.com/the-last-word/photon-transfer-curves-a-unifying-tool-for-sensor-noise/
  • https://www.msss.com/http/near_cal/robinson_stuff/ltc.html



Physical setup

The detector may be characterized in situ in a beamline or ex situ (bench testing).

For out-of-vacuum bench testing, thermoelectric cooling (TEC) must never be turned on, because it will cause ambient moisture to condense onto the sensor. The condensation will contaminate the sensor, obscure the signal, and cause stresses that can physically damage the sensor. Although TEC should never be turned on, it is ideal to circulate a coolant through the cooling lines to maintain a stable sensor temperature while capturing images.



Data collection

The set of data described here will provide all of the information necessary to analyze the detector performance metrics.

Gather the following equipment/resources:

  • The camera to be tested
  • A means to isolate the sensor from all light. E.g., a Nikon lens with a cap for bench testing.
  • Light source. Ideally, this would be a flat-field illumination source whose intensity does not drift over time. However, this may be difficult in practice especially for the X-ray illumination data, as the beam is a small size, and its intensity oscillates over a period of ~1 day.


Frame rate stack: Close any shutters, and ensure that the camera sensor is shielded from any light. Set the camera to the lowest-possible exposure time (burst-capture). Capture 200-500 frames at full sensor ROI. Repeat at smaller ROIs.

Dark stack: Close any shutters, and ensure that the camera sensor is shielded from any light. Capture replicates of dark frames at exposure times spanning across the camera’s capability (e.g., 10 us to 300 s). For bench tests, only capture these images at room temperature (or include a few temperatures close to room temperature). However, in the beamline, it is good to capture stacks between room temperature and the lowest operating temperature. Note, the full ROI frame rate stack can be used as part of the lowest exposure time datapoint for the dark stack as well.


Visible light stack: Set up a light source with fixed, continuous intensity (no flickering or drifting). In a first pass, try capturing images at a few exposure times and pinpoint the point at which saturation occurs. In a second pass, divide the time between minimum exposure and saturation exposure into 20-30 steps. Add ~10 additional steps for exposure times past saturation. Then collect replicates of images at each exposure time in this series.


At all exposure times, capture at least 2 replicates. More replicates will improve statistics. Other than the burst-frame stack, ensure that there is sufficient time between images for the sensor to cool down to its temperature setpoint.


X-ray stack:



Detector metrics and analysis

To avoid confounding metrics, the order of data processing should follow the hierarchy of measurement dependencies. Generally, most electronic performance metrics should be characterized first, as the characterization of many quantum performance metrics depends on correction of the electronic performance metrics.

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

import copy
from typing import Dict, List, Optional, Sequence
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm


def create_figure(
    number_rows: int = 1,
    number_columns: int = 1,
    subplot_size: Sequence[float] = (3.25, 3.25),
    margin_subplot_size: Sequence[float] = (0.6, 0.6),
    margin_figure_size: Sequence[float] = (0.8, 1.2, 0.4, 0.6),
    figure_title: str = "",
):
    """
    Initialize a figure with deterministic subplot geometry for publication plots and animations.

    This utility intentionally avoids plt.tight_layout() because automatic layout adjustments can
    produce frame-to-frame subplot drift in animations. Instead, absolute subplot and margin sizes
    are used to keep axis boxes spatially stable across renders.

    Parameters
    ----------
    number_rows : int
        Number of subplot rows.
    number_columns : int
        Number of subplot columns.
    subplot_size : sequence of float
        (width_in, height_in) for each subplot panel.
    margin_subplot_size : sequence of float
        (horizontal_gap_in, vertical_gap_in) between adjacent subplots.
    margin_figure_size : sequence of float
        (left_in, right_in, top_in, bottom_in) figure margins in inches.
    figure_title : str
        Optional figure title.

    Returns
    -------
    tuple
        (fig, axs) where axs is always returned as a 2D NumPy array for consistent indexing.
    """
    figure_width = (
        number_columns * subplot_size[0]
        + margin_figure_size[0]
        + margin_figure_size[1]
        + (number_columns - 1) * margin_subplot_size[0]
    )
    figure_height = (
        number_rows * subplot_size[1]
        + margin_figure_size[2]
        + margin_figure_size[3]
        + (number_rows - 1) * margin_subplot_size[1]
    )

    fig, axs = plt.subplots(
        nrows=number_rows,
        ncols=number_columns,
        figsize=(figure_width, figure_height),
        edgecolor=(0, 0, 0, 0),
        linewidth=3,
        dpi=150,
    )

    # Convert inch margins/gaps to figure-relative fractions for subplots_adjust().
    fig.subplots_adjust(
        left=margin_figure_size[0] / figure_width,
        right=1 - (margin_figure_size[1] / figure_width),
        top=1 - (margin_figure_size[2] / figure_height),
        bottom=margin_figure_size[3] / figure_height,
        wspace=margin_subplot_size[0] / subplot_size[0],
        hspace=margin_subplot_size[1] / subplot_size[1],
    )

    # Ensure axs is always 2D for predictable indexing: axs[row, col].
    if number_rows == 1 and number_columns == 1:
        axs = np.array([[axs]])
    elif number_rows == 1 and number_columns > 1:
        axs = np.asarray(axs).reshape(1, number_columns)
    elif number_rows > 1 and number_columns == 1:
        axs = np.asarray(axs).reshape(number_rows, 1)

    fig.suptitle(figure_title, color=(0, 0, 0, 1), fontname="Calibri", size=12)
    return fig, axs


def format_subplot_1d(
    ax,
    plots: Optional[List[Dict]] = None,
    plots_yright: Optional[List[Dict]] = None,
    plots_errorbar: Optional[List[Dict]] = None,
    plots_errorbar_yright: Optional[List[Dict]] = None,
    axvlines: Optional[List[Dict]] = None,
    axhlines: Optional[List[Dict]] = None,
    x_axis_label: Optional[Dict] = None,
    y_axis_label: Optional[Dict] = None,
    yright_axis_label: Optional[Dict] = None,
    xscale: str = "linear",
    yscale: str = "linear",
    yrightscale: str = "linear",
    xlim: Optional[Dict] = None,
    ylim: Optional[Dict] = None,
    yrightlim: Optional[Dict] = None,
    border_colors: Sequence = ((0, 0, 0, 1), (0, 0, 0, 1), (0, 0, 0, 1), (0, 0, 0, 1)),
):
    """
    Format a single 1D subplot with optional primary and secondary y-axis data.

    All plotting inputs are dictionaries of Matplotlib kwargs plus required coordinate fields.
    For example, each element in `plots` should include keys `x` and `y`, plus optional style keys.

    Parameters
    ----------
    ax : matplotlib.axes.Axes
        Target axis.
    plots, plots_yright, plots_errorbar, plots_errorbar_yright : list of dict, optional
        Plot specification dictionaries for left and right y axes.
    axvlines, axhlines : list of dict, optional
        Matplotlib axvline/axhline keyword dictionaries.
    x_axis_label, y_axis_label, yright_axis_label : dict, optional
        Label dictionaries. Use keys:
        - x_axis_label: {'xlabel': 'Text', ...matplotlib kwargs...}
        - y_axis_label: {'ylabel': 'Text', ...matplotlib kwargs...}
        - yright_axis_label: {'ylabel': 'Text', ...matplotlib kwargs...}
    xscale, yscale, yrightscale : str
        Axis scales, typically 'linear' or 'log'.
    xlim, ylim, yrightlim : dict, optional
        Axis limit dictionaries forwarded to set_xlim/set_ylim.
    border_colors : sequence
        RGBA colors for [left, right, top, bottom] spines.

    Returns
    -------
    dict
        Dictionary with references to plotted axes: {'ax': ax, 'axright': axright or None}.
    """
    plots = [] if plots is None else plots
    plots_yright = [] if plots_yright is None else plots_yright
    plots_errorbar = [] if plots_errorbar is None else plots_errorbar
    plots_errorbar_yright = [] if plots_errorbar_yright is None else plots_errorbar_yright
    axvlines = [] if axvlines is None else axvlines
    axhlines = [] if axhlines is None else axhlines

    xlim = {"left": None, "right": None, "auto": False} if xlim is None else xlim
    ylim = {"bottom": None, "top": None, "auto": False} if ylim is None else ylim
    yrightlim = (
        {"bottom": None, "top": None, "auto": False} if yrightlim is None else yrightlim
    )

    yright_conditions = bool(plots_yright or plots_errorbar_yright)
    axright = ax.twinx() if yright_conditions else None

    plot_settings_default = {
        "marker": "o",
        "markersize": 0,
        "markerfacecolor": (0, 0, 0, 1),
        "markeredgecolor": (0, 0, 0, 1),
        "markeredgewidth": 1,
        "color": (0, 0, 0, 1),
        "linestyle": "solid",
        "linewidth": 2,
    }

    def render_plots(target_ax, line_specs, errorbar_specs):
        for spec_group in [line_specs, errorbar_specs]:
            plot_colors = cm.rainbow(np.linspace(1, 0, len(spec_group))) if len(spec_group) > 0 else []
            if len(spec_group) == 1:
                plot_colors = [(0, 0, 0, 1)]

            for idx, user_settings in enumerate(spec_group):
                settings = {
                    **plot_settings_default,
                    "markerfacecolor": plot_colors[idx],
                    "color": plot_colors[idx],
                    **user_settings,
                }
                x = settings.pop("x")
                y = settings.pop("y")
                target_ax.plot(x, y, **settings)

    render_plots(ax, plots, plots_errorbar)
    if yright_conditions and axright is not None:
        render_plots(axright, plots_yright, plots_errorbar_yright)

    line_settings_default = {
        **plot_settings_default,
        "color": (0.5, 0.5, 0.5, 1),
        "linestyle": "dotted",
    }

    for line_type in [axvlines, axhlines]:
        plot_colors = cm.rainbow(np.linspace(1, 0, len(line_type))) if len(line_type) > 0 else []
        if len(line_type) == 1:
            plot_colors = [(0.5, 0.5, 0.5, 1)]

        for idx, line_settings_user in enumerate(line_type):
            line_settings = {
                **line_settings_default,
                "markerfacecolor": plot_colors[idx],
                "color": plot_colors[idx],
                **line_settings_user,
            }
            if line_type is axvlines:
                ax.axvline(**line_settings)
            else:
                ax.axhline(**line_settings)

    if x_axis_label:
        x_axis_label_copy = copy.deepcopy(x_axis_label)
        xlabel = x_axis_label_copy.pop("xlabel")
        ax.set_xlabel(xlabel, **x_axis_label_copy)
    if y_axis_label:
        y_axis_label_copy = copy.deepcopy(y_axis_label)
        ylabel = y_axis_label_copy.pop("ylabel")
        ax.set_ylabel(ylabel, **y_axis_label_copy)
    if yright_axis_label and axright is not None:
        yright_axis_label_copy = copy.deepcopy(yright_axis_label)
        yrightlabel = yright_axis_label_copy.pop("ylabel")
        axright.set_ylabel(yrightlabel, **yright_axis_label_copy)

    ax.set_xscale(xscale)
    ax.set_yscale(yscale)
    if yright_conditions and axright is not None:
        axright.set_yscale(yrightscale)

    ax.set_xlim(**xlim)
    ax.set_ylim(**ylim)
    if yright_conditions and axright is not None:
        axright.set_ylim(**yrightlim)

    def border_formatting(target_ax, axes_tick_params=(3, 0)):
        for idx, border in enumerate(["left", "right", "top", "bottom"]):
            target_ax.spines[border].set_linewidth(2)
            target_ax.spines[border].set_color(border_colors[idx])
        for idx_axis, axis in enumerate(["x", "y"]):
            target_ax.tick_params(
                axis=axis,
                colors=border_colors[axes_tick_params[idx_axis]],
                width=2,
            )
        target_ax.tick_params(axis="both", which="major", labelsize=10)

    border_formatting(ax)
    if yright_conditions and axright is not None:
        border_formatting(axright, axes_tick_params=(3, 1))

    return {"ax": ax, "axright": axright}
```



Electronic performance



Frame rate

Physics:

Calculation: Take the burst-frame dark stack. Use the camera’s internal FPGA clock to extract timestamps and calculate time difference between consecutive frames. Check for frame dropping, in which the interval suddenly increases because the computer’s USB bus could not keep up and dropped a frame; eliminate such data points. From the good datapoints, caluclate a mean/median time difference. The frame rate is equal to 1 divided by the mean time difference. Repeat these calculations for all ROIs and plot frame rate vs. ROI.

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

import numpy as np
from typing import Any, Dict, Optional, Sequence


def compute_frame_intervals_from_timestamps(
    timestamps_s: Sequence[float],
    assume_monotonic: bool = True,
    sort_if_needed: bool = False,
    dropped_frame_zscore_threshold: float = 5.0,
    dropped_frame_absolute_factor: float = 1.5,
    expected_interval_s: Optional[float] = None,
) -> Dict[str, Any]:
    """
    Compute frame-to-frame timing statistics from camera timestamps and flag dropped-frame events.

    This function is designed for the burst-frame dark stack workflow described in the notebook text.
    It takes a 1D sequence of timestamps (typically from a camera FPGA clock, converted to seconds),
    computes frame intervals, and returns robust summary metrics. Potential dropped frames are detected
    with two complementary criteria:
      1) Robust outlier criterion using a MAD-based z-score of the frame intervals.
      2) Engineering criterion that an interval is larger than an expected cadence by a
         configurable multiplicative factor (default 1.5x).

    Parameters
    ----------
    timestamps_s : Sequence[float]
        Monotonic frame timestamps in seconds (or any consistent time unit, as long as intervals are
        interpreted in the same unit).
    assume_monotonic : bool, optional
        If True, require non-decreasing timestamp order unless sort_if_needed=True.
    sort_if_needed : bool, optional
        If True and timestamps are not monotonic, timestamps are sorted before interval computation.
    dropped_frame_zscore_threshold : float, optional
        Threshold on robust z-score of frame intervals for outlier detection.
    dropped_frame_absolute_factor : float, optional
        Intervals larger than dropped_frame_absolute_factor * expected interval are flagged.
    expected_interval_s : float, optional
        Known nominal frame interval (e.g., from camera settings). If omitted, median interval is used.

    Returns
    -------
    Dict[str, Any]
        Dictionary containing:
        - interval arrays (all and filtered)
        - boolean mask for dropped intervals
        - robust cadence estimates
        - frame rate estimates before/after filtering
        - warning messages for diagnostics

    Notes
    -----
    The final frame-rate estimate for characterization should usually use the filtered intervals
    (dropped-frame events removed), while retaining the unfiltered estimate for traceability.
    """
    timestamps = np.asarray(timestamps_s, dtype=float)
    if timestamps.ndim != 1 or timestamps.size < 2:
        raise ValueError("timestamps_s must be a 1D sequence with at least 2 elements.")

    if assume_monotonic and np.any(np.diff(timestamps) < 0):
        if sort_if_needed:
            timestamps = np.sort(timestamps)
        else:
            raise ValueError(
                "timestamps are not monotonic. Set sort_if_needed=True or pass monotonic input."
            )

    intervals = np.diff(timestamps)
    if np.any(intervals <= 0):
        raise ValueError(
            "Non-positive intervals found. Check timestamp source, time unit conversion, and ordering."
        )

    median_interval = float(np.median(intervals))
    if expected_interval_s is None:
        expected_interval_s = median_interval

    # Robust z-score based on median absolute deviation (MAD).
    mad = float(np.median(np.abs(intervals - median_interval)))
    if mad == 0.0:
        robust_z = np.zeros_like(intervals)
    else:
        robust_z = 0.6745 * (intervals - median_interval) / mad

    dropped_by_z = robust_z > dropped_frame_zscore_threshold
    dropped_by_abs = intervals > (dropped_frame_absolute_factor * expected_interval_s)
    dropped_mask = dropped_by_z | dropped_by_abs

    intervals_kept = intervals[~dropped_mask]
    warnings = []
    if intervals_kept.size == 0:
        warnings.append(
            "All intervals were flagged as dropped; returning unfiltered values for rate estimates."
        )
        intervals_kept = intervals.copy()

    result = {
        "n_frames": int(timestamps.size),
        "n_intervals": int(intervals.size),
        "n_dropped_intervals": int(np.sum(dropped_mask)),
        "fraction_dropped": float(np.mean(dropped_mask)),
        "intervals_s": intervals,
        "intervals_s_filtered": intervals_kept,
        "dropped_mask": dropped_mask,
        "expected_interval_s": float(expected_interval_s),
        "median_interval_s": median_interval,
        "mean_interval_s": float(np.mean(intervals)),
        "median_interval_s_filtered": float(np.median(intervals_kept)),
        "mean_interval_s_filtered": float(np.mean(intervals_kept)),
        "frame_rate_hz_raw_mean": float(1.0 / np.mean(intervals)),
        "frame_rate_hz_raw_median": float(1.0 / np.median(intervals)),
        "frame_rate_hz_filtered_mean": float(1.0 / np.mean(intervals_kept)),
        "frame_rate_hz_filtered_median": float(1.0 / np.median(intervals_kept)),
        "warnings": warnings,
    }
    return result


def extract_timestamps_from_xarray_like(
    frame_stack: Any,
    time_coord_name: str = "timestamp_s",
) -> np.ndarray:
    """
    Extract frame timestamps from an xarray-like object.

    This helper intentionally keeps the interface flexible so you can pass either:
      - an xarray.DataArray with a time coordinate, or
      - an xarray.Dataset variable with a matching coordinate.

    Parameters
    ----------
    frame_stack : Any
        xarray.DataArray or object exposing `.coords[time_coord_name]`.
    time_coord_name : str, optional
        Name of the coordinate containing timestamps in seconds.

    Returns
    -------
    np.ndarray
        1D timestamp array in seconds.
    """
    if not hasattr(frame_stack, "coords") or time_coord_name not in frame_stack.coords:
        raise ValueError(
            f"Input object does not contain coordinate '{time_coord_name}'."
        )
    timestamps = np.asarray(frame_stack.coords[time_coord_name].values, dtype=float)
    if timestamps.ndim != 1:
        raise ValueError(
            f"Coordinate '{time_coord_name}' must be 1D; got shape {timestamps.shape}."
        )
    return timestamps


def summarize_frame_rate_from_stack(
    frame_stack: Any,
    time_coord_name: str = "timestamp_s",
    **kwargs: Any,
) -> Dict[str, Any]:
    """
    Convenience wrapper: read timestamps from a loaded stack and compute frame-rate diagnostics.

    This function is useful once your data are loaded into xarray and include a timestamp coordinate.
    Any additional keyword arguments are forwarded to compute_frame_intervals_from_timestamps.
    """
    timestamps = extract_timestamps_from_xarray_like(
        frame_stack=frame_stack,
        time_coord_name=time_coord_name,
    )
    return compute_frame_intervals_from_timestamps(timestamps_s=timestamps, **kwargs)
```



Generate the photon transfer curve (PTC)

For each set of replicates at a single exposure time in the dark and visible light stacks, calculate the following:

  • Temporal mean image: pixel-by-pixel mean of the replicates. \(\bar{I}_{x,y} = \frac{1}{N} \sum_{i=1}^{N} I_{i,x,y}\)
  • Temporal variance image: pixel-by-pixel variance of the replicates. \(\sigma^2_{x,y} = \frac{1}{N-1} \sum_{i=1}^{N} (I_{i,x,y} - \bar{I}_{x,y})^2\)
  • Temporal mean raw signal (\(S_{\text{raw}}\)): spatial average (across all pixels) of the temporal mean image
  • Temporal variance signal (\(\sigma^2_{\text{temporal}}\)): spatial average of temporal variance image

Plot \(S_{\text{raw}}\) vs. exposure time for the dark stack. Extrapolate to 0 exposure time. This is the temporal mean dark/bias signal (\(S_{\text{dark}}\)). Note that other sources such as Teledyne estimate the bias signal by taking an image at the shortest exposure time.

Subtract the temporal mean dark/bias signal from list of temporal mean raw signals from the visible light stack (\(S_{\text{DN}} = S_{\text{raw}} - S_{\text{dark}}\)). This is the temporal mean signal (\(S_{\text{DN}}\)).

Generate a log-log plot of the temporal variance signals (\(\sigma_{\text{temporal}}\)) vs. temporal mean signals (\(S_{\text{DN}}\)) for the visible light stack.




Note that the statistical implmenetation of the PTC can vary across the detector community.


The approach used by Teledyne

  • captures only two images per exposure time with flat-field illumination
  • calculates a difference between these images to remove isolate the temporal noise from the incident light signal and spatial pixel-to-pixel signal variations that do not depend on time
  • calculates the spatial variance across pixels of the difference image
  • calculates the noise variance by dividing the spatial variance by 2 to correct for the statistical effects of subtraction.

While this approach priorizes speed and low data volume, it relies heavily on perfect uniformity of the flat-field illumination while calculating the spatial variance across pixels in the difference image.


The statistical uncertainty can be reduced by capturing three or more replicate images per exposure time. James R. Janesick calculated multiple difference images at each exposure time to combat the drift in sensor temperature and light source intensity that occurred during the long read-out times for CCD cameras in the late 1900s. As readout speeds are much faster in current cameras, the approach used at SST-1 RSoXS calculates a simple variance for each pixel across time. In addition to improved statistics,

If desired, a comparison across the difference-image approach and pixel-by-pixel temporal variance approach can be performed for data that was collected with reliably flat-field illlumination.

Code

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

from typing import Dict, List, Optional, Sequence, Tuple

import matplotlib.pyplot as plt
import numpy as np


def _validate_stack_3d(image_stack: np.ndarray, min_frames: int = 2) -> np.ndarray:
    """
    Validate and normalize a 3D image stack to float ndarray form.

    Parameters
    ----------
    image_stack : np.ndarray
        Input image stack with expected shape (n_frames, n_rows, n_cols).
    min_frames : int
        Minimum number of frames required by downstream statistics.

    Returns
    -------
    np.ndarray
        Float view/copy of the validated stack.

    Raises
    ------
    ValueError
        If input is not 3D or does not satisfy the minimum frame requirement.
    """
    arr = np.asarray(image_stack)
    if arr.ndim != 3:
        raise ValueError(
            f"image_stack must be 3D [frame, row, col]; received ndim={arr.ndim}."
        )
    if arr.shape[0] < min_frames:
        raise ValueError(
            f"image_stack requires at least {min_frames} frames; got {arr.shape[0]}."
        )
    return arr.astype(float, copy=False)


def calculate_pixel_metrics(image_stack: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculate temporal per-pixel mean and variance for a replicate stack.

    Parameters
    ----------
    image_stack : np.ndarray
        3D NumPy array with shape (n_frames, n_rows, n_cols).

    Returns
    -------
    Tuple[np.ndarray, np.ndarray]
        mean_image_2d, var_image_2d where var uses Bessel correction (ddof=1).

    Notes
    -----
    This N>2 temporal-variance implementation is the statistically robust workflow described in your
    notebook text. Because variance is computed through time for each individual pixel, fixed spatial
    pattern terms are not folded into the temporal variance estimate at this stage.
    """
    stack = _validate_stack_3d(image_stack, min_frames=2)
    mean_image_2d = np.mean(stack, axis=0)
    var_image_2d = np.var(stack, axis=0, ddof=1)
    return mean_image_2d, var_image_2d


def reduce_to_roi_scalars(
    mean_image_2d: np.ndarray,
    var_image_2d: np.ndarray,
    roi: Optional[Tuple[int, int, int, int]] = None,
) -> Tuple[float, float]:
    """
    Reduce 2D mean/variance maps to scalar ROI summaries.

    Parameters
    ----------
    mean_image_2d : np.ndarray
        2D temporal mean map.
    var_image_2d : np.ndarray
        2D temporal variance map.
    roi : tuple, optional
        ROI in (row_start, row_end, col_start, col_end) format. If None, uses full frame.

    Returns
    -------
    Tuple[float, float]
        raw_scalar_mean, raw_scalar_variance (both in DN units).
    """
    mean_map = np.asarray(mean_image_2d, dtype=float)
    var_map = np.asarray(var_image_2d, dtype=float)
    if mean_map.shape != var_map.shape:
        raise ValueError("mean_image_2d and var_image_2d must have the same shape.")

    if roi is None:
        mean_roi = mean_map
        var_roi = var_map
    else:
        r_start, r_end, c_start, c_end = roi
        mean_roi = mean_map[r_start:r_end, c_start:c_end]
        var_roi = var_map[r_start:r_end, c_start:c_end]

    raw_scalar_mean = float(np.mean(mean_roi))
    raw_scalar_variance = float(np.mean(var_roi))
    return raw_scalar_mean, raw_scalar_variance


def extrapolate_dark_bias(
    dark_exposure_times: Sequence[float],
    dark_raw_means: Sequence[float],
    return_fit_details: bool = False,
) -> Dict[str, float] | float:
    """
    Extrapolate dark/bias baseline at t=0 from dark mean signal vs exposure time.

    Parameters
    ----------
    dark_exposure_times : sequence of float
        Exposure times (seconds) for dark measurements.
    dark_raw_means : sequence of float
        ROI mean dark signals (DN) corresponding to each exposure time.
    return_fit_details : bool, optional
        If True, returns fit slope/intercept and basic diagnostics; otherwise returns intercept only.

    Returns
    -------
    float or Dict[str, float]
        Either S_dark intercept (DN) or fit dictionary containing intercept/slope/r2.
    """
    t = np.asarray(dark_exposure_times, dtype=float)
    y = np.asarray(dark_raw_means, dtype=float)

    if t.size != y.size or t.size < 2:
        raise ValueError(
            "dark_exposure_times and dark_raw_means must have same length and at least 2 points."
        )

    # Fit y = a*t + b using ordinary least squares.
    A = np.vstack([t, np.ones_like(t)]).T
    slope, intercept = np.linalg.lstsq(A, y, rcond=None)[0]
    y_hat = slope * t + intercept

    ss_res = float(np.sum((y - y_hat) ** 2))
    ss_tot = float(np.sum((y - np.mean(y)) ** 2))
    r2 = float(1.0 - ss_res / ss_tot) if ss_tot > 0 else float("nan")

    if return_fit_details:
        return {
            "dark_slope_dn_per_s": float(slope),
            "dark_intercept_dn": float(intercept),
            "dark_fit_r2": r2,
        }
    return float(intercept)


def subtract_dark_baseline(light_raw_means: Sequence[float], s_dark: float) -> np.ndarray:
    """
    Subtract extrapolated dark baseline from raw light means.

    Parameters
    ----------
    light_raw_means : sequence of float
        Raw ROI means from light-stack exposure series (DN).
    s_dark : float
        Dark/bias baseline at t=0 (DN).

    Returns
    -------
    np.ndarray
        Baseline-corrected signal array S_DN.

    Notes
    -----
    Negative values can occur in very low-signal conditions after subtraction; downstream analysis
    functions filter non-positive points where logarithms are required.
    """
    return np.asarray(light_raw_means, dtype=float) - float(s_dark)


def generate_log_log_plot(
    s_dn: np.ndarray,
    temporal_variances: np.ndarray,
    y_mode: str = "noise",
    ax: Optional[plt.Axes] = None,
) -> plt.Axes:
    """
    Plot PTC on log-log axes with explicit Y-axis mode.

    Parameters
    ----------
    s_dn : np.ndarray
        Baseline-corrected mean signal values (DN).
    temporal_variances : np.ndarray
        Temporal variance values (DN^2).
    y_mode : str
        'noise' plots sigma vs signal (shot-noise slope target ~0.5).
        'variance' plots sigma^2 vs signal (shot-noise slope target ~1.0).
    ax : matplotlib.axes.Axes, optional
        Existing axis for plotting; if omitted, a new figure/axis is created.

    Returns
    -------
    matplotlib.axes.Axes
        Axis containing the plot.
    """
    s = np.asarray(s_dn, dtype=float)
    v = np.asarray(temporal_variances, dtype=float)
    valid = (s > 0) & (v > 0)
    s = s[valid]
    v = v[valid]
    if s.size < 2:
        raise ValueError("Need at least two positive points to create log-log plot.")

    if ax is None:
        _, ax = plt.subplots(figsize=(8, 6))

    if y_mode == "noise":
        y = np.sqrt(v)
        ylabel = "Temporal RMS Noise ($\\sigma_S$) [DN]"
        title = "Photon Transfer Curve (Noise vs Signal)"
    elif y_mode == "variance":
        y = v
        ylabel = "Temporal Variance ($\\sigma_S^2$) [DN$^2$]"
        title = "Photon Transfer Curve (Variance vs Signal)"
    else:
        raise ValueError("y_mode must be either 'noise' or 'variance'.")

    ax.loglog(s, y, "o-", label="Measured sensor profile")
    ax.set_xlabel("True Mean Signal ($S_{DN}$)")
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    ax.grid(True, which="both", ls="--", alpha=0.5)
    ax.legend()
    return ax


def isolate_shot_noise_regime(
    s_dn: np.ndarray,
    temporal_variances: np.ndarray,
    y_mode: str = "noise",
    slope_tolerance: float = 0.15,
    min_points: int = 4,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Isolate shot-noise regime using local slope on log-log PTC coordinates.

    Parameters
    ----------
    s_dn : np.ndarray
        Baseline-corrected mean signal values (DN).
    temporal_variances : np.ndarray
        Temporal variance values (DN^2).
    y_mode : str
        'noise' expects target slope 0.5, 'variance' expects target slope 1.0.
    slope_tolerance : float
        Allowed absolute deviation around the target slope.
    min_points : int
        Minimum acceptable number of selected points.

    Returns
    -------
    Tuple[np.ndarray, np.ndarray, np.ndarray]
        s_dn_shot, var_shot, shot_indices_in_sorted_valid_arrays.
    """
    s = np.asarray(s_dn, dtype=float)
    v = np.asarray(temporal_variances, dtype=float)
    valid = (s > 0) & (v > 0)
    s = s[valid]
    v = v[valid]

    if s.size < max(3, min_points):
        raise ValueError("Insufficient valid positive points to isolate shot-noise regime.")

    sort_idx = np.argsort(s)
    s_sorted = s[sort_idx]
    v_sorted = v[sort_idx]

    log_s = np.log10(s_sorted)
    if y_mode == "noise":
        log_y = np.log10(np.sqrt(v_sorted))
        target_slope = 0.5
    elif y_mode == "variance":
        log_y = np.log10(v_sorted)
        target_slope = 1.0
    else:
        raise ValueError("y_mode must be either 'noise' or 'variance'.")

    local_slopes = np.gradient(log_y, log_s)
    shot_mask = np.abs(local_slopes - target_slope) <= slope_tolerance
    shot_indices = np.where(shot_mask)[0]

    if shot_indices.size < min_points:
        raise ValueError(
            "Could not isolate a robust shot-noise region. Adjust tolerance or inspect data quality."
        )

    return s_sorted[shot_indices], v_sorted[shot_indices], shot_indices


def calculate_conversion_gain(s_dn_shot: np.ndarray, var_shot: np.ndarray) -> Dict[str, float]:
    """
    Calculate conversion gain K (e-/DN) from shot-noise regime using origin-constrained fit.

    Parameters
    ----------
    s_dn_shot : np.ndarray
        Shot-noise-region signal values (DN).
    var_shot : np.ndarray
        Shot-noise-region temporal variances (DN^2).

    Returns
    -------
    Dict[str, float]
        Fit outputs including slope (1/K), gain K, and residual RMSE in DN^2.
    """
    x = np.asarray(s_dn_shot, dtype=float)
    y = np.asarray(var_shot, dtype=float)
    valid = (x > 0) & (y > 0)
    x = x[valid]
    y = y[valid]

    if x.size < 2:
        raise ValueError("Need at least two valid shot-noise points for gain estimation.")

    denom = float(np.sum(x ** 2))
    if denom <= 0.0:
        raise ValueError("Degenerate input: sum(S_DN^2) is not positive.")

    slope = float(np.sum(x * y) / denom)  # slope = 1/K in variance-vs-signal representation
    if slope <= 0.0:
        raise ValueError("Non-positive slope encountered; gain cannot be resolved.")

    y_fit = slope * x
    rmse = float(np.sqrt(np.mean((y - y_fit) ** 2)))
    gain_k = 1.0 / slope

    return {
        "linear_linear_slope_1_over_K": slope,
        "conversion_gain_K_e_per_DN": float(gain_k),
        "shot_fit_rmse_dn2": rmse,
    }
```



Conversion gain (K)

Physics: This is the conversion factor in electrons per DN (raw digital units) that converts the ditital signal to the number of electrons. The gain is calculated as \(K = \frac{S_{\text{DN}}}{\sigma^2_{\text{temporal}}}\), which is derived from the Poisson physics in the shot-noise regime.


Calculation: Visually or algorithmically define the boundaries of the Poisson shot noise regime in the PTC where the slope of the log-log plot is 1/2. Plot the shot-noise regime portion of the (\(\sigma^2_{\text{temporal}}\)) vs. (\(S_{\text{DN}}\)) data on a linear-linear scale. The y-intercept is through the origin. Perform a linear regression to find the slope with the fixed y = 0 intercept. K = 1/slope.

In theory, K also could be calculated as the intercept of the log-log representation of the shot-noise regime (\(log(\sigma) = \frac{1}{2} log(S_{DN}) + log(\frac{1}{\sqrt{K}})\)). However, performing the linear regression on the linear-linear data is more robust and will have less uncertainty. The linear regression assumes that uncertainties are homeoscedastic; i.e., all datapoints have equal uncertainty. Taking a logarithm would stretch the uncertainties at low signals while compressing uncertainties at higher signals, and more weight would be placed on the low-signal data points. Similarly, calculating the logarithm compress the dynamic range, which can make the analysis sensitive to rounding errors. The regression on the linear-linear data also eliminates a degree of freedom by fixing the intercept to 0 and makes the analysis more stable.

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

from typing import Dict, List, Optional, Tuple

import numpy as np


def run_characterization_pipeline(
    raw_dark_stacks: List[np.ndarray],
    dark_times_s: List[float],
    raw_light_stacks: List[np.ndarray],
    roi: Optional[Tuple[int, int, int, int]] = None,
    ptc_y_mode: str = "noise",
    shot_slope_tolerance: float = 0.15,
) -> Dict[str, object]:
    """
    Execute a high-level PTC and gain-estimation workflow from dark and light replicate stacks.

    Parameters
    ----------
    raw_dark_stacks : list of np.ndarray
        Dark stacks, one 3D array per exposure condition.
    dark_times_s : list of float
        Exposure times corresponding to raw_dark_stacks.
    raw_light_stacks : list of np.ndarray
        Light stacks, one 3D array per exposure condition.
    roi : tuple, optional
        (row_start, row_end, col_start, col_end). If None, full frame is used.
    ptc_y_mode : str
        'noise' or 'variance' for log-log PTC representation.
    shot_slope_tolerance : float
        Tolerance used for shot-noise slope-based point selection.

    Returns
    -------
    Dict[str, object]
        Comprehensive dictionary containing reduced arrays, fit results, and diagnostics.

    Notes
    -----
    This function does not load data from disk. It expects arrays already loaded into memory (for
    example from xarray containers). That separation keeps the analysis code deterministic and testable.
    """
    if len(raw_dark_stacks) != len(dark_times_s):
        raise ValueError("raw_dark_stacks and dark_times_s must have the same length.")
    if len(raw_light_stacks) < 2:
        raise ValueError("At least two light exposure steps are recommended for PTC analysis.")

    # 1) Reduce dark stacks to ROI scalar means/variances and estimate t=0 dark baseline.
    dark_means_dn: List[float] = []
    dark_vars_dn2: List[float] = []
    for d_stack in raw_dark_stacks:
        m_2d, v_2d = calculate_pixel_metrics(d_stack)
        s_mean, s_var = reduce_to_roi_scalars(m_2d, v_2d, roi=roi)
        dark_means_dn.append(s_mean)
        dark_vars_dn2.append(s_var)

    dark_fit = extrapolate_dark_bias(
        dark_exposure_times=dark_times_s,
        dark_raw_means=dark_means_dn,
        return_fit_details=True,
    )
    s_dark_dn = float(dark_fit["dark_intercept_dn"]

    )
    # 2) Reduce light stacks to ROI scalar means/variances.
    light_raw_means_dn: List[float] = []
    light_vars_dn2: List[float] = []
    for l_stack in raw_light_stacks:
        m_2d, v_2d = calculate_pixel_metrics(l_stack)
        s_mean, s_var = reduce_to_roi_scalars(m_2d, v_2d, roi=roi)
        light_raw_means_dn.append(s_mean)
        light_vars_dn2.append(s_var)

    # 3) Subtract dark baseline to construct true mean signal S_DN.
    s_dn = subtract_dark_baseline(light_raw_means=light_raw_means_dn, s_dark=s_dark_dn)
    temporal_variances = np.asarray(light_vars_dn2, dtype=float)

    # 4) Build PTC plot using explicit y-mode semantics (noise slope 0.5 vs variance slope 1.0).
    ax = generate_log_log_plot(
        s_dn=s_dn,
        temporal_variances=temporal_variances,
        y_mode=ptc_y_mode,
    )
    _ = ax  # Keep explicit reference for notebook workflows that inspect returned axes.


    # 5) Isolate shot-noise window from local log-log slope.
    s_dn_shot, var_shot, shot_idx = isolate_shot_noise_regime(
        s_dn=s_dn,
        temporal_variances=temporal_variances,
        y_mode=ptc_y_mode,
        slope_tolerance=shot_slope_tolerance,
    )

    # 6) Compute origin-constrained linear fit in variance-vs-signal space to extract gain.
    gain_results = calculate_conversion_gain(s_dn_shot=s_dn_shot, var_shot=var_shot)

    return {
        "dark_fit": dark_fit,
        "s_dark_dn": s_dark_dn,
        "dark_means_dn": np.asarray(dark_means_dn, dtype=float),
        "dark_vars_dn2": np.asarray(dark_vars_dn2, dtype=float),
        "light_raw_means_dn": np.asarray(light_raw_means_dn, dtype=float),
        "light_vars_dn2": temporal_variances,
        "s_dn": np.asarray(s_dn, dtype=float),
        "shot_indices_sorted_valid": np.asarray(shot_idx, dtype=int),
        "s_dn_shot": np.asarray(s_dn_shot, dtype=float),
        "var_shot_dn2": np.asarray(var_shot, dtype=float),
        "gain_results": gain_results,
    }
```



Fixed pattern noise

Physics:

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

import numpy as np
from typing import Any, Dict, List, Optional, Sequence, Tuple


def _validate_stack_3d(image_stack: np.ndarray, min_frames: int = 2) -> np.ndarray:
    """
    Internal validator for image stacks with shape (frames, rows, cols).
    """
    arr = np.asarray(image_stack)
    if arr.ndim != 3:
        raise ValueError(
            f"image_stack must be 3D with shape (frames, rows, cols); got ndim={arr.ndim}."
        )
    if arr.shape[0] < min_frames:
        raise ValueError(
            f"image_stack must contain at least {min_frames} frames; got {arr.shape[0]}."
        )
    return arr.astype(float, copy=False)


def compute_mean_variance_maps(image_stack: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute temporal mean and temporal variance maps for one stack.

    Parameters
    ----------
    image_stack : np.ndarray
        3D array (n_frames, n_rows, n_cols) acquired at identical settings.

    Returns
    -------
    Tuple[np.ndarray, np.ndarray]
        mean_map, var_map where:
        - mean_map is temporal mean per pixel.
        - var_map is Bessel-corrected temporal variance per pixel (ddof=1).

    Notes
    -----
    The temporal variance map naturally excludes static fixed-pattern structure because each pixel is
    compared only to itself through time. This is useful in your N>2 replicate workflow.
    """
    stack = _validate_stack_3d(image_stack, min_frames=2)
    mean_map = np.mean(stack, axis=0)
    var_map = np.var(stack, axis=0, ddof=1)
    return mean_map, var_map


def estimate_fpn_from_maps(
    mean_map: np.ndarray,
    temporal_var_map: np.ndarray,
    roi: Optional[Tuple[int, int, int, int]] = None,
) -> Dict[str, float]:
    """
    Estimate fixed-pattern-noise (FPN) variance by separating spatial vs temporal contributions.

    Parameters
    ----------
    mean_map : np.ndarray
        Temporal mean image (2D).
    temporal_var_map : np.ndarray
        Temporal variance image (2D) computed from the same stack.
    roi : tuple, optional
        (row_start, row_end, col_start, col_end). If None, use full frame.

    Returns
    -------
    Dict[str, float]
        Spatial and temporal summary metrics in DN units:
        - spatial_variance_total_dn2
        - temporal_variance_mean_dn2
        - fpn_variance_dn2
        - fpn_rms_dn

    Explanation
    -----------
    In a flat-field step, the spatial variance of the temporal mean image contains both true spatial
    structure and residual temporal terms. A practical FPN estimate is:
        FPN_variance ~= Var_spatial(mean_map) - Mean_temporal_variance
    clipped at zero to avoid non-physical negative values from finite-sample noise.
    """
    mean_map = np.asarray(mean_map, dtype=float)
    temporal_var_map = np.asarray(temporal_var_map, dtype=float)
    if mean_map.shape != temporal_var_map.shape:
        raise ValueError(
            "mean_map and temporal_var_map must have the same shape."
        )

    if roi is None:
        r0, r1, c0, c1 = 0, mean_map.shape[0], 0, mean_map.shape[1]
    else:
        r0, r1, c0, c1 = roi

    mean_roi = mean_map[r0:r1, c0:c1]
    tvar_roi = temporal_var_map[r0:r1, c0:c1]

    spatial_variance_total = float(np.var(mean_roi, ddof=1))
    temporal_variance_mean = float(np.mean(tvar_roi))
    fpn_variance = max(spatial_variance_total - temporal_variance_mean, 0.0)

    return {
        "spatial_variance_total_dn2": spatial_variance_total,
        "temporal_variance_mean_dn2": temporal_variance_mean,
        "fpn_variance_dn2": float(fpn_variance),
        "fpn_rms_dn": float(np.sqrt(fpn_variance)),
    }


def estimate_offset_and_prnu_maps(
    dark_stack: np.ndarray,
    flat_stack: np.ndarray,
) -> Dict[str, np.ndarray]:
    """
    Estimate spatial offset map and pixel-response non-uniformity (PRNU) map.

    Parameters
    ----------
    dark_stack : np.ndarray
        Dark replicates at a given exposure setting, shape (frames, rows, cols).
    flat_stack : np.ndarray
        Flat-field replicates at same exposure setting, shape (frames, rows, cols).

    Returns
    -------
    Dict[str, np.ndarray]
        - offset_map_dn: spatial offset/bias map from temporal mean dark image.
        - signal_map_dn: dark-corrected flat-field signal map.
        - prnu_map_fraction: normalized gain non-uniformity map, (signal/mean_signal - 1).

    Notes
    -----
    PRNU is one operational measure of spatial pixel gain noise. For robust PRNU estimation, use a
    flat-field level within linear response (not near noise floor, not near saturation).
    """
    dark_mean, _ = compute_mean_variance_maps(dark_stack)
    flat_mean, _ = compute_mean_variance_maps(flat_stack)

    signal_map = flat_mean - dark_mean
    mean_signal = float(np.mean(signal_map))
    if mean_signal <= 0:
        raise ValueError(
            "Mean dark-corrected flat signal is non-positive; cannot compute PRNU."
        )
    prnu_map = (signal_map / mean_signal) - 1.0

    return {
        "offset_map_dn": dark_mean,
        "signal_map_dn": signal_map,
        "prnu_map_fraction": prnu_map,
    }


def compute_pair_difference_variance(
    image_stack: np.ndarray,
    roi: Optional[Tuple[int, int, int, int]] = None,
) -> Dict[str, float]:
    """
    Compute temporal variance from pairwise frame differencing (Janesick/Teledyne style).

    Parameters
    ----------
    image_stack : np.ndarray
        3D stack for one exposure level. Frames are paired in order: (0,1), (2,3), ...
    roi : tuple, optional
        (row_start, row_end, col_start, col_end). If None, use full frame.

    Returns
    -------
    Dict[str, float]
        - pair_count
        - variance_diff_image_dn2
        - temporal_variance_single_frame_dn2 (variance_diff/2)

    Notes
    -----
    This is useful for direct comparison with pair-difference methods in the literature. If an odd
    number of frames is supplied, the last frame is ignored.
    """
    stack = _validate_stack_3d(image_stack, min_frames=2)
    n_pairs = stack.shape[0] // 2
    if n_pairs < 1:
        raise ValueError("At least two frames are required to form one pair.")

    a = stack[0 : 2 * n_pairs : 2]
    b = stack[1 : 2 * n_pairs : 2]
    diffs = a - b

    if roi is None:
        diff_roi = diffs
    else:
        r0, r1, c0, c1 = roi
        diff_roi = diffs[:, r0:r1, c0:c1]

    variance_diff = float(np.var(diff_roi, ddof=1))
    temporal_variance_single = variance_diff / 2.0

    return {
        "pair_count": int(n_pairs),
        "variance_diff_image_dn2": variance_diff,
        "temporal_variance_single_frame_dn2": temporal_variance_single,
    }


def summarize_fpn_over_exposure_series(
    stacks_by_exposure: Sequence[np.ndarray],
    exposure_times_s: Sequence[float],
    roi: Optional[Tuple[int, int, int, int]] = None,
) -> List[Dict[str, float]]:
    """
    Summarize FPN metrics versus exposure time for a full stack series.

    Parameters
    ----------
    stacks_by_exposure : sequence of np.ndarray
        Each element is a 3D stack captured at one exposure time.
    exposure_times_s : sequence of float
        Exposure times matching stacks_by_exposure.
    roi : tuple, optional
        Region of interest for metric reduction.

    Returns
    -------
    List[Dict[str, float]]
        One dictionary per exposure level with FPN and temporal variance summary metrics.
    """
    if len(stacks_by_exposure) != len(exposure_times_s):
        raise ValueError(
            "stacks_by_exposure and exposure_times_s must have the same length."
        )

    rows: List[Dict[str, float]] = []
    for t_s, stack in zip(exposure_times_s, stacks_by_exposure):
        mean_map, var_map = compute_mean_variance_maps(stack)
        fpn = estimate_fpn_from_maps(mean_map, var_map, roi=roi)
        rows.append({"exposure_time_s": float(t_s), **fpn})
    return rows
```



Read noise and dark current rate

Physics: Dark current is produced when thermal energy present inside the sensor generates electron-hole pairs even in the absence of any incident photons on the sensor. Dark current can be minimized by cooling the sensor. Read noise is generated while converting electrons to digital signal.

Calculation: Calculate a mean temporal signal and mean temporal variance for each set of replicates at a particular exposure time and temperature. Plot the variance vs. exposure time for every temperature set. Perform a global intercept linear regression where all lines are forced to have the same Y-intercept; this allows for a more robust calculation. The Y-intercept is the read variance, and the slope is the dark current rate. The dark noise is the dark rate multiplied by the exposure time.

Code
```{python}
#| code-fold: true
## Insert code here
from __future__ import annotations

import numpy as np
from typing import Dict, List, Optional, Sequence, Tuple


def fit_global_intercept_dark_model(
    exposure_times_by_temperature: Sequence[Sequence[float]],
    variances_by_temperature: Sequence[Sequence[float]],
) -> Dict[str, object]:
    """
    Fit a global-intercept linear model for dark variance vs exposure time.

    Model
    -----
    For each temperature index j and data point i:
        variance_ij = b0 + m_j * time_ij + eps_ij

    where:
      - b0 is a single, shared intercept across all temperatures (read variance in DN^2),
      - m_j is temperature-specific slope (dark variance growth rate in DN^2/s).

    Parameters
    ----------
    exposure_times_by_temperature : sequence of sequence of float
        Per-temperature exposure time arrays in seconds.
    variances_by_temperature : sequence of sequence of float
        Per-temperature temporal variance arrays in DN^2.

    Returns
    -------
    Dict[str, object]
        - read_variance_dn2: shared intercept b0
        - dark_variance_rate_dn2_per_s_by_temperature: list of m_j
        - fitted_values_by_temperature: list of fitted arrays
        - residual_std_dn2: pooled residual standard deviation
        - design_matrix_rank: least-squares matrix rank

    Notes
    -----
    This directly implements the "global intercept" strategy described in your notebook text, giving
    a robust read-noise estimate that is less sensitive to per-temperature slope differences.
    """
    if len(exposure_times_by_temperature) != len(variances_by_temperature):
        raise ValueError(
            "exposure_times_by_temperature and variances_by_temperature must have equal length."
        )
    if len(exposure_times_by_temperature) < 1:
        raise ValueError("At least one temperature series is required.")

    n_groups = len(exposure_times_by_temperature)
    y_blocks: List[np.ndarray] = []
    X_blocks: List[np.ndarray] = []

    for j, (times_j, vars_j) in enumerate(
        zip(exposure_times_by_temperature, variances_by_temperature)
    ):
        t = np.asarray(times_j, dtype=float)
        v = np.asarray(vars_j, dtype=float)
        if t.size != v.size:
            raise ValueError(f"Group {j}: time and variance arrays must have same length.")
        if t.size < 2:
            raise ValueError(f"Group {j}: need at least 2 points for a slope.")

        # Build block design matrix: one shared intercept + one slope column per group.
        # For rows from group j, only slope column j is populated with time values.
        block = np.zeros((t.size, 1 + n_groups), dtype=float)
        block[:, 0] = 1.0
        block[:, 1 + j] = t

        X_blocks.append(block)
        y_blocks.append(v)

    X = np.vstack(X_blocks)
    y = np.concatenate(y_blocks)

    beta, residuals, rank, _ = np.linalg.lstsq(X, y, rcond=None)
    intercept = float(beta[0])
    slopes = [float(s) for s in beta[1:]]

    fitted_values_by_group: List[np.ndarray] = []
    for j, times_j in enumerate(exposure_times_by_temperature):
        t = np.asarray(times_j, dtype=float)
        fitted_values_by_group.append(intercept + slopes[j] * t)

    if residuals.size > 0:
        # Least squares residuals returns SSE for overdetermined systems.
        dof = max(X.shape[0] - rank, 1)
        residual_std = float(np.sqrt(residuals[0] / dof))
    else:
        residual_std = float("nan")

    return {
        "read_variance_dn2": intercept,
        "dark_variance_rate_dn2_per_s_by_temperature": slopes,
        "fitted_values_by_temperature": fitted_values_by_group,
        "residual_std_dn2": residual_std,
        "design_matrix_rank": int(rank),
    }


def convert_dark_model_to_electron_units(
    dark_model_result: Dict[str, object],
    gain_e_per_dn: float,
) -> Dict[str, object]:
    """
    Convert dark-model outputs from DN units to electron units.

    Parameters
    ----------
    dark_model_result : dict
        Output from fit_global_intercept_dark_model.
    gain_e_per_dn : float
        Conversion gain K in e-/DN.

    Returns
    -------
    Dict[str, object]
        Extended dictionary with electron-scale metrics:
        - read_variance_e2
        - read_noise_e_rms
        - dark_current_rate_e_per_s_by_temperature

    Notes
    -----
    Variance transforms with K^2, while RMS transforms with K.
    """
    if gain_e_per_dn <= 0:
        raise ValueError("gain_e_per_dn must be positive.")

    read_var_dn2 = float(dark_model_result["read_variance_dn2"])
    dark_rate_dn2_per_s = [
        float(x) for x in dark_model_result["dark_variance_rate_dn2_per_s_by_temperature"]
    ]

    read_var_e2 = read_var_dn2 * (gain_e_per_dn ** 2)
    read_noise_e = np.sqrt(max(read_var_e2, 0.0))
    dark_rate_e_per_s = [r * (gain_e_per_dn ** 2) for r in dark_rate_dn2_per_s]

    return {
        **dark_model_result,
        "read_variance_e2": float(read_var_e2),
        "read_noise_e_rms": float(read_noise_e),
        "dark_current_rate_e_per_s_by_temperature": [float(x) for x in dark_rate_e_per_s],
    }


def estimate_dark_noise_vs_exposure(
    dark_current_rate_e_per_s: float,
    exposure_times_s: Sequence[float],
) -> Dict[str, np.ndarray]:
    """
    Compute dark noise metrics at specified exposure times.

    Parameters
    ----------
    dark_current_rate_e_per_s : float
        Dark current expressed as variance growth rate in e^2/s.
    exposure_times_s : sequence of float
        Exposure times in seconds.

    Returns
    -------
    Dict[str, np.ndarray]
        - dark_variance_e2
        - dark_noise_e_rms

    Notes
    -----
    If dark current follows Poisson statistics, variance grows linearly with time and
    RMS dark noise grows as square-root of exposure time.
    """
    t = np.asarray(exposure_times_s, dtype=float)
    if np.any(t < 0):
        raise ValueError("Exposure times must be non-negative.")
    var_e2 = dark_current_rate_e_per_s * t
    noise_e = np.sqrt(np.clip(var_e2, 0.0, None))
    return {
        "dark_variance_e2": var_e2,
        "dark_noise_e_rms": noise_e,
    }


def analyze_linearity_full_well_dynamic_range(
    exposure_times_s: Sequence[float],
    mean_signal_dn: Sequence[float],
    read_noise_dn_rms: float,
    gain_e_per_dn: Optional[float] = None,
    linearity_tolerance_percent: float = 1.0,
) -> Dict[str, object]:
    """
    Analyze linearity, full-well estimate, and dynamic range from an exposure sweep.

    Parameters
    ----------
    exposure_times_s : sequence of float
        Exposure times in seconds for light-stack mean signal values.
    mean_signal_dn : sequence of float
        Mean signal values in DN (typically dark-corrected S_DN).
    read_noise_dn_rms : float
        Read noise RMS in DN.
    gain_e_per_dn : float, optional
        Conversion gain for electron-unit outputs.
    linearity_tolerance_percent : float, optional
        Maximum absolute percent deviation allowed for the "linear operating region."

    Returns
    -------
    Dict[str, object]
        Dictionary with linear fit, deviation profile, estimated linear window, full-well, and
        dynamic range in both ratio and dB.

    Strategy
    --------
    1) Fit a provisional linear model in a conservative pre-saturation region (5% to 70% of max).
    2) Compute percent deviation from model across all points.
    3) Define linear points as those within tolerance.
    4) Report full-well from maximum observed signal (or optionally user-refined cutoff).
    """
    t = np.asarray(exposure_times_s, dtype=float)
    s = np.asarray(mean_signal_dn, dtype=float)

    if t.size != s.size or t.size < 3:
        raise ValueError(
            "exposure_times_s and mean_signal_dn must have equal length and at least 3 points."
        )
    if np.any(t <= 0):
        raise ValueError("Exposure times must be positive for linearity analysis.")
    if read_noise_dn_rms <= 0:
        raise ValueError("read_noise_dn_rms must be positive.")

    s_max = float(np.max(s))
    if s_max <= 0:
        raise ValueError("mean_signal_dn must contain positive values.")

    # Build a conservative calibration window away from floor and saturation.
    fit_mask = (s >= 0.05 * s_max) & (s <= 0.70 * s_max)
    if np.sum(fit_mask) < 2:
        # Fallback: use first half of points if the heuristic window is too small.
        idx = np.argsort(t)
        take = max(2, t.size // 2)
        fit_mask = np.zeros_like(t, dtype=bool)
        fit_mask[idx[:take]] = True

    # Fit y = a*t + b in the pre-saturation region.
    A = np.vstack([t[fit_mask], np.ones(np.sum(fit_mask))]).T
    a, b = np.linalg.lstsq(A, s[fit_mask], rcond=None)[0]

    predicted = a * t + b
    denom = np.maximum(np.abs(predicted), 1e-12)
    pct_dev = 100.0 * (s - predicted) / denom

    linear_mask = np.abs(pct_dev) <= linearity_tolerance_percent

    # Full-well estimate from observed maximum signal (in DN).
    full_well_dn = s_max
    dynamic_range_ratio_dn = full_well_dn / read_noise_dn_rms
    dynamic_range_db = 20.0 * np.log10(dynamic_range_ratio_dn)

    result: Dict[str, object] = {
        "linearity_fit_slope_dn_per_s": float(a),
        "linearity_fit_intercept_dn": float(b),
        "predicted_signal_dn": predicted,
        "percent_deviation_from_linear_fit": pct_dev,
        "linearity_mask": linear_mask,
        "linearity_tolerance_percent": float(linearity_tolerance_percent),
        "full_well_dn": float(full_well_dn),
        "dynamic_range_ratio_dn": float(dynamic_range_ratio_dn),
        "dynamic_range_db": float(dynamic_range_db),
    }

    if np.any(linear_mask):
        linear_times = t[linear_mask]
        result["linear_region_time_start_s"] = float(np.min(linear_times))
        result["linear_region_time_end_s"] = float(np.max(linear_times))
    else:
        result["linear_region_time_start_s"] = float("nan")
        result["linear_region_time_end_s"] = float("nan")

    if gain_e_per_dn is not None:
        if gain_e_per_dn <= 0:
            raise ValueError("gain_e_per_dn must be positive when provided.")
        full_well_e = full_well_dn * gain_e_per_dn
        read_noise_e = read_noise_dn_rms * gain_e_per_dn
        dynamic_range_ratio_e = full_well_e / read_noise_e
        dynamic_range_db_e = 20.0 * np.log10(dynamic_range_ratio_e)
        result.update({
            "full_well_e": float(full_well_e),
            "read_noise_e_rms": float(read_noise_e),
            "dynamic_range_ratio_e": float(dynamic_range_ratio_e),
            "dynamic_range_db_e": float(dynamic_range_db_e),
        })

    return result
```

Point spread function



Quantum performance