Defective pixels: mask a persistent patch, not only a lone pixel; sub-lattice ask reads the R contrast

Two defects that together turned a tetragonal small-molecule crystal (open-arm cuhf2, published
P4/nmm) into P222.

1. HotPixelFinder masked a persistent strong pixel only when it stood alone or in a pair ("a
   larger patch is a feature of the scattering"). On cuhf2 a ring of ~20 pixels reads 1e5-3e5
   counts on EVERY frame (dead centre) - stationary in the lab, so no reflection of the rotating
   crystal. Unmasked, the reflections crossing it ((2,6,+-9)) merged to 3.5e6 / 6.7e5 against 4e4
   for the strongest real reflection. The final merge's outlier test removes them (12
   equivalents), but the space-group search's P1-like merge has no equivalents to judge by, and
   the h<->k operator read CC 0.18 / R 0.28 instead of ~0.99 / 0.013. The persistence and strength
   tests are unchanged; only the isolation requirement is dropped.

2. The sub-lattice ask kept every metric two-fold that passes min_operator_cc. On a pseudo-cubic
   cell two false cubic two-folds still correlate at 0.31-0.33, LePage closes them with the
   genuine ones into the full cubic metric, and the tetragonal class in between is never asked.
   An operator now also has to pass the R contrast every promotion's added operators are held to
   (min_operator_r_contrast against random_pairing_r / global_best_operator_r, read only where the
   search reads it): the false ones sit at 0.2, the genuine ones at 1.0.

cuhf2: P222 -> P422, R_meas 2.9% -> 2.6%, ISa 45.5 -> 52.3 (overnight-cint had found P422 by the
sub-lattice path; rc174-all lost it). Prescan survey, masked pixels base -> new: unchanged on
lyso_x06da_ref/5keV/atten_wedge, thau, insu, cytc, myob_split, 9qw8, aspirin20, HEPES, YAG,
metformin, nidppe; cuhf2 3 -> 42, 5reo 2 -> 6, lcystine25 0 -> 9 (two clusters of noisy pixels at
3-7 counts/frame on a zero background). Targeted battery (sm2rest-fix vs rc174all-ctl / smt-all):
the 15 protein / sub-lattice sets (lyso x3, thau, insu, cytc, myob_split, 5reo, 9qw8, 6iu6, 6iu8,
6iu9, 6z8o, 8t7r, 9ea5) identical to the 4th digit; SM sets identical except cuhf2 (above) and
lcystine25, SHELXL R1 .112 -> .151: the nine masked pixels flip the fulls-only per-frame scale
smoother on that polycrystalline sweep from "settled after 29 iterations" to "stopped settling"
(the base run reproduces bit-identically) - a sensitivity of that smoother, reported to the
scaling work, not of the mask.

Test: HotPixelFinder_PersistentPixelNotBragg gains a 3x3 patch that must be masked. The
device-vs-host test now waits for each queued frame before overwriting the shared device buffer;
it remains flaky at ~4/25 on the base binary as well (a borderline pixel or two differ), which is
pre-existing and not addressed here.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB

rugnux/HotPixels.cpp

@@ -335,16 +335,13 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp
             }
     });
 
