2. SwissSPAD2 Data Processing#

Apply detector-specific processing to a real SwissSPAD2 acquisition with PyFLI.

The processing pipeline preserves a raw reference acquisition and separately loads a processed dataset using explicitly configured corrections. This makes it possible to compare detector output before and after processing.

This notebook walks through:

  • configuring the SwissSPAD2 processing pipeline,

  • loading an unprocessed reference acquisition,

  • applying pile-up correction when requested,

  • applying background subtraction and hot-pixel correction to HDF5 folder acquisitions when requested,

  • applying temporal alignment and folding when requested,

  • comparing raw and processed intensity images,

  • comparing raw and processed temporal responses,

  • applying an optional PyFLI mask, and

  • reviewing the complete processing metadata.

import os
from pathlib import Path
from pprint import pprint

import matplotlib.pyplot as plt
import numpy as np

from pyfli.data_cc import DataPreprocessing
from pyfli.io import Detector

2.1. Configure processing#

Machine-specific paths and acquisition-specific processing choices are read from environment variables.

Required variables:

  • PYFLI_SS2_DATA_PATH

  • PYFLI_SS2_BIT_SIZE

Processing switches default to disabled. At least one processing operation must be explicitly enabled.

Supported switches:

  • PYFLI_SS2_PILE_UP

  • PYFLI_SS2_FOLD

  • PYFLI_SS2_SUB_BG

  • PYFLI_SS2_HOT_PIXEL

  • PYFLI_SS2_MAKE_HP_MAP

Optional paths:

  • PYFLI_SS2_IRF_PATH

  • PYFLI_SS2_BG_PATH

  • PYFLI_SS2_MASK_PATH

  • PYFLI_SS2_HP_PATH

Optional acquisition and folding parameters:

  • PYFLI_SS2_GATE_COUNT

  • PYFLI_SS2_FOLD_REPETITIONS

  • PYFLI_SS2_PERIOD_BINS

  • PYFLI_SS2_PHASE_SHIFT

  • PYFLI_SS2_DETECTOR_FREQ_MHZ

  • PYFLI_SS2_LASER_FREQ_MHZ

  • PYFLI_SS2_THRESHOLD_SIGMA

  • PYFLI_SS2_FOLDER_MODE

Background subtraction and detector-derived hot-pixel correction currently apply to HDF5 folder acquisitions. Native BIN acquisitions support the shared pile-up and temporal-folding pipeline.

def require_env(name: str) -> str:
    value = os.environ.get(name)

    if value is None or not value.strip():
        raise RuntimeError(
            f"Environment variable '{name}' must be set before running this notebook."
        )

    return value.strip()


def require_path(name: str) -> Path:
    path = Path(require_env(name)).expanduser().resolve()

    if not path.exists():
        raise FileNotFoundError(
            f"Path configured by '{name}' does not exist: {path}"
        )

    return path


def optional_path(name: str) -> Path | None:
    value = os.environ.get(name)

    if value is None or not value.strip():
        return None

    path = Path(value.strip()).expanduser().resolve()

    if not path.exists():
        raise FileNotFoundError(
            f"Path configured by '{name}' does not exist: {path}"
        )

    return path


def env_bool(
    name: str,
    default: bool = False,
) -> bool:
    value = os.environ.get(name)

    if value is None or not value.strip():
        return default

    normalized = value.strip().lower()

    true_values = {
        "1",
        "true",
        "yes",
        "on",
    }

    false_values = {
        "0",
        "false",
        "no",
        "off",
    }

    if normalized in true_values:
        return True

    if normalized in false_values:
        return False

    raise ValueError(
        f"Environment variable '{name}' must be one of "
        f"{sorted(true_values | false_values)}, got '{value}'."
    )


def optional_int(name: str) -> int | None:
    value = os.environ.get(name)

    if value is None or not value.strip():
        return None

    return int(value.strip())


