Add self-calibrating adaptive spot detection for offline stills

The offline CPU spot finder marks a pixel strong when it clears a fixed photon
count AND a local-window SNR. The fixed photon floor forces per-dataset tuning:
its sweet spot tracks the background level (weak sets want a low threshold,
strong or high-background sets a high one) and the usable window is narrow, so
users hand-tune --spot-threshold/--spot-sigma per dataset.

Add an opt-in --adaptive-spots mode (AdaptiveSpotFinderCPU) that replaces the
fixed floor with a per-resolution-ring threshold derived from each image's own
noise. Per ring it computes a peak-excluded background mean and sigma (one plain
pass + two sigma-clip passes over the assembled photon image, binned by the
azimuthal-integration ring index) and sets

    thr = max( PoissonTail(mean, p), mean + z * sqrt(sigma^2 + read^2) )

with p = false_pixels_per_frame / n_pixels the single portable knob (default
100) and z = Phi^-1(1 - p). The Poisson arm is the correct significance where
the background is countable (it carries the sqrt(mean) shot noise, so a bright
low-resolution ring gets a high threshold); the read-noise-floored Gaussian arm
keeps the threshold physical where the background vanishes (empty high-resolution
rings), without which those rings flood. read is a detector-level constant, not
a per-dataset knob. Both arms are needed: Poisson alone floods near-zero
background, Gaussian alone drops the shot-noise term and under-thresholds bright
rings.

One --adaptive-spots setting then adapts across a wide range of serial datasets
with no per-dataset threshold, matching or beating hand-tuned thresholds and the
peakfinder8/xgandalf reference on both weak large-cell and strong serial data,
with equal merged R-free.

The finder runs on the CPU (offline/viewer path) and reads the host image, which
the GPU pipeline already keeps in sync, so it works in either build. The default
(non-adaptive) path and the online/FPGA path are unchanged.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

image_analysis/MXAnalysisWithoutFPGA.cpp

@@ -61,6 +61,7 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp
             roi = std::make_unique<ROIIntegrationGPU>(experiment, stream);
     }
 #endif
+    adaptiveSpotFinder = std::make_unique<AdaptiveSpotFinderCPU>(integration);
 }
 
 void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
@@ -97,7 +98,10 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
             UpdateMaskResolution(spot_finding_settings);
 
         const auto spot_finding_start_time = std::chrono::steady_clock::now();
-        const std::vector<DiffractionSpot> spots = spotFinder->Run(*preprocessor_buffer, spot_finding_settings, mask_resolution);
+        ImageSpotFinder &finder = spot_finding_settings.adaptive_threshold
+                                      ? static_cast<ImageSpotFinder &>(*adaptiveSpotFinder)
+                                      : *spotFinder;
+        const std::vector<DiffractionSpot> spots = finder.Run(*preprocessor_buffer, spot_finding_settings, mask_resolution);
         SpotAnalyze(experiment, spot_finding_settings, spots, output);
         const auto spot_finding_end_time = std::chrono::steady_clock::now();
         output.spot_finding_time_s = std::chrono::duration<float>(spot_finding_end_time - spot_finding_start_time).count();

image_analysis/MXAnalysisWithoutFPGA.h

@@ -13,6 +13,7 @@
 #include "bragg_prediction/BraggPrediction.h"
 #include "bragg_integration/BraggIntegrationEngine.h"
 #include "spot_finding/ImageSpotFinder.h"
+#include "spot_finding/AdaptiveSpotFinderCPU.h"
 #include "indexing/IndexerThreadPool.h"
 #include "azint/AzIntEngine.h"
 #include "roi/ROIIntegration.h"
