diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 813782e4..35180df5 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -554,9 +554,10 @@ void RotationScaleMerge::RefineDecay(int n_groups) { // noise / aggregate structure (not real per-reflection decay) would amplify into a high-res // mis-correction and worsen the symmetry-equivalent scatter that sets ISa. An opt-in correction must // never add noise where its systematic is absent (a low-dose reference set), so the slope is - // CROSS-VALIDATED: fit on even-frame fulls and score the held-out odd-frame equivalent chi^2, and vice - // versa. Real decay generalizes across the split (both held-out scores improve); a noise slope does - // not. Apply the full-data slope only if the mean held-out chi^2 strictly improves; else it is a no-op. + // CROSS-VALIDATED: fit on even-frame fulls and score the held-out odd-frame equivalent agreement (a + // sigma-independent Rmeas-like metric), and vice versa. Real decay generalizes across the split (both + // held-out scores improve); a noise slope does not. Apply the full-data slope only if the mean held-out + // agreement strictly improves; else it is a no-op. if (fulls.empty() || n_frames <= 0) return; const double fcenter = 0.5 * n_frames; @@ -593,9 +594,14 @@ void RotationScaleMerge::RefineDecay(int n_groups) { // strong damage (slope = 2 dB/dframe, so +-1.0 admits total relative-B up to ~n_frames/2 A^2). return (var > 0.0) ? std::clamp((Sxy - Sx * Sy / Sw) / var, -1.0, 1.0) : 0.0; }; - // Mean studentized squared deviation over the subset {frame&1 == parity} when `slope` is applied, - // scored against that subset's own reference (no leakage). Lower = tighter equivalents. - auto subset_chi2 = [&](int parity, double slope) -> double { + // Rmeas-like (sigma-INDEPENDENT) disagreement of the equivalents over the subset {frame&1 == parity} + // when `slope` is applied, scored against that subset's own reference (no leakage): sum|Is-Iref|/sum|Iref|. + // Scoring on the studentized deviation instead would let the slope "improve" the held-out chi^2 by + // reshaping sigma - it scales sigma up on late / high-angle frames, down-weighting them - without + // actually tightening the intensities; a fractional metric cannot be gamed that way. This mirrors the + // absorption / modulation gate in ApplyCellSurface (which was switched off studentized chi^2 for the + // same reason); the reference is still the inverse-variance mean. + auto subset_disagreement = [&](int parity, double slope) -> double { std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); for (const auto &o : fulls) { if (!usable(o) || (o.frame & 1) != parity) continue; @@ -604,17 +610,15 @@ void RotationScaleMerge::RefineDecay(int n_groups) { const double w = 1.0 / (sc * sc); sw[o.group] += w; swI[o.group] += w * Is; } - double sum = 0.0; size_t n = 0; + double num = 0.0, den = 0.0; for (const auto &o : fulls) { if (!usable(o) || (o.frame & 1) != parity || sw[o.group] <= 0.0) continue; const double Iref = swI[o.group] / sw[o.group]; const double Is = static_cast(o.I) * o.corr * decay_factor(o, slope); - const double sc = static_cast(o.sigma) * o.corr * decay_factor(o, slope); - if (!std::isfinite(Iref) || !(sc > 0.0)) continue; - const double dd = (Is - Iref) / sc; - sum += dd * dd; ++n; + if (!std::isfinite(Iref) || Iref <= 0.0) continue; + num += std::abs(Is - Iref); den += Iref; } - return n > 0 ? sum / static_cast(n) : 0.0; + return den > 0.0 ? num / den : 0.0; }; // A meaningful radiation-damage relative-B accrues several A^2 across a dataset. Below a physical floor @@ -631,8 +635,8 @@ void RotationScaleMerge::RefineDecay(int n_groups) { } // And require the slope to cross-validate by a clear margin: fit on even frames, score the held-out // odd equivalents (and vice versa). Real damage generalizes; a fluke does not. - const double base = subset_chi2(1, 0.0) + subset_chi2(0, 0.0); - const double gain = base - (subset_chi2(1, fit_slope(0)) + subset_chi2(0, fit_slope(1))); + const double base = subset_disagreement(1, 0.0) + subset_disagreement(0, 0.0); + const double gain = base - (subset_disagreement(1, fit_slope(0)) + subset_disagreement(0, fit_slope(1))); if (!(gain > CV_MIN_RELATIVE_GAIN * base)) { logger.Info("Decay correction: not cross-validated (dB = {:.2f} A^2, held-out gain {:.1f}%, skipped)", total_delta_B, 100.0 * gain / std::max(base, 1e-30));