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; } }