rugnux: handle ice rings in --scale as the full pipeline does

--scale did none of the ice handling the run that wrote the _process.h5 had
done, so re-scaling a stored dataset silently produced a different - and
flatteringly more complete - answer than the pipeline it was meant to
reproduce. Three separate gaps:

  * --detect-ice-rings was accepted and ignored. The --scale block returns
    before the line that applies it.
  * Reflections were never flagged as sitting on an ice ring, so the per-image
    scale fit included them. The flag is not stored per reflection, so it has
    to be recomputed from the resolution.
  * RotationScaleMerge was constructed with the ice half-width hardcoded to
    zero. That is what turns a resolution into a ring index, so every ice test
    inside the merge was a no-op whatever was passed to it.

The CC1/2 ring test that decides which rings to drop moves into
FindDecorrelatedIceRings, shared with the full pipeline so both reach the same
verdict on the same data, and --scale now re-merges with the mask the way the
pipeline does. The stills branch re-runs only the merge: the scaling has
already been applied to the reflections and repeating it would compound it.

Measured on a rotation dataset with three decorrelated rings, --scale went
from 8765 unique / 36.3% completeness / R-meas 18.5% / <I/sig> 1.1 to
7638 / 31.6% / 18.0% / 1.3, against the full pipeline's 7692 / 31.8% / 17.9% /
1.3 - the reported completeness had been inflated by reflections the pipeline
drops. The full pipeline is bit-identical across the refactor.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

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

image_analysis/scale_merge/IceRingMask.cpp

@@ -0,0 +1,42 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <algorithm>
+#include <cmath>
+
+#include "IceRingMask.h"
+#include "../../common/Definitions.h"        // ICE_RING_RES_A
+#include "../../common/CorrelationCoefficient.h"
+
+std::vector<char> FindDecorrelatedIceRings(const std::vector<MergedReflection> &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<char> 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;
+}

image_analysis/scale_merge/IceRingMask.h

@@ -0,0 +1,20 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+
+#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<char> FindDecorrelatedIceRings(const std::vector<MergedReflection> &merged,
+                                           float half_width_q_recipA, Logger &logger);

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<char> 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);

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<MergedReflection> 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<char> 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<char> &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<int64_t>(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<char> &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<int64_t>(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<double> 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<double> 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)",