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

rugnux: add detector-plane modulation (flat-field) correction surface

Browse Source 
Rotation scaling gained a third cross-validated correction surface alongside
decay and absorption: a smooth multiplicative factor over the predicted detector
position (px, py) where each full lands. It absorbs detector-response and
geometric flat-field systematics that vary across the detector plane and inflate
R-meas; symmetry-equivalents of one reflection land at different detector
positions as the crystal rotates, over-determining the surface. Because it lives
in the detector frame (not the rotation), the same correction concept applies to
stills.

The surface fit/cross-validate/apply machinery shared with the absorption surface
is factored into ApplyCellSurface; RefineAbsorption and the new RefineModulation
just build their per-full cell assignment (goniometer-frame vs detector-frame)
and call it. The cross-validation now scores a sigma-INDEPENDENT, R-meas-like
fractional agreement of the held-out equivalents rather than a studentized chi^2,
so a surface can no longer pass CV by reshaping sigma without tightening the
intensities - this removes an over-fit regression on mis-indexed data and hardens
the absorption surface too.

On the /data/rotation_test rotation battery (A/B, modulation off vs on): 17 of 22
processed crystals improve R-meas, 0 regress, e.g. lysoC 23.2->16.3%, lyso_2
47.4->28.5%, EcwtAL500 53.2->27.6%, EP_cs_01-17 CC1/2 60->93%; CC1/2 held or
improved everywhere and the anomalous S signal preserved/improved (lyso_ref
SD_MET 4.37->4.46, lysoC 3.07->3.25). On by default with the other surfaces.

