From 0ed50d91c923a8ecbbbc8a1491e112abf77722c8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 23 Jul 2026 17:23:19 +0200 Subject: [PATCH] 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) --- image_analysis/MXAnalysisWithoutFPGA.cpp | 6 +- image_analysis/MXAnalysisWithoutFPGA.h | 5 + .../spot_finding/AdaptiveSpotFinderCPU.cpp | 198 ++++++++++++++++++ .../spot_finding/AdaptiveSpotFinderCPU.h | 44 ++++ image_analysis/spot_finding/CMakeLists.txt | 2 + .../spot_finding/SpotFindingSettings.h | 8 + rugnux/RugnuxCommandLine.cpp | 2 + rugnux/rugnux_cli.cpp | 19 ++ 8 files changed, 283 insertions(+), 1 deletion(-) create mode 100644 image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp create mode 100644 image_analysis/spot_finding/AdaptiveSpotFinderCPU.h diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 537c3493..bd256424 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -61,6 +61,7 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp roi = std::make_unique(experiment, stream); } #endif + adaptiveSpotFinder = std::make_unique(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 spots = spotFinder->Run(*preprocessor_buffer, spot_finding_settings, mask_resolution); + ImageSpotFinder &finder = spot_finding_settings.adaptive_threshold + ? static_cast(*adaptiveSpotFinder) + : *spotFinder; + const std::vector 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(spot_finding_end_time - spot_finding_start_time).count(); diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 16cd0560..1b40f2f4 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/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 azint; std::unique_ptr roi; std::unique_ptr 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 adaptiveSpotFinder; IndexAndRefine &indexer; std::unique_ptr prediction; std::unique_ptr bragg_engine; diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp new file mode 100644 index 00000000..ef4eece4 --- /dev/null +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -0,0 +1,198 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#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(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(k + 1); +} + +} // namespace + +AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping) + : ImageSpotFinder(static_cast(in_mapping.GetWidth()), + static_cast(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(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(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(m); + ring_sigma[b] = static_cast(std::sqrt(var)); + } + } +} + +std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask) { + const auto &pixel_to_bin = mapping.GetPixelToBin(); + const size_t nbins = ring_sum.size(); + const size_t npix = static_cast(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(n_total); + p = std::min(std::max(p, 1e-9), 0.1); + const float z = static_cast(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(p), static_cast(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(g_mean), static_cast(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); +} diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h new file mode 100644 index 00000000..c28f7290 --- /dev/null +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -0,0 +1,44 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#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 ring_sum; + std::vector ring_sum2; + std::vector ring_cnt; + std::vector ring_mean; + std::vector ring_sigma; + std::vector ring_thr; + + void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k); + +public: + explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); + std::vector Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask) override; +}; diff --git a/image_analysis/spot_finding/CMakeLists.txt b/image_analysis/spot_finding/CMakeLists.txt index 44bd6ca1..0a94ca7d 100644 --- a/image_analysis/spot_finding/CMakeLists.txt +++ b/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 diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 126d6442..b69b7072 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/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; }; diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 4d8450a6..0b6b6bfa 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/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())); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index eafc7967..6d526ce6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -65,6 +65,8 @@ void print_usage() { std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; std::cout << " --min-pix-per-spot 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 Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; std::cout << " --spot-low-resolution 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 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(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(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(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)