diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 89c938d7..ffd6caf3 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -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) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 17285937..1d881e70 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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((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(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 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(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(sum / cnt); + } + } else { + for (int o = 0; o < n_frames; ++o) mos_smooth[o] = static_cast(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(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(g_partial[f]); - o.mosaicity_deg = static_cast(mosaicity_deg); + o.mosaicity_deg = (f < static_cast(mos_smooth.size()) && std::isfinite(mos_smooth[f])) + ? mos_smooth[f] : static_cast(mosaicity_deg); if (std::isfinite(cc[f])) { o.image_scale_cc = static_cast(cc[f]); o.image_scale_cc_n = cc_n[f]; } else { o.image_scale_cc.reset(); o.image_scale_cc_n.reset(); } } else { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index f3d7cdb6..8fafb1d0 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -117,6 +117,10 @@ private: // Set by FitPerFrameG: which frames were fitted this call (so corr/G is updated only there). std::vector 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 mos_smooth; + // Working per-group arrays (sized to the current group count; reused). std::vector group_h, group_k, group_l; @@ -143,6 +147,13 @@ private: const std::vector &frame_scaled) const; void SmoothG(std::vector &obs, std::vector &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