Add soft per-spot quality weighting for adaptive spot detection

Add --soft-weight (implies --adaptive-spots): give every detected spot a
continuous quality weight in (0,1] and keep the highest-weight spots rather than
the brightest, so a deliberately loose detector self-cleans -- bright ice / salt
/ jet blobs and single-pixel noise no longer evict faint clean Bragg spots from
the max-spots cut.

The weight is a product of dimensionless gates (AdaptiveSpotFinderCPU::ApplyWeights,
computed against the per-ring background the adaptive finder already builds): a
logistic ramp in the spot's SNR and a soft size band (rises from one pixel,
plateaus, falls for oversized ice/salt/streak blobs). It carries on
DiffractionSpot -> SpotToSave and is consumed by FilterSpotsByCount, which ranks
by {non-ice, weight, intensity} when requested and by intensity otherwise, so the
classic and FPGA paths are unchanged.

Honest result: on the serial-stills battery this is index-rate-NEUTRAL. The
weighted ranking only changes the outcome when the spot count exceeds the
max-spots cap and the weight disagrees with intensity in a way that affects
indexing; the adaptive detectors already produce clean spot lists and the weak
sets sit under the cap, so re-ranking is a wash there (and a wash, not a
regression, on the one set that floods). Its intended benefit -- robustness to
ice/jet-contaminated frames and to a loosened detector -- is not exercised by
this battery; kept opt-in as the substrate for that.

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

common/DiffractionSpot.cpp

@@ -85,6 +85,7 @@ std::optional<SpotToSave> DiffractionSpot::Export(const DiffractionGeometry &geo
         .lattice = -1,
         .image = image_num,
         .d_A = d,
+        .weight = weight,
         .ice_ring = false,
         .indexed = false
     };

common/DiffractionSpot.h

@@ -15,7 +15,10 @@ class DiffractionSpot {
     int64_t pixel_count = 0;
     int64_t photons = 0;     // total photon count
     int64_t max_photons = INT64_MIN; // maximum number of counts per pixel in the spot
+    float weight = 1.0f;     // soft quality weight in (0,1]; 1 = "no opinion" (default)
 public:
+    void SetWeight(float w) { weight = w; }
+    float Weight() const { return weight; }
     DiffractionSpot() = default;
     DiffractionSpot(uint32_t col, uint32_t line, int64_t photons);
     DiffractionSpot(const SpotToSave &save);

common/SpotToSave.h

@@ -16,6 +16,7 @@ struct SpotToSave {
     int64_t h = 0, k = 0, l= 0;
     float d_A = 0.0;
     float dist_ewald_sphere = 0.0;
+    float weight = 1.0f;   // soft quality weight in (0,1]; used to rank the kept spots (see --soft-weight)
     bool ice_ring = false;
     bool indexed = false;
 

image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp

@@ -199,7 +199,34 @@ std::vector<DiffractionSpot> AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB
         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);
+    auto spots = ExtractSpots(image, settings, res_mask);
+    if (settings.soft_weight)
+        ApplyWeights(spots);
+    return spots;
+}
+
+void AdaptiveSpotFinderCPU::ApplyWeights(std::vector<DiffractionSpot> &spots) const {
+    const auto &pixel_to_bin = mapping.GetPixelToBin();
+    const size_t nbins = ring_mean.size();
+    const float READ = 1.0f;
+    for (auto &s : spots) {
+        const Coord c = s.RawCoord();   // flux-weighted centroid (col, row)
+        const int col = std::min(std::max(static_cast<int>(std::lround(c.x)), 0), width - 1);
+        const int row = std::min(std::max(static_cast<int>(std::lround(c.y)), 0), height - 1);
+        const uint16_t b = pixel_to_bin[static_cast<size_t>(row) * width + col];
+        const float mu = (b < nbins) ? ring_mean[b] : 0.0f;
+        const double N = std::max<int64_t>(s.PixelCount(), 1);
+        const double tot = std::max<int64_t>(s.Count(), 0);
+        const double signal = tot - N * mu;
+        const double noise = std::sqrt(std::max(1.0, tot + N * static_cast<double>(READ) * READ));
+        const double snr = signal / noise;
+        // Dimensionless gates (sigma, pixels): high SNR -> keep; a reasonable pixel count -> keep, while
+        // 1-pixel noise (rising edge) and oversized ice/salt/streak blobs (falling edge) -> ~0.
+        const float w_snr = 1.0f / (1.0f + std::exp(-static_cast<float>(snr - 4.0) / 1.5f));
+        const float w_size = (1.0f / (1.0f + std::exp(-(static_cast<float>(N) - 1.5f) / 0.7f)))
+                           * (1.0f / (1.0f + std::exp(-(40.0f - static_cast<float>(N)) / 8.0f)));
+        s.SetWeight(std::min(std::max(w_snr * w_size, 0.0f), 1.0f));
+    }
 }
 
 // Threshold-free variant. Build a noise-normalised image z = (I - ring_mean)/sqrt(ring_sigma^2+READ^2)
@@ -311,5 +338,7 @@ 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.soft_weight)
+        ApplyWeights(out);
     return out;
 }

image_analysis/spot_finding/AdaptiveSpotFinderCPU.h

