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

Remove threshold-free persistence spot-detection variant

Browse Source 
Drops --persistence-spots and AdaptiveSpotFinderCPU::RunPersistence (the 0-D
topological-persistence detector added in 5a33b0743). It was a research variant
that never beat the hard-threshold adaptive detector on a CC1/2 basis and is a
GPU dead-end (global candidate sort + union-find), so it is not a production
path. The hard-threshold --adaptive-spots detector is unaffected.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-25 12:41:44 +02:00
co-authored by Claude Opus 4.8
parent ca7cbe206a
commit ecf79af018
4 changed files with 0 additions and 142 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
-120
View File
@@ -77,7 +77,6 @@ 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
@@ -120,9 +119,6 @@ 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;
@@ -279,119 +275,3 @@ void AdaptiveSpotFinderCPU::ApplyWeights(std::vector<DiffractionSpot> &spots) co
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)
// (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);
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
-5
View File
@@ -34,13 +34,8 @@ 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);
// 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.
image_analysis/spot_finding/SpotFindingSettings.h
-7
View File
@@ -32,13 +32,6 @@ struct SpotFindingSettings {
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;
// 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
rugnux/rugnux_cli.cpp
-10
View File
@@ -67,7 +67,6 @@ 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 << " --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;
@@ -152,7 +151,6 @@ enum {
OPT_MIN_PIX_PER_SPOT,
OPT_ADAPTIVE_SPOTS,
OPT_SPOT_FALSE_PIXELS,
OPT_PERSISTENCE_SPOTS,
OPT_SOFT_WEIGHT,
OPT_LOCAL_SNR,
OPT_MIN_INDEXED_FRACTION,
@@ -268,7 +266,6 @@ 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},
{"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},
@@ -538,7 +535,6 @@ 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 soft_weight_flag = false;
float local_snr = 0.0f;
std::optional<float> min_indexed_fraction;
@@ -781,11 +777,6 @@ 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_SOFT_WEIGHT:
adaptive_spots = true;
soft_weight_flag = true;
@@ -1543,7 +1534,6 @@ 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;
spot_settings.soft_weight = soft_weight_flag;
spot_settings.local_snr = local_snr;
if (d_min_spot_finding > 0.0f)