From 5a33b07435403af5f33af9149474ee54270171a9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 23 Jul 2026 17:38:18 +0200 Subject: [PATCH] Add threshold-free persistence variant of adaptive spot detection Add --persistence-spots, a second parameter-free detector alongside --adaptive-spots. Instead of a hard per-ring threshold it builds the noise-normalised image z = (I - ring_mean) / sqrt(ring_sigma^2 + read^2) (same per-ring background as the hard variant) and scores every intensity maximum by its 0-D topological persistence: sweeping the height from high to low, each maximum is born and, when its basin meets a taller one at a saddle, dies with persistence = birth - saddle, in sigma. A lone noise spike merges into the background almost immediately (persistence ~1 sigma); a real peak stands many sigma proud. Emitting maxima whose persistence clears the same z(E) significance bar needs no photon threshold and no min-pix, and it deblends touching peaks (each keeps its own maximum). Implemented with the same union-find idiom as the connected-component labeller. On serial stills this auto-adapts with no per-dataset tuning like --adaptive-spots, finding fewer but cleaner (deblended) spots; the hard-threshold variant remains more sensitive on the very weakest data. Both share the per-ring background and read-noise floor. comp_of is allocated lazily so the default and hard-adaptive paths pay nothing. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../spot_finding/AdaptiveSpotFinderCPU.cpp | 117 ++++++++++++++++++ .../spot_finding/AdaptiveSpotFinderCPU.h | 5 + .../spot_finding/SpotFindingSettings.h | 7 ++ rugnux/rugnux_cli.cpp | 10 ++ 4 files changed, 139 insertions(+) diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index ef4eece4..06479019 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include "AdaptiveSpotFinderCPU.h" @@ -76,6 +77,7 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping & ring_mean.assign(nbins, 0.0f); ring_sigma.assign(nbins, 0.0f); ring_thr.assign(nbins, 0.0f); + // comp_of is allocated lazily (only the persistence variant needs it). } // Accumulate per-ring mean/variance from the raw (photon) image. clip_k <= 0 -> use every valid @@ -118,6 +120,9 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask) { + if (settings.spot_persistence) + return RunPersistence(image, settings, res_mask); + const auto &pixel_to_bin = mapping.GetPixelToBin(); const size_t nbins = ring_sum.size(); const size_t npix = static_cast(width) * height; @@ -196,3 +201,115 @@ std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB // --- Stage D: connected components + resolution mask + min/max-pix (shared with classic path) --- return ExtractSpots(image, settings, res_mask); } + +// Threshold-free variant. Build a noise-normalised image z = (I - ring_mean)/sqrt(ring_sigma^2+READ^2) +// (same per-ring background and read-noise floor as the hard variant), then score every intensity +// maximum by its 0-D topological persistence: sweep the height from high to low, each maximum is +// "born" and, when its basin meets a taller one at a saddle, "dies" with persistence = birth - saddle +// (in sigma units). A lone noise spike merges into the sea almost immediately (persistence ~1); a real +// peak stands many sigma proud. Emitting maxima with persistence >= z(E) needs no photon threshold and +// no min-pix, and deblends touching peaks (each keeps its own maximum). Union-find, same idiom as the +// classic connected-component labeller. This is the offline/rugnux "soft" alternative to the hard cut. +std::vector AdaptiveSpotFinderCPU::RunPersistence(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; + if (comp_of.size() != npix) + comp_of.assign(npix, -1); + + AccumulateRings(image, 0.0f); + AccumulateRings(image, 3.0f); + AccumulateRings(image, 3.0f); + + 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 = std::min(std::max(E / static_cast(n_total), 1e-9), 0.1); + const float z = static_cast(NormalQuantile(1.0 - p)); + const float READ = 1.0f; + const float PERS_THR = z; // a maximum must stand z sigmas above its saddle to be a spot + const float Z_FLOOR = 2.0f; // loose landscape floor: a compute bound, not a detection threshold + const float g_mean = static_cast(g_sum / n_total); + const float g_sigma = static_cast(std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * (double)g_mean))); + + auto mu_of = [&](uint16_t b) { return ring_cnt[b] < MIN_RING_PIXELS ? g_mean : ring_mean[b]; }; + auto se_of = [&](uint16_t b) { + const float s = ring_cnt[b] < MIN_RING_PIXELS ? g_sigma : ring_sigma[b]; + return std::sqrt(s * s + READ * READ); + }; + + // Candidate pixels: everything above a loose noise-normalised floor (in-resolution, not masked). + struct Cand { float z; int32_t pxl; }; + std::vector cand; + for (size_t pxl = 0; pxl < npix; ++pxl) { + if (res_mask[pxl]) continue; + const int32_t v = image[pxl]; + if (v == INT32_MIN) continue; + const uint16_t b = pixel_to_bin[pxl]; + if (b >= nbins) continue; + const float zz = (v == INT32_MAX) ? 1.0e6f : (v - mu_of(b)) / se_of(b); + if (zz > Z_FLOOR) cand.push_back({zz, static_cast(pxl)}); + } + if (cand.empty()) + return {}; + std::sort(cand.begin(), cand.end(), [](const Cand &a, const Cand &b) { return a.z > b.z; }); + + // Union-find over candidates, processed highest first. parent/birth/pers are per-component. + std::vector parent; + std::vector birth, pers; + parent.reserve(cand.size()); birth.reserve(cand.size()); pers.reserve(cand.size()); + auto find = [&](int32_t c) { while (parent[c] != c) { parent[c] = parent[parent[c]]; c = parent[c]; } return c; }; + + for (const auto &cd : cand) { + const int32_t pxl = cd.pxl; + const int32_t col = pxl % width, row = pxl / width; + int32_t roots[8]; int nr = 0; + for (int dr = -1; dr <= 1; ++dr) for (int dc = -1; dc <= 1; ++dc) { + if (dr == 0 && dc == 0) continue; + const int rr = row + dr, cc = col + dc; + if (rr < 0 || rr >= height || cc < 0 || cc >= width) continue; + const int32_t np = rr * width + cc; + if (comp_of[np] < 0) continue; // neighbour not yet processed (lower z) + const int32_t r = find(comp_of[np]); + bool dup = false; + for (int i = 0; i < nr; ++i) if (roots[i] == r) dup = true; + if (!dup && nr < 8) roots[nr++] = r; + } + if (nr == 0) { // new maximum born + const int32_t c = static_cast(parent.size()); + parent.push_back(c); birth.push_back(cd.z); pers.push_back(1.0e9f); + comp_of[pxl] = c; + } else { + int32_t tall = roots[0]; + for (int i = 1; i < nr; ++i) if (birth[roots[i]] > birth[tall]) tall = roots[i]; + for (int i = 0; i < nr; ++i) + if (roots[i] != tall) { pers[roots[i]] = birth[roots[i]] - cd.z; parent[roots[i]] = tall; } + comp_of[pxl] = tall; + } + } + for (size_t c = 0; c < parent.size(); ++c) + if (parent[c] == static_cast(c)) pers[c] = birth[c] - Z_FLOOR; // survivors + + // One spot per surviving maximum whose persistence clears the significance bar. + std::unordered_map spots; + for (const auto &cd : cand) { + const int32_t root = find(comp_of[cd.pxl]); + if (pers[root] < PERS_THR) continue; + const int32_t pxl = cd.pxl; + const int64_t v = (image[pxl] == INT32_MAX) ? 65535 : image[pxl]; + spots[root].AddPixel(pxl % width, pxl / width, v); + } + + for (const auto &cd : cand) comp_of[cd.pxl] = -1; // reset for the next frame (touched pixels only) + + std::vector out; + out.reserve(spots.size()); + for (auto &kv : spots) out.push_back(kv.second); + return out; +} diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index c28f7290..42424f53 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -33,8 +33,13 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { std::vector ring_mean; std::vector ring_sigma; std::vector ring_thr; + std::vector comp_of; // per-pixel component id for the persistence variant (-1 = unset) void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k); + // Threshold-free variant: 0-D topological persistence on the noise-normalised image. + std::vector RunPersistence(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask); public: explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index b69b7072..43852c70 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -31,4 +31,11 @@ struct SpotFindingSettings { // per frame (the threshold's operating point), ~100 for a multi-megapixel detector. bool adaptive_threshold = false; float false_pixels_per_frame = 100.0f; + + // Threshold-free variant of the adaptive detector (implies adaptive_threshold): instead of a hard + // per-ring cut, score each intensity maximum by its topological persistence (how many sigma it + // stands above the saddle joining it to higher ground) on the noise-normalised image. Persistence + // 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; }; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 6d526ce6..b43984d0 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -67,6 +67,7 @@ void print_usage() { 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 << " --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 << " --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; @@ -148,6 +149,7 @@ enum { OPT_MIN_PIX_PER_SPOT, OPT_ADAPTIVE_SPOTS, OPT_SPOT_FALSE_PIXELS, + OPT_PERSISTENCE_SPOTS, OPT_SPOT_RESOLUTION, OPT_SPOT_LOW_RESOLUTION, OPT_MAX_SPOTS, @@ -260,6 +262,7 @@ static option long_options[] = { {"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}, + {"persistence-spots", no_argument, nullptr, OPT_PERSISTENCE_SPOTS}, {"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}, @@ -526,6 +529,7 @@ int main(int argc, char **argv) { int64_t min_pix_per_spot = 2; bool adaptive_spots = false; float false_pixels_per_frame = 100.0f; + bool persistence_spots = false; bool refine_bfactor = false; std::string ref_mtz; std::string ref_column; @@ -765,6 +769,11 @@ int main(int argc, char **argv) { adaptive_spots = true; logger.Info("Adaptive spot detection: expected false pixels/frame set to {:.0f}", false_pixels_per_frame); break; + case OPT_PERSISTENCE_SPOTS: + adaptive_spots = true; + persistence_spots = true; + logger.Info("Threshold-free (topological-persistence) spot detection enabled"); + 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); @@ -1506,6 +1515,7 @@ int main(int argc, char **argv) { 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; + spot_settings.spot_persistence = persistence_spots; if (d_min_spot_finding > 0.0f) spot_settings.high_resolution_limit = d_min_spot_finding; if (d_max_spot_finding > 0.0f)