The dynamic cutoff in C++, with SMAP's nearest-rank quantiles

NOTES.md

@@ -433,16 +433,12 @@ M1 Pro, 46,005 frames of 200x200, spline fit, ROI 13: **35.9 s** (1283 frames/s,
 
 Findings worth keeping:
 
-* the C++ single-pass maximum search beats `scipy.ndimage.maximum_filter` 5.2x
-  (short-circuit evaluation rejects most pixels after one comparison)
-* the 2-D DoG kernel has **rank 2**, so it is not separable -- subtracting the
-  1-D kernels is wrong.  Subtracting the 2-D kernels and doing one convolution
-  is exact and, below radius 4, faster than the fused separable passes
-  (`_SEPARABLE_FROM_RADIUS`)
+* the C++ single-pass maximum search beats `scipy.ndimage.maximum_filter` 5.2x (short-circuit evaluation rejects most pixels after one comparison)
+* the 2-D DoG kernel has **rank 2**, so it is not separable -- subtracting the 1-D kernels is wrong.  Subtracting the 2-D kernels and doing one convolution is exact and, below radius 4, faster than the fused separable passes (`_SEPARABLE_FROM_RADIUS`)
 * SMAP's DoG window is `max(ceil(6*sigma-1), 3)`: radius 3 for sigma 1.2
 * threading gives ~4x on filtering and maxima, ~6.5x on fitting
-* the per-frame dynamic-cutoff loop is the last serial part (~2.7 s, 80% of it
-  numpy dispatch overhead).  It could be vectorised with one lexsort per chunk.
+* the dynamic cutoff ran as a Python loop over frames (~2.7 s of the detection's 12.9 s, 80% of it numpy dispatch overhead).  It needs three order statistics, not a sort: `csrc/cutoff.hpp` finds them with nth_element, frames over threads, **9.3 -> 1.1 ms** per 200 frames.  A lexsort over the block, the obvious vectorisation, is slower than the loop it replaced.
+* its quantiles are SMAP's nearest rank (`myquantilefast`); smappy had `np.quantile`'s interpolation only because it is numpy's default.  The cutoff rises slightly, dropping 0.01-0.14 % of the candidates, those at the threshold.  Note that clang fuses `a * b - c` into an FMA by default, so a kernel that has to match a reference bit for bit must prevent it (`cutoff.hpp`).
 
 ### Rendering
 

csrc/cutoff.hpp

@@ -0,0 +1,50 @@
+// The dynamic cutoff of every frame of a block: median + factor * the slope of
+// the 20-80 % quantile range of that frame's local maxima -- SMAP's
+// getdynamiccutoff, and `smappy.detect.DynamicCutoff.__call__`, bit for bit.
+//
+// A quantile is the ceil(n p)-th smallest value, as SMAP's myquantilefast
+// takes it: no interpolation, and three order statistics are not a sort, so
+// nth_element finds them in linear time, each search starting where the last
+// one left off.
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <vector>
+
+namespace smappy {
+
+// Cutoffs for the segments values[starts[i], starts[i+1]), i in [first, last).
+// A segment with fewer than `min_count` maxima gets NaN: below ten the rule is
+// a mean (`DynamicCutoff.__call__`), which the caller keeps.
+inline void segment_cutoffs(const float* values, const int64_t* starts,
+                            long long first, long long last, double factor,
+                            long long min_count, float* out) {
+    const double levels[3] = {0.2, 0.5, 0.8};
+    std::vector<float> v;
+    for (long long s = first; s < last; ++s) {
+        const long long n = starts[s + 1] - starts[s];
+        if (n < min_count || n <= 0) {
+            out[s] = std::numeric_limits<float>::quiet_NaN();
+            continue;
+        }
+        v.assign(values + starts[s], values + starts[s + 1]);
+        double q[3];
+        long long from = 0;
+        for (int k = 0; k < 3; ++k) {
+            const long long rank = static_cast<long long>(std::ceil(n * levels[k])) - 1;
+            std::nth_element(v.begin() + from, v.begin() + rank, v.end());
+            q[k] = v[rank];
+            from = rank;
+        }
+        // 0.8 - 0.2 as MATLAB divides by it, 0.6000000000000001 and not 0.6
+        const double slope = (q[2] - q[0]) / (0.8 - 0.2);
+        // apart from the sum: fused into an FMA it would round once, not twice
+        volatile double rise = slope * factor;
+        out[s] = static_cast<float>(q[1] + rise);
+    }
+}
+
+}  // namespace smappy

csrc/fit.cpp

@@ -12,6 +12,7 @@
 #include <thread>
 #include <vector>
 
+#include "cutoff.hpp"
 #include "filters.hpp"
 #include "lm.hpp"
 #include "maxima.hpp"
@@ -339,6 +340,30 @@ py::tuple find_maxima(const Array& images, float threshold, int n_threads) {
                           std::move(value));
 }
 
+py::array_t<float> segment_cutoffs(
+        const py::array_t<float, py::array::c_style | py::array::forcecast>& values,
+        const py::array_t<int64_t, py::array::c_style | py::array::forcecast>& starts,
+        double factor, long long min_count, int n_threads) {
+    if (values.ndim() != 1 || starts.ndim() != 1 || starts.shape(0) < 1)
+        throw std::invalid_argument("values and starts must be 1-D, starts non-empty");
+    const long long segments = starts.shape(0) - 1;
+    const int64_t* s = starts.data();
+    for (long long i = 0; i < segments; ++i)
+        if (s[i] < 0 || s[i] > s[i + 1] || s[i + 1] > values.shape(0))
+            throw std::invalid_argument("starts must ascend within values");
+    py::array_t<float> out(segments);
+    float* o = out.mutable_data();
+    const float* v = values.data();
+    {
+        py::gil_scoped_release release;
+        smappy::parallel_ranges(segments, smappy::resolve_threads(n_threads, segments),
+                                [&](long long begin, long long end, int) {
+            smappy::segment_cutoffs(v, s, begin, end, factor, min_count, o);
+        });
+    }
+    return out;
+}
+
 }  // namespace
 
 PYBIND11_MODULE(_fit3d, m) {
@@ -362,6 +387,10 @@ PYBIND11_MODULE(_fit3d, m) {
           py::arg("threshold") = -std::numeric_limits<float>::infinity(),
           py::arg("n_threads") = 0,
           "Strict 3x3 local maxima of an image block -> (frame, y, x, value).");
+    m.def("segment_cutoffs", &segment_cutoffs, py::arg("values"), py::arg("starts"),
+          py::arg("factor"), py::arg("min_count") = 10, py::arg("n_threads") = 0,
+          "Dynamic cutoff of each segment values[starts[i]:starts[i+1]], as float32; "
+          "NaN below min_count.");
 
     m.def("fit_cspline", &fit_cspline, py::arg("rois"), py::arg("coeff"),
           py::arg("z_start"), py::arg("iterations") = 50, py::arg("n_threads") = 0,

docs/examples/fit_and_measure.chain.yaml

@@ -6,7 +6,7 @@ description: Gaussian fit, localization statistics
 steps:
   - plugin: Localize/Gaussian 2D
     label: fit
-    version: "2"
+    version: "3"
     values:
       camera.conversion: 0.5           # null (or leaving it out) reads it from the file
       camera.offset: 100.0

src/smappy/detect.py

@@ -17,6 +17,7 @@
 
 from __future__ import annotations
 
+import math
 from dataclasses import dataclass
 from typing import Optional, Protocol
 
@@ -143,9 +144,19 @@ def __str__(self) -> str:
 class DynamicCutoff:
     """Threshold derived from the distribution of local-maximum intensities.
 
-    Follows SMAP: the cutoff is the median of the maxima plus ``factor`` times
-    the slope of their 20-80 % quantile range.  Computed per frame, so it
-    adapts to changing background without needing the whole movie.
+    SMAP's ``getdynamiccutoff``: the cutoff is the median of the maxima plus
+    ``factor`` times the slope of their 20-80 % quantile range.  Computed per
+    frame, so it adapts to changing background without needing the whole
+    movie.
+
+    A quantile is the ceil(n p)-th smallest maximum, as SMAP's
+    ``myquantilefast`` takes it.  Up to 0.3.1 it was ``np.quantile``'s linear
+    interpolation -- numpy's default, and nothing records choosing it.  The
+    two differ in most frames, by 0.04 % of the cutoff at 1500 maxima a frame
+    and about 1 % at 140, and nearest rank is the higher -- it drops 0.01-0.14 %
+    of the candidates, the ones at the threshold, and adds none.  Matching SMAP
+    is the point of the port, and the kernel that computes it in C++ is the
+    simpler for it (no interpolation to reproduce to the last bit).
     """
 
     factor: float = 1.7
@@ -156,9 +167,44 @@ def __call__(self, maxima_values: np.ndarray) -> float:
             return 0.0
         if n < 10:
             return float(np.mean(maxima_values) * self.factor)
-        q20, q50, q80 = np.quantile(maxima_values, (0.2, 0.5, 0.8))
-        slope = (q80 - q20) / 0.6
-        return float(q50 + slope * self.factor)
+        ranks = [math.ceil(n * p) - 1 for p in (0.2, 0.5, 0.8)]
+        q20, q50, q80 = (float(q) for q in np.partition(maxima_values, ranks)[ranks])
+        # 0.8 - 0.2 is what SMAP divides by, and 0.6000000000000001 in double
+        slope = (q80 - q20) / (0.8 - 0.2)
+        return q50 + slope * self.factor
+
+    def thresholds(self, values: np.ndarray, starts: np.ndarray,
+                   n_threads: int = 0) -> np.ndarray:
+        """The cutoff of every segment ``values[starts[i]:starts[i+1]]`` at once.
+
+        The same numbers as calling the cutoff once per segment, bit for bit,
+        returned as float32 because that is what ``values > cutoff`` compared
+        in: a Python float against a float32 array is rounded to float32 first.
+
+        Called per frame, the cutoff was a Python loop in detection -- three
+        quantiles a frame, 80 % of it numpy's dispatch: 9.3 ms for 200
+        frames of 1500 maxima.  One lexsort over the block, the obvious
+        vectorisation, is slower than that (52 ms, a comparison sort over
+        300 k values), and sorting a padded frames x maxima array along its
+        rows is 5.2 ms, most of it the sort.  A cutoff needs three order
+        statistics, not a sort: `csrc/cutoff.hpp` finds them with nth_element,
+        frames split over threads, in 1.1 ms.  Without the extension each
+        segment goes through `__call__`, which is the definition anyway.
+        """
+        values = np.ascontiguousarray(values, dtype=np.float32)
+        starts = np.ascontiguousarray(starts, dtype=np.int64)
+        counts = np.diff(starts)
+        kernel = getattr(_fit3d, "segment_cutoffs", None)
+        if kernel is not None:
+            out = kernel(values, starts, float(self.factor), 10, n_threads)
+            # below ten maxima the rule is a mean, which the kernel leaves NaN
+            rest = np.flatnonzero(np.isnan(out) & (counts > 0))
+        else:
+            out = np.full(counts.size, np.nan, np.float32)
+            rest = np.flatnonzero(counts > 0)
+        for i in rest:
+            out[i] = self(values[starts[i]:starts[i + 1]])
+        return out
 
     def __str__(self) -> str:
         return f"dynamic(factor={self.factor:g})"
@@ -265,6 +311,18 @@ def find_candidates(filtered: np.ndarray, cutoff,
     if split is not None:
         axis, position = split
         side = (rows if axis == 0 else cols) >= position
+    if isinstance(cutoff, DynamicCutoff):
+        # each half is still sorted by frame, so its frames are segments too
+        for part in ((slice(None),) if side is None else (side, ~side)):
+            index = np.arange(frames.size)[part]
+            starts = np.searchsorted(frames[index], np.arange(n_frames + 1))
+            limit = np.repeat(cutoff.thresholds(values[index], starts, n_threads),
+                              np.diff(starts))
+            keep[index] = values[index] > limit
+        return Candidates(frame=frames[keep].astype(np.int64),
+                          x=cols[keep].astype(np.int32),
+                          y=rows[keep].astype(np.int32),
+                          value=values[keep])
     for start, stop in zip(bounds[:-1], bounds[1:]):
         if stop <= start:
             continue

src/smappy/plugins/docs/Localize/Gaussian 2D 2C.md

@@ -1,5 +1,5 @@
 ---
-version: "3"
+version: "4"
 covers: [smappy.dualfit.combine_peaks, smappy.dualfit.build_link, smappy.dualfit.cut_paired_rois, smappy.dualfit.paired_to_localizations, smappy.dualfit.DualChannelEngine, smappy.psf.GlobalGaussianPSF, smappy.detect.find_candidates, smappy.plugins.fit.DualGaussianFit, smappy.plugins.fit.calibration_blocks, smappy.plugins.fit.finish_localizations, smappy.calibrate.transform.register_channels]
 ---
 

src/smappy/plugins/docs/Localize/Gaussian 2D.md

@@ -1,5 +1,5 @@
 ---
-version: "2"
+version: "3"
 covers: [smappy.pipeline.LocalizationEngine, smappy.detect.PeakFinder, smappy.detect.DynamicCutoff, smappy.psf.GaussianPSF, smappy.camera.to_photons, smappy.rawframes.RawFrameKeeper]
 ---
 
@@ -82,7 +82,9 @@ cutoff is set above it,
 $$\mathrm{cutoff} = \mathrm{median} + f \cdot \frac{q_{80} - q_{20}}{0.6} ,$$
 
 with $q_{20}$ and $q_{80}$ the 20th and 80th percentiles of the maxima of that
-frame, and $f$ the *cutoff value* (1.7).  The fraction is a robust measure of
+frame, and $f$ the *cutoff value* (1.7).  A percentile is taken as SMAP takes
+it: the $p$-th of $n$ maxima is the $\lceil np \rceil$-th smallest, with no
+interpolation between neighbours.  The fraction is a robust measure of
 the spread of the noise, so the cutoff is "$f$ noise widths above the typical
 maximum" whatever the background.  The **absolute** cutoff is a fixed
 number, in photons of the filtered image.  Lower values find dimmer

src/smappy/plugins/docs/Localize/Spline 3D 2C.md

@@ -1,5 +1,5 @@
 ---
-version: "5"
+version: "6"
 covers: [smappy.psf.GlobalSplinePSF, smappy.dualfit.combine_peaks, smappy.dualfit.build_link, smappy.dualfit.cut_paired_rois, smappy.dualfit.paired_to_localizations, smappy.dualfit.DualChannelEngine, smappy.calibrate.dual.build_dual_calibration, smappy.calibrate.dual.fit_dual_transform, smappy.calibrate.core.positive_pair_models, smappy.calibrate.dual.load_dual_color_calibration, smappy.plugins.fit.DualSplineFit, smappy.plugins.fit.finish_localizations]
 ---
 

src/smappy/plugins/docs/Localize/Spline 3D.md

@@ -1,5 +1,5 @@
 ---
-version: "2"
+version: "3"
 covers: [smappy.psf.SplinePSF, smappy.io.calibration.load_spline_calibration, smappy.io.calibration.SplineCalibration, smappy.io.calibration.evaluate_spline, smappy.io.calibration.warn_on_em_mismatch, smappy.roi.cut_rois, smappy.calibrate.core.build_calibration]
 ---
 

src/smappy/plugins/fit.py

@@ -934,7 +934,8 @@ class GaussianFit(_FitPlugin):
     description = "Detect and fit with a free-width Gaussian PSF: x, y, photons, sigma."
     Settings = GaussianFitSettings
     # 2: the camera's read noise in the likelihood (1 e- without EM gain)
-    version = "2"
+    # 3: the dynamic cutoff's quantiles nearest-rank, as SMAP takes them
+    version = "3"
     params = {f"fit.{k}": v for k, v in FIT_PARAMS.items()}
 
     def model(self, settings, camera):
@@ -947,7 +948,8 @@ class SplineFit(_FitPlugin):
                    "calibration: adds z.")
     Settings = SplineFitSettings
     # 2: the camera's read noise in the likelihood (1 e- without EM gain)
-    version = "2"
+    # 3: the dynamic cutoff's quantiles nearest-rank, as SMAP takes them
+    version = "3"
     params = GaussianFit.params
 
     def model(self, settings, camera):
@@ -975,7 +977,8 @@ class DualSplineFit(_FitPlugin):
     #    the background (the version is the calibration's as much as the fit's)
     # 5: the camera's read noise in the likelihood; the calibration's photon
     #    shares from the light before clipping
-    version = "5"
+    # 6: the dynamic cutoff's quantiles nearest-rank, as SMAP takes them
+    version = "6"
     Settings = DualSplineFitSettings
     params = {**GaussianFit.params, **finish_params()}
 
@@ -1285,7 +1288,8 @@ class DualGaussianFit(_FitPlugin):
                    "registration can be measured from the movie itself.")
     # 2: linked photons are the total over both halves, not the main half's
     # 3: the camera's read noise in the likelihood (1 e- without EM gain)
-    version = "3"
+    # 4: the dynamic cutoff's quantiles nearest-rank, as SMAP takes them
+    version = "4"
     Settings = DualGaussianFitSettings
     params = {**GaussianFit.params, **finish_params(), **REGISTRATION_PARAMS}
 

tests/test_detect_roi.py

@@ -112,6 +112,58 @@ def test_dynamic_cutoff_scales_with_the_distribution():
     assert DynamicCutoff(3.0)(rng.normal(10, 1, 1000)) > weak
 
 
+def test_the_dynamic_cutoff_takes_smaps_nearest_rank_quantiles():
+    """20 maxima: the 4th, 10th and 16th smallest, as SMAP's myquantilefast
+    picks them.  np.quantile would interpolate to 0.88, 1.35 and 2.42."""
+    values = np.array([0.4, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5,
+                       1.6, 1.8, 2.0, 2.3, 2.9, 3.5, 41.0, 95.0], np.float32)
+    q20, q50, q80 = 0.8, 1.3, 2.3
+    expected = q50 + 1.7 * (q80 - q20) / 0.6
+    assert DynamicCutoff(1.7)(values) == pytest.approx(expected, abs=1e-5)
+    starts = np.array([0, values.size])
+    assert DynamicCutoff(1.7).thresholds(values, starts)[0] == pytest.approx(expected,
+                                                                              abs=1e-5)
+
+
+@pytest.mark.parametrize("extension", [True, False])
+def test_the_cutoff_of_a_block_is_the_cutoff_of_each_frame_bit_for_bit(
+        extension, monkeypatch):
+    """Empty frames, frames under ten maxima, ties, and values far from 1."""
+    if not extension:
+        monkeypatch.setattr("smappy.detect._fit3d", None)
+    rng = np.random.default_rng(3)
+    cutoff = DynamicCutoff(1.7)
+    for trial in range(60):
+        counts = rng.integers(0, 2000 if trial % 2 else 25, rng.integers(1, 40))
+        starts = np.concatenate([[0], np.cumsum(counts)])
+        values = (rng.gamma(2, 3, starts[-1]) * 10.0 ** rng.integers(-3, 4))
+        if trial % 3 == 0:
+            values = np.round(values, 1)
+        values = values.astype(np.float32)
+        block = cutoff.thresholds(values, starts)
+        for i, (a, b) in enumerate(zip(starts[:-1], starts[1:])):
+            if b > a:
+                assert block[i] == np.float32(cutoff(values[a:b]))
+
+
+@pytest.mark.parametrize("split", [None, (1, 30), (0, 17)])
+def test_detection_keeps_what_the_per_frame_loop_keeps(split):
+    """A cutoff that is not a DynamicCutoff still takes the loop: the reference."""
+    rng = np.random.default_rng(4)
+    images = rng.normal(20, 3, (20, 48, 60)).astype(np.float32)
+    for frame in images:
+        frame += _spots(rng.uniform(5, 55, (4, 2)), shape=(48, 60), background=0.0)
+    images[3] = 0.0                                  # a frame with no maxima
+    images[5] = _spots([(20, 20), (40, 30)], shape=(48, 60))  # one with two
+    filtered = DoGFilter(1.2)(images)
+    cutoff = DynamicCutoff(1.7)
+    fast = find_candidates(filtered, cutoff, split=split)
+    loop = find_candidates(filtered, lambda v: cutoff(v), split=split)
+    assert len(fast) > 0
+    for name in ("frame", "x", "y", "value"):
+        assert np.array_equal(getattr(fast, name), getattr(loop, name))
+
+
 @pytest.mark.parametrize("n_threads", [1, 2, 8])
 def test_threading_does_not_change_detection(n_threads):
     """Threaded filtering and peak finding must be bit-identical to serial."""