def optional_float(name: str) -> float | None:
    value = os.environ.get(name)

    if value is None or not value.strip():
        return None

    return float(value.strip())


DATA_PATH = require_path(
    "PYFLI_SS2_DATA_PATH"
)

BIT_SIZE = int(
    require_env(
        "PYFLI_SS2_BIT_SIZE"
    )
)

IRF_PATH = optional_path(
    "PYFLI_SS2_IRF_PATH"
)

BG_PATH = optional_path(
    "PYFLI_SS2_BG_PATH"
)

MASK_PATH = optional_path(
    "PYFLI_SS2_MASK_PATH"
)

HP_PATH = optional_path(
    "PYFLI_SS2_HP_PATH"
)

PILE_UP = env_bool(
    "PYFLI_SS2_PILE_UP",
    default=False,
)

FOLD = env_bool(
    "PYFLI_SS2_FOLD",
    default=False,
)

SUB_BG = env_bool(
    "PYFLI_SS2_SUB_BG",
    default=False,
)

HOT_PIXEL = env_bool(
    "PYFLI_SS2_HOT_PIXEL",
    default=False,
)

MAKE_HP_MAP = env_bool(
    "PYFLI_SS2_MAKE_HP_MAP",
    default=True,
)

EXPECTED_GATE_COUNT = optional_int(
    "PYFLI_SS2_GATE_COUNT"
)

FOLD_REPETITIONS = optional_int(
    "PYFLI_SS2_FOLD_REPETITIONS"
)

PERIOD_BINS = optional_int(
    "PYFLI_SS2_PERIOD_BINS"
)

PHASE_SHIFT = optional_int(
    "PYFLI_SS2_PHASE_SHIFT"
)

DETECTOR_FREQUENCY_MHZ = optional_float(
    "PYFLI_SS2_DETECTOR_FREQ_MHZ"
)

LASER_FREQUENCY_MHZ = optional_float(
    "PYFLI_SS2_LASER_FREQ_MHZ"
)

threshold_value = optional_float(
    "PYFLI_SS2_THRESHOLD_SIGMA"
)

THRESHOLD_SIGMA = (
    5.0
    if threshold_value is None
    else threshold_value
)

FOLDER_MODE = os.environ.get(
    "PYFLI_SS2_FOLDER_MODE",
    "sum",
).strip().lower()

if BIT_SIZE < 1:
    raise ValueError(
        f"PYFLI_SS2_BIT_SIZE must be positive, got {BIT_SIZE}."
    )

if THRESHOLD_SIGMA <= 0:
    raise ValueError(
        f"PYFLI_SS2_THRESHOLD_SIGMA must be positive, got {THRESHOLD_SIGMA}."
    )

if FOLDER_MODE not in {
    "sum",
    "mean",
}:
    raise ValueError(
        "PYFLI_SS2_FOLDER_MODE must be 'sum' or 'mean', "
        f"got '{FOLDER_MODE}'."
    )

if not any(
    (
        PILE_UP,
        FOLD,
        SUB_BG,
        HOT_PIXEL,
    )
):
    raise RuntimeError(
        "At least one SwissSPAD2 processing operation must be explicitly enabled "
        "through PYFLI_SS2_PILE_UP, PYFLI_SS2_FOLD, PYFLI_SS2_SUB_BG, "
        "or PYFLI_SS2_HOT_PIXEL."
    )

if SUB_BG and BG_PATH is None:
    raise RuntimeError(
        "PYFLI_SS2_BG_PATH must be configured when PYFLI_SS2_SUB_BG is enabled."
    )

if HOT_PIXEL and MAKE_HP_MAP and BG_PATH is None:
    raise RuntimeError(
        "PYFLI_SS2_BG_PATH must be configured when automatic hot-pixel "
        "mapping is enabled."
    )

