From 7c5bedfd74e8ae20dd3e3144d70d50ae0dfcb029 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 23 Jul 2026 19:53:55 +0200 Subject: [PATCH] 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) --- common/DiffractionSpot.cpp | 1 + common/DiffractionSpot.h | 3 ++ common/SpotToSave.h | 1 + .../spot_finding/AdaptiveSpotFinderCPU.cpp | 31 ++++++++++++++++++- .../spot_finding/AdaptiveSpotFinderCPU.h | 4 +++ .../spot_finding/SpotFindingSettings.h | 7 +++++ image_analysis/spot_finding/SpotUtils.cpp | 16 ++++++---- image_analysis/spot_finding/SpotUtils.h | 2 +- rugnux/rugnux_cli.cpp | 10 ++++++ 9 files changed, 67 insertions(+), 8 deletions(-) diff --git a/common/DiffractionSpot.cpp b/common/DiffractionSpot.cpp index 5e6a8857..9d65d668 100644 --- a/common/DiffractionSpot.cpp +++ b/common/DiffractionSpot.cpp @@ -85,6 +85,7 @@ std::optional DiffractionSpot::Export(const DiffractionGeometry &geo .lattice = -1, .image = image_num, .d_A = d, + .weight = weight, .ice_ring = false, .indexed = false }; diff --git a/common/DiffractionSpot.h b/common/DiffractionSpot.h index 2614b4c9..43c00ea1 100644 --- a/common/DiffractionSpot.h +++ b/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); diff --git a/common/SpotToSave.h b/common/SpotToSave.h index d866acbc..491f7a96 100644 --- a/common/SpotToSave.h +++ b/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; diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index 06479019..5ce1b687 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -199,7 +199,34 @@ std::vector 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 &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(std::lround(c.x)), 0), width - 1); + const int row = std::min(std::max(static_cast(std::lround(c.y)), 0), height - 1); + const uint16_t b = pixel_to_bin[static_cast(row) * width + col]; + const float mu = (b < nbins) ? ring_mean[b] : 0.0f; + const double N = std::max(s.PixelCount(), 1); + const double tot = std::max(s.Count(), 0); + const double signal = tot - N * mu; + const double noise = std::sqrt(std::max(1.0, tot + N * static_cast(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(snr - 4.0) / 1.5f)); + const float w_size = (1.0f / (1.0f + std::exp(-(static_cast(N) - 1.5f) / 0.7f))) + * (1.0f / (1.0f + std::exp(-(40.0f - static_cast(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 AdaptiveSpotFinderCPU::RunPersistence(const ImagePr std::vector out; out.reserve(spots.size()); for (auto &kv : spots) out.push_back(kv.second); + 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 c5e046d9..039653a4 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -41,6 +41,10 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { std::vector RunPersistence(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &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 &spots) const; public: explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 43852c70..970df0ae 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/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; }; diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index d272fb3f..da082741 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -37,15 +37,19 @@ void MarkIceRings(std::vector &spots, float tolerance_q_recipA) { } } -void FilterSpotsByCount(std::vector &input, int64_t count) { +void FilterSpotsByCount(std::vector &input, int64_t count, bool by_weight) { size_t output_size = std::min(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; } diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index 9ba458be..f85441d5 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -18,7 +18,7 @@ void CountSpots(DataMessage &msg, void MarkIceRings(std::vector &spots, float tolerance_q_recipA); -void FilterSpotsByCount(std::vector &input, int64_t count); +void FilterSpotsByCount(std::vector &input, int64_t count, bool by_weight = false); void FilterSpuriousHighResolutionSpots(std::vector &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) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index b43984d0..73d3f2fd 100644 --- a/rugnux/rugnux_cli.cpp +++ b/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 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 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; @@ -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(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)