diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index e5430ad5..e891604d 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -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 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; @@ -279,119 +275,3 @@ void AdaptiveSpotFinderCPU::ApplyWeights(std::vector &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 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); - if (settings.local_snr > 0.0f) - FilterByLocalSNR(image, out, settings.local_snr); - 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 a62dd792..0faf7590 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -34,13 +34,8 @@ 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); // 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. diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index f79c5354..4d1dd099 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -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 diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 624cd654..182d41b2 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -67,7 +67,6 @@ 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 << " --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 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 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 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)