From ca7cbe206a9b64239e1b31fcbecfb81a1f041bd0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 24 Jul 2026 11:12:22 +0200 Subject: [PATCH] Add opt-in local-SNR spot gate and acceptance-fraction knob (serial stills) Two opt-in tools for weak serial-stills tuning; both default-off, so the default pipeline is bit-identical (verified: a serial-stills reference run reproduces HEAD's 7.85% indexing rate exactly). --local-snr (AdaptiveSpotFinderCPU::FilterByLocalSNR): after the loose per-ring adaptive threshold builds connected-component spots, drop any spot that does not stand this many sigmas above its OWN LOCAL background (robust median/MAD of a square annulus), not just the azimuthal ring mean. On structured-background (XFEL) frames the ring mean underestimates the local diffuse level in some sectors, so the ring threshold floods; a real Bragg peak still stands many local sigmas proud. Validated on XFEL stills to separate real peaks from flood at the pixel level (real median local-SNR ~70 vs flood ~2.6; SNR>=5 keeps ~99.8% of real peaks, ~14% of flood). GPU-portable (a per-spot local reduction). NOTE: on the current serial-stills battery it is index-rate/CC1/2 neutral -- the flood that survives as CC clusters overlaps weak-real spots, and only lattice-fit separates those -- but it is the correct tool for genuinely floody data (ice/jet/loosened detector) and the right substrate for the online FPGA path. --min-indexed-fraction : exposes the previously hardcoded 0.20 minimum indexed-spot fraction (AnalyzeIndexing) as a per-run setting. Lowering it admits weaker/sparser crystals; on flooded XFEL data the extra lattices are spurious (pair with --min-image-cc to gate them), on clean synchrotron data there are no marginal frames so it is a no-op -- useful as a gating-experiment primitive. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/IndexingSettings.cpp | 9 ++++ common/IndexingSettings.h | 6 +++ image_analysis/indexing/AnalyzeIndexing.cpp | 3 +- image_analysis/indexing/AnalyzeIndexing.h | 2 - .../spot_finding/AdaptiveSpotFinderCPU.cpp | 53 +++++++++++++++++++ .../spot_finding/AdaptiveSpotFinderCPU.h | 5 ++ .../spot_finding/SpotFindingSettings.h | 10 ++++ rugnux/RugnuxCommandLine.cpp | 2 + rugnux/rugnux_cli.cpp | 20 +++++++ 9 files changed, 107 insertions(+), 3 deletions(-) diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 2de46bcb..1400f85c 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -28,6 +28,15 @@ int64_t IndexingSettings::GetViableCellMinSpots() const { return viable_cell_min_spots; } +IndexingSettings &IndexingSettings::MinIndexedSpotFraction(float input) { + min_indexed_spot_fraction = input; + return *this; +} + +float IndexingSettings::GetMinIndexedSpotFraction() const { + return min_indexed_spot_fraction; +} + IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) { switch (input) { case IndexingAlgorithmEnum::Auto: diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 258e15ba..2872fc6d 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -24,6 +24,10 @@ class IndexingSettings { static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree int64_t indexing_threads = 4; int64_t viable_cell_min_spots = 9; + // Minimum fraction of the in-resolution spots a candidate lattice must index to be accepted. + // Lowering it admits weaker/sparser crystals (more real ones on flooded XFEL frames, but also + // more spurious lattices that a downstream merge-consistency gate must remove). + float min_indexed_spot_fraction = 0.20f; int64_t max_extra_lattices = 2; @@ -39,6 +43,7 @@ public: IndexingSettings(); IndexingSettings& ViableCellMinSpots(int64_t input); + IndexingSettings& MinIndexedSpotFraction(float input); IndexingSettings& Algorithm(IndexingAlgorithmEnum input); IndexingSettings& FFT_MaxUnitCell_A(float input); IndexingSettings& FFT_MinUnitCell_A(float input); @@ -58,6 +63,7 @@ public: IndexingSettings& MaxExtraLattices(int64_t input); [[nodiscard]] int64_t GetViableCellMinSpots() const; + [[nodiscard]] float GetMinIndexedSpotFraction() const; [[nodiscard]] IndexingAlgorithmEnum GetAlgorithm() const; [[nodiscard]] GeomRefinementAlgorithmEnum GetGeomRefinementAlgorithm() const; [[nodiscard]] float GetFFT_MaxUnitCell_A() const; diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 1a4be08d..bc1661b8 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -375,7 +375,8 @@ bool AnalyzeIndexing(DataMessage &message, int64_t indexing_lattice_count = 0; bool outcome = false; - if (nspots_indexed >= viable_cell_min_spots && nspots_indexed >= std::lround(min_percentage_spots * nspots_ref)) { + const float min_frac = experiment.GetIndexingSettings().GetMinIndexedSpotFraction(); + if (nspots_indexed >= viable_cell_min_spots && nspots_indexed >= std::lround(min_frac * nspots_ref)) { auto uc = latt.GetUnitCell(); if (ok(uc.a) && ok(uc.b) && ok(uc.c) && ok(uc.alpha) && ok(uc.beta) && ok(uc.gamma)) { message.indexing_result = true; diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index a2b347be..1e612343 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -7,8 +7,6 @@ #include "../../common/DiffractionExperiment.h" #include "../../common/JFJochMessages.h" -constexpr static float min_percentage_spots = 0.20f; - bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt, diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index 5ce1b687..e5430ad5 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -200,11 +200,62 @@ std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB // --- Stage D: connected components + resolution mask + min/max-pix (shared with classic path) --- auto spots = ExtractSpots(image, settings, res_mask); + if (settings.local_snr > 0.0f) + FilterByLocalSNR(image, spots, settings.local_snr); if (settings.soft_weight) ApplyWeights(spots); return spots; } +// Reject spots that do not stand out against their LOCAL background. The loose per-ring threshold +// keeps ~100% of real Bragg peaks but, on structured-background (XFEL) frames, also floods with +// spurious pixels: the azimuthal ring mean underestimates the diffuse level in some sectors, so a +// locally-high background pixel clears it. A real peak stands many sigmas above the background in +// its IMMEDIATE neighbourhood, a flood pixel does not. For each spot the background mean and scatter +// are measured from a square annulus around its centroid (robust median / MAD, so a neighbouring +// peak in the annulus cannot bias it), and the spot is kept only if its integrated signal exceeds +// k local sigmas. k is in sigma units -- self-calibrating, no photon threshold. +void AdaptiveSpotFinderCPU::FilterByLocalSNR(const ImagePreprocessorBuffer &image, + std::vector &spots, float k) const { + constexpr int RIN = 3; // half-width of the excluded core (7x7) + constexpr int ROUT = 6; // half-width of the background annulus (13x13) + std::vector bg; + bg.reserve((2 * ROUT + 1) * (2 * ROUT + 1)); + std::vector kept; + kept.reserve(spots.size()); + for (const auto &s : spots) { + const Coord c = s.RawCoord(); + const int col = static_cast(std::lround(c.x)); + const int row = static_cast(std::lround(c.y)); + bg.clear(); + for (int dr = -ROUT; dr <= ROUT; ++dr) { + const int rr = row + dr; + if (rr < 0 || rr >= height) continue; + for (int dc = -ROUT; dc <= ROUT; ++dc) { + if (std::max(std::abs(dr), std::abs(dc)) <= RIN) continue; // skip the peak core + const int cc = col + dc; + if (cc < 0 || cc >= width) continue; + const int32_t v = image[static_cast(rr) * width + cc]; + if (v == INT32_MIN || v == INT32_MAX) continue; // masked / saturated + bg.push_back(static_cast(v)); + } + } + if (bg.size() < 8) { kept.push_back(s); continue; } // too few bg pixels to judge + const size_t mid = bg.size() / 2; + std::nth_element(bg.begin(), bg.begin() + mid, bg.end()); + const float bg_med = bg[mid]; + for (auto &v : bg) v = std::fabs(v - bg_med); + std::nth_element(bg.begin(), bg.begin() + mid, bg.end()); + const float sigma = std::max(1.4826f * bg[mid], 1.0f); + const double npix = static_cast(std::max(s.PixelCount(), 1)); + const double signal = static_cast(s.Count()) - bg_med * npix; + const double snr = signal / (sigma * std::sqrt(npix)); + if (snr >= static_cast(k)) + kept.push_back(s); + } + spots.swap(kept); +} + void AdaptiveSpotFinderCPU::ApplyWeights(std::vector &spots) const { const auto &pixel_to_bin = mapping.GetPixelToBin(); const size_t nbins = ring_mean.size(); @@ -338,6 +389,8 @@ std::vector AdaptiveSpotFinderCPU::RunPersistence(const ImagePr std::vector out; out.reserve(spots.size()); for (auto &kv : spots) out.push_back(kv.second); + if (settings.local_snr > 0.0f) + FilterByLocalSNR(image, out, settings.local_snr); if (settings.soft_weight) ApplyWeights(out); return out; diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index 039653a4..a62dd792 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -45,6 +45,11 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { // consumed downstream to keep the best spots rather than the brightest. Needs the ring background, // so it must run after AccumulateRings. void ApplyWeights(std::vector &spots) const; + // Drop spots that do not stand k sigmas above their own LOCAL (annulus) background. The + // discriminator the per-ring threshold lacks on structured-background frames; see the header + // comment on SpotFindingSettings::local_snr. + void FilterByLocalSNR(const ImagePreprocessorBuffer &image, + std::vector &spots, float k) const; public: explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 970df0ae..f79c5354 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -45,4 +45,14 @@ struct SpotFindingSettings { // deliberately loose detector self-clean -- bright ice/salt/jet blobs and single-pixel noise no // longer evict faint clean Bragg spots. See AdaptiveSpotFinderCPU::ApplyWeights. bool soft_weight = false; + + // Local-background SNR gate (implies adaptive detection): after connected-component spots are + // built from the loose per-ring threshold, reject any spot whose integrated signal does not + // clear this many sigmas above its OWN LOCAL background (measured from a robust annulus around + // it), not just the azimuthal ring mean. On structured-background (XFEL) frames the ring mean + // underestimates the local diffuse level in some sectors, so the ring threshold floods with + // spurious pixels; a real Bragg peak still stands many local sigmas proud, so this recovers the + // clean spot list a global threshold cannot. In sigma units -> self-calibrating, no photon + // threshold. 0 disables. See AdaptiveSpotFinderCPU::FilterByLocalSNR. + float local_snr = 0.0f; }; diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 0b6b6bfa..7c258ddb 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -87,6 +87,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, add("--spot-threshold", std::to_string(sf.photon_count_threshold)); if (sf.adaptive_threshold) add("--spot-false-pixels", num(sf.false_pixels_per_frame)); + if (sf.local_snr > 0.0f) + add("--local-snr", num(sf.local_snr)); 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 73d3f2fd..624cd654 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -69,6 +69,8 @@ void print_usage() { std::cout << " --spot-false-pixels Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl; std::cout << " --persistence-spots Threshold-free variant of --adaptive-spots: score each intensity maximum by its topological persistence (no hard cut, no min-pix)" << std::endl; std::cout << " --soft-weight With --adaptive-spots: weight each spot by SNR + a soft size band and keep the highest-quality (not brightest) spots, so bright ice/salt/noise cannot crowd out faint Bragg" << std::endl; + std::cout << " --local-snr With --adaptive-spots: keep only spots standing this many sigma above their own LOCAL background (annulus), removing structured-background (XFEL) flood the ring threshold lets through (e.g. 5; implies --adaptive-spots)" << std::endl; + std::cout << " --min-indexed-fraction Minimum fraction of in-resolution spots a lattice must index to be accepted (default 0.20); lower to admit weaker/sparser crystals (pair with --min-image-cc to gate the extra spurious ones)" << 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; @@ -152,6 +154,8 @@ enum { OPT_SPOT_FALSE_PIXELS, OPT_PERSISTENCE_SPOTS, OPT_SOFT_WEIGHT, + OPT_LOCAL_SNR, + OPT_MIN_INDEXED_FRACTION, OPT_SPOT_RESOLUTION, OPT_SPOT_LOW_RESOLUTION, OPT_MAX_SPOTS, @@ -266,6 +270,8 @@ static option long_options[] = { {"spot-false-pixels", required_argument, nullptr, OPT_SPOT_FALSE_PIXELS}, {"persistence-spots", no_argument, nullptr, OPT_PERSISTENCE_SPOTS}, {"soft-weight", no_argument, nullptr, OPT_SOFT_WEIGHT}, + {"local-snr", required_argument, nullptr, OPT_LOCAL_SNR}, + {"min-indexed-fraction", required_argument, nullptr, OPT_MIN_INDEXED_FRACTION}, {"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}, @@ -534,6 +540,8 @@ int main(int argc, char **argv) { float false_pixels_per_frame = 100.0f; bool persistence_spots = false; bool soft_weight_flag = false; + float local_snr = 0.0f; + std::optional min_indexed_fraction; bool refine_bfactor = false; std::string ref_mtz; std::string ref_column; @@ -783,6 +791,15 @@ int main(int argc, char **argv) { soft_weight_flag = true; logger.Info("Soft per-spot weighting enabled (keep highest-quality spots, not brightest)"); break; + case OPT_LOCAL_SNR: + local_snr = parse_number_arg(optarg, "--local-snr", logger, 0.0f); + adaptive_spots = true; + logger.Info("Local-background SNR gate enabled: keep spots >= {:.1f} sigma above their local background", local_snr); + break; + case OPT_MIN_INDEXED_FRACTION: + min_indexed_fraction = parse_number_arg(optarg, "--min-indexed-fraction", logger, 0.0f); + logger.Info("Minimum indexed-spot fraction for acceptance set to {:.2f}", min_indexed_fraction.value()); + 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); @@ -1421,6 +1438,8 @@ int main(int argc, char **argv) { if (rotation_indexing_range.has_value()) indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value()); indexing_settings.GeomRefinementAlgorithm(refinement_algorithm); + if (min_indexed_fraction.has_value()) + indexing_settings.MinIndexedSpotFraction(min_indexed_fraction.value()); experiment.ImportIndexingSettings(indexing_settings); // --detect-ice-rings[=on|off] overrides the value carried in from the dataset (HDF5MetadataSource @@ -1526,6 +1545,7 @@ int main(int argc, char **argv) { spot_settings.false_pixels_per_frame = false_pixels_per_frame; spot_settings.spot_persistence = persistence_spots; spot_settings.soft_weight = soft_weight_flag; + spot_settings.local_snr = local_snr; if (d_min_spot_finding > 0.0f) spot_settings.high_resolution_limit = d_min_spot_finding; if (d_max_spot_finding > 0.0f)