@@ -41,6 +41,10 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder {
     std::vector<DiffractionSpot> RunPersistence(const ImagePreprocessorBuffer &image,
                                                 const SpotFindingSettings &settings,
                                                 const std::vector<bool> &res_mask);
+    // Assign each spot a soft quality weight (SNR against the per-ring background x a soft size band),
+    // 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;
 
 public:
     explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);

image_analysis/spot_finding/SpotFindingSettings.h

@@ -38,4 +38,11 @@ struct SpotFindingSettings {
     // is a graded per-spot significance and needs no min-pix (a lone noise spike has ~1 sigma
     // persistence; a real peak much more). See AdaptiveSpotFinderCPU::RunPersistence.
     bool spot_persistence = false;
+
+    // Soft per-spot weighting (implies adaptive detection): assign every detected spot a continuous
+    // quality weight from its SNR (against the per-ring background) and a soft size band (too few or
+    // too many pixels -> low), then keep the highest-WEIGHT spots rather than the brightest. Lets a
+    // 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;
 };

image_analysis/spot_finding/SpotUtils.cpp

@@ -37,15 +37,19 @@ void MarkIceRings(std::vector<SpotToSave> &spots, float tolerance_q_recipA) {
     }
 }
 
-void FilterSpotsByCount(std::vector<SpotToSave> &input, int64_t count) {
+void FilterSpotsByCount(std::vector<SpotToSave> &input, int64_t count, bool by_weight) {
     size_t output_size = std::min<size_t>(input.size(), count);
 
     std::ranges::partial_sort(input, input.begin() + output_size,
                               std::ranges::less{}, // comparator on the projected key
-                              [](const SpotToSave &s) {
-                                  // projection: key to compare by
-                                  return std::tuple{s.ice_ring, -s.intensity};
-                                  // false < true → non-ice first; negate intensity → higher first
+                              [by_weight](const SpotToSave &s) {
+                                  // projection: key to compare by. non-ice first (false < true), then
+                                  // by soft quality weight (higher first) when requested -- so a loose
+                                  // detector's bright junk cannot evict faint clean Bragg -- else by
+                                  // raw intensity. Intensity is the tie-breaker under the weight.
+                                  if (by_weight)
+                                      return std::tuple{s.ice_ring, -s.weight, -s.intensity};
+                                  return std::tuple{s.ice_ring, 0.0f, -s.intensity};
                               });
     input.resize(output_size);
 }
@@ -152,7 +156,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment,
 
     output.resolution_estimate = GetResolution(spots_out);
 
-    FilterSpotsByCount(spots_out, experiment.GetMaxSpotCount());
+    FilterSpotsByCount(spots_out, experiment.GetMaxSpotCount(), spot_finding_settings.soft_weight);
 
     output.spots = spots_out;
 }

image_analysis/spot_finding/SpotUtils.h

@@ -18,7 +18,7 @@ void CountSpots(DataMessage &msg,
 
 void MarkIceRings(std::vector<SpotToSave> &spots, float tolerance_q_recipA);
 
-void FilterSpotsByCount(std::vector<SpotToSave> &input, int64_t count);
+void FilterSpotsByCount(std::vector<SpotToSave> &input, int64_t count, bool by_weight = false);
 
 void FilterSpuriousHighResolutionSpots(std::vector<SpotToSave> &spots, float threshold);
 // Ignore high res. spots if there is a gap in (1/d) between two spots of dist_threshold (default: 0.25 A^-1)

rugnux/rugnux_cli.cpp

@@ -68,6 +68,7 @@ void print_usage() {
     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 << "   --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 << "   --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;
@@ -150,6 +151,7 @@ enum {
     OPT_ADAPTIVE_SPOTS,
     OPT_SPOT_FALSE_PIXELS,
     OPT_PERSISTENCE_SPOTS,
+    OPT_SOFT_WEIGHT,
     OPT_SPOT_RESOLUTION,
     OPT_SPOT_LOW_RESOLUTION,
     OPT_MAX_SPOTS,
@@ -263,6 +265,7 @@ static option long_options[] = {
     {"adaptive-spots", no_argument, nullptr, OPT_ADAPTIVE_SPOTS},
     {"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},
     {"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},
@@ -530,6 +533,7 @@ int main(int argc, char **argv) {
     bool adaptive_spots = false;
     float false_pixels_per_frame = 100.0f;
     bool persistence_spots = false;
+    bool soft_weight_flag = false;
     bool refine_bfactor = false;
     std::string ref_mtz;
     std::string ref_column;
@@ -774,6 +778,11 @@ int main(int argc, char **argv) {
             persistence_spots = true;
             logger.Info("Threshold-free (topological-persistence) spot detection enabled");
             break;
+        case OPT_SOFT_WEIGHT:
+            adaptive_spots = true;
+            soft_weight_flag = true;
+            logger.Info("Soft per-spot weighting enabled (keep highest-quality spots, not brightest)");
+            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);
@@ -1516,6 +1525,7 @@ int main(int argc, char **argv) {
     spot_settings.adaptive_threshold = adaptive_spots;
     spot_settings.false_pixels_per_frame = false_pixels_per_frame;
     spot_settings.spot_persistence = persistence_spots;
+    spot_settings.soft_weight = soft_weight_flag;
     if (d_min_spot_finding > 0.0f)
         spot_settings.high_resolution_limit = d_min_spot_finding;
     if (d_max_spot_finding > 0.0f)