From 697ffc9a182259545d959a8d8371aa196c97b85e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 8 Oct 2026 06:53:40 +0200 Subject: [PATCH] ResolutionCutoff: a bin's CC1/2 is the median over groups of its reflections The automatic cut fitted its logistic to one Pearson CC1/2 per bin, which has no bound on what one pair can do to it. In a weak shell a single reflection with a huge intensity in both halves - an artefact the halves reproduce, or one observation with a sigma fifty times the shell's - carries nearly all of the covariance. Each bin's reflections are now dealt into up to seven groups by their indices (groups of at least ten reflections), and the bin reads the median of the groups' CC1/2. On a clean shell every group estimates the same CC1/2, so the median is the Pearson value; an outlier takes one group with it and the median does not follow. Only the cut decision changes; the reported shell table is unchanged. Found on a weak split myoglobin crystal: one reflection carried 50-99 % of the covariance in most fit bins past 2.6 A (e.g. 43647 in a bin with median I 1.3), reading them at CC1/2 0.93-0.98 where the median of groups reads 0.2. The cut swung between 1.21 and 2.03 A with --prepass-fraction; it is now 1.85 A at 0.5 and 1.89 A at 1 (XDS 1.51 A, rc.174 1.96 A). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- .../scale_merge/ResolutionCutoff.cpp | 58 +++++++++++++------ tests/MergeScaleTest.cpp | 24 ++++++++ 2 files changed, 65 insertions(+), 17 deletions(-) diff --git a/image_analysis/scale_merge/ResolutionCutoff.cpp b/image_analysis/scale_merge/ResolutionCutoff.cpp index 35c2946ca..9f7b5aa62 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.cpp +++ b/image_analysis/scale_merge/ResolutionCutoff.cpp @@ -4,6 +4,7 @@ #include "ResolutionCutoff.h" #include +#include #include #include @@ -15,6 +16,7 @@ namespace { constexpr int N_FIT_BINS = 25; constexpr int MIN_BIN_COUNT = 10; constexpr int MIN_FIT_BINS = 5; // need at least this many usable bins for a fit + constexpr int MAX_CC_GROUPS = 7; // a bin's CC1/2 is the median over this many groups (odd) constexpr int EXTEND_BINS_PAST_FALLOFF = 2; // bins kept beyond the first sub-target bin constexpr double SHELLS_FOR_EXTENSION = 10.0; // "+1 shell" = one 10-shell width in s (report-independent) @@ -80,12 +82,13 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector 0.0f && std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1]); + }; double s_lo = std::numeric_limits::max(), s_hi = 0.0; double d_data_min = std::numeric_limits::max(); for (const auto &m : merged) { - if (m.on_ice_ring) - continue; - if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) + if (!usable(m)) continue; const double s = 1.0 / (static_cast(m.d) * m.d); s_lo = std::min(s_lo, s); @@ -99,28 +102,49 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector bin_cc(N_FIT_BINS); - std::vector bin_n(N_FIT_BINS, 0); - for (const auto &m : merged) { - if (m.on_ice_ring) - continue; - if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) - continue; + auto bin_of = [&](const MergedReflection &m) { const double s = 1.0 / (static_cast(m.d) * m.d); - int b = static_cast((s - s_lo) / bin_w); - b = std::clamp(b, 0, N_FIT_BINS - 1); - bin_cc[b].Add(m.I_half[0], m.I_half[1], m.cc_weight); - ++bin_n[b]; + return std::clamp(static_cast((s - s_lo) / bin_w), 0, N_FIT_BINS - 1); + }; + std::vector bin_n(N_FIT_BINS, 0); + for (const auto &m : merged) + if (usable(m)) ++bin_n[bin_of(m)]; + + // A bin's CC1/2 is the median of the Pearson CC1/2 of a few groups its reflections are dealt into by + // their indices, not one Pearson over the whole bin. A Pearson CC has no bound on + // what one pair can do to it: in a weak shell a single reflection with a huge intensity in both + // halves - an artefact the halves reproduce, or one observation with a sigma fifty times the shell's + // - carries nearly all of the covariance and reads the shell at CC1/2 0.98 where the rest of it is + // noise, and the cut runs on through it. Here such a reflection takes one group with it, and the + // median does not follow. On a shell without one, every group estimates the same CC1/2 and the + // median is the Pearson value. Groups of fewer than MIN_BIN_COUNT reflections are not formed: a bin + // too small for three groups is one group, the plain Pearson CC1/2. + std::vector bin_groups(N_FIT_BINS); + for (int b = 0; b < N_FIT_BINS; ++b) { + const int g = std::min(MAX_CC_GROUPS, std::max(1, bin_n[b] / MIN_BIN_COUNT)); + bin_groups[b] = g % 2 == 1 ? g : g - 1; + } + std::vector> group_cc(N_FIT_BINS); + for (const auto &m : merged) { + if (!usable(m)) continue; + const int b = bin_of(m); + const uint32_t hash = static_cast(m.h) * 73856093u ^ static_cast(m.k) * 19349663u + ^ static_cast(m.l) * 83492791u; + group_cc[b][hash % bin_groups[b]].Add(m.I_half[0], m.I_half[1], m.cc_weight); } // Usable bins (enough counts), in ascending-s order, with their bin-centre s. std::vector s_bin, cc_bin; for (int b = 0; b < N_FIT_BINS; ++b) { if (bin_n[b] < MIN_BIN_COUNT) continue; - const double cc = bin_cc[b].GetCC(); - if (!std::isfinite(cc)) continue; + std::vector cc; + for (int g = 0; g < bin_groups[b]; ++g) + if (std::isfinite(group_cc[b][g].GetCC())) + cc.push_back(group_cc[b][g].GetCC()); + if (cc.empty()) continue; + std::nth_element(cc.begin(), cc.begin() + cc.size() / 2, cc.end()); s_bin.push_back(s_lo + (b + 0.5) * bin_w); - cc_bin.push_back(cc); + cc_bin.push_back(cc[cc.size() / 2]); } if (static_cast(s_bin.size()) < MIN_FIT_BINS) { result.note = "too few usable CC1/2 bins"; diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 6dcb3c6dc..74c3819f0 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -316,6 +316,7 @@ namespace { std::normal_distribution g(0.0, sigma); for (int j = 0; j < n; ++j) { MergedReflection m; + m.h = static_cast(v.size()); // every reflection its own index m.d = static_cast(1.0 / std::sqrt(s_from + (j + 0.5) * (s_to - s_from) / n)); const double a = g(rng), b = g(rng); m.I_half[0] = static_cast(a); @@ -345,6 +346,29 @@ TEST_CASE("ResolutionCutoff_CleanFallOff") { CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15)); } +// Past the fall-off, a few single reflections read hundreds of times the shell's noise in both halves - +// an artefact the half-sets reproduce. Each carries nearly all of its shell's covariance, so a Pearson +// CC1/2 reads those shells as signal and the cut runs on through them; one reflection per shell must not +// decide where the signal stops. +TEST_CASE("ResolutionCutoff_SingleReflectionsDoNotHoldTheCutOpen") { + Logger logger("test"); + std::mt19937 rng(12345); + std::vector merged; + AddBand(merged, rng, 0.01, 0.25, 480, true); + AddBand(merged, rng, 0.25, 0.51, 520, false); + for (const double s : {0.26, 0.28, 0.30, 0.32}) { + MergedReflection m; + m.h = static_cast(merged.size()); + m.d = static_cast(1.0 / std::sqrt(s)); + m.I_half[0] = m.I_half[1] = 300.0f; + merged.push_back(m); + } + + const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger); + REQUIRE(rc.d_fit); + CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15)); +} + // A sweep with a long stretch where the crystal barely diffracted: the reflections measured only there // are noise scaled up by 1/G, sigmas with them, and there are many of them at every resolution. Counted // by the information they carry they must not hide where the well-measured reflections stop.