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 <sigma> (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 <f>: 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) <noreply@anthropic.com>

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:

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;

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;

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,

image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp

@@ -200,11 +200,62 @@ std::vector<DiffractionSpot> 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<DiffractionSpot> &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<float> bg;
+    bg.reserve((2 * ROUT + 1) * (2 * ROUT + 1));
+    std::vector<DiffractionSpot> kept;
+    kept.reserve(spots.size());
+    for (const auto &s : spots) {
+        const Coord c = s.RawCoord();
+        const int col = static_cast<int>(std::lround(c.x));
+        const int row = static_cast<int>(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<size_t>(rr) * width + cc];
+                if (v == INT32_MIN || v == INT32_MAX) continue;             // masked / saturated
+                bg.push_back(static_cast<float>(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<double>(std::max<int64_t>(s.PixelCount(), 1));
+        const double signal = static_cast<double>(s.Count()) - bg_med * npix;
+        const double snr = signal / (sigma * std::sqrt(npix));
+        if (snr >= static_cast<double>(k))
+            kept.push_back(s);
+    }
+    spots.swap(kept);
+}
+
 void AdaptiveSpotFinderCPU::ApplyWeights(std::vector<DiffractionSpot> &spots) const {
     const auto &pixel_to_bin = mapping.GetPixelToBin();
     const size_t nbins = ring_mean.size();
@@ -338,6 +389,8 @@ std::vector<DiffractionSpot> AdaptiveSpotFinderCPU::RunPersistence(const ImagePr
     std::vector<DiffractionSpot> 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;

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<DiffractionSpot> &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<DiffractionSpot> &spots, float k) const;
 
 public:
     explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);

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;
 };

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()));
 

rugnux/rugnux_cli.cpp

@@ -69,6 +69,8 @@ void print_usage() {
     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 << "   --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 <sigma>              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 <f>       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 <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;
@@ -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<float> 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<float>(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<float>(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<float>(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)