diff --git a/NOTES.md b/NOTES.md index 3536a3a..ab3f5f6 100644 --- a/NOTES.md +++ b/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 diff --git a/csrc/cutoff.hpp b/csrc/cutoff.hpp new file mode 100644 index 0000000..5cddfcf --- /dev/null +++ b/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 +#include +#include +#include +#include + +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 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::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(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(q[1] + rise); + } +} + +} // namespace smappy diff --git a/csrc/fit.cpp b/csrc/fit.cpp index 64f083d..8104841 100644 --- a/csrc/fit.cpp +++ b/csrc/fit.cpp @@ -12,6 +12,7 @@ #include #include +#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 segment_cutoffs( + const py::array_t& values, + const py::array_t& 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 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::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, diff --git a/docs/examples/fit_and_measure.chain.yaml b/docs/examples/fit_and_measure.chain.yaml index 0379fe0..90980a0 100644 --- a/docs/examples/fit_and_measure.chain.yaml +++ b/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 diff --git a/src/smappy/detect.py b/src/smappy/detect.py index 9ffec13..7887ad8 100644 --- a/src/smappy/detect.py +++ b/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 diff --git a/src/smappy/plugins/docs/Localize/Gaussian 2D 2C.md b/src/smappy/plugins/docs/Localize/Gaussian 2D 2C.md index 1030962..f235dfb 100644 --- a/src/smappy/plugins/docs/Localize/Gaussian 2D 2C.md +++ b/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] --- diff --git a/src/smappy/plugins/docs/Localize/Gaussian 2D.md b/src/smappy/plugins/docs/Localize/Gaussian 2D.md index 9a123f7..209c3a9 100644 --- a/src/smappy/plugins/docs/Localize/Gaussian 2D.md +++ b/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 diff --git a/src/smappy/plugins/docs/Localize/Spline 3D 2C.md b/src/smappy/plugins/docs/Localize/Spline 3D 2C.md index c656e1f..362d0c1 100644 --- a/src/smappy/plugins/docs/Localize/Spline 3D 2C.md +++ b/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] --- diff --git a/src/smappy/plugins/docs/Localize/Spline 3D.md b/src/smappy/plugins/docs/Localize/Spline 3D.md index 5da7f1e..c072509 100644 --- a/src/smappy/plugins/docs/Localize/Spline 3D.md +++ b/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] --- diff --git a/src/smappy/plugins/fit.py b/src/smappy/plugins/fit.py index 25288d2..17212b7 100644 --- a/src/smappy/plugins/fit.py +++ b/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} diff --git a/tests/test_detect_roi.py b/tests/test_detect_roi.py index 87ff59c..e2d8862 100644 --- a/tests/test_detect_roi.py +++ b/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."""