if HOT_PIXEL and not MAKE_HP_MAP and HP_PATH is None:
    raise RuntimeError(
        "PYFLI_SS2_HP_PATH must be configured when hot-pixel correction is "
        "enabled with PYFLI_SS2_MAKE_HP_MAP disabled."
    )

if (
    DETECTOR_FREQUENCY_MHZ is None
) != (
    LASER_FREQUENCY_MHZ is None
):
    raise RuntimeError(
        "PYFLI_SS2_DETECTOR_FREQ_MHZ and PYFLI_SS2_LASER_FREQ_MHZ "
        "must either both be configured or both be omitted."
    )

fold_parameters_present = any(
    value is not None
    for value in (
        FOLD_REPETITIONS,
        PERIOD_BINS,
        PHASE_SHIFT,
        DETECTOR_FREQUENCY_MHZ,
        LASER_FREQUENCY_MHZ,
    )
)

if fold_parameters_present and not FOLD:
    raise RuntimeError(
        "Folding parameters were configured while PYFLI_SS2_FOLD is disabled."
    )

if EXPECTED_GATE_COUNT is not None and EXPECTED_GATE_COUNT < 1:
    raise ValueError(
        f"PYFLI_SS2_GATE_COUNT must be positive, got {EXPECTED_GATE_COUNT}."
    )

if FOLD_REPETITIONS is not None and FOLD_REPETITIONS < 2:
    raise ValueError(
        "PYFLI_SS2_FOLD_REPETITIONS must be at least 2."
    )

if PERIOD_BINS is not None and PERIOD_BINS < 2:
    raise ValueError(
        "PYFLI_SS2_PERIOD_BINS must be at least 2."
    )

print(f"Data path: {DATA_PATH}")
print(f"Bit size: {BIT_SIZE}")
print(f"Pile-up correction: {PILE_UP}")
print(f"Temporal folding: {FOLD}")
print(f"Background subtraction: {SUB_BG}")
print(f"Hot-pixel correction: {HOT_PIXEL}")
print(f"Automatic hot-pixel map: {MAKE_HP_MAP}")

if IRF_PATH is not None:
    print(f"IRF path: {IRF_PATH}")

if BG_PATH is not None:
    print(f"Background path: {BG_PATH}")

if MASK_PATH is not None:
    print(f"Mask path: {MASK_PATH}")

if HP_PATH is not None:
    print(f"Hot-pixel mask path: {HP_PATH}")

2.2. Load an unprocessed reference#

The same acquisition is first loaded without detector corrections. The resulting cube is retained as a reference for direct comparison with the processed result.

raw_config = {
    "input_format": "auto",
    "bit_depth": BIT_SIZE,
    "pile_up": False,
    "fold": False,
    "hdf5_folder_mode": FOLDER_MODE,
}

if EXPECTED_GATE_COUNT is not None:
    raw_config[
        "ss2_expected_gate_count"
    ] = EXPECTED_GATE_COUNT

raw_loader = Detector(
    data_path=str(DATA_PATH),
    bit_size=BIT_SIZE,
)

raw_dataset = raw_loader.SS2(
    name="SwissSPAD2_raw",
    sub_bg=False,
    pile_up=False,
    hot_pixel=False,
    make_hp_map=False,
    config=raw_config,
)

raw_decay = np.asarray(
    raw_dataset["raw_data"]["decay"]
)

raw_processing = raw_dataset[
    "metadata"
][
    "processing"
]

input_format = raw_processing[
    "input_format"
]

if raw_decay.ndim != 3:
    raise RuntimeError(
        f"Raw SwissSPAD2 data must be three-dimensional, got {raw_decay.shape}."
    )

if not np.all(
    np.isfinite(raw_decay)
):
    raise RuntimeError(
        "Raw SwissSPAD2 data contains non-finite values."
    )

print(f"Input format: {input_format}")
print(f"Raw shape: {raw_decay.shape}")
print(f"Raw dtype: {raw_decay.dtype}")

folder_processing_requested = (
    SUB_BG
    or HOT_PIXEL
)

if folder_processing_requested:
    if input_format != "hdf5":
        raise RuntimeError(
            "Background subtraction and hot-pixel correction require "
            "SwissSPAD2 HDF5 folder input in the current detector pipeline. "
            f"The acquisition resolved to '{input_format}'."
        )

    if not DATA_PATH.is_dir():
        raise RuntimeError(
            "Background subtraction and hot-pixel correction require "
            "an HDF5 acquisition directory rather than a single HDF5 file."
        )

    if BG_PATH is not None and not BG_PATH.is_dir():
        raise RuntimeError(
            "PYFLI_SS2_BG_PATH must be an HDF5 directory when folder "
            "processing is enabled."
        )

if (
    input_format == "hdf5"
    and DATA_PATH.is_dir()
    and IRF_PATH is not None
    and not IRF_PATH.is_dir()
):
    raise RuntimeError(
        "PYFLI_SS2_IRF_PATH must be an HDF5 directory when the data "
        "acquisition is an HDF5 directory."
    )

2.3. Apply detector processing#

The processing configuration is assembled only from explicitly supplied acquisition parameters. PyFLI applies detector corrections before optional temporal folding and records the processing choices in the returned metadata.

processed_config = {
    "input_format": "auto",
    "bit_depth": BIT_SIZE,
    "pile_up": PILE_UP,
    "fold": FOLD,
    "hdf5_folder_mode": FOLDER_MODE,
}

if EXPECTED_GATE_COUNT is not None:
    processed_config[
        "ss2_expected_gate_count"
    ] = EXPECTED_GATE_COUNT

if FOLD_REPETITIONS is not None:
    processed_config[
        "fold_repetitions"
    ] = FOLD_REPETITIONS

if PERIOD_BINS is not None:
    processed_config[
        "period_bins"
    ] = PERIOD_BINS

if PHASE_SHIFT is not None:
    processed_config[
        "phase_shift"
    ] = PHASE_SHIFT

if DETECTOR_FREQUENCY_MHZ is not None:
    processed_config[
        "detector_frequency_mhz"
    ] = DETECTOR_FREQUENCY_MHZ

if LASER_FREQUENCY_MHZ is not None:
    processed_config[
        "laser_frequency_mhz"
    ] = LASER_FREQUENCY_MHZ

processed_loader = Detector(
    data_path=str(DATA_PATH),
    irf_path=(
        str(IRF_PATH)
        if IRF_PATH is not None
        else None
    ),
    bg_path=(
        str(BG_PATH)
        if BG_PATH is not None
        else None
    ),
    mask_path=(
        str(MASK_PATH)
        if MASK_PATH is not None
        else None
    ),
    hp_path=(
        str(HP_PATH)
        if HP_PATH is not None
        else None
    ),
    bit_size=BIT_SIZE,
)

processed_dataset = processed_loader.SS2(
    name="SwissSPAD2_processed",
    sub_bg=SUB_BG,
    pile_up=PILE_UP,
    hot_pixel=HOT_PIXEL,
    make_hp_map=MAKE_HP_MAP,
    threshold_sigma=THRESHOLD_SIGMA,
    config=processed_config,
)

processed_decay = np.asarray(
    processed_dataset["raw_data"]["decay"]
)

processed_irf = processed_dataset[
    "raw_data"
][
    "irf"
]

processed_background = processed_dataset[
    "raw_data"
][
    "background"
]

processed_mask = processed_dataset[
    "raw_data"
][
    "mask"
]

processing = processed_dataset[
    "metadata"
][
    "processing"
]

if processed_decay.ndim != 3:
    raise RuntimeError(
        "Processed SwissSPAD2 data must be three-dimensional, "
        f"got {processed_decay.shape}."
    )

if not np.all(
    np.isfinite(processed_decay)
):
    raise RuntimeError(
        "Processed SwissSPAD2 data contains non-finite values."
    )