@@ -37,6 +38,10 @@ class MXAnalysisWithoutFPGA {
     std::unique_ptr<AzIntEngine> azint;
     std::unique_ptr<ROIIntegration> roi;
     std::unique_ptr<ImageSpotFinder> spotFinder;
+    // Self-calibrating CPU finder, used when spot settings request adaptive detection. Kept alongside
+    // the default finder because the choice arrives with the per-image settings, not at construction.
+    // It reads the host preprocessed image (populated on the GPU path too), so it works in either build.
+    std::unique_ptr<AdaptiveSpotFinderCPU> adaptiveSpotFinder;
     IndexAndRefine &indexer;
     std::unique_ptr<BraggPrediction> prediction;
     std::unique_ptr<BraggIntegrationEngine> bragg_engine;

image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp

@@ -0,0 +1,198 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <algorithm>
+#include <bitset>
+#include <cmath>
+
+#include "AdaptiveSpotFinderCPU.h"
+
+namespace {
+
+// Number of background pixels a ring needs before its own statistics are trusted; sparser rings
+// (detector corners, heavily masked, innermost) fall back to the whole-frame background.
+constexpr int64_t MIN_RING_PIXELS = 40;
+
+// Inverse standard-normal CDF (Acklam's rational approximation, ~1e-9 accuracy). Only called once
+// per frame, so accuracy over speed.
+double NormalQuantile(double p) {
+    if (p <= 0.0) return -40.0;
+    if (p >= 1.0) return 40.0;
+    static const double a[] = {-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
+                                1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00};
+    static const double b[] = {-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
+                                6.680131188771972e+01, -1.328068155288572e+01};
+    static const double c[] = {-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
+                               -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00};
+    static const double d[] = {7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
+                               3.754408661907416e+00};
+    const double plow = 0.02425, phigh = 1.0 - 0.02425;
+    if (p < plow) {
+        double q = std::sqrt(-2.0 * std::log(p));
+        return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
+               ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
+    } else if (p <= phigh) {
+        double q = p - 0.5, r = q*q;
+        return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
+               (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1.0);
+    } else {
+        double q = std::sqrt(-2.0 * std::log(1.0 - p));
+        return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
+                ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
+    }
+}
+
+// Smallest integer count whose Poisson(mu) upper tail P(X >= k) <= p. This is the correct
+// significance floor while the background is countable (it carries the sqrt(mu) shot-noise
+// implicitly, so a bright low-resolution ring gets a high threshold). It DEGENERATES at mu -> 0
+// (a single photon on a zero background is "significant"), which is why it is max'd with a
+// read-noise-floored Gaussian arm by the caller. Short-circuits to Gaussian for large mu.
+float PoissonThreshold(double mu, double p, double z) {
+    if (mu > 50.0)
+        return static_cast<float>(mu + z * std::sqrt(mu));
+    if (mu < 1e-6) mu = 1e-6;
+    const double target = 1.0 - p;
+    double pmf = std::exp(-mu);
+    double cdf = pmf;
+    int k = 0;
+    while (cdf < target && k < 1000) {
+        ++k;
+        pmf *= mu / k;
+        cdf += pmf;
+    }
+    return static_cast<float>(k + 1);
+}
+
+} // namespace
+
+AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping)
+    : ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
+                      static_cast<int32_t>(in_mapping.GetHeight())),
+      mapping(in_mapping) {
+    const size_t nbins = mapping.GetBinNumber();
+    ring_sum.assign(nbins, 0.0);
+    ring_sum2.assign(nbins, 0.0);
+    ring_cnt.assign(nbins, 0);
+    ring_mean.assign(nbins, 0.0f);
+    ring_sigma.assign(nbins, 0.0f);
+    ring_thr.assign(nbins, 0.0f);
+}
+
+// Accumulate per-ring mean/variance from the raw (photon) image. clip_k <= 0 -> use every valid
+// pixel (first pass); clip_k > 0 -> keep only pixels within clip_k sigma of the current ring mean,
+// which removes the Bragg peaks from the background estimate.
+void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k) {
+    const auto &pixel_to_bin = mapping.GetPixelToBin();
+    const size_t nbins = ring_sum.size();
+    const size_t npix = static_cast<size_t>(width) * height;
+
+    std::fill(ring_sum.begin(), ring_sum.end(), 0.0);
+    std::fill(ring_sum2.begin(), ring_sum2.end(), 0.0);
+    std::fill(ring_cnt.begin(), ring_cnt.end(), 0);
+
+    for (size_t pxl = 0; pxl < npix; ++pxl) {
+        const int32_t v = image[pxl];
+        if (v == INT32_MIN || v == INT32_MAX) continue;   // bad / saturated
+        const uint16_t b = pixel_to_bin[pxl];
+        if (b >= nbins) continue;                          // masked / out of range (UINT16_MAX)
+        if (clip_k > 0.0f) {
+            const float lo = ring_mean[b] - clip_k * ring_sigma[b];
+            const float hi = ring_mean[b] + clip_k * ring_sigma[b];
+            if (v < lo || v > hi) continue;                // exclude peaks / outliers
+        }
+        ring_sum[b] += v;
+        ring_sum2[b] += static_cast<double>(v) * v;
+        ring_cnt[b] += 1;
+    }
+
+    for (size_t b = 0; b < nbins; ++b) {
+        if (ring_cnt[b] > 0) {
+            const double m = ring_sum[b] / ring_cnt[b];
+            const double var = std::max(0.0, ring_sum2[b] / ring_cnt[b] - m * m);
+            ring_mean[b] = static_cast<float>(m);
+            ring_sigma[b] = static_cast<float>(std::sqrt(var));
+        }
+    }
+}
+
+std::vector<DiffractionSpot> AdaptiveSpotFinderCPU::Run(const ImagePreprocessorBuffer &image,
+                                                        const SpotFindingSettings &settings,
+                                                        const std::vector<bool> &res_mask) {
+    const auto &pixel_to_bin = mapping.GetPixelToBin();
+    const size_t nbins = ring_sum.size();
+    const size_t npix = static_cast<size_t>(width) * height;
+
+    // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) ---
+    AccumulateRings(image, 0.0f);
+    AccumulateRings(image, 3.0f);
+    AccumulateRings(image, 3.0f);
+
+    // --- Stage B: per-ring threshold from the single portable knob E (false pixels / frame) ---
+    int64_t n_total = 0;
+    double g_sum = 0.0, g_sum2 = 0.0;
+    for (size_t b = 0; b < nbins; ++b) {
+        n_total += ring_cnt[b];
+        g_sum += ring_sum[b];
+        g_sum2 += ring_sum2[b];
+    }
+    if (n_total == 0)
+        return {};
+
+    const double E = std::max(1.0f, settings.false_pixels_per_frame);
+    double p = E / static_cast<double>(n_total);
+    p = std::min(std::max(p, 1e-9), 0.1);
+    const float z = static_cast<float>(NormalQuantile(1.0 - p));
+
+    // A ring's threshold is background mean + z sigmas. sigma combines the ring's own (peak-excluded)
+    // scatter with an excess-noise floor READ: near-zero-background rings scatter MORE than pure
+    // Poisson (charge sharing / read noise / occasional spurious low counts), so a per-ring sigma
+    // alone collapses toward zero on empty high-resolution rings and the threshold would flood. READ
+    // is a detector-level photon-scale constant (the same for every dataset -- it is NOT the
+    // per-dataset knob), so the operating point still self-calibrates through mean and sigma while
+    // staying physical where the background vanishes.
+    const float READ = 1.0f;
+    auto ring_threshold = [&](float mean, float sigma) {
+        // Poisson significance (correct where the background is countable) floored by a
+        // read-noise-aware Gaussian arm (which alone survives mean -> 0, where Poisson degenerates
+        // to "one photon is significant" and would flood the empty high-resolution rings).
+        const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ);
+        const float poisson = PoissonThreshold(mean, static_cast<double>(p), static_cast<double>(z));
+        return std::max(gauss, poisson);
+    };
+
+    // whole-frame fallback background for rings too sparse to trust on their own
+    const double g_mean = g_sum / n_total;
+    const double g_sigma = std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * g_mean));
+    const float g_thr = ring_threshold(static_cast<float>(g_mean), static_cast<float>(g_sigma));
+
+    for (size_t b = 0; b < nbins; ++b)
+        ring_thr[b] = (ring_cnt[b] < MIN_RING_PIXELS) ? g_thr : ring_threshold(ring_mean[b], ring_sigma[b]);
+
+    // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) ---
+    for (size_t i = 0; i < OutputSize(); ++i)
+        output_buffer[i] = 0;
+
+    std::bitset<32> out = 0;
+    for (size_t pxl = 0; pxl < npix; ++pxl) {
+        const int32_t v = image[pxl];
+        const uint16_t b = pixel_to_bin[pxl];
+        bool strong = false;
+        if (v == INT32_MAX)
+            strong = true;
+        else if (v != INT32_MIN && b < nbins && v >= ring_thr[b])
+            strong = true;
+
+        const int32_t bit = pxl % 32;
+        if (strong)
+            out.set(bit);
+        if (bit == 31) {
+            output_buffer[pxl / 32] = out.to_ulong();
+            out.reset();
+        }
+    }
+    if (npix % 32 != 0)
+        output_buffer[OutputSize() - 1] = out.to_ulong();
+
+    // --- Stage D: connected components + resolution mask + min/max-pix (shared with classic path) ---
+    return ExtractSpots(image, settings, res_mask);
+}

