Bombcell integration and template metrics refactoring

[Julie-Fabre]

Update 26/2/26

After reviewing a lot of detected peaks, we removed smoothing and refined the detection as follows:

1. Use the specified prominence threshold to detect peaks. If multiple are found, the most prominent is selected as the main extremum.
2. If no peaks are found at the initial threshold, the threshold is halved and detection is attempted again. If multiple "peaks" are found at the half threshold, the most prominent is selected as the main extremum.
3. If still no peaks are found, a last resort method is used: use the global maximum in the search window.

Update 12/2/26

This PR ports bombcell-style unit classification to SpikeInterface and refactors/extends template metrics extensively.

Metrics

Template-metrics refactoring

Template metrics have been extended and refactored to align with the bombcell intefration. The new get_trough_and_peak_idx() function that uses scipy.signal.find_peaks()and finds all peaks (positive extrema) and troughs (negative extrema). Optionally, templates can be smoothed with a Savitzky-Golay filter.

Important

The smoothing is off by default

New metrics have been added and other renamed. See new docs here - https://spikeinterface--4306.org.readthedocs.build/en/4306/modules/metrics/template_metrics.html
Main additions are waveform_ratios, waveform_widths, and waveform_baseline_flatness

Quality metrics changes

Modified SNR implementation to match bombcell one as follows:

• use the median template (when available)
• use peak_sign="both" by default, which should handle also non-somatic spikes correctly

Bombcell labeling

Added an how-to to show the different labelign options available in SI: https://spikeinterface--4306.org.readthedocs.build/en/4306/how_to/auto_label_units.html

• New in spikeinterface.curation:

import spikeinterface.curation as sc

thresholds = sc.bombcell_get_default_thresholds()
unit_labels = sc.bombcell_classify_units(
    quality_metrics,
    thresholds=thresholds,
    classify_non_somatic=True,
)
summary = sc.get_classification_summary(unit_type)

Units get classified as noise → mua → good based on successive threshold checks. Optional non_soma category for non-somatic waveforms.

Important

The function has been refactored (both implementations are still there at the moment) to internally reuse the threshold_metrics_label_units, which makes the bombcell logic much clearer

Important

The bombcell dictionary has been refactored to contain different secions: "noise"/"mua"/"non-somatic". This makes it much easier to read and potentially to modify. We think this is better, so we don't have to rely on the hard-coded lists of metrics to distribute them in the function

Plots

• Added plots for classification summaries, metric histograms with threshold lines, waveform overlays by category, and UpSet plots.

from spikeinterface.widgets import (
    plot_unit_labels,  # former WaveformOverlay
    plot_metric_histograms,  # these are not per se connected to classification (so they are part of metrics.py
    plot_bombcell_labels_upset  # this is bombcell specific
)

plot_unit_labels(analyzer, unit_labels)
plot_metric_histograms(analyzer, thresholds=thresholds)
plot_bobcell_upset(analyzer, unit_labels, thresholds)

or a wrapper for all plots:

plots = plot_bombcell_unit_labeling_all(
    sorting_analyzer,
    unit_labels,
    thresholds=thresholds,            # optional, uses defaults
    include_upset=True,
)

Original PR body

This PR ports bombcell-style unit classification to SpikeInterface.

Template metrics

• Rewrote peak/trough detection with a new get_trough_and_peak_idx() function that uses scipy.signal.find_peaks(). Since SpikeInterface stores templates based on raw data rather than the heavily smoothed templates used in template matching, the waveforms can be noisy—so you can optionally apply Savitzky-Golay smoothing before detection. The function returns dicts for troughs, peaks before, and peaks after, each containing indices, values, prominences, and widths.

from spikeinterface.postprocessing import get_trough_and_peak_idx

troughs, peaks_before, peaks_after = get_trough_and_peak_idx(
    templates,
    sampling_frequency,
    smooth=True,
    min_thresh_detect_peaks_troughs=0.4,
)

• New metrics: peak_before_to_trough_ratio, peak_after_to_trough_ratio, waveform_baseline_flatness, peak_before_width, trough_width, main_peak_to_trough_ratio.

• Renamed peak_to_valley to peak_to_trough_duration.

analyzer.compute("template_metrics", metric_names=[
    "peak_before_to_trough_ratio",
    "waveform_baseline_flatness",
    "trough_width",
])

Quality metrics

• Added snr_bombcell—peak amplitude over baseline MAD.

analyzer.compute("quality_metrics", metric_names=["snr_bombcell"])

• amplitude_cutoff now has parameters for controlling the histogram fitting:

analyzer.compute("quality_metrics", metric_names=["amplitude_cutoff"], qm_params={
    "amplitude_cutoff": {
        "num_histogram_bins": 100,
        "histogram_smoothing_value": 3,
    }
})

Unit classification

• New in spikeinterface.curation:

import spikeinterface.comparison as sc

thresholds = sc.bombcell_get_default_thresholds()
unit_type, unit_type_string = sc.bombcell_classify_units(
    quality_metrics,
    thresholds=thresholds,
    classify_non_somatic=True,
)
summary = sc.get_classification_summary(unit_type, unit_type_string)

Units get classified as NOISE → MUA → GOOD based on successive threshold checks. Optional NON_SOMA category for non-somatic waveforms.

Plots

• Added plots for classification summaries, metric histograms with threshold lines, waveform overlays by category, and UpSet plots.

from spikeinterface.widgets import (
    plot_unit_classification,
    plot_classification_histograms,
    plot_waveform_overlay,
    plot_upset,
)

plot_unit_classification(analyzer, unit_type, unit_type_string)
plot_classification_histograms(quality_metrics, thresholds=thresholds)
plot_waveform_overlay(analyzer, unit_type, unit_type_string)
plot_upset(quality_metrics, unit_type, unit_type_string)

or a wrapper for all plots:

plots = plot_unit_classification_all(
    sorting_analyzer,
    unit_type,
    unit_type_string,
    quality_metrics=quality_metrics,  # optional, will try to get from analyzer
    thresholds=thresholds,            # optional, uses defaults
    split_non_somatic=False,
    include_upset=True,
)

[alejoe91]

@chrishalcrow @samuelgarcia @Julie-Fabre this is ready to merge!

I updated the docs, added a "new tool" entry in README and docs, and also added a reference also in the citation page.

I also changed the default column name returned by the bombcell_label_units function (bombcell_label) and also unitrefine_label_units (unitrefine_label). This way we can simply merge the returned dataframes to the metrics