Merge branch 'rc186-ice-int' into rc186-merge

docs/ACKNOWLEDGEMENT.md

@@ -472,6 +472,14 @@ hands at I/sigma 3, reported beside SigAno and CCanom — is phenix.xtriage's: P
 indicators in protein crystallography" (2005), Acta Cryst. D61, 1437-1448
 [doi:10.1107/S0907444905023589](https://doi.org/10.1107/S0907444905023589).
 
+**Ice in the merged intensities** — the merged reflections on the hexagonal-ice rings are compared with
+the off-ring reflections beside them in resolution, and the rings' excesses are fitted together in
+proportion to hexagonal ice's line strengths; the table of those strengths, and the idea of reading ice
+off a merge's own intensities at the ring positions, are phenix.xtriage's: P. H. Zwart, R. W.
+Grosse-Kunstleve and P. D. Adams, "Xtriage and Fest: automatic assessment of X-ray data and substructure
+structure factor estimation" (2005), CCP4 Newsletter on Protein Crystallography 43 (no DOI). rugnux's
+comparison is local in resolution, where xtriage's is a z-score of a shell against the whole Wilson curve.
+
 ### Physical corrections and calibration
 
 **Sensor absorption at oblique incidence, and the flight path** — the angle-dependent quantum

docs/CPU_DATA_ANALYSIS_INTEGRATION.md

@@ -487,6 +487,8 @@ A twin law is only partly an ambiguity where the lattice is pseudo-merohedral: a
 
 ### 10.10 Ice rings at the scale and merge stages
 
+Whatever the gate of §3.3 found, every written merge is also read for ice in its own intensities: for each hexagonal-ice ring the merged reflections in the ring's core ($|q - q_{\mathrm{ring}}| < 0.015$ Å$^{-1}$) and in its flanks (out to $0.03$ Å$^{-1}$) are compared with a Wilson line fitted to the off-ring reflections within $0.15$ Å$^{-1}$ of the ring only, as the ratio of their mean intensity to the line and that ratio's excess in standard errors of the band's own scatter. A span so short in resolution cancels anisotropy and the curvature of the Wilson plot, which is what a z-score of a shell against a whole-dataset Wilson curve cannot do, and on the battery that z-score (phenix.xtriage's ice test) flagged broad features of the Wilson plot on some thirty crystals with flat radial profiles, while this one did not. The rings are then fitted together as one ring system - an additive intensity in proportion to hexagonal ice's line strengths - twice: the cores against the line, which is the size of the bias, and the cores against their own flanks, which decides. A ring is narrower than the flank offset, so it leaves a contrast between the two whichever way the integration bent it (a ring under the spots lifts the cores; a ring in the background annulus beside them depresses the cores and lifts the flanks - both are seen on the battery), while a feature of the Wilson plot broader than a ring moves core and flanks alike and cancels. The merge is read as carrying ice where the contrast stands at 5 sigma, either sign, over at least two rings (`ICE_INTENSITY_Z`). Where it fires on a run whose images showed no ring, the report warns and names the shells; the rings' reflections are not set aside after the fact.
+
 Where the gate of §3.3 has found ice, reflections falling within $\pm w$ in $q = 2\pi/d$ (§1.2) of a hexagonal-ice band ($w=0.03$ Å$^{-1}$ offline, about the measured ring half-width) are marked. Marked reflections are **excluded where a model is fitted** — the per-frame scale $G$, the per-image correlation, and the $P1$ merge the space-group search runs on — because ice contamination is a *positive bias*, not extra scatter, and a least-squares scale absorbs it into $G$ and into the error-model $b$, where it damages every other reflection on the same frame. They are also left out of the **resolution-cut fit** (§13.4), which is the one consumer that reads the *merged* reflections rather than the observations: an ice-flagged observation marks its merged reflection, on the rotation merge's own accumulator (host and device alike) as well as on the stills one. The stills merge keeps them in the final merge. The **rotation merge decides ring by ring**, because what a ring does to its reflections differs between rings and between runs: it can bias them (up on the ring, where its powder intensity is under the spot; down on the flanks, where it is in the background ring), only add scatter, or leave them alone. Deleting every band unconditionally was implemented once and removed: on a crystal whose band was harmless it removed better-than-average reflections and 6 to 8 points of completeness in the affected shell for no measurable gain. So each ring is split into its **core** ($|q-q_\mathrm{ring}|<0.01$ Å$^{-1}$) and its two **flanks** (out to $w$), and every band is compared with the off-ring merged reflections within $0.02$ Å$^{-2}$ of the ring in $1/d^2$. Both are read against a log-linear **Wilson fit** of the mean intensity over those off-ring reflections (the fit behind the reported Wilson $B$; over all off-ring reflections where the window cannot be fitted), as $t=(I-I_\mathrm{exp})/\sqrt{\sigma^2+I_\mathrm{exp}^2}$, which has unit scatter whether a reflection is signal- or noise-dominated. A band whose mean $t$ departs from its neighbours' by more than 3 standard errors is **rejected** — a bias, which no larger $\sigma$ would remove; one whose variance of $t$ exceeds the neighbours' by more than 3 standard errors of the log ratio (taken from the samples' own kurtosis: Wilson intensities are far from Gaussian) is **down-weighted**, its excess variance added to each reflection's own; the rest is **merged** as measured. A band of fewer than 20 reflections, or a ring with fewer than 50 neighbours, is merged. The decision is applied to the merge's accumulators, so the written reflections, the shell table and the half-set $\mathrm{CC}$ all see the same thing, and each band's verdict is logged and reported (`ICE_RING_BAND_nn`, see the results report). It complements the radial background correction (§9), which removes most of the bias at its source: on 58 rotation crystals with ice and a model, it lowered the model's R-free against merging every band on 26 and raised it on one, by 0.003.
 
 ---

docs/RUGNUX_REPORT.md

@@ -740,13 +740,42 @@ against the ice-free flanks beside them (ice in large crystallites); 1 is no ice
 **`ICE_RINGS_DETECTED=`** is `TRUE` where either passes its gate (1.5 and 2.0), and
 **`ICE_REFLECTIONS_ON_RINGS_PCT=`** is then the share of the integrated reflections set aside from
 scaling. **`ICE_RING_BAND_nn=`** (rotation only) is one line per band of a ring - its core, within