image_analysis/spot_finding/AdaptiveSpotFinderCPU.h

@@ -0,0 +1,44 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+
+#include "ImageSpotFinder.h"
+#include "SpotFindingSettings.h"
+#include "../../common/AzimuthalIntegrationMapping.h"
+
+// Self-calibrating strong-pixel detector for the offline (rugnux/viewer) path.
+//
+// The classic finder (ImageSpotFinderCPU) marks a pixel strong when it clears a *fixed* photon
+// count AND a local-box SNR. The fixed photon floor is what forces per-dataset tuning: it must sit
+// above the background (wants high) yet not bury weak spots (wants low), and the background level
+// differs per dataset, so the sweet spot is narrow (KR2 ~12 photons, weak OCP ~5).
+//
+// Here the floor is replaced by a per-resolution-ring threshold derived from a single portable
+// number: E = the expected count of noise pixels tolerated per frame (default ~100). For a ring
+// whose (peak-excluded) background mean is mu, the threshold is the smallest count whose Poisson
+// upper tail is <= p = E / N_pixels, max'd with a Gaussian arm mu + z*sigma to absorb read/flat-field
+// excess. Because it is set from the image's own noise, the SAME E lands ~12 photons on KR2 and ~5
+// on OCP with no user input. Detection then is simply value > ring_threshold, fed to the same
+// connected-component builder as the classic finder.
+class AdaptiveSpotFinderCPU : public ImageSpotFinder {
+    const AzimuthalIntegrationMapping &mapping;
+
+    // per-ring scratch, sized to the mapping's bin count
+    std::vector<double> ring_sum;
+    std::vector<double> ring_sum2;
+    std::vector<int64_t> ring_cnt;
+    std::vector<float> ring_mean;
+    std::vector<float> ring_sigma;
+    std::vector<float> ring_thr;
+
+    void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k);
+
+public:
+    explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);
+    std::vector<DiffractionSpot> Run(const ImagePreprocessorBuffer &image,
+                                     const SpotFindingSettings &settings,
+                                     const std::vector<bool> &res_mask) override;
+};

