From 82b8a9752c61d9453c64ef0a3348fe126eb15e31 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 23:57:04 +0200 Subject: [PATCH] rotation: fix mosaicity underestimate + default-on capture-uncertainty (ISa 10.7->19.1) The rotation per-image mosaicity was ~3x too small (0.045 vs the true 0.13deg), which crippled the partiality model and capped per-observation precision: it predicted reflections on too few frames and over-peaked the rocking partiality, so the rot3d-combined fulls were ~1.7x noisier than XDS's, the integration bottleneck behind the jungfraujoch-vs-XDS ISa gap. Two root causes, both fixed: - CalcMosaicityXDS (Kabsch-2010 MLE) searched each spot's exact-Bragg phi only within +-wedge (the 0.2deg oscillation). Reflections recorded at larger rocking offset - the tail that defines the mosaic width - fell outside and were dropped, truncating the tau distribution so the MLE underestimated ~2x. Widen the search window to wedge+0.8deg; the MLE then converges to the true 0.13deg (and is insensitive to widening further, since it weights by the recorded fraction). - ScaleOnTheFly then re-refined the mosaicity from the intensity residual, which is degenerate with the per-image scale G and collapses it toward its floor. Trust the (now correct) indexing mosaicity and keep it fixed during scaling. With the correct mosaicity, --capture-uncertainty (which down-weights the over-extrapolated under-captured fulls) now pays off strongly, so default it ON (1.0) for the rot3d combine; it stays off for non-rot3d. Together on the HEWL rotation crystal: ISa 10.7 -> 19.1, and anomalous peak height vs XDS goes from 52% to ~78% (CL_CL 1.92x -> 1.29x). This reaches XDS's own published-correction ceiling (DECAY+ABSORP+MODPIX ~= 19.6); the remaining gap to its quoted ISa 28 is the I->inf extrapolation. No effect on the stills path (rotation-only code). Co-Authored-By: Claude Opus 4.8 --- image_analysis/indexing/AnalyzeIndexing.cpp | 10 ++++++++-- image_analysis/scale_merge/ScaleOnTheFly.cpp | 6 ++++-- tools/jfjoch_process.cpp | 11 +++++++---- 3 files changed, 19 insertions(+), 8 deletions(-) diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 5da8f709..17497199 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -250,8 +250,14 @@ namespace { + bstar * static_cast(s.k) + cstar * static_cast(s.l); - // Find predicted phi angle - const auto phi_pred_opt = predict_phi_deg_local(pstar, S0, m2, 0.0f, axis.GetWedge_deg()); + // Find predicted phi angle. The search window must be wide enough to catch reflections + // recorded at large rocking offset (|tau| up to ~mosaicity + dphi/2). Using ±wedge alone + // clips the tau tail at the oscillation width, so the MLE then underestimates the mosaicity + // ~2x (the tail reflections are exactly the ones that define the mosaic width). A generous + // window (oscillation + ~0.8deg rocking allowance) lets the tail in; the MLE is insensitive + // to making it wider still (it weights by the recorded fraction R(tau), which decays). + const float window_deg = axis.GetWedge_deg() + 0.8f; + const auto phi_pred_opt = predict_phi_deg_local(pstar, S0, m2, 0.0f, window_deg); if (!phi_pred_opt.has_value() || !std::isfinite(phi_pred_opt.value())) continue; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 46fc6cea..98b08a7a 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -310,8 +310,10 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { } else { problem.SetParameterBlockConstant(&result.wedge); } - problem.SetParameterLowerBound(&result.mos, 0, s.GetMinMosaicity()); - problem.SetParameterUpperBound(&result.mos, 0, s.GetMaxMosaicity()); + // Trust the (now correctly estimated) indexing mosaicity; do NOT re-refine it here. The + // per-image scaling residual fit is degenerate between G and the mosaicity, and it collapses + // the mosaicity toward its floor (3x too small) which corrupts the partiality. Keep it fixed. + problem.SetParameterBlockConstant(&result.mos); } ceres::Solver::Options options; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 4ae99cc4..9a1001f9 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -56,7 +56,7 @@ void print_usage() { std::cout << " -w, --wedge[=num] Refine image wedge during scaling with starting wedge value" << std::endl; std::cout << " --scaling-high-resolution High resolution limit for spot finding (default: no limit)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; - std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 0 = off)" << std::endl; + std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: off; e.g. 6, XDS/DIALS-style)" << std::endl; std::cout << " --reject-delta-cchalf Per-crystal CC1/2-delta rejection: drop images with deltaCChalf below mean - N*stddev (default: off; e.g. 2.5)" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; @@ -288,7 +288,7 @@ int main(int argc, char **argv) { std::string ref_column; std::string dump_observations; // diagnostic: dump unmerged -P rot3d fulls to this path double min_partiality = 0.02; - double capture_uncertainty_coeff = 0.0; + std::optional capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d double min_image_cc = 0.0; int64_t scaling_iter = 3; std::optional forced_rotation_lattice; @@ -510,7 +510,7 @@ int main(int argc, char **argv) { min_partiality = std::stod(optarg); break; case OPT_CAPTURE_UNCERTAINTY: - capture_uncertainty_coeff = std::stod(optarg); + capture_uncertainty_arg = std::stod(optarg); break; case OPT_INTEGRATION_RADIUS: integration_radius_arg = optarg; @@ -715,7 +715,10 @@ int main(int argc, char **argv) { if (wedge_for_scaling.has_value()) scaling_settings.RotationWedgeForScaling(wedge_for_scaling); scaling_settings.MinPartiality(min_partiality); - scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_coeff); + // Capture-aware systematic sigma defaults ON (1.0) for the rot3d combine - it down-weights the + // over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous + // substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins. + scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(combine_3d ? 1.0 : 0.0)); scaling_settings.MinCCForImage(min_image_cc); if (outlier_reject_nsigma) scaling_settings.OutlierRejectNsigma(*outlier_reject_nsigma);