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.