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

Rotation scaling: do not trust a per-frame scale that has collapsed toward zero

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-27 22:43:58 +02:00
co-authored by Claude Opus 5
parent 0cd8cb7ba3
commit 7040987125
1 changed files with 37 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+37 -1
View File
@@ -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<double> &g, int
const int n = static_cast<int>(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<double> 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;
}
}