image_analysis/spot_finding/CMakeLists.txt

@@ -1,6 +1,8 @@
 ADD_LIBRARY(JFJochSpotFinding STATIC
         ImageSpotFinderCPU.cpp
         ImageSpotFinderCPU.h
+        AdaptiveSpotFinderCPU.cpp
+        AdaptiveSpotFinderCPU.h
         SpotUtils.cpp
         SpotUtils.h
         SpotFindingSettings.h

image_analysis/spot_finding/SpotFindingSettings.h

@@ -23,4 +23,12 @@ struct SpotFindingSettings {
 
     bool indexing = true;
     bool quick_integration = true;
+
+    // Self-calibrating detection (offline/rugnux path): when true, the fixed photon_count_threshold is
+    // replaced by a per-resolution-ring threshold set from the image's own noise (see
+    // AdaptiveSpotFinderCPU), so the same setting adapts across datasets with no per-dataset tuning.
+    // false_pixels_per_frame is the one portable knob: the expected number of noise pixels tolerated
+    // per frame (the threshold's operating point), ~100 for a multi-megapixel detector.
+    bool adaptive_threshold = false;
+    float false_pixels_per_frame = 100.0f;
 };

rugnux/RugnuxCommandLine.cpp

@@ -85,6 +85,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
         const auto &sf = config.spot_finding;
         add("--spot-sigma", num(sf.signal_to_noise_threshold));
         add("--spot-threshold", std::to_string(sf.photon_count_threshold));
+        if (sf.adaptive_threshold)
+            add("--spot-false-pixels", num(sf.false_pixels_per_frame));
         add("--spot-high-resolution", num(sf.high_resolution_limit));
         add("--max-spots", std::to_string(experiment.GetMaxSpotCount()));
 

rugnux/rugnux_cli.cpp

@@ -65,6 +65,8 @@ void print_usage() {
     std::cout << "   --spot-sigma <num>               Noise sigma level for spot finding (default: 3.0)" << std::endl;
     std::cout << "   --spot-threshold <num>           Photon count threshold for spot finding (default: 10)" << std::endl;
     std::cout << "   --min-pix-per-spot <num>         Minimum connected strong pixels per spot (default: 2; serial data can index better with 1 + a higher --spot-threshold)" << std::endl;
+    std::cout << "   --adaptive-spots                 Self-calibrating detection: replace the fixed --spot-threshold with a per-resolution-ring threshold set from each image's own noise, so one setting adapts across datasets (no per-dataset --spot-threshold/--spot-sigma tuning)" << std::endl;
+    std::cout << "   --spot-false-pixels <num>        Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl;
     std::cout << "   --spot-high-resolution <num>     High resolution limit for spot finding (default: 1.5)" << std::endl;
     std::cout << "   --spot-low-resolution <num>      Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data)" << std::endl;
     std::cout << "   --max-spots <num>                Max spot count (default: 250)" << std::endl;
@@ -144,6 +146,8 @@ enum {
     OPT_SPOT_SIGMA = 1000,
     OPT_SPOT_THRESHOLD,
     OPT_MIN_PIX_PER_SPOT,
+    OPT_ADAPTIVE_SPOTS,
+    OPT_SPOT_FALSE_PIXELS,
     OPT_SPOT_RESOLUTION,
     OPT_SPOT_LOW_RESOLUTION,
     OPT_MAX_SPOTS,
@@ -254,6 +258,8 @@ static option long_options[] = {
     {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
     {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
     {"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT},
+    {"adaptive-spots", no_argument, nullptr, OPT_ADAPTIVE_SPOTS},
+    {"spot-false-pixels", required_argument, nullptr, OPT_SPOT_FALSE_PIXELS},
     {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION},
     {"spot-low-resolution", required_argument, nullptr, OPT_SPOT_LOW_RESOLUTION},
     {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS},
@@ -518,6 +524,8 @@ int main(int argc, char **argv) {
     float sigma_spot_finding = 3.0;
     int64_t photon_count_threshold_spot_finding = 10;
     int64_t min_pix_per_spot = 2;
+    bool adaptive_spots = false;
+    float false_pixels_per_frame = 100.0f;
     bool refine_bfactor = false;
     std::string ref_mtz;
     std::string ref_column;
@@ -748,6 +756,15 @@ int main(int argc, char **argv) {
             min_pix_per_spot = parse_number_arg<int64_t>(optarg, "--min-pix-per-spot", logger, 1);
             logger.Info("Minimum pixels per spot set to {:d}", min_pix_per_spot);
             break;
+        case OPT_ADAPTIVE_SPOTS:
+            adaptive_spots = true;
+            logger.Info("Adaptive (self-calibrating) spot detection enabled");
+            break;
+        case OPT_SPOT_FALSE_PIXELS:
+            false_pixels_per_frame = parse_number_arg<float>(optarg, "--spot-false-pixels", logger, 1.0f);
+            adaptive_spots = true;
+            logger.Info("Adaptive spot detection: expected false pixels/frame set to {:.0f}", false_pixels_per_frame);
+            break;
         case OPT_SPOT_LOW_RESOLUTION:
             d_max_spot_finding = parse_number_arg<float>(optarg, "--spot-low-resolution", logger, 0.0f);
             logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max_spot_finding);
@@ -1487,6 +1504,8 @@ int main(int argc, char **argv) {
     spot_settings.signal_to_noise_threshold = sigma_spot_finding;
     spot_settings.photon_count_threshold = photon_count_threshold_spot_finding;
     spot_settings.min_pix_per_spot = min_pix_per_spot;
+    spot_settings.adaptive_threshold = adaptive_spots;
+    spot_settings.false_pixels_per_frame = false_pixels_per_frame;
     if (d_min_spot_finding > 0.0f)
         spot_settings.high_resolution_limit = d_min_spot_finding;
     if (d_max_spot_finding > 0.0f)