-    // Masked: a persistent pixel standing alone - on its own or as one of a pair; a larger patch is
-    // a feature of the scattering, not a defect - whose mean is STRONG_RATIO times its ring's and
-    // whose excess is above the Poisson bound. A weaker one cannot make an outlier.
-    const auto persistent_neighbours = [&](size_t x, size_t y) {
-        int count = 0;
-        for (size_t yy = (y > 0 ? y - 1 : 0); yy <= std::min(height - 1, y + 1); yy++)
-            for (size_t xx = (x > 0 ? x - 1 : 0); xx <= std::min(width - 1, x + 1); xx++)
-                if ((xx != x || yy != y) && persistent[yy * width + xx]) count++;
-        return count;
-    };
+    // Masked: a persistent pixel whose mean is STRONG_RATIO times its ring's and whose excess is above
+    // the Poisson bound - a weaker one cannot make an outlier. Alone or in a patch: a patch lit on frame
+    // after frame is stationary in the lab and so no reflection of a rotating crystal, whatever made it,
+    // and every reflection that crosses it is integrated on top of it. Measured: a ring of about twenty
+    // pixels at 1e5-3e5 counts on every frame (a dead centre in the middle) merged the reflections
+    // crossing it to 90 times the strongest real one; the space-group search, which merges without the
+    // equivalents an outlier test needs, read the crystal's four-fold off them as unrelated.
     ParallelChunks(static_cast<int>(height), nthreads, [&](int y0, int y1) {
         for (size_t y = y0; y < static_cast<size_t>(y1); y++)
             for (size_t x = 0; x < width; x++) {
@@ -355,17 +352,9 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp
                     continue;
                 }
                 if (!persistent[i]) continue;
-                const int neighbours = persistent_neighbours(x, y);
-                bool isolated = neighbours == 0;
-                if (neighbours == 1)
-                    for (size_t yy = (y > 0 ? y - 1 : 0); yy <= std::min(height - 1, y + 1); yy++)
-                        for (size_t xx = (x > 0 ? x - 1 : 0); xx <= std::min(width - 1, x + 1); xx++)
-                            if ((xx != x || yy != y) && persistent[yy * width + xx])
-                                isolated = persistent_neighbours(xx, yy) == 1;
                 const double mean_level = static_cast<double>(sum_level(i)) / n_valid(i);
                 const double excess = static_cast<double>(sum_value[i] - sum_level(i)) / n_valid(i);
-                if (isolated
-                    && static_cast<double>(sum_value[i]) >= STRONG_RATIO * static_cast<double>(sum_level(i))
+                if (static_cast<double>(sum_value[i]) >= STRONG_RATIO * static_cast<double>(sum_level(i))
                     && excess >= LIT_NSIGMA * std::sqrt(std::max(mean_level, 0.0)) + LIT_OFFSET)
                     ret.mask[i] = 1;
             }

rugnux/HotPixels.h

@@ -27,9 +27,10 @@
 //     frames (measured, on the ring's pixels that are not candidates), and the number of lit frames
 //     that fewer than 0.01 of all the detector's pixels would reach by chance is the binomial bound.
 //
-// A persistent pixel is masked only if it stands alone (by itself or in a pair - a larger patch is a
-// feature of the scattering, not a defect), reads on average at least ten times its ring, and its
-// mean excess is above the Poisson bound: a weaker one cannot make an outlier. A pixel holding the
+// A persistent pixel is masked if it reads on average at least ten times its ring and its mean excess
+// is above the Poisson bound: a weaker one cannot make an outlier. Alone or as part of a patch - a
+// patch lit through the whole sweep is stationary in the lab, so it is not a reflection of the
+// rotating crystal, and a reflection crossing it reads the patch. A pixel holding the
 // detector's error value on most frames is masked with it: every frame already treats it as invalid,
 // but the static mask did not know. The detector's outermost row and column get no rule of their
 // own - a hot pixel there is caught like any other.

rugnux/Rugnux.cpp

@@ -7886,12 +7886,30 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
                         const gemmi::Mat33 from_primitive = to_primitive.inverse();
                         // An operator counts as real on the same bound Stage A calls one present at -
                         // SearchSpaceGroupOptions::min_operator_cc, applied by OperatorCorrelation over
-                        // the same reflection population, so nothing new is calibrated here.
+                        // the same reflection population - AND on the R contrast every promotion's added
+                        // operators are held to (SearchSpaceGroupOptions::min_operator_r_contrast,
+                        // against the same two ends the search above measured), so nothing new is
+                        // calibrated here. The correlation alone cannot tell a pseudo-symmetric
+                        // structure's near-operators from real ones: measured on a pseudo-cubic
+                        // tetragonal crystal, two of the false cubic two-folds correlate at 0.31-0.33 and
+                        // pass it, LePage closes them with the genuine ones into the whole cubic metric,
+                        // and the tetragonal class between the two is never asked about. Their R reads
+                        // 0.39 against 0.013 for the genuine operators and 0.49 for unrelated
+                        // reflections - a contrast of 0.2 against the genuine operators' 1.0.
+                        const double r_random = sg_search.random_pairing_r;
+                        const double r_best = sg_search.global_best_operator_r;
+                        // Read only where the search itself reads it: a best operator close to the
+                        // random end leaves no scale to measure on (max_reference_over_random).
+                        const bool contrast_known = std::isfinite(r_random) && std::isfinite(r_best)
+                                                    && r_best <= sg_opts.max_reference_over_random * r_random;
                         const auto keep_operator = [&](const gemmi::Mat33 &m_primitive) {
                             const gemmi::Mat33 m_indexed =
                                     from_primitive.multiply(m_primitive).multiply(to_primitive);
                             const auto score = OperatorCorrelation(sm.merged, m_indexed, sg_opts);
-                            return score.has_value() && score->present;
+                            if (!score.has_value() || !score->present)
+                                return false;
+                            return !contrast_known
+                                   || (r_random - score->r_stat) / (r_random - r_best) >= sg_opts.min_operator_r_contrast;
                         };
                         const auto sub = LePageLattice(primitive, LATTICE_MAX_OBLIQUITY_DEG, keep_operator);
                         if (sub && static_cast<size_t>(sub->n_operators) > sg_search.point_group_order

tests/HotPixelFinderTest.cpp

@@ -25,7 +25,8 @@ namespace {
 // high on every frame, one only 20 high - persistent, but too weak to make an outlier - one that holds
 // the error value on every frame, and one crossed by a genuine reflection. The last sits close to the rotation axis, where one reflection stays on a pixel for a
 // long stretch of rotation - here nine consecutive sampled frames, more than the chance bound allows
-// but fewer than the one-reflection bound at that zeta. Only the first and the third are masked.
+// but fewer than the one-reflection bound at that zeta. Only the first and the third are masked - and
+// a 3x3 patch reading high on every frame, which is stationary in the lab and so no reflection either.
 TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     DiffractionExperiment x(DetDECTRIS(W, H, "Test detector", ""));
     x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(10.0f);
@@ -36,6 +37,7 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
 
     constexpr int HOT_X = 200, HOT_Y = 200, WARM_X = 190, WARM_Y = 60, ERR_X = 60, ERR_Y = 190;
     constexpr int BRAGG_X = 40, BRAGG_Y = 115;
+    constexpr int PATCH_X = 150, PATCH_Y = 225;
     std::mt19937 rng(1);
     // The powder ring is a smooth radial profile, as a real one is: 45 counts over the background at
     // 82 px, 4 px sigma.
@@ -50,6 +52,9 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
         frame[I(HOT_X, HOT_Y)] += 60;
         frame[I(WARM_X, WARM_Y)] += 20;
         frame[I(ERR_X, ERR_Y)] = INT32_MIN;
+        for (int dy = -1; dy <= 1; dy++)
+            for (int dx = -1; dx <= 1; dx++)
+                frame[I(PATCH_X + dx, PATCH_Y + dy)] += 60;
         if (f >= 20 && f < 29)
             frame[I(BRAGG_X, BRAGG_Y)] += 500;
         finder.AddImage(frame.data(), scratch);
@@ -61,7 +66,10 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     CHECK(result.mask[I(ERR_X, ERR_Y)] == 2);
     CHECK(result.mask[I(WARM_X, WARM_Y)] == 0);
     CHECK(result.mask[I(BRAGG_X, BRAGG_Y)] == 0);
-    CHECK(result.hot == 1);
+    for (int dy = -1; dy <= 1; dy++)
+        for (int dx = -1; dx <= 1; dx++)
+            CHECK(result.mask[I(PATCH_X + dx, PATCH_Y + dy)] == 1);
+    CHECK(result.hot == 10);
     CHECK(result.error == 1);
 }
 
@@ -122,6 +130,8 @@ TEST_CASE("HotPixelFinder_DeviceMatchesHost", "[HotPixelFinder]") {
         REQUIRE(cudaMemcpy(device_image, image.data(), image.size() * sizeof(int32_t), cudaMemcpyHostToDevice)
                 == cudaSuccess);
         device.AddDeviceImage(device_image, frame);
+        // The frame is queued, not processed: the next copy must not overwrite it before its kernels ran.
+        REQUIRE(cudaStreamSynchronize(*frame.stream) == cudaSuccess);
     }
 
     const auto expected = host.GetMask(OSC_DEG, SPACING_DEG, 4);