diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 338b1290..459a17e1 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -13,6 +13,8 @@ ADD_LIBRARY(JFJochScaleMerge RotationScaleMerge.h ResolutionCutoff.cpp ResolutionCutoff.h + IceRingMask.cpp + IceRingMask.h HKLKey.cpp HKLKey.h RfreeFlags.cpp diff --git a/image_analysis/scale_merge/IceRingMask.cpp b/image_analysis/scale_merge/IceRingMask.cpp new file mode 100644 index 00000000..eb68aefb --- /dev/null +++ b/image_analysis/scale_merge/IceRingMask.cpp @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "IceRingMask.h" +#include "../../common/Definitions.h" // ICE_RING_RES_A +#include "../../common/CorrelationCoefficient.h" + +std::vector FindDecorrelatedIceRings(const std::vector &merged, + float half_width_q_recipA, Logger &logger) { + if (merged.empty()) + return {}; + + constexpr float two_pi = 6.283185307f; + const float w = half_width_q_recipA; + std::vector mask(ICE_RING_RES_A.size(), 0); + + for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i) { + const float q_ring = two_pi / ICE_RING_RES_A[i]; + CorrelationCoefficient ring, shoulder; + size_t n_ring = 0, n_shoulder = 0; + for (const auto &m : merged) { + if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) + continue; + const float dq = std::fabs(two_pi / m.d - q_ring); + if (dq < w) { ring.Add(m.I_half[0], m.I_half[1]); ++n_ring; } + else if (dq < 3.0f * w) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; } + } + if (n_ring >= 20 && n_shoulder >= 20 && shoulder.GetCC() > 0.5 + && ring.GetCC() < shoulder.GetCC() - 0.05) { + mask[i] = 1; + logger.Info("Ice-ring mask: {:.2f} A ring CC1/2 {:.3f} << shoulders {:.3f}; masked from merge", + ICE_RING_RES_A[i], ring.GetCC(), shoulder.GetCC()); + } + } + + if (std::none_of(mask.begin(), mask.end(), [](char c) { return c != 0; })) + return {}; + return mask; +} diff --git a/image_analysis/scale_merge/IceRingMask.h b/image_analysis/scale_merge/IceRingMask.h new file mode 100644 index 00000000..5fdaa604 --- /dev/null +++ b/image_analysis/scale_merge/IceRingMask.h @@ -0,0 +1,20 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../../common/Logger.h" +#include "../../common/Reflection.h" // MergedReflection + +// Hexagonal-ice rings whose merged half-set CC1/2 has collapsed well below the resolution shoulders +// either side: ice has decorrelated them and their Bragg intensity is unrecoverable, so the caller +// should re-merge with them dropped. A weak or absent ring tracks its neighbours and is not +// returned, which is what keeps completeness on a clean crystal untouched. +// +// Returns a mask indexed like ICE_RING_RES_A, or an empty vector when no ring qualifies - so the +// caller can test emptiness rather than scanning. Shared by the full pipeline and the offline +// --scale path, which have to reach the same verdict on the same data. +std::vector FindDecorrelatedIceRings(const std::vector &merged, + float half_width_q_recipA, Logger &logger); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 5ef45192..fe837de5 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -34,6 +34,7 @@ #include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/IceRingMask.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" #include "../image_analysis/scale_merge/ScalingResult.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" @@ -1652,28 +1653,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // re-merge. Weak/absent rings track their neighbours and stay, so completeness on clean crystals // is untouched. Uses the ring band (+/- ice width in q=2pi/d) vs the shoulders either side. if (experiment_.IsDetectIceRings() && !sm.merged.empty()) { - constexpr float two_pi = 6.283185307f; - const float w = config_.spot_finding.ice_ring_width_Q_recipA; - std::vector mask(ICE_RING_RES_A.size(), 0); - for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i) { - const float q_ring = two_pi / ICE_RING_RES_A[i]; - CorrelationCoefficient ring, shoulder; - size_t n_ring = 0, n_shoulder = 0; - for (const auto &m : sm.merged) { - if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) - continue; - const float dq = std::fabs(two_pi / m.d - q_ring); - if (dq < w) { ring.Add(m.I_half[0], m.I_half[1]); ++n_ring; } - else if (dq < 3.0f * w) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; } - } - if (n_ring >= 20 && n_shoulder >= 20 && shoulder.GetCC() > 0.5 - && ring.GetCC() < shoulder.GetCC() - 0.05) { - mask[i] = 1; - logger.Info("Ice-ring mask: {:.2f} A ring CC1/2 {:.3f} << shoulders {:.3f}; masked from merge", - ICE_RING_RES_A[i], ring.GetCC(), shoulder.GetCC()); - } - } - if (std::any_of(mask.begin(), mask.end(), [](char c) { return c != 0; })) { + auto mask = FindDecorrelatedIceRings(sm.merged, config_.spot_finding.ice_ring_width_Q_recipA, + logger); + if (!mask.empty()) { masked_ice_rings = std::move(mask); const auto final_sg = experiment_.GetGemmiSpaceGroup(); sm = scale_and_merge(final_sg ? final_sg->short_name() : "P1", false); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 5f75fe6b..387dcaa5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -31,6 +31,7 @@ #include "../image_analysis/scale_merge/StillsPartialityRefine.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/IceRingMask.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "Rugnux.h" @@ -1112,6 +1113,10 @@ static int RunRugnux(int argc, char **argv) { experiment.SpaceGroupNumber(space_group_number); if (fixed_reference_unit_cell.has_value()) experiment.SetUnitCell(fixed_reference_unit_cell); + // --detect-ice-rings, applied here as well as on the full path below: this block returns + // before that one runs, so without it the flag is silently ignored by --scale. + if (detect_ice_rings.has_value()) + experiment.DetectIceRings(detect_ice_rings.value()); // A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling. IndexingSettings indexing_settings; @@ -1155,10 +1160,34 @@ static int RunRugnux(int argc, char **argv) { logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); experiment.ImagesPerTrigger(refl_stats.n_images); + // Ice-ring handling, as the full pipeline does it (Rugnux.cpp): flag reflections on a + // hexagonal-ice powder ring so scaling skips them while the merge keeps them. The flag is not + // stored per reflection, so it has to be recomputed here from the resolution just assigned - + // otherwise --scale re-scales a dataset the writing run had scaled without those reflections, + // and the per-image scales come out of a different fit than the ones in the file. + const float ice_width = SpotFindingSettings().ice_ring_width_Q_recipA; + if (experiment.IsDetectIceRings()) { + size_t total = 0, flagged = 0; + for (auto &outcome : reflections) { + for (auto &r : outcome.reflections) { + ++total; + r.on_ice_ring = IsOnIceRing(r.d, ice_width); + if (r.on_ice_ring) + ++flagged; + } + } + logger.Info("Ice-ring handling: flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); " + "excluded from scaling, kept for merging", flagged, total, ice_width); + } + const auto scale_start = std::chrono::steady_clock::now(); std::vector merged_reflections; MergeStatistics merged_statistics; double error_model_isa = 0.0; + // Ice rings dropped from the merge because their CC1/2 collapsed. Decided from a first merge + // and applied by a second, as the full pipeline does - flagging alone only keeps the ice + // reflections out of the SCALE fit, which on its own costs a little and buys nothing. + std::vector masked_ice_rings; // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> // merge, including the default-on decay + absorption correction surfaces. It does not support @@ -1172,13 +1201,25 @@ static int RunRugnux(int argc, char **argv) { throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Rotation scaling/merging (RotationScaleMerge) does not support reference " "scaling or wedge refinement"); + // The ice half-width has to be the real one: it is what turns a reflection's resolution + // into a ring index, so a zero here makes every ice test inside the merge a no-op. RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), - scaling_iter, 0.0f, nthreads, logger); + scaling_iter, ice_width, nthreads, logger); rsm.Ingest(); - auto r = rsm.Run(false); - merged_reflections = std::move(r.merged); - merged_statistics = std::move(r.statistics); - error_model_isa = r.isa; + // Ingest() is separate from Run() precisely so the merge can be repeated; the ice-ring + // mask below needs a first merge before it can be decided. + auto run = [&](const std::vector &masked) { + auto r = rsm.Run(false, masked); + merged_reflections = std::move(r.merged); + merged_statistics = std::move(r.statistics); + error_model_isa = r.isa; + }; + run({}); + if (experiment.IsDetectIceRings()) { + masked_ice_rings = FindDecorrelatedIceRings(merged_reflections, ice_width, logger); + if (!masked_ice_rings.empty()) + run(masked_ice_rings); + } } else { // Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports // CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image @@ -1194,35 +1235,48 @@ static int RunRugnux(int argc, char **argv) { const double mean_tilt = refiner.Run(reflections, nthreads); logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); } - MergeOnTheFly merge_engine(experiment); - merge_engine.ReferenceCell(experiment.GetUnitCell()); - // --min-image-cc has to hold for the merge itself, not only for the reported statistics. - merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); - // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full - // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity - // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor). - merge_engine.RefineErrorModel(reflections); - if (merge_engine.ErrorModelActive()) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), - merge_engine.ErrorModelB(), - merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, - merge_engine.ErrorModelChi2()); - for (size_t i = 0; i < reflections.size(); ++i) - merge_engine.AddImage(reflections[i], static_cast(i)); - merged_reflections = merge_engine.ExportReflections(); + // The merge alone, repeatable for the ice-ring mask below. The scaling above is NOT redone: + // it has already been applied to `reflections`, and running it twice would compound the + // correction. + auto merge = [&](const std::vector &masked) { + MergeOnTheFly merge_engine(experiment); + merge_engine.ReferenceCell(experiment.GetUnitCell()); + // --min-image-cc has to hold for the merge itself, not only for the reported statistics. + merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); + if (!masked.empty()) + merge_engine.MaskIceRings(masked, ice_width); + // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full + // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity + // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor). + merge_engine.RefineErrorModel(reflections); + if (merge_engine.ErrorModelActive()) + logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), + merge_engine.ErrorModelB(), + merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, + merge_engine.ErrorModelChi2()); + for (size_t i = 0; i < reflections.size(); ++i) + merge_engine.AddImage(reflections[i], static_cast(i)); + merged_reflections = merge_engine.ExportReflections(); - // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual - // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells - // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.) - const auto &cut_ss = experiment.GetScalingSettings(); - // The offline --scale path re-scales a stored _process.h5 and is never a P1 search merge. - const std::optional effective_d_min = ApplyResolutionCutoff( - merged_reflections, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(), - cut_ss.GetResolutionCCTarget(), /*for_search=*/false, logger); + // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual + // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells + // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.) + const auto &cut_ss = experiment.GetScalingSettings(); + // The offline --scale path re-scales a stored _process.h5 and is never a P1 search merge. + const std::optional effective_d_min = ApplyResolutionCutoff( + merged_reflections, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(), + cut_ss.GetResolutionCCTarget(), /*for_search=*/false, logger); - merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, - effective_d_min); - error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, + effective_d_min); + error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + }; + merge({}); + if (experiment.IsDetectIceRings()) { + masked_ice_rings = FindDecorrelatedIceRings(merged_reflections, ice_width, logger); + if (!masked_ice_rings.empty()) + merge(masked_ice_rings); + } } logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)",