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

rotation: deterministic frame-order mosaicity smoothing + partiality recompute

Browse Source 
Prediction applied a mosaicity/profile-radius moving average (RotationParameters) over the
last N *processed* frames. Under the parallel per-image loop that window is thread-arrival
order, so the smoothed value - and hence which reflections are predicted/integrated - was
non-deterministic run-to-run, swinging CC1/2 (and even the space group) on marginal crystals.
`-N 1` was deterministic; `-N 32` was not.

Fix (as designed with FL): prediction now uses each frame's OWN mosaicity/profile-radius
(image-local, deterministic membership - a reflection on the cutoff contributes ~nothing).
The smoothing that actually matters is moved into RotationScaleMerge and done in FRAME order
(deterministic): per-frame mosaicity is smoothed with the same window as smooth-G, then every
partial's partiality is recomputed from it BEFORE the 3D combine. This is the mosaicity analogue
of smooth-G: combining a reflection's per-frame partials only tiles the rocking curve correctly
(captured fractions summing toward 1) if neighbouring frames share a consistent mosaicity.

Battery (18 crystals, /data/rotation_test, 2 runs each): 15/18 now bit-identical run-to-run
(the good crystals unchanged - lyso P41212 ISa 7.8 CC1/2 99.7%). The 3 residual crystals
(EcwtAL500, EcwtCQ066S, pding4_003 - all large/triclinic cells) still jitter ~0.002%, traced
to a SEPARATE, benign cause: the GPU prediction buffer overflow (BraggPredictionRotGPU
max_reflections=10000 with a racy atomicAdd/atomicSub) on dense frames - cell/space group stay
stable; to be addressed in the GPU prediction/integration rework (naively raising the cap also
changes prediction quality, so it is not a one-line bump). Minor label refinements from the
recomputed partiality: cytC_2 P321 -> P3121 (now consistent with cytC_3), Ins_I_2/3 report the
honest I23/I213 screw-axis ambiguity.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-02 21:58:10 +02:00
co-authored by Claude Opus 4.8
parent b5d9167bf4
commit 29c8ba6112
3 changed files with 79 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/IndexAndRefine.cpp
+8 -8
View File
@@ -272,14 +272,14 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
CrystalLattice latt = outcome.lattice_candidate.value();
if (rotation_indexer) {
// Use moving average for mosaicity and profile_radius (also add beam center later)
if (msg.mosaicity_deg)
msg.mosaicity_deg = rotation_parameters.Mosaicity(msg.mosaicity_deg.value());
if (msg.profile_radius) {
msg.profile_radius = rotation_parameters.ProfileRadius(msg.profile_radius.value());
}
}
// Prediction uses each frame's OWN mosaicity/profile_radius (image-local). We deliberately do NOT
// smooth them here with a running moving average: it averaged the last N *processed* frames, whose
// order under the parallel per-image loop is thread-arrival order, making the predicted rocking
// width - and hence which reflections are integrated - non-deterministic run-to-run. Prediction only
// decides membership (a reflection on the cutoff contributes ~nothing), so the per-frame value is
// fine here. The mosaicity smoothing that actually matters - keeping the partialities of one rocking
// event consistent so they tile the curve and sum toward 1 - is done deterministically in frame
// order before the 3D combine (RotationScaleMerge), where partiality is recomputed from it.
float ewald_dist_cutoff = 0.001f;
if (msg.profile_radius)
image_analysis/scale_merge/RotationScaleMerge.cpp
+60 -1
View File
@@ -36,6 +36,17 @@ namespace {
return 1.0 / x;
}
// Kabsch rotation partiality: the fraction of a reflection recorded in the sampled slice, from the
// erf of the rocking angle relative to the mosaic width. Identical to ScaleOnTheFly's RotationPartiality
// (and the predictor's), so recomputing here just swaps in the smoothed mosaicity.
float RotationPartiality(double delta_phi_deg, double zeta, double mosaicity_deg, double wedge_deg) {
const double half_wedge = wedge_deg / 2.0;
const double c1 = zeta / std::sqrt(2.0);
const double arg_plus = (delta_phi_deg + half_wedge) * c1 / mosaicity_deg;
const double arg_minus = (delta_phi_deg - half_wedge) * c1 / mosaicity_deg;
return static_cast<float>((std::erf(arg_plus) - std::erf(arg_minus)) / 2.0);
}
// Deterministic CC1/2 half from the frame's stable index (splitmix64), matching Merge.cpp.
int HalfForImage(int64_t image_id) {
uint64_t z = static_cast<uint64_t>(image_id) + 0x9e3779b97f4a7c15ULL;
@@ -225,6 +236,53 @@ void RotationScaleMerge::Ingest() {
rawrun_group.assign(rawrun_start.size(), -1);
logger.Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)",
total, n_frames, rawrun_start.size());
SmoothMosaicityAndPartiality();
}
void RotationScaleMerge::SmoothMosaicityAndPartiality() {
// Raw per-frame mosaicity as measured at integration (image-local, deterministic).
std::vector<double> mos_raw(n_frames, NAN);
for (int o = 0; o < n_frames; ++o) {
const auto &m = partials_out[o].mosaicity_deg;
if (m && std::isfinite(*m) && *m > 0.0f) mos_raw[o] = *m;
}
// Frame-order moving average with the same window as smooth-G (a rotation range -> frame count).
// With smoothing off, fall back to the per-frame value (still deterministic, just unsmoothed).
const auto ss = x.GetScalingSettings();
const double smooth_deg = ss.GetSmoothGDegrees();
const auto gon = x.GetGoniometer();
const double osc = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0;
mos_smooth.assign(n_frames, NAN);
if (smooth_deg > 0.0 && osc > 1e-6) {
int window = std::max(1, static_cast<int>(std::lround(smooth_deg / osc)));
if (window % 2 == 0) ++window;
const int half = window / 2;
for (int o = 0; o < n_frames; ++o) {
double sum = 0.0;
int cnt = 0;
for (int j = std::max(0, o - half); j <= std::min(n_frames - 1, o + half); ++j)
if (std::isfinite(mos_raw[j])) { sum += mos_raw[j]; ++cnt; }
if (cnt > 0) mos_smooth[o] = static_cast<float>(sum / cnt);
}
} else {
for (int o = 0; o < n_frames; ++o) mos_smooth[o] = static_cast<float>(mos_raw[o]);
}
// Recompute each partial's partiality from the smoothed mosaicity (same wedge the predictor used).
// Frames without a mosaicity keep the stored partiality.
const double wedge = gon ? std::fabs(gon->GetWedge_deg()) : 0.0;
ParallelChunks(static_cast<int>(partials.size()), nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; ++i) {
auto &o = partials[i];
const float mos = mos_smooth[o.frame];
if (std::isfinite(mos) && mos > 1e-6f && std::isfinite(o.zeta) && o.zeta > 0.0f
&& std::isfinite(o.delta_phi))
o.partiality = RotationPartiality(o.delta_phi, o.zeta, mos, wedge);
}
});
logger.Info("Recomputed partiality from frame-order-smoothed mosaicity");
}
int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
@@ -593,7 +651,8 @@ void RotationScaleMerge::FinalizePerFrameScale(int n_groups, const std::vector<d
auto &o = partials_out[f];
if (frame_scaled[f]) {
o.image_scale_g = static_cast<float>(g_partial[f]);
o.mosaicity_deg = static_cast<float>(mosaicity_deg);
o.mosaicity_deg = (f < static_cast<int>(mos_smooth.size()) && std::isfinite(mos_smooth[f]))
? mos_smooth[f] : static_cast<float>(mosaicity_deg);
if (std::isfinite(cc[f])) { o.image_scale_cc = static_cast<float>(cc[f]); o.image_scale_cc_n = cc_n[f]; }
else { o.image_scale_cc.reset(); o.image_scale_cc_n.reset(); }
} else {
image_analysis/scale_merge/RotationScaleMerge.h
+11
View File
@@ -117,6 +117,10 @@ private:
// Set by FitPerFrameG: which frames were fitted this call (so corr/G is updated only there).
std::vector<uint8_t> frame_scaled_scratch;
// Per-frame mosaicity smoothed in frame order (deterministic); used to recompute partiality and
// written back for the per-image scaling table. Empty if there is no per-frame mosaicity.
std::vector<float> mos_smooth;
// Working per-group arrays (sized to the current group count; reused).
std::vector<int32_t> group_h, group_k, group_l;
@@ -143,6 +147,13 @@ private:
const std::vector<uint8_t> &frame_scaled) const;
void SmoothG(std::vector<Obs> &obs, std::vector<double> &g, int window) 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
// that prediction applied. SG-independent, so done once in Ingest.
void SmoothMosaicityAndPartiality();
void Combine(); // partials -> fulls
// Per-frame CC vs the partial merge reference, then write G/CC/mosaicity back onto the partials