-0.01 Å$^{-1}$ in $q$ of the ring, and its low- and high-$q$ flanks out to the flag width: the ring's
-$d$ (Å), the band, what the merge did with its reflections (`MERGE`, `DOWN_WEIGHT`, `REJECT`), how
+0.01 Å$^{-1}$ in $q$ of the ring, and on each side an inner flank out to 0.02 and an outer one out to the
+flag width (a ring's shoulder reaches into the inner flank, so it is judged apart from the outer): the ring's
+$d$ (Å), the band (`CORE`, `LOW_FLANK`, `LOW_OUTER_FLANK`, ...), what the merge did with its reflections (`MERGE`, `DOWN_WEIGHT`, `REJECT`), how
 many merged reflections it holds, their bias against the off-ring reflections around the ring in
 standard errors, and the ratio of their variance to those reflections' (both on intensities
 normalised by a local Wilson fit). A band is rejected where it is biased beyond 3 standard errors,
-down-weighted where it is only noisier beyond 3, and merged otherwise. The summary has an `Ice` line, and where the merge found ice an `ICE_RINGS` warning prompts a
-look at the ice-ring shells (one warning line, whichever of the two measurements raised it). The measurement is described in
+down-weighted where it is only noisier beyond 3, and merged otherwise.
+
+Ice is also looked for in the merged intensities themselves, on every written merge and whether or
+not the images showed a ring - a ring too weak or too diffuse for the radial profile still leaves the
+reflections on it too strong:
+
+```
+ICE_INTENSITY_Z= +8.4
+ICE_INTENSITY_BIAS_Z= +8.2
+ICE_INTENSITY_RINGS= 6
+ICE_INTENSITY_RING_01= 3.895 182 1.03 +1.2 1.02
+ICE_INTENSITY_RING_05= 2.249 551 1.31 +5.8 0.90
+```
+
+**`ICE_INTENSITY_RING_nn=`** is one line per hexagonal-ice ring the merge reaches: the ring's $d$ (Å),
+the merged reflections in its core (within 0.015 Å$^{-1}$ in $q$), their mean intensity over a Wilson
+line fitted to the off-ring reflections within 0.15 Å$^{-1}$ of the ring, that excess in standard
+errors of the core's own scatter, and the same mean over the ring's flanks (0.015-0.03 Å$^{-1}$). The
+comparison is local in resolution, so anisotropy and the curvature of the Wilson plot cancel, and a small
+cell with few reflections per band is not read as ice. **`ICE_INTENSITY_Z=`** fits the cores against
+their flanks over all rings as one hexagonal-ice ring system, in proportion to the known line strengths,
+in sigma: a ring is narrower than the flank offset, while a feature of the Wilson plot broader than one
+moves core and flanks alike and cancels, so a strong shell that happens to sit on one ring position does
+not pass as ice. Its sign says which way the integration bent the ring - positive where the ring under
+the spots lifted the cores, negative where the ring in the background annulus beside them depressed
+the cores and lifted the flanks. **`ICE_INTENSITY_BIAS_Z=`** fits the cores against the line instead,
+the size of the bias. The summary has an `Ice` line carrying both readings, and an `ICE_RINGS` warning
+prompts a look at the ice-ring shells where the merge found ice on the images, or where the intensities
+show the ring system at 5 sigma or more (either sign) over at least two rings although the images did
+not (one warning line, whichever measurement raised it). The measurement is described in
 [CPU/GPU data analysis ▸ Resolution and ice-ring handling](CPU_DATA_ANALYSIS_IMAGE.md#33-resolution-and-ice-ring-handling).
 
 ## Overloads and overlapping spots

image_analysis/scale_merge/CMakeLists.txt

@@ -38,6 +38,8 @@ ADD_LIBRARY(JFJochScaleMerge
         ScalingResult.h
         ScalingResult.cpp
         WilsonOutliers.cpp
+        IceRingIntensity.cpp
+        IceRingIntensity.h
         WilsonOutliers.h
         ErrorModel.cpp
         ErrorModel.h

image_analysis/scale_merge/IceRingIntensity.cpp

@@ -0,0 +1,130 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "IceRingIntensity.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "../../common/Definitions.h"
+#include "../../common/JFJochMath.h"
+
+namespace {
+    // Hexagonal ice's line strengths at the rings a merge commonly reaches; the rest of ICE_RING_RES_A
+    // are weak. Following phenix.xtriage's ice-ring table (Zwart, Grosse-Kunstleve & Adams, CCP4
+    // Newsletter 43, 2005).
+    constexpr std::array<float, 11> ICE_REL = {1.000f, 0.750f, 0.530f, 0.170f, 0.390f, 0.300f,
+                                               0.040f, 0.180f, 0.030f, 0.020f, 0.020f};
+
+    struct Point { double s2, I; };
+}
+
+bool WilsonLineThrough(const std::vector<MergedReflection> &reflections, double &intercept, double &slope) {
+    std::vector<Point> base;
+    base.reserve(reflections.size());
+    for (const auto &r : reflections)
+        if (std::isfinite(r.I) && r.d > 0.0f)
+            base.push_back({1.0 / (static_cast<double>(r.d) * r.d), r.I});
+    std::sort(base.begin(), base.end(), [](const Point &p, const Point &q) { return p.s2 < q.s2; });
+    constexpr int NSHELL = 8;
+    std::vector<double> x, y, w;
+    for (int i = 0; i < NSHELL; ++i) {
+        const size_t lo = base.size() * i / NSHELL, hi = base.size() * (i + 1) / NSHELL;
+        if (hi - lo < 5) continue;
+        double sx = 0, sI = 0, sI2 = 0;
+        for (size_t j = lo; j < hi; ++j) { sx += base[j].s2; sI += base[j].I; sI2 += base[j].I * base[j].I; }
+        const double n = static_cast<double>(hi - lo), m = sI / n, var = std::max(sI2 / n - m * m, 1e-30);
+        if (m <= 0.0) continue;
+        x.push_back(sx / n);
+        y.push_back(std::log(m));
+        w.push_back(n * m * m / var);     // 1 / var(ln mean)
+    }
+    if (x.size() < 3) return false;
+    // Weighted least squares on the shell points.
+    double sw = 0, swx = 0, swy = 0, swxx = 0, swxy = 0;
+    for (size_t i = 0; i < x.size(); ++i) {
+        sw += w[i]; swx += w[i] * x[i]; swy += w[i] * y[i]; swxx += w[i] * x[i] * x[i]; swxy += w[i] * x[i] * y[i];
+    }
+    const double det = sw * swxx - swx * swx;
+    if (!(std::fabs(det) > 0.0)) return false;
+    slope = (sw * swxy - swx * swy) / det;
+    intercept = (swy - slope * swx) / sw;
+    return true;
+}
+
+IceRingIntensity MeasureIceRingIntensity(const std::vector<MergedReflection> &merged) {
+    IceRingIntensity out;
+    std::vector<double> q;
+    q.reserve(merged.size());
+    for (const auto &r : merged)
+        q.push_back(std::isfinite(r.I) && r.d > 0.0f ? 2.0 * PI / r.d : NAN);
+    std::vector<char> on_band(merged.size(), 0);
+    for (size_t i = 0; i < merged.size(); ++i)
+        if (std::isfinite(q[i]))
+            for (const float d : ICE_RING_RES_A)
+                if (std::fabs(q[i] - 2.0 * PI / d) < ICE_INTENSITY_BAND_Q) { on_band[i] = 1; break; }
+
+    // Across-ring fits: y_k against x_k = rel_k / <I>_k, weighted by 1 / se_k^2, y_k the core's excess
+    // over the line and then the core's excess over the flanks.
+    double sxx = 0.0, sxy = 0.0, cxx = 0.0, cxy = 0.0;
+    for (size_t k = 0; k < ICE_REL.size(); ++k) {
+        const double qk = 2.0 * PI / ICE_RING_RES_A[k];
+        std::vector<MergedReflection> base;
+        int below = 0, above = 0;
+        for (size_t i = 0; i < merged.size(); ++i) {
+            if (!std::isfinite(q[i]) || on_band[i] || std::fabs(q[i] - qk) >= ICE_INTENSITY_WINDOW_Q) continue;
+            base.push_back(merged[i]);
+            (q[i] < qk ? below : above)++;
+        }
+        if (static_cast<int>(base.size()) < ICE_INTENSITY_MIN_BASE || below < 20 || above < 20) continue;
+        double a, b;
+        if (!WilsonLineThrough(base, a, b)) continue;
+
+        // The core and the flanks, each reflection as its ratio to the line.
+        struct Band { double sum = 0.0, sum2 = 0.0, sum_exp = 0.0; int n = 0; } core, flank;
+        for (size_t i = 0; i < merged.size(); ++i) {
+            if (!std::isfinite(q[i]) || std::fabs(q[i] - qk) >= ICE_INTENSITY_BAND_Q) continue;
+            const double I_exp = std::exp(a + b / (static_cast<double>(merged[i].d) * merged[i].d));
+            const double ratio = merged[i].I / I_exp;
+            Band &band = std::fabs(q[i] - qk) < ICE_INTENSITY_CORE_Q ? core : flank;
+            band.sum += ratio; band.sum2 += ratio * ratio; band.sum_exp += I_exp; ++band.n;
+        }
+        if (core.n < ICE_INTENSITY_MIN_CORE || flank.n < ICE_INTENSITY_MIN_CORE) continue;
+        auto mean_se = [](const Band &band, double &mean, double &se) {
+            mean = band.sum / band.n;
+            const double var = (band.sum2 - band.sum * band.sum / band.n) / (band.n - 1);
+            se = std::sqrt(std::max(var, 0.0) / band.n);
+        };
+        double mean, se, mean_f, se_f;
+        mean_se(core, mean, se);
+        mean_se(flank, mean_f, se_f);
+        if (!(se > 0.0) || !(se_f > 0.0)) continue;
+        IceRingIntensityRing ring;
+        ring.ring_d = ICE_RING_RES_A[k];
+        ring.rel = ICE_REL[k];
+        ring.n_core = core.n;
+        ring.n_flank = flank.n;
+        ring.n_base = static_cast<int>(base.size());
+        ring.ratio = mean;
+        ring.z = (mean - 1.0) / se;
+        ring.ratio_flank = mean_f;
+        ring.z_flank = (mean_f - 1.0) / se_f;
+        ring.mean_expected = core.sum_exp / core.n;
+        out.rings.push_back(ring);
+        if (ring.mean_expected > 0.0) {
+            const double x = ring.rel / ring.mean_expected;
+            sxx += x * x / (se * se);
+            sxy += x * (mean - 1.0) / (se * se);
+            cxx += x * x / (se * se + se_f * se_f);
+            cxy += x * (mean - mean_f) / (se * se + se_f * se_f);
+            ++out.n_rings;
+        }
+    }
+    if (sxx > 0.0) {
+        out.amplitude = sxy / sxx;
+        out.amplitude_z = out.amplitude * std::sqrt(sxx);
+        out.contrast = cxy / cxx;
+        out.contrast_z = out.contrast * std::sqrt(cxx);
+    }
+    return out;
+}

image_analysis/scale_merge/IceRingIntensity.h

@@ -0,0 +1,67 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cmath>
+#include <vector>
+
+#include "../../common/Reflection.h"
+
+// Ice seen in the merged intensities, whatever the images said.
+//
+// A powder ring that the image-side gate did not catch - too weak for the radial profile, or diffuse -
+// still leaves its mark on the integrated reflections: those on the ring come out too strong, since
+// the ring under a spot is not what the background annulus beside it measures. The mark is narrow in
+// resolution, so each ring's core (|q - q_ring| < ICE_INTENSITY_CORE_Q) and its flanks (out to
+// ICE_INTENSITY_BAND_Q) are compared with a Wilson line fitted to the off-ring reflections within
+// ICE_INTENSITY_WINDOW_Q of it only: anisotropy and the curvature of the Wilson plot cancel over so
+// short a span, which is what a z-score against a whole-dataset Wilson curve (xtriage's ice test)
+// cannot do, and the standard error is the band's own scatter, so a small cell with few reflections per
+// band is not read as ice. The rings are then fitted together as one hypothesis - an additive ring
+// intensity in proportion to the known ice line strengths - twice: the cores against the line, which
+// is the size of the bias, and the cores against their own flanks, which is what decides. A feature of
+// the Wilson plot broader than a ring moves core and flank alike and cancels in the second; a ring,
+// narrower than the flank offset, does not, whichever way the integration bent it (a ring under the
+// spot lifts the core; a ring in the background annulus beside it can as well depress the core and lift
+// the flanks), so the decision is on the contrast's size, not its sign.
+struct IceRingIntensityRing {
+    float ring_d = 0.0f;       // the ice ring (A)
+    float rel = 0.0f;          // its relative strength in hexagonal ice
+    int n_core = 0;            // merged reflections in the ring's core
+    int n_flank = 0;           // and in its two flanks together
+    int n_base = 0;            // off-ring reflections the Wilson line was fitted to
+    double ratio = 1.0;        // mean I / I_expected over the core
+    double z = 0.0;            // (ratio - 1) / its standard error
+    double ratio_flank = 1.0;  // the same over the flanks
+    double z_flank = 0.0;
+    double mean_expected = 0;  // the fitted intensity at the core, in the merge's units
+};
+
+struct IceRingIntensity {
+    std::vector<IceRingIntensityRing> rings;   // the rings the merge reaches with enough reflections
+    // The ice hypothesis across the rings, ratio_k - 1 = amplitude * rel_k / <I>_k: the cores against
+    // the line (the bias) and the cores against their flanks (the contrast, which decides).
+    double amplitude = 0.0;
+    double amplitude_z = 0.0;
+    double contrast = 0.0;
+    double contrast_z = 0.0;
+    int n_rings = 0;                            // rings that entered the fits
+    [[nodiscard]] bool Detected() const { return n_rings >= 2 && std::fabs(contrast_z) >= ICE_INTENSITY_Z; }
+    static constexpr double ICE_INTENSITY_Z = 5.0;
+};
+
+constexpr double ICE_INTENSITY_CORE_Q = 0.015;    // the core, |q - q_ring| in A^-1
+constexpr double ICE_INTENSITY_BAND_Q = 0.03;     // reflections this close to any ring are not baseline
+constexpr double ICE_INTENSITY_WINDOW_Q = 0.15;   // the baseline is fitted within this of the ring
+constexpr int ICE_INTENSITY_MIN_CORE = 20;
+constexpr int ICE_INTENSITY_MIN_BASE = 60;
+
+// `merged` needs I and d.
+IceRingIntensity MeasureIceRingIntensity(const std::vector<MergedReflection> &merged);
+
+// ln <I> against 1/d^2 through eight equal-count shells of `reflections`, each weighted by its own
+// precision, with no I/sigma cut on the shells (a weak dataset has none that would pass one, and its
+// rings still need a line to be read against). False where fewer than three shells have a positive
+// mean.
+bool WilsonLineThrough(const std::vector<MergedReflection> &reflections, double &intercept, double &slope);

image_analysis/scale_merge/IceRingVerdict.cpp

@@ -7,7 +7,7 @@
 
 #include "../../common/Definitions.h"
 #include "../../common/JFJochMath.h"
-#include "../bragg_integration/CalcISigma.h"
+#include "IceRingIntensity.h"
 
 namespace {
     struct Moments {
@@ -29,6 +29,16 @@ namespace {
     double OneOverD2(float d) { return 1.0 / (static_cast<double>(d) * d); }
 }
 
+const char *IceRingZoneName(int zone) {
+    switch (zone) {
+        case -2: return "low-q outer flank";
+        case -1: return "low-q flank";
+        case 1: return "high-q flank";
+        case 2: return "high-q outer flank";
+        default: return "core";
+    }
+}
+
 const char *IceRingActionName(IceRingAction a) {
     switch (a) {
         case IceRingAction::DownWeight: return "down-weight";
@@ -45,7 +55,6 @@ IceRingDecision DecideIceRings(const std::vector<MergedReflection> &merged) {
     for (const auto &r : merged)
         if (!r.on_ice_ring)
             off_ring.push_back(r);
-    const GlobalWilsonB global = CalcGlobalWilsonB(off_ring);
 
     // The ring each flagged reflection belongs to: the nearest in q.
     std::vector<int> ring_of(merged.size(), -1);
@@ -72,29 +81,55 @@ IceRingDecision DecideIceRings(const std::vector<MergedReflection> &merged) {
         for (const auto &r : off_ring)
             if (r.d > 0.0f && std::fabs(OneOverD2(r.d) - s2_ring) < ICE_VERDICT_WINDOW_S2)
                 window.push_back(r);
-        const GlobalWilsonB local = CalcGlobalWilsonB(window);
-        const GlobalWilsonB &wilson = local.n_shells >= 3 ? local : global;
+        // The line the bands are read against: through the window's own shells, with no I/sigma cut
+        // (a weak low-resolution dataset has no shell that would pass one, and its rings - the ones
+        // that reach it - still need judging); flat at the window's mean where even that cannot be
+        // fitted.
+        double intercept = 0.0, slope = 0.0;
+        bool have_line = WilsonLineThrough(window, intercept, slope);
+        if (!have_line && !window.empty()) {
+            double sum = 0.0;
+            for (const auto &r : window) sum += r.I;
+            if (sum > 0.0) {
+                intercept = std::log(sum / static_cast<double>(window.size()));
+                have_line = true;
+            }
+        }
 
         auto t_of = [&](const MergedReflection &r) {
-            const double I_exp = std::exp(wilson.intercept + wilson.slope * OneOverD2(r.d));
+            const double I_exp = std::exp(intercept + slope * OneOverD2(r.d));
             return (r.I - I_exp) / std::sqrt(static_cast<double>(r.sigma) * r.sigma + I_exp * I_exp);
         };
 
         Moments nb;
-        if (wilson.n_shells >= 3)
+        if (have_line)
             for (const auto &r : window)
                 nb.Add(t_of(r));
 
-        Moments zone[3];
-        for (size_t i = 0; i < merged.size(); ++i)
-            if (ring_of[i] == k && nb.n > 0)
-                zone[dq_of[i] < -ICE_VERDICT_CORE_Q ? 0 : (dq_of[i] > ICE_VERDICT_CORE_Q ? 2 : 1)].Add(t_of(merged[i]));
+        // Each band's reflections are counted whether or not there is a baseline to judge them
+        // against, so a band merged for want of one is still reported as such.
+        // Five bands: the core and, on each side, an inner and an outer flank. A ring's shoulder reaches
+        // past the core into the inner flank; judged together with the outer flank it was diluted into a
+        // merge (a 5x shoulder was kept on one crystal), so each flank is judged on its own.
+        auto zone_of = [](double dq) {
+            return dq < -ICE_VERDICT_FLANK_Q ? 0 : dq < -ICE_VERDICT_CORE_Q ? 1
+                 : dq <= ICE_VERDICT_CORE_Q ? 2 : dq <= ICE_VERDICT_FLANK_Q ? 3 : 4;
+        };
+        Moments zone[5];
+        int in_zone[5] = {0, 0, 0, 0, 0};
+        for (size_t i = 0; i < merged.size(); ++i) {
+            if (ring_of[i] != k) continue;
+            const int z = zone_of(dq_of[i]);
+            ++in_zone[z];
+            if (nb.n > 0)
+                zone[z].Add(t_of(merged[i]));
+        }
 
-        IceRingZoneVerdict v[3];
-        for (int z = 0; z < 3; ++z) {
+        IceRingZoneVerdict v[5];
+        for (int z = 0; z < 5; ++z) {
             v[z].ring_d = ICE_RING_RES_A[k];
-            v[z].zone = z - 1;
-            v[z].n = zone[z].n;
+            v[z].zone = z - 2;
+            v[z].n = in_zone[z];
             v[z].n_neighbours = nb.n;
             if (zone[z].n < ICE_VERDICT_MIN_ZONE || nb.n < ICE_VERDICT_MIN_NEIGHBOURS)
                 continue;
@@ -110,15 +145,15 @@ IceRingDecision DecideIceRings(const std::vector<MergedReflection> &merged) {
 
         for (size_t i = 0; i < merged.size(); ++i) {
             if (ring_of[i] != k) continue;
-            const IceRingZoneVerdict &zv = v[dq_of[i] < -ICE_VERDICT_CORE_Q ? 0 : (dq_of[i] > ICE_VERDICT_CORE_Q ? 2 : 1)];
+            const IceRingZoneVerdict &zv = v[zone_of(dq_of[i])];
             if (zv.action == IceRingAction::Reject)
                 out.weight_factor[i] = 0.0;
             else if (zv.action == IceRingAction::DownWeight) {
                 // The band's excess variance of t, in the reflection's own units.
                 const auto &r = merged[i];
-                const double I_exp = std::exp(wilson.intercept + wilson.slope * OneOverD2(r.d));
+                const double I_exp = std::exp(intercept + slope * OneOverD2(r.d));
                 const double var = static_cast<double>(r.sigma) * r.sigma;
-                const double excess = (zone[zv.zone + 1].Var() - nb.Var()) * (var + I_exp * I_exp);
+                const double excess = (zone[zv.zone + 2].Var() - nb.Var()) * (var + I_exp * I_exp);
                 out.weight_factor[i] = var / (var + excess);
             }
         }

image_analysis/scale_merge/IceRingVerdict.h

@@ -14,10 +14,11 @@
 // noisier than the crystal's own reflections - or untouched, where the ring is weak. Neither CC1/2 nor
 // a scatter test sees a bias: the ring is the same in both half-sets. So each band is compared with the
 // off-ring reflections around it in resolution, through the Wilson curve: the expected intensity comes
-// from a log-linear Wilson fit (CalcGlobalWilsonB) over the off-ring reflections within
-// ICE_VERDICT_WINDOW_S2 of the ring, or over all of them where that window cannot be fitted, and each
+// from a log-linear Wilson fit (WilsonLineThrough, no I/sigma cut on its shells) over the off-ring
+// reflections within ICE_VERDICT_WINDOW_S2 of the ring, flat at their mean where that cannot be fitted, and each
 // reflection is read as t = (I - I_exp) / sqrt(sigma^2 + I_exp^2) - unit scatter whether the reflection
-// is signal- or noise-dominated. Per band (the core of the ring and each flank):
+// is signal- or noise-dominated. Per band (the core of the ring, and an inner and an outer flank on each
+// side - a ring's shoulder reaches into the inner flank, and judged with the outer one it was diluted):
 //   * the band's mean t differs from its neighbours' by more than ICE_VERDICT_SIGMA standard errors:
 //     biased - REJECT, since a larger sigma would not remove a bias;
 //   * otherwise its variance of t exceeds the neighbours' by more than ICE_VERDICT_SIGMA standard errors
@@ -29,7 +30,7 @@ enum class IceRingAction { Merge, DownWeight, Reject };
 
 struct IceRingZoneVerdict {
     float ring_d = 0.0f;        // the ice ring (A)
-    int zone = 0;               // -1 low-q flank, 0 core, +1 high-q flank
+    int zone = 0;               // -2 low-q outer flank, -1 low-q flank, 0 core, +1 high-q flank, +2 high-q outer flank
     int n = 0;                  // merged reflections in the band
     int n_neighbours = 0;       // off-ring reflections they are compared with
     double bias_z = 0.0;        // (mean t band - mean t neighbours) / its standard error
@@ -47,6 +48,7 @@ struct IceRingDecision {
 
 constexpr double ICE_VERDICT_SIGMA = 3.0;
 constexpr double ICE_VERDICT_CORE_Q = 0.01;        // |q - q_ring| below this is the core (A^-1)
+constexpr double ICE_VERDICT_FLANK_Q = 0.02;       // ... below this the inner flank, beyond it the outer
 constexpr double ICE_VERDICT_WINDOW_S2 = 0.02;     // neighbours within this of the ring in 1/d^2 (A^-2)
 constexpr int ICE_VERDICT_MIN_ZONE = 20;
 constexpr int ICE_VERDICT_MIN_NEIGHBOURS = 50;
@@ -55,3 +57,4 @@ constexpr int ICE_VERDICT_MIN_NEIGHBOURS = 50;
 IceRingDecision DecideIceRings(const std::vector<MergedReflection> &merged);
 
 const char *IceRingActionName(IceRingAction a);
+const char *IceRingZoneName(int zone);

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -6272,7 +6272,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
                 for (const auto &z : ice.zones)
                     logger->Info("Ice ring {:.3f} A {}: {} reflections against {} off the ring, bias {:+.1f} sigma, "
                                  "variance x{:.2f} ({:+.1f} sigma) - {}", z.ring_d,
-                                 z.zone < 0 ? "low-q flank" : (z.zone > 0 ? "high-q flank" : "core"), z.n,
+                                 IceRingZoneName(z.zone), z.n,
                                  z.n_neighbours, z.bias_z, z.var_ratio, z.var_ratio_z, IceRingActionName(z.action));
                 result.ice_ring_verdicts = ice.zones;
             }

rugnux/ResultReport.cpp

@@ -2276,7 +2276,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                 const auto &z = result.ice_ring_verdicts[i];
                 Add(s, KeyText(fmt::format("ICE_RING_BAND_{:02}", i + 1).c_str(),
                                fmt::format("{:.3f} {} {} {} {:+.1f} {:.2f}", z.ring_d,
-                                           z.zone < 0 ? "LOW_FLANK" : (z.zone > 0 ? "HIGH_FLANK" : "CORE"),
+                                           z.zone == -2 ? "LOW_OUTER_FLANK" : z.zone == -1 ? "LOW_FLANK" : z.zone == 0 ? "CORE"
+                                           : z.zone == 1 ? "HIGH_FLANK" : "HIGH_OUTER_FLANK",
                                            z.action == IceRingAction::Reject ? "REJECT"
                                            : (z.action == IceRingAction::DownWeight ? "DOWN_WEIGHT" : "MERGE"),
                                            z.n, z.bias_z, z.var_ratio)));
@@ -2305,6 +2306,52 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                          " flanks beside them (ice in large crystallites); 1 is no ice on either.", true));
             Add(s, Blank());
         }
+        // ---- ice as the merged intensities show it, whatever the images said
+        const auto &ii = result.ice_intensity;
+        if (!ii.rings.empty()) {
+            Add(s, KeyReal("ICE_INTENSITY_Z", ii.contrast_z, "{:+.1f}"));
+            Add(s, KeyReal("ICE_INTENSITY_BIAS_Z", ii.amplitude_z, "{:+.1f}"));
+            Add(s, KeyInt("ICE_INTENSITY_RINGS", ii.n_rings));
+            for (size_t i = 0; i < ii.rings.size(); ++i) {
+                const auto &r = ii.rings[i];
+                Add(s, KeyText(fmt::format("ICE_INTENSITY_RING_{:02}", i + 1).c_str(),
+                               fmt::format("{:.3f} {} {:.2f} {:+.1f} {:.2f}", r.ring_d, r.n_core, r.ratio, r.z, r.ratio_flank)));
+            }
+            Add(s, Blank());
+            // The ring whose core stands furthest from its flanks, in standard errors.
+            const auto strongest = std::max_element(ii.rings.begin(), ii.rings.end(),
+                                                    [](const IceRingIntensityRing &a, const IceRingIntensityRing &b) {
+                                                        return std::fabs(a.z - a.z_flank) < std::fabs(b.z - b.z_flank); });
+            const bool above = strongest->ratio > strongest->ratio_flank;
+            Add(s, Prose(ii.Detected()
+                    ? fmt::format(" Ice in the intensities: DETECTED. The merged reflections on the hexagonal-ice rings\n"
+                                  " stand {} the reflections beside them in resolution as a ring would, {:.1f} sigma\n"
+                                  " over {} rings (the strongest, {:.2f} A: core {:.2f}, flanks {:.2f} of the Wilson line).{}",
+                                  above ? "above" : "below", std::fabs(ii.contrast_z), ii.n_rings, strongest->ring_d,
+                                  strongest->ratio, strongest->ratio_flank,
+                                  result.ice_rings_detected
+                                      ? " The ring handling above has treated what it could."
+                                      : " The images' radial profile and spots showed no ring, so nothing was set\n"
+                                        " aside for it: check the ice-ring shells.")
+                    : fmt::format(" Ice in the intensities: no ring system. The ring cores stand {:+.1f} sigma from the\n"
+                                  " reflections beside them over {} rings, inside the 5 sigma a ring system would show.",
+                                  ii.contrast_z, ii.n_rings)));
+            Add(s, Prose(" ICE_INTENSITY_RING_nn: ring d (A), merged reflections in its core, their mean intensity\n"
+                         " over a Wilson line fitted to the off-ring reflections within 0.15 A^-1 of the ring, that\n"
+                         " excess in standard errors, and the same mean over the ring's flanks. ICE_INTENSITY_Z fits\n"
+                         " core against flanks over all rings as one hexagonal-ice ring system (a ring is narrower\n"
+                         " than the flank offset; a broad feature of the Wilson plot moves both alike), in sigma;\n"
+                         " ICE_INTENSITY_BIAS_Z fits the cores against the line, the size of the bias.", true));
+            const bool ice_warned = std::any_of(doc.warnings.begin(), doc.warnings.end(),
+                                                [](const ReportWarning &w) { return w.code == std::string(PathologyCode::ICE_RINGS); });
+            if (ii.Detected() && !ice_warned)
+                Warn(doc, PathologyCode::ICE_RINGS, fmt::format(
+                        "Ice rings: the merged intensities on the hexagonal-ice rings stand {} their neighbours as a ring "
+                        "would ({:.1f} sigma; at {:.2f} A the core is {:.2f} and the flanks {:.2f} of the Wilson line) "
+                        "although the images showed no ring - check the ice-ring shells", above ? "above" : "below",
+                        std::fabs(ii.contrast_z), strongest->ring_d, strongest->ratio, strongest->ratio_flank));
+            Add(s, Blank());
+        }
 
         // ---- overloads and overlapping spots: two ways the detector or the pattern took reflections
         //      the crystal gave, each with what to change at the next collection
@@ -3643,14 +3690,21 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                                       pw.rings_q_recipA.size() - pw.ice_ring_count, 100.0f * pw.non_ice_spot_fraction,
                                       result.powder_excluded_from_indexing ? "; left out of indexing" : ""));
             }
-            if (result.ice_ring_score || result.ice_spot_ratio)
-                checked_row("Ice", {PathologyCode::ICE_RINGS}, result.ice_rings_detected
+            if (result.ice_ring_score || result.ice_spot_ratio || !result.ice_intensity.rings.empty()) {
+                const auto &ii = result.ice_intensity;
+                const std::string in_intensities = ii.rings.empty() ? std::string()
+                        : ii.Detected() ? fmt::format("; the merged intensities on the rings stand {:.1f} sigma {} their "
+                                                      "neighbours", std::fabs(ii.contrast_z),
+                                                      ii.contrast_z > 0 ? "above" : "below")
+                                        : fmt::format("; no ring system in the merged intensities (core-flank contrast "
+                                                      "{:+.1f} sigma)", ii.contrast_z);
+                checked_row("Ice", {PathologyCode::ICE_RINGS}, (result.ice_rings_detected
                         ? fmt::format("rings DETECTED: the {:.1f}% of the reflections on them were left out of "
                                       "scaling{}", 100.0 * result.ice_flagged_fraction,
                                       result.ice_ring_verdicts.empty() ? ", kept in the merge"
                                                                        : ", merged ring by ring")
-                        : std::string("no rings detected"));
-            else
+                        : std::string("no rings detected on the images")) + in_intensities);
+            } else
                 untested("Ice", {PathologyCode::ICE_RINGS});
             if (result.harmonic.detected)
                 checked_row("Beam", {PathologyCode::HARMONIC_CONTAMINATION},

rugnux/Rugnux.h

@@ -35,6 +35,7 @@
 #include "../image_analysis/scale_merge/StillsPartialityRefine.h"   // StillsPartialityRefine::RunWidth
 #include "../image_analysis/scale_merge/TranslationalNCS.h"   // TranslationalNCSResult
 #include "../image_analysis/scale_merge/IceRingVerdict.h"   // IceRingZoneVerdict
+#include "../image_analysis/scale_merge/IceRingIntensity.h"
 #include "../image_analysis/indexing/HarmonicContamination.h"  // HarmonicContaminationResult
 #include "../image_analysis/geom_refinement/CellEsd.h"   // CellEsd
 #include "../image_analysis/structure_refinement/ModelValidation.h"                  // ModelValidationResult
@@ -847,6 +848,10 @@ struct ProcessResult {
     // What the rotation merge did with each ice-ring band (DecideIceRings): merged, down-weighted or
     // rejected. Empty on stills and where nothing was flagged.
     std::vector<IceRingZoneVerdict> ice_ring_verdicts;
+    // Ice as the merged intensities show it (MeasureIceRingIntensity), measured on every written
+    // merge whether or not the gate above found any: the ring cores against the off-ring reflections
+    // beside them, and the ice hypothesis fitted across the rings.
+    IceRingIntensity ice_intensity;
     // The reindex operators a rotation run's cell and space group admit (" ; "-joined), where no
     // reference or model fixed the choice. Within one sweep the choice is arbitrary but consistent;
     // it matters only against other data. Empty where there is none, or where it was resolved.

rugnux/RugnuxScaleMerge.cpp

@@ -1958,8 +1958,16 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
         result.segment_decay_model = r.segment_decay_model;
         if (for_search)
             result.search_scaling_converged = r.scaling_converged;
-        else
+        else {
             result.ice_ring_verdicts = r.ice_ring_verdicts;
+            result.ice_intensity = MeasureIceRingIntensity(r.merged);
+            if (result.ice_intensity.Detected())
+                logger.Info("Ice in the merged intensities ({}): the ring cores stand {} the reflections beside "
+                            "them as a ring would, {:.1f} sigma over {} rings{}", label,
+                            result.ice_intensity.contrast_z > 0 ? "above" : "below",
+                            std::fabs(result.ice_intensity.contrast_z), result.ice_intensity.n_rings,
+                            result.ice_rings_detected ? "" : " - the images' profile and spots showed none");
+        }
         Note(fmt::format("merge in {}: {} unique", label, r.merged.size()));
         return ScaleMergeResult{std::move(r.merged), std::move(r.statistics),
                                 r.cc_half_before_corrections,

tests/CMakeLists.txt

@@ -213,6 +213,7 @@ ADD_EXECUTABLE(jfjoch_test
         SphericalHarmonicSurfaceTest.cpp
         WilsonOutliersTest.cpp
         IceRingVerdictTest.cpp
+        IceRingIntensityTest.cpp
         PairSortGPUTest.cpp
         ReindexAmbiguityTest.cpp
         CrystalSettingTest.cpp

tests/IceRingIntensityTest.cpp

@@ -0,0 +1,81 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+#include <random>
+
+#include "../common/Definitions.h"
+#include "../common/JFJochMath.h"
+#include "../image_analysis/scale_merge/IceRingIntensity.h"
+
+namespace {
+    // A Wilson crystal (B = 20 A^2) to 1.8 A, acentric intensities with 10% plus constant noise, with
+    // an additive powder ring of strength `ring(d)` (in units of the local mean) on the core of
+    // every hexagonal-ice ring - or, with `whole_band`, over core and flanks alike.
+    template <class F>
+    std::vector<MergedReflection> Crystal(F ring, bool whole_band = false) {
+        std::mt19937 rng(4321);
+        std::exponential_distribution<double> wilson(1.0);
+        std::normal_distribution<double> noise(0.0, 1.0);
+        std::uniform_real_distribution<double> s2(1.0 / (8.0 * 8.0), 1.0 / (1.8 * 1.8));
+        std::vector<MergedReflection> out;
+        for (int i = 0; i < 40000; ++i) {
+            MergedReflection r;
+            r.d = static_cast<float>(1.0 / std::sqrt(s2(rng)));
+            const double I_mean = 1000.0 * std::exp(-0.5 * 20.0 / (r.d * r.d));
+            double I_true = I_mean * wilson(rng);
+            const double q = 2.0 * PI / r.d;
+            for (const float ice_d : ICE_RING_RES_A)
+                if (std::fabs(q - 2.0 * PI / ice_d) < (whole_band ? ICE_INTENSITY_BAND_Q : ICE_INTENSITY_CORE_Q))
+                    I_true += ring(ice_d) * I_mean;
+            r.sigma = static_cast<float>(std::sqrt(0.01 * I_true * I_true + 25.0));
+            r.I = static_cast<float>(I_true + r.sigma * noise(rng));
+            out.push_back(r);
+        }
+        return out;
+    }
+}
+
+TEST_CASE("IceRingIntensity_clean_crystal_is_not_ice") {
+    const auto m = MeasureIceRingIntensity(Crystal([](float) { return 0.0; }));
+    REQUIRE(m.n_rings >= 5);
+    CHECK(!m.Detected());
+    CHECK(std::fabs(m.amplitude_z) < 3.0);
+    CHECK(std::fabs(m.contrast_z) < 3.0);
+    for (const auto &r : m.rings)
+        CHECK(std::fabs(r.z) < 4.0);
+}
+
+TEST_CASE("IceRingIntensity_hexagonal_ring_system_is_detected") {
+    // Ice in proportion to its line strengths: the 3.90 A ring adds 100% of the local mean.
+    const auto m = MeasureIceRingIntensity(Crystal([](float d) {
+        return d > 3.8f ? 1.0 : d > 3.6f ? 0.75 : d > 3.4f ? 0.53 : d > 2.2f && d < 2.3f ? 0.39 : 0.0;
+    }));
+    REQUIRE(m.Detected());
+    CHECK(m.amplitude_z > 8.0);
+    CHECK(m.contrast_z > 5.0);
+    for (const auto &r : m.rings)
+        if (r.ring_d > 3.8f)
+            CHECK(r.ratio == Catch::Approx(2.0).margin(0.2));
+}
+
+TEST_CASE("IceRingIntensity_one_strong_shell_is_not_ice") {
+    // A crystal whose transform is strong over a shell at one ring position, wider than the ring: the
+    // cores stand out against the Wilson line but not against their flanks.
+    const auto m = MeasureIceRingIntensity(Crystal([](float d) { return d > 2.6f && d < 2.7f ? 1.0 : 0.0; }, true));
+    for (const auto &r : m.rings)
+        if (r.ring_d > 2.6f && r.ring_d < 2.7f) {
+            CHECK(r.z > 5.0);
+            CHECK(r.ratio_flank == Catch::Approx(r.ratio).margin(0.2));
+        }
+    CHECK(!m.Detected());
+}
+
+TEST_CASE("IceRingIntensity_depressed_cores_are_a_ring_too") {
+    // A ring in the background annulus beside the spots takes intensity out of the cores instead.
+    const auto m = MeasureIceRingIntensity(Crystal([](float d) {
+        return d > 3.8f ? -0.5 : d > 3.6f ? -0.4 : d > 3.4f ? -0.3 : d > 2.2f && d < 2.3f ? -0.25 : 0.0;
+    }));
+    REQUIRE(m.Detected());
+    CHECK(m.contrast_z < -5.0);
+}

tests/IceRingVerdictTest.cpp

@@ -13,19 +13,19 @@ namespace {
     // reflections within the flag width of an ice ring are flagged, and `ring` modifies the ones on
     // the 2.667 A ring: dq is their offset from it in q.
     template <class F>
-    std::vector<MergedReflection> Crystal(F ring) {
+    std::vector<MergedReflection> Crystal(F ring, double d_min = 2.0, double noise = 1.0) {
         std::mt19937 rng(1234);
         std::exponential_distribution<double> wilson(1.0);
-        std::normal_distribution<double> noise(0.0, 1.0);
-        std::uniform_real_distribution<double> s2(1.0 / (6.0 * 6.0), 1.0 / (2.0 * 2.0));
+        std::normal_distribution<double> gauss(0.0, 1.0);
+        std::uniform_real_distribution<double> s2(1.0 / (6.0 * 6.0), 1.0 / (d_min * d_min));
         std::vector<MergedReflection> out;
         for (int i = 0; i < 20000; ++i) {
             MergedReflection r;
             r.d = static_cast<float>(1.0 / std::sqrt(s2(rng)));
             const double I_mean = 1000.0 * std::exp(-0.5 * 20.0 / (r.d * r.d));
             const double I_true = I_mean * wilson(rng);
-            r.sigma = static_cast<float>(std::sqrt(0.01 * I_true * I_true + 25.0));
-            r.I = static_cast<float>(I_true + r.sigma * noise(rng));
+            r.sigma = static_cast<float>(noise * std::sqrt(0.01 * I_true * I_true + 25.0));
+            r.I = static_cast<float>(I_true + r.sigma * gauss(rng));
             r.on_ice_ring = IsOnIceRing(r.d, 0.03f);
             const double dq = 2.0 * PI / r.d - 2.0 * PI / 2.667;
             if (std::fabs(dq) < 0.03)
@@ -60,8 +60,8 @@ TEST_CASE("IceRingVerdict_biased_core_is_rejected") {
     });
     const auto d = DecideIceRings(merged);
     CHECK(ActionOf(d, 2.667f, 0) == IceRingAction::Reject);
-    CHECK(ActionOf(d, 2.667f, -1) == IceRingAction::Reject);
-    CHECK(ActionOf(d, 2.667f, 1) == IceRingAction::Reject);
+    for (int zone : {-2, -1, 1, 2})
+        CHECK(ActionOf(d, 2.667f, zone) == IceRingAction::Reject);
     CHECK(ActionOf(d, 2.249f, 0) == IceRingAction::Merge);
     size_t rejected = 0;
     for (double f : d.weight_factor)
@@ -83,3 +83,17 @@ TEST_CASE("IceRingVerdict_noisy_unbiased_ring_is_down_weighted") {
             CHECK(d.weight_factor[i] > 0.0);
         }
 }
+
+TEST_CASE("IceRingVerdict_rings_near_the_resolution_limit_get_a_verdict") {
+    // A weak 3.1 A dataset, every shell at I/sigma below 1: the 3.44 A ring sits just inside the limit
+    // and the 3.66 A one close to it. Each ring with flagged reflections must be judged, not merged
+    // silently for want of a Wilson fit.
+    const auto merged = Crystal([](MergedReflection &, double, double, std::mt19937 &) {}, 3.1, 4.0);
+    const auto d = DecideIceRings(merged);
+    for (float ring_d : {3.895f, 3.661f, 3.438f}) {
+        bool judged = false;
+        for (const auto &z : d.zones)
+            judged |= z.ring_d == ring_d && z.zone == 0 && z.n_neighbours >= ICE_VERDICT_MIN_NEIGHBOURS;
+        CHECK(judged);
+    }
+}