Also: --dump-observations now writes the per-full centroid frame and predicted
px/py columns, so the dumped combined fulls are a complete unmerged export.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-16 13:24:10 +02:00
co-authored by Claude Opus 4.8
parent eae7ea90dc
commit 818ec8fcbf
4 changed files with 85 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+59 -21
View File
@@ -171,6 +171,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
// so a no-op when their systematic is absent). Decoupled from the stills-only -B / RefineB flag.
refine_decay_b = s.GetCorrectionSurfaces();
absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0;
modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0;
scaling_iter = std::max(1, scaling_iterations);
resolution_cutoff_method = s.GetResolutionCutoff();
resolution_cc_target = s.GetResolutionCCTarget();
@@ -661,7 +662,6 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) {
return;
constexpr int NB = 8;
const int ncell = NB * NB;
std::vector<int32_t> cell(fulls.size(), -1);
for (size_t i = 0; i < fulls.size(); ++i) {
const Obs &o = fulls[i];
@@ -673,14 +673,46 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) {
const int iy = std::clamp(static_cast<int>((u.y + 1.0f) * 0.5f * NB), 0, NB - 1);
cell[i] = ix * NB + iy;
}
ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Absorption (goniometer-frame 8x8)");
}
void RotationScaleMerge::RefineModulation(int n_iter, int n_groups) {
// Detector-plane modulation (flat-field): a smooth multiplicative factor over the position where a
// reflection lands on the detector (predicted px, py). It absorbs detector-response and geometric
// systematics that vary across the detector plane; symmetry-equivalents of one reflection land at
// different positions as the crystal rotates, over-determining the surface. Because it is a property of
// the detector + beam (not the rotation), the same correction concept transfers to stills.
if (fulls.empty())
return;
float pxmin = std::numeric_limits<float>::infinity(), pxmax = -pxmin, pymin = pxmin, pymax = -pxmin;
for (const Obs &o : fulls) {
if (!std::isfinite(o.px) || !std::isfinite(o.py)) continue;
pxmin = std::min(pxmin, o.px); pxmax = std::max(pxmax, o.px);
pymin = std::min(pymin, o.py); pymax = std::max(pymax, o.py);
}
if (!(pxmax > pxmin) || !(pymax > pymin))
return;
constexpr int NB = 16;
const float sx = NB / (pxmax - pxmin), sy = NB / (pymax - pymin);
std::vector<int32_t> cell(fulls.size(), -1);
for (size_t i = 0; i < fulls.size(); ++i) {
const Obs &o = fulls[i];
if (!std::isfinite(o.px) || !std::isfinite(o.py)) continue;
const int ix = std::clamp(static_cast<int>((o.px - pxmin) * sx), 0, NB - 1);
const int iy = std::clamp(static_cast<int>((o.py - pymin) * sy), 0, NB - 1);
cell[i] = ix * NB + iy;
}
ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Modulation (detector-frame 16x16)");
}
void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int ncell, int n_iter,
int n_groups, const char *name) {
auto usable = [&](const Obs &o) {
return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality;
};
// Fit the per-cell absorption factor over the subset {frame&1 == parity} (parity < 0 = all fulls),
// n_iter alternating rounds against that subset's own reference (Tikhonov pull to 1, gauge-fixed to a
// den-weighted geometric mean of 1 so it never drifts the overall scale). Returns the accumulated
// per-cell factor A[cell].
// Fit the per-cell factor over the subset {frame&1 == parity} (parity < 0 = all fulls), n_iter
// alternating rounds against that subset's own reference (Tikhonov pull to 1, gauge-fixed to a
// den-weighted geometric mean of 1 so it never drifts the overall scale). Returns the per-cell factor.
auto fit_surface = [&](int parity) -> std::vector<double> {
std::vector<double> A(ncell, 1.0);
for (int it = 0; it < n_iter; ++it) {
@@ -725,39 +757,40 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) {
const double sc = static_cast<double>(o.sigma) * o.corr * a, w = 1.0 / (sc * sc);
sw[o.group] += w; swI[o.group] += w * Is;
}
double sum = 0.0; size_t n = 0;
// Rmeas-like (sigma-INDEPENDENT) agreement of the held-out equivalents: sum|Is - Iref| / sum|Iref|.
// Scoring on the studentized deviation instead lets a surface "improve" the held-out chi^2 by
// reshaping sigma (via corr) without tightening the actual intensities - which on mis-indexed / bad
// data passes cross-validation yet worsens Rmeas. A fractional metric cannot be gamed that way.
double num = 0.0, den = 0.0;
for (size_t i = 0; i < fulls.size(); ++i) {
const Obs &o = fulls[i];
if (!usable(o) || cell[i] < 0 || (o.frame & 1) != parity || sw[o.group] <= 0.0) continue;
const double a = A[cell[i]], Is = static_cast<double>(o.I) * o.corr * a;
const double sc = static_cast<double>(o.sigma) * o.corr * a;
const double Iref = swI[o.group] / sw[o.group];
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<double>(n) : 0.0;
return den > 0.0 ? num / den : 0.0;
};
// Cross-validate: fit the surface on even frames and score the held-out odd equivalents (and vice
// versa). A real absorption surface generalizes; an over-fit one (few obs / cell) does not. Apply only
// if the held-out equivalent chi^2 improves by a clear margin - so the correction can never worsen the
// data or act on marginal noise.
// versa). A real surface generalizes; an over-fit one (few obs / cell) does not. Apply only if the
// held-out equivalent chi^2 improves by a clear margin - so the correction can never worsen the data.
const std::vector<double> ident(ncell, 1.0);
const std::vector<double> A_even = fit_surface(0), A_odd = fit_surface(1);
const double base = score(1, ident) + score(0, ident);
const double gain = base - (score(1, A_even) + score(0, A_odd));
if (!(gain > CV_MIN_RELATIVE_GAIN * base)) {
logger.Info("Absorption correction: not cross-validated (held-out gain {:.1f}%, skipped)",
100.0 * gain / std::max(base, 1e-30));
logger.Info("{} correction: not cross-validated (held-out gain {:.1f}%, skipped)",
name, 100.0 * gain / std::max(base, 1e-30));
return;
}
const std::vector<double> A = fit_surface(-1);
for (size_t i = 0; i < fulls.size(); ++i)
if (cell[i] >= 0)
fulls[i].corr = static_cast<float>(fulls[i].corr * A[cell[i]]);
logger.Info("Absorption correction: goniometer-frame {}x{} surface (cross-validated, held-out gain {:.1f}%)",
NB, NB, 100.0 * gain / std::max(base, 1e-30));
logger.Info("{} correction: cross-validated, held-out gain {:.1f}%",
name, 100.0 * gain / std::max(base, 1e-30));
}
void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &g, int window,
@@ -909,7 +942,7 @@ void RotationScaleMerge::Combine() {
if (dump != nullptr)
*dump << full.h << ' ' << full.k << ' ' << full.l << ' ' << full.I << ' '
<< full.sigma << ' ' << d << ' ' << n_frames_event << ' ' << sum_partiality << ' '
<< static_cast<int>(peak_frame) << '\n';
<< static_cast<int>(peak_frame) << ' ' << full.px << ' ' << full.py << '\n';
}
return ev.size();
};
@@ -922,7 +955,9 @@ void RotationScaleMerge::Combine() {
std::ofstream dump;
if (!observation_dump_path.empty()) {
dump.open(observation_dump_path);
dump << "# h k l I sigma d n_frames captured_fraction\n";
// Columns are enough to rescale/merge the combined fulls in an external program: peak_frame
// is a per-full batch (ZD) and px/py the predicted detector position (XD/YD).
dump << "# h k l I sigma d n_frames captured_fraction peak_frame px py\n";
}
fulls.reserve(n_run);
for (int r = 0; r < n_run; ++r)
@@ -1600,11 +1635,14 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search,
// the merged reference. Cheap host loops over the downloaded fulls; applied only on the final
// in-symmetry merge, never the P1 space-group search pass (corrections there add risk and can
// perturb the symmetry determination). ---
[[maybe_unused]] const bool corrections = !for_search && (refine_decay_b || absorption_iter > 0);
[[maybe_unused]] const bool corrections =
!for_search && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0);
if (!for_search && refine_decay_b)
RefineDecay(n_groups);
if (!for_search && absorption_iter > 0)
RefineAbsorption(absorption_iter, n_groups);
if (!for_search && modulation_iter > 0)
RefineModulation(modulation_iter, n_groups);
#ifdef JFJOCH_USE_CUDA
// The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push
// the corrected corr back to the device so the merge reads it.