From ffbf38d2abd8ec6d93bc3c016129457b5b25618a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 30 Jul 2026 17:51:44 +0200 Subject: [PATCH] Rotation scaling: guard the per-frame scales whatever else is switched on The protection against a per-frame scale collapsing toward zero lived inside ComputeSmoothGWindow, so it only existed when smooth-G did: --smooth-g=0, a dataset whose oscillation width is unknown, and any caller that never sets a smoothing range - the viewer among them - merged with no guard at all. A collapsed G multiplies that frame's intensities by 1/G and its sigmas by the same factor, so nothing downstream can see it; the merge's n-sigma cut scales with the number that is wrong. Pull it out into ReplaceCollapsedScales, called unconditionally right after the partial scaling loop, and let the smooth-G window assume what it now guarantees instead of computing its own median and floor. The fulls guard built its median from every frame including those never fitted - those sit at the combine's corr = 1, so a run with many unfitted frames dragged the median toward 1 and the floor with it. It also reported the absolute amplification where the message says "below the run median". Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 116 +++++++++++++----- .../scale_merge/RotationScaleMerge.h | 6 + 2 files changed, 92 insertions(+), 30 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 87162f0f..ffc7281e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1049,43 +1049,72 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int 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; - + // No collapsed-scale handling here: ReplaceCollapsedScales has already run on g, so every fitted + // entry is credible. That matters because the window average is a GEOMETRIC mean - a scale + // collapsing toward zero would not just corrupt its own frame, its logarithm would drag the whole + // window down, and where the only fitted frames in a window ARE the collapsed ones the mean would + // average the fault with itself and make it permanent. 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] >= g_floor && g[j] > 0.0) { + if (frame_scaled_scratch[j] && std::isfinite(g[j]) && 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; } } +// Replace any per-frame scale that fitted far below the run's median with that median, and report the +// per-frame ratio by which the frame's corr has to be rescaled to match (1 where nothing changed). +// Returns false when there was nothing to replace. +// +// A collapsed G multiplies the frame's intensities by 1/G AND its sigmas by the same factor, so no +// downstream test can recognise it: the merge's n-sigma outlier cut scales with the very number that is +// wrong. Only a total collapse gives itself away, by overflowing corr to infinity. +bool RotationScaleMerge::ReplaceCollapsedScales(const std::vector &fitted_mask, + std::vector &g, + std::vector &apply, + std::vector &ratio) const { + const int n = static_cast(g.size()); + apply.assign(n, 0); + ratio.assign(n, 1.0); + + std::vector fitted; + fitted.reserve(n); + for (int f = 0; f < n; ++f) + if (fitted_mask[f] && std::isfinite(g[f]) && g[f] > 0.0) + fitted.push_back(g[f]); + if (fitted.empty()) + return false; + + const size_t mid = fitted.size() / 2; + std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); + const double g_typ = fitted[mid]; + const double g_floor = g_typ * MIN_CREDIBLE_SCALE_RATIO; + + int n_replaced = 0; + double worst = 1.0; + for (int f = 0; f < n; ++f) { + if (!fitted_mask[f] || !std::isfinite(g[f]) || g[f] <= 0.0 || g[f] >= g_floor) + continue; + // corr = rlp / (partiality * G), so swapping G for g_typ scales this frame's corr by G / g_typ. + apply[f] = 1; + ratio[f] = g[f] / g_typ; + worst = std::max(worst, g_typ / g[f]); + g[f] = g_typ; + ++n_replaced; + } + if (n_replaced > 0) + logger.Warning("Replaced the fitted scale of {} frame(s) that came out more than {:.0f}x below " + "the run median (worst {:.0f}x) - a scale that small amplifies the frame's " + "intensities, and its sigmas with them, by the same factor", + n_replaced, 1.0 / MIN_CREDIBLE_SCALE_RATIO, worst); + return n_replaced > 0; +} + void RotationScaleMerge::SmoothG(std::vector &obs, std::vector &g, int window) const { const int n = static_cast(g.size()); std::vector g_smooth; @@ -1103,10 +1132,14 @@ void RotationScaleMerge::SmoothG(std::vector &obs, std::vector &g, } bool RotationScaleMerge::RejectCollapsedFullScales() { - // The Unity model leaves corr = 1/G, and every full of a frame carries that frame's G. + // The Unity model leaves corr = 1/G, and every full of a frame carries that frame's G. A full whose + // frame was NOT fitted still carries the 1.0 the combine gave it, on both the host and the device + // path - those are not measurements of anything and must stay out of the median, or a run with many + // unfitted frames drags the median toward 1 and the floor with it. std::vector g_frame(n_frames, NAN); for (const auto &o : fulls) - if (o.frame >= 0 && o.frame < n_frames && std::isfinite(o.corr) && o.corr > 0.0f) + if (o.frame >= 0 && o.frame < n_frames && std::isfinite(o.corr) && o.corr > 0.0f + && o.corr != 1.0f) g_frame[o.frame] = 1.0 / static_cast(o.corr); std::vector fitted; @@ -1118,7 +1151,8 @@ bool RotationScaleMerge::RejectCollapsedFullScales() { return false; const size_t mid = fitted.size() / 2; std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); - const double g_floor = fitted[mid] * MIN_CREDIBLE_SCALE_RATIO; + const double fitted_median = fitted[mid]; + const double g_floor = fitted_median * MIN_CREDIBLE_SCALE_RATIO; int n_rejected = 0; double worst = 1.0; @@ -1128,14 +1162,14 @@ bool RotationScaleMerge::RejectCollapsedFullScales() { const double gf = g_frame[o.frame]; if (!std::isfinite(gf) || gf >= g_floor) continue; - worst = std::max(worst, 1.0 / gf); + worst = std::max(worst, fitted_median / gf); o.corr = 1.0f; // as if the frame's scale had never been fitted ++n_rejected; } if (n_rejected > 0) logger.Warning("Rejected the fitted scale of {} full(s) whose frame scaled to less than {:.0f}x " "below the run median (worst {:.0f}x amplification); those frames merge unscaled", - n_rejected, 1.0 / MIN_CREDIBLE_SCALE_RATIO, worst / fitted[mid]); + n_rejected, 1.0 / MIN_CREDIBLE_SCALE_RATIO, worst); return n_rejected > 0; } @@ -1953,6 +1987,28 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, } const std::vector partial_scaled = frame_scaled_scratch; + // --- 1b. Guard the per-frame partial scales. Unconditional: the smooth-G window below used to be + // where this lived, so --smooth-g 0, a dataset with no oscillation width, and every caller + // that never sets a smoothing range at all (the viewer among them) had no guard whatsoever. + { + std::vector apply; + std::vector ratio; + if (ReplaceCollapsedScales(frame_scaled_scratch, g_partial, apply, ratio)) { + bool applied_on_gpu = false; +#ifdef JFJOCH_USE_CUDA + if (gpu_active_) { + // Same kernel smooth-G uses: multiply the resident corr by a per-frame ratio in place. + gpu_->SmoothCorr(apply.data(), ratio.data()); + applied_on_gpu = true; + } +#endif + if (!applied_on_gpu) + for (auto &o : partials) + if (apply[o.frame] && std::isfinite(o.corr)) + o.corr = static_cast(o.corr * ratio[o.frame]); + } + } + // --- 2. Smooth G across frames (XDS DELPHI-like) before the combine. --- const auto s = x.GetScalingSettings(); const double smooth_g_deg = s.GetSmoothGDegrees(); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 413863d9..12ffa0ff 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -216,6 +216,12 @@ private: void ComputeSmoothGWindow(const std::vector &g, int window, std::vector &g_smooth) const; + // Replace any fitted per-frame scale that collapsed far below the run median with that median, + // reporting the per-frame corr ratio the caller has to apply. See the .cpp for why nothing + // downstream can catch a collapsed scale on its own. + bool ReplaceCollapsedScales(const std::vector &fitted_mask, std::vector &g, + std::vector &apply, std::vector &ratio) const; + // Smooth per-frame mosaicity in frame order and recompute each partial's partiality from it, so the // per-frame partials of one rocking event tile the curve consistently (they sum toward 1) before the // 3D combine. Deterministic (frame order); replaces the old arrival-order mosaicity moving average