From 7040987125f18d906f08017db4624b97c40fa2a7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 27 Jul 2026 22:43:58 +0200 Subject: [PATCH] Rotation scaling: do not trust a per-frame scale that has collapsed toward zero The per-frame scale enters every intensity as 1/G, and SolveScaleIRLS floors G at zero and nothing else. A frame whose fit is not determined by its data can return G ~ 0.002 against a run median of 0.865, and every observation it carries is then multiplied by ~500 - sigma by the identical factor, which is why no sigma-based outlier test can see it and why this looked for a long time like a partiality problem. (The 1/partiality path is in fact guarded: min_captured_fraction floors it at 0.7 by default on rotation.) The window smoothing that should have absorbed such a frame instead made it permanent. It averages log G over a window, so a scale collapsing toward zero does not merely corrupt its own frame - its logarithm drags the whole window down. Worse, where a run has a stretch of frames too sparse to fit at all, the only FITTED frames in a window can be the collapsed ones, and the geometric mean then averages the fault with itself. Measured on a multi-lattice dataset: frames 816 and 818 fitted G = 0.0023 and 0.0014 with every neighbour from 800 to 839 unfitted, so smoothing set G = 0.0018 across the whole neighbourhood - a 546x amplification. About 500 observations of 152000 (0.66%) then carried 99% of sum(I^2), and the merged CC1/2 read 17.2% where the same data with the classic finder read 93.7%. Treat a fitted scale far below the run's median as what it is - an undetermined scale, exactly like the too-few-reflections case the code already handles - rather than as a successful fit. Such frames no longer contribute to the smoothing mean, and a frame whose own scale is not credible takes the neighbourhood's, or the run's typical scale when the neighbourhood holds nothing credible either. The bound is a RATIO to the run's own median because the rotation per-frame G is not gauge-fixed: G and the group means have an exact global multiplicative degeneracy, and the fitted median drifts over 0.745-1.358 across the battery. An absolute floor would reject everything in a run that drifted low. MIN_CREDIBLE_SCALE_RATIO = 0.02 was chosen from measurement over 12 crystals in the default configuration, where the smallest legitimate min(G)/median(G) is 0.070; the failing case sat at 0.0017. It is 3.5x below anything real and 12x above the failure. Effect on the intensity tail of the failing case: max I 10224 -> 438, and the top 1000 observations' share of sum(I^2) 0.990 -> 0.421 (the classic-finder reference is 0.632, so the tail is now cleaner than the run this was compared against). Rotation battery, 33 crystals in the default configuration: ZERO crystals differ - no space group, CC1/2, high-shell CC or ISa change anywhere. The guard fires only on the pathology. It does NOT rescue that dataset: with the amplification gone its CC1/2 is 26.2% and R_meas 49.2% against the classic finder's 93.7% and 27.8%. Adaptive detection degrades those intensities for a second, independent reason that is still open. This commit removes a latent hazard for any run with a sparse stretch of frames; it is not the fix for that dataset. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 38 ++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 290bef99..c08cbf3b 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -31,6 +31,17 @@ namespace { // implementation is numerically identical - see the comments there for the details. constexpr size_t MIN_REFLECTIONS = 20; // per-frame scale needs at least this many constexpr double SCALE_ROBUST_K = 3.0; // Cauchy loss scale (sigma units) for the per-frame G fit + + // A fitted per-frame scale below this fraction of the run's median is not a measurement of anything: + // it would mean the frame received 2% of its neighbours' dose while still producing indexable spots. + // It matters because the scale enters as 1/G, so a G that collapses toward zero multiplies every + // intensity on that frame without bound - and sigma by the identical factor, which is why no + // sigma-based outlier test can see it. Measured over 12 rotation crystals in the default + // configuration, the smallest LEGITIMATE min(G)/median(G) is 0.070; the case this bound exists for + // sat at 0.0017 (a 546x amplification that put 99% of sum(I^2) into 0.66% of the observations and + // took the merged CC1/2 from 93.7% to 17.2%). 0.02 is 3.5x below anything real and 12x above the + // failure. + constexpr double MIN_CREDIBLE_SCALE_RATIO = 0.02; constexpr float MAX_FRAME_GAP = 2.0f; // a rocking event is a run of frames no more apart than this constexpr double CHI2_1_MEDIAN = 0.454936; // A post-scale-fulls correction surface (decay / absorption) is applied only if its held-out @@ -1034,16 +1045,41 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int const int n = static_cast(g.size()); const int half = window / 2; g_smooth.assign(n, NAN); + + // The run's typical scale, and the floor below which a fitted scale is not credible + // (MIN_CREDIBLE_SCALE_RATIO). The window average is a GEOMETRIC mean, so a scale collapsing toward + // zero does not just corrupt its own frame - its logarithm drags the whole window down. Where a run + // has a stretch of frames too sparse to fit, the only fitted frames in a window can BE the collapsed + // ones, and the mean then averages the fault with itself and makes it permanent. + std::vector fitted; + fitted.reserve(n); + for (int j = 0; j < n; ++j) + if (frame_scaled_scratch[j] && std::isfinite(g[j]) && g[j] > 0.0) + fitted.push_back(g[j]); + double g_typ = NAN; + if (!fitted.empty()) { + const size_t mid = fitted.size() / 2; + std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); + g_typ = fitted[mid]; + } + const double g_floor = std::isfinite(g_typ) ? g_typ * MIN_CREDIBLE_SCALE_RATIO : 0.0; + for (int o = 0; o < n; ++o) { double sum_log = 0.0; int count = 0; for (int j = std::max(0, o - half); j <= std::min(n - 1, o + half); ++j) { - if (frame_scaled_scratch[j] && std::isfinite(g[j]) && g[j] > 0.0) { + if (frame_scaled_scratch[j] && std::isfinite(g[j]) && g[j] >= g_floor && g[j] > 0.0) { sum_log += std::log(g[j]); ++count; } } if (count > 0) g_smooth[o] = std::exp(sum_log / count); + // A frame whose own scale is not credible takes the neighbourhood's, or - when the neighbourhood + // holds nothing credible either - the run's typical scale. Callers divide g by g_smooth, so this + // replaces the collapsed scale rather than layering another correction on top of it. + if (frame_scaled_scratch[o] && std::isfinite(g[o]) && g[o] > 0.0 && g[o] < g_floor + && !std::isfinite(g_smooth[o])) + g_smooth[o] = g_typ; } }