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

Rotation scaling: guard the per-frame scales whatever else is switched on

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-30 17:51:44 +02:00
co-authored by Claude Opus 5
parent a770f0857b
commit ffbf38d2ab
2 changed files with 92 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+86 -30
View File
@@ -1049,43 +1049,72 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &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<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;
// 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<uint8_t> &fitted_mask,
std::vector<double> &g,
std::vector<uint8_t> &apply,
std::vector<double> &ratio) const {
const int n = static_cast<int>(g.size());
apply.assign(n, 0);
ratio.assign(n, 1.0);
std::vector<double> 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> &obs, std::vector<double> &g, int window) const {
const int n = static_cast<int>(g.size());
std::vector<double> g_smooth;
@@ -1103,10 +1132,14 @@ void RotationScaleMerge::SmoothG(std::vector<Obs> &obs, std::vector<double> &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<double> 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<double>(o.corr);
std::vector<double> 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<uint8_t> 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<uint8_t> apply;
std::vector<double> 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<float>(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();
image_analysis/scale_merge/RotationScaleMerge.h
+6
View File
@@ -216,6 +216,12 @@ private:
void ComputeSmoothGWindow(const std::vector<double> &g, int window,
std::vector<double> &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<uint8_t> &fitted_mask, std::vector<double> &g,
std::vector<uint8_t> &apply, std::vector<double> &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