print(f"Processed shape: {processed_decay.shape}")
print(f"Processed dtype: {processed_decay.dtype}")
print(f"Pile-up applied: {processing['pile_up']}")
print(f"Temporal folding applied: {processing['fold']}")
print(f"Background subtraction applied: {processing['sub_bg']}")
print(f"Hot-pixel correction applied: {processing['hot_pixel']}")

2.4. Compare integrated intensity#

The raw and processed acquisitions are displayed independently because detector corrections can change the count scale.

raw_intensity = np.sum(
    raw_decay,
    axis=-1,
    dtype=np.float64,
)

processed_intensity = np.sum(
    processed_decay,
    axis=-1,
    dtype=np.float64,
)

fig, axes = plt.subplots(
    1,
    2,
    figsize=(12, 5),
)

raw_image = axes[0].imshow(
    raw_intensity,
    cmap="turbo",
    origin="upper",
)

axes[0].set_title("Raw intensity")
axes[0].set_xlabel("X pixel")
axes[0].set_ylabel("Y pixel")

fig.colorbar(
    raw_image,
    ax=axes[0],
    fraction=0.046,
    pad=0.04,
    label="Integrated counts",
)

processed_image = axes[1].imshow(
    processed_intensity,
    cmap="turbo",
    origin="upper",
)

axes[1].set_title("Processed intensity")
axes[1].set_xlabel("X pixel")
axes[1].set_ylabel("Y pixel")

fig.colorbar(
    processed_image,
    ax=axes[1],
    fraction=0.046,
    pad=0.04,
    label="Integrated counts",
)

fig.tight_layout()
plt.show()

2.5. Compare temporal responses#

The spatially integrated raw and processed temporal responses are shown separately so changes in gate count or count scale introduced by folding and detector corrections remain explicit.

raw_trace = np.sum(
    raw_decay,
    axis=(0, 1),
    dtype=np.float64,
)

processed_trace = np.sum(
    processed_decay,
    axis=(0, 1),
    dtype=np.float64,
)

fig, axes = plt.subplots(
    2,
    1,
    figsize=(10, 7),
)

axes[0].plot(
    np.arange(raw_trace.size),
    raw_trace,
)

axes[0].set_title("Raw temporal response")
axes[0].set_xlabel("Gate index")
axes[0].set_ylabel("Integrated counts")
axes[0].grid(alpha=0.25)

axes[1].plot(
    np.arange(processed_trace.size),
    processed_trace,
)

axes[1].set_title("Processed temporal response")
axes[1].set_xlabel("Gate index")
axes[1].set_ylabel("Integrated counts")
axes[1].grid(alpha=0.25)

fig.tight_layout()
plt.show()

2.6. Apply the dataset mask#

When a mask is supplied to Detector, PyFLI packages it with the acquisition. DataPreprocessing can then apply the same spatial mask to the processed decay cube before analysis.

if processed_mask is not None:
    mask = np.asarray(
        processed_mask,
        dtype=bool,
    )

    if mask.shape != processed_decay.shape[:2]:
        raise RuntimeError(
            f"Mask shape {mask.shape} does not match processed spatial shape "
            f"{processed_decay.shape[:2]}."
        )

    (
        masked_decay,
    ) = DataPreprocessing(
        processed_decay
    ).apply_mask(
        mask=mask
    )

    masked_intensity = np.sum(
        masked_decay,
        axis=-1,
        dtype=np.float64,
    )

    fig, axes = plt.subplots(
        1,
        2,
        figsize=(11, 5),
    )

    axes[0].imshow(
        mask,
        cmap="gray",
        origin="upper",
    )

    axes[0].set_title("Analysis mask")
    axes[0].set_xlabel("X pixel")
    axes[0].set_ylabel("Y pixel")

    masked_image = axes[1].imshow(
        masked_intensity,
        cmap="turbo",
        origin="upper",
    )

    axes[1].set_title("Masked processed intensity")
    axes[1].set_xlabel("X pixel")
    axes[1].set_ylabel("Y pixel")

    fig.colorbar(
        masked_image,
        ax=axes[1],
        fraction=0.046,
        pad=0.04,
        label="Integrated counts",
    )

    fig.tight_layout()
    plt.show()

else:
    masked_decay = processed_decay

    print(
        "No analysis mask was supplied with this acquisition. "
        "The processed decay cube is retained without spatial masking."
    )

2.7. Inspect the IRF#

When an instrument response acquisition is supplied, it is loaded through the same detector-specific SwissSPAD2 path. The integrated IRF and decay traces are peak-normalized here only for visual comparison of their temporal positions.

if processed_irf is not None:
    processed_irf = np.asarray(
        processed_irf
    )

    if processed_irf.shape != processed_decay.shape:
        raise RuntimeError(
            f"Processed IRF shape {processed_irf.shape} does not match "
            f"processed decay shape {processed_decay.shape}."
        )

    irf_trace = np.sum(
        processed_irf,
        axis=(0, 1),
        dtype=np.float64,
    )

    decay_trace = np.sum(
        processed_decay,
        axis=(0, 1),
        dtype=np.float64,
    )

    irf_peak = np.max(
        irf_trace
    )

    decay_peak = np.max(
        decay_trace
    )

    if irf_peak <= 0:
        raise RuntimeError(
            "The integrated IRF has no positive signal."
        )

    if decay_peak <= 0:
        raise RuntimeError(
            "The integrated processed decay has no positive signal."
        )

    normalized_irf = irf_trace / irf_peak
    normalized_decay = decay_trace / decay_peak

    fig, ax = plt.subplots(figsize=(9, 4))

    ax.plot(
        np.arange(normalized_decay.size),
        normalized_decay,
        label="Decay",
    )

    ax.plot(
        np.arange(normalized_irf.size),
        normalized_irf,
        label="IRF",
    )

    ax.set_title("Normalized decay and IRF")
    ax.set_xlabel("Gate index")
    ax.set_ylabel("Normalized intensity")
    ax.legend()
    ax.grid(alpha=0.25)

    fig.tight_layout()
    plt.show()

else:
    print(
        "No IRF acquisition was supplied for this processing run."
    )

2.8. Inspect the processed background#

When background subtraction is enabled for an HDF5 folder acquisition, PyFLI packages the background used by the processing pipeline with the returned dataset.

if processed_background is not None:
    processed_background = np.asarray(
        processed_background
    )

    if processed_background.shape != processed_decay.shape:
        raise RuntimeError(
            f"Background shape {processed_background.shape} does not match "
            f"processed decay shape {processed_decay.shape}."
        )

    background_intensity = np.sum(
        processed_background,
        axis=-1,
        dtype=np.float64,
    )

    fig, ax = plt.subplots(figsize=(7, 6))

    image = ax.imshow(
        background_intensity,
        cmap="turbo",
        origin="upper",
    )

    fig.colorbar(
        image,
        ax=ax,
        label="Integrated background counts",
    )

    ax.set_title("SwissSPAD2 background")
    ax.set_xlabel("X pixel")
    ax.set_ylabel("Y pixel")

    fig.tight_layout()
    plt.show()

else:
    print(
        "No processed background cube is packaged with this acquisition."
    )

2.9. Review processing metadata#

The processing metadata records the resolved input format, detector settings, correction state, folding information, and underlying SPAD reader metadata. This provides a reproducible record of how the returned decay cube was produced.

pprint(
    processing,
    sort_dicts=False,
)

2.10. Summary#

The SwissSPAD2 acquisition has been processed through the detector-specific PyFLI pipeline and is now available in the standard PyFLI dataset structure.

The processed decay, optional IRF, optional background, and optional mask can be passed to downstream PyFLI lifetime-analysis workflows without detector-specific file handling.