mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

Add threshold-free persistence variant of adaptive spot detection

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-23 17:38:18 +02:00
co-authored by Claude Opus 4.8
parent 0ed50d91c9
commit 5a33b07435
4 changed files with 139 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
+117
View File
@@ -4,6 +4,7 @@
#include <algorithm>
#include <bitset>
#include <cmath>
#include <unordered_map>
#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<DiffractionSpot> AdaptiveSpotFinderCPU::Run(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings,
const std::vector<bool> &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<size_t>(width) * height;
@@ -196,3 +201,115 @@ std::vector<DiffractionSpot> 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<DiffractionSpot> AdaptiveSpotFinderCPU::RunPersistence(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings,
const std::vector<bool> &res_mask) {
const auto &pixel_to_bin = mapping.GetPixelToBin();
const size_t nbins = ring_sum.size();
const size_t npix = static_cast<size_t>(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<double>(n_total), 1e-9), 0.1);
const float z = static_cast<float>(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<float>(g_sum / n_total);
const float g_sigma = static_cast<float>(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> 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<int32_t>(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<int32_t> parent;
std::vector<float> 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<int32_t>(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<int32_t>(c)) pers[c] = birth[c] - Z_FLOOR; // survivors
// One spot per surviving maximum whose persistence clears the significance bar.
std::unordered_map<int32_t, DiffractionSpot> 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<DiffractionSpot> out;
out.reserve(spots.size());
for (auto &kv : spots) out.push_back(kv.second);
return out;
}
image_analysis/spot_finding/AdaptiveSpotFinderCPU.h
+5
View File
@@ -33,8 +33,13 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder {
std::vector<float> ring_mean;
std::vector<float> ring_sigma;
std::vector<float> ring_thr;
std::vector<int32_t> 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<DiffractionSpot> RunPersistence(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings,
const std::vector<bool> &res_mask);
public:
explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);
image_analysis/spot_finding/SpotFindingSettings.h
+7
View File
@@ -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;
};
rugnux/rugnux_cli.cpp
+10
View File
@@ -67,6 +67,7 @@ void print_usage() {
std::cout << " --min-pix-per-spot <num> 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 <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 << " --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;
@@ -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<float>(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)