merge: mask ice rings whose CC1/2 collapses below their resolution neighbours

An ice-contaminated ring decorrelates its reflections, so its merged half-set
CC1/2 falls well below the Bragg trend of its resolution shoulders (e.g.
EP_cs_01-17 2.25A ring CC1/2 0.12 vs shoulders 0.89). After the final merge,
compute per-hexagonal-ice-ring CC1/2 in the ring band (+/- ice width in q=2pi/d)
against the shoulders either side, mask any ring more than 0.05 below its
shoulders (reliable shoulder CC1/2 > 0.5, >=20 reflections each) and re-merge
without them. Weak/absent rings track their neighbours and stay, so clean
crystals are untouched.

Data-driven and robust on /data/rotation_test: 16/18 crystals mask nothing and
are unchanged (every clean crystal, plus icy-but-fine EP_cs_02-10); EP_cs_01-17
masks 6 rings, CC1/2 6.9% -> 28%; CQ066/pding4_003 mask one marginal ring each.
A well-scaled CC1/2 test, replacing the fragile per-image azimuthal cutoffs.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

common/Definitions.h

@@ -70,3 +70,16 @@ inline bool IsOnIceRing(float d_A, float half_width_q_recipA) {
             return true;
     return false;
 }
+
+// Index into ICE_RING_RES_A of the hexagonal-ice ring resolution d sits on (within half_width in
+// q = 2*pi/d), or -1 if none. Used to look up that ring's per-image strength for the merge gate.
+inline int IceRingIndex(float d_A, float half_width_q_recipA) {
+    if (!(d_A > 0.0f))
+        return -1;
+    constexpr float two_pi = 6.283185307f;
+    const float q = two_pi / d_A;
+    for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i)
+        if (std::fabs(q - two_pi / ICE_RING_RES_A[i]) < half_width_q_recipA)
+            return static_cast<int>(i);
+    return -1;
+}

image_analysis/scale_merge/Merge.cpp

@@ -15,6 +15,7 @@
 
 #include "../../common/CorrelationCoefficient.h"
 #include "../../common/ResolutionShells.h"
+#include "../../common/Definitions.h"
 #include "HKLKey.h"
 
 namespace {
@@ -54,6 +55,13 @@ MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional<UnitCell> &cell)
     return *this;
 }
 
+bool MergeOnTheFly::IsMaskedRing(const Reflection &r) const {
+    if (masked_ice_rings.empty())
+        return false;
+    const int ring = IceRingIndex(r.d, mask_ice_half_width_q);
+    return ring >= 0 && ring < static_cast<int>(masked_ice_rings.size()) && masked_ice_rings[ring];
+}
+
 void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id, bool cc_mask) {
     std::unique_lock ul(merged_mutex);
 
@@ -72,6 +80,8 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id
             continue;
         if (exclude_ice_rings && r.on_ice_ring)
             continue;
+        if (IsMaskedRing(r))
+            continue;
         if (r.partiality < min_partiality)
             continue;
 
@@ -168,6 +178,8 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &outc
                 continue;
             if (exclude_ice_rings && r.on_ice_ring)
                 continue;
+            if (IsMaskedRing(r))
+                continue;
             if (r.partiality < min_partiality)
                 continue;
             const float I_corr = r.I * r.image_scale_corr;

image_analysis/scale_merge/Merge.h

@@ -73,6 +73,11 @@ class MergeOnTheFly {
     // see the ice-contaminated intensities. The final in-symmetry merge keeps them (for completeness).
     bool exclude_ice_rings = false;
 
+    // Ice rings (indices into ICE_RING_RES_A) to exclude from this merge because their merged CC1/2
+    // collapsed relative to their resolution neighbours - ice decorrelated them. Empty = none masked.
+    std::vector<char> masked_ice_rings;
+    float mask_ice_half_width_q = 0.03f;
+
     HKLKeyGenerator generator;
 
     std::map<uint64_t, MergeAccum> accumulator;
@@ -103,10 +108,14 @@ class MergeOnTheFly {
     size_t reject_count = 0;
 
     bool Mask(const IntegrationOutcome &outcome, bool cc_mask);
+    [[nodiscard]] bool IsMaskedRing(const Reflection &r) const;
 public:
     MergeOnTheFly(const DiffractionExperiment &x);
     MergeOnTheFly& ReferenceCell(const std::optional<UnitCell> &cell);
     MergeOnTheFly& ExcludeIceRings(bool input) { exclude_ice_rings = input; return *this; }
+    MergeOnTheFly& MaskIceRings(std::vector<char> masked, float half_width_q) {
+        masked_ice_rings = std::move(masked); mask_ice_half_width_q = half_width_q; return *this;
+    }
 
     // Fit the global error model from the spread of symmetry-equivalent observations.
     // Call once before merging; AddImage then applies it.

process/JFJochProcess.cpp

@@ -35,6 +35,7 @@
 #include "../image_analysis/scale_merge/HKLKey.h"
 #include "../image_analysis/WriteReflections.h"
 #include "../common/Definitions.h"
+#include "../common/CorrelationCoefficient.h"
 #include <array>
 #include <map>
 #include <Eigen/Dense>
@@ -566,6 +567,9 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) {
         // Smoothing it more than once would compound the correction, so do it only on the first
         // pass. The no-reference path recomputes G from scratch each pass and re-smooths correctly.
         bool reference_g_smoothed = false;
+        // Ice rings masked from the merge by the CC1/2 test after the final merge; empty until then,
+        // at which point a re-merge applies them.
+        std::vector<char> masked_ice_rings;
         auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult {
             // ScaleOnTheFly self-scaling is only for the no-reference path; with a reference each
             // image is already scaled against it during the per-image pass, so we merge directly.
@@ -620,6 +624,7 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) {
             // intensities and the error model, so the space-group search sees clean data; the final
             // in-symmetry merge keeps them so completeness is not lost.
             merge_engine.ExcludeIceRings(for_search);
+            merge_engine.MaskIceRings(masked_ice_rings, config_.spot_finding.ice_ring_width_Q_recipA);
             if (result.consensus_cell.has_value())
                 merge_engine.ReferenceCell(*result.consensus_cell);
             merge_engine.RefineErrorModel(merge_input);
@@ -670,6 +675,39 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) {
             }
         }
 
+        // Ice-ring CC1/2 mask: a hexagonal-ice ring whose merged half-set CC1/2 collapses well below its
+        // resolution shoulders has been decorrelated by ice (its Bragg is unrecoverable) - drop it and
+        // 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; })) {
+                masked_ice_rings = std::move(mask);
+                const auto final_sg = experiment_.GetGemmiSpaceGroup();
+                sm = scale_and_merge(final_sg ? final_sg->short_name() : "P1", false);
+            }
+        }
+
         const auto twin_sg_number = experiment_.GetSpaceGroupNumber();
         const gemmi::SpaceGroup *twin_sg = twin_sg_number
             ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr;