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

Remove --soft-weight and --local-snr spot-finder options

Browse Source 
Both were opt-in adaptive-spot refinements that did not help. Soft per-spot
weighting was index-rate neutral across the battery (re-ranking only bites when
spots exceed the max-spot cap, which weak serial data does not reach). The
local-SNR gate was neutral on index rate and degraded merged CC1/2 on flooded
XFEL data. Drops the flags, ApplyWeights/FilterByLocalSNR, the per-spot weight
field, and the by-weight FilterSpotsByCount branch (now strongest-first only).
--adaptive-spots itself is unchanged.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-25 17:31:35 +02:00
co-authored by Claude Opus 4.8
parent f72b4484e2
commit 20bbcb1cd3
10 changed files with 7 additions and 143 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
+1 -79
View File
@@ -195,83 +195,5 @@ std::vector<DiffractionSpot> AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB
output_buffer[OutputSize() - 1] = out.to_ulong();
// --- Stage D: connected components + resolution mask + min/max-pix (shared with classic path) ---
auto spots = ExtractSpots(image, settings, res_mask);
if (settings.local_snr > 0.0f)
FilterByLocalSNR(image, spots, settings.local_snr);
if (settings.soft_weight)
ApplyWeights(spots);
return spots;
}
// Reject spots that do not stand out against their LOCAL background. The loose per-ring threshold
// keeps ~100% of real Bragg peaks but, on structured-background (XFEL) frames, also floods with
// spurious pixels: the azimuthal ring mean underestimates the diffuse level in some sectors, so a
// locally-high background pixel clears it. A real peak stands many sigmas above the background in
// its IMMEDIATE neighbourhood, a flood pixel does not. For each spot the background mean and scatter
// are measured from a square annulus around its centroid (robust median / MAD, so a neighbouring
// peak in the annulus cannot bias it), and the spot is kept only if its integrated signal exceeds
// k local sigmas. k is in sigma units -- self-calibrating, no photon threshold.
void AdaptiveSpotFinderCPU::FilterByLocalSNR(const ImagePreprocessorBuffer &image,
std::vector<DiffractionSpot> &spots, float k) const {
constexpr int RIN = 3; // half-width of the excluded core (7x7)
constexpr int ROUT = 6; // half-width of the background annulus (13x13)
std::vector<float> bg;
bg.reserve((2 * ROUT + 1) * (2 * ROUT + 1));
std::vector<DiffractionSpot> kept;
kept.reserve(spots.size());
for (const auto &s : spots) {
const Coord c = s.RawCoord();
const int col = static_cast<int>(std::lround(c.x));
const int row = static_cast<int>(std::lround(c.y));
bg.clear();
for (int dr = -ROUT; dr <= ROUT; ++dr) {
const int rr = row + dr;
if (rr < 0 || rr >= height) continue;
for (int dc = -ROUT; dc <= ROUT; ++dc) {
if (std::max(std::abs(dr), std::abs(dc)) <= RIN) continue; // skip the peak core
const int cc = col + dc;
if (cc < 0 || cc >= width) continue;
const int32_t v = image[static_cast<size_t>(rr) * width + cc];
if (v == INT32_MIN || v == INT32_MAX) continue; // masked / saturated
bg.push_back(static_cast<float>(v));
}
}
if (bg.size() < 8) { kept.push_back(s); continue; } // too few bg pixels to judge
const size_t mid = bg.size() / 2;
std::nth_element(bg.begin(), bg.begin() + mid, bg.end());
const float bg_med = bg[mid];
for (auto &v : bg) v = std::fabs(v - bg_med);
std::nth_element(bg.begin(), bg.begin() + mid, bg.end());
const float sigma = std::max(1.4826f * bg[mid], 1.0f);
const double npix = static_cast<double>(std::max<int64_t>(s.PixelCount(), 1));
const double signal = static_cast<double>(s.Count()) - bg_med * npix;
const double snr = signal / (sigma * std::sqrt(npix));
if (snr >= static_cast<double>(k))
kept.push_back(s);
}
spots.swap(kept);
}
void AdaptiveSpotFinderCPU::ApplyWeights(std::vector<DiffractionSpot> &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<int>(std::lround(c.x)), 0), width - 1);
const int row = std::min(std::max(static_cast<int>(std::lround(c.y)), 0), height - 1);
const uint16_t b = pixel_to_bin[static_cast<size_t>(row) * width + col];
const float mu = (b < nbins) ? ring_mean[b] : 0.0f;
const double N = std::max<int64_t>(s.PixelCount(), 1);
const double tot = std::max<int64_t>(s.Count(), 0);
const double signal = tot - N * mu;
const double noise = std::sqrt(std::max(1.0, tot + N * static_cast<double>(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<float>(snr - 4.0) / 1.5f));
const float w_size = (1.0f / (1.0f + std::exp(-(static_cast<float>(N) - 1.5f) / 0.7f)))
* (1.0f / (1.0f + std::exp(-(40.0f - static_cast<float>(N)) / 8.0f)));
s.SetWeight(std::min(std::max(w_snr * w_size, 0.0f), 1.0f));
}
return ExtractSpots(image, settings, res_mask);
}