Correction surfaces: the resolution gauge counts observations

The detector-frame surfaces are held to a mean of one within each resolution
shell, because a factor that is a function of d alone is a null direction of
their fit. That mean was weighted by 1/sigma^2, and sigma is smallest where
little light arrives: on a detector part-shadowed by the pin the dark cells
carry the shell's mean, and their large factors go onto the lit cells as small
ones. How much of a ring is shadowed varies with resolution, so this bends the
merged intensities along resolution: REFMAC's fitted scale against the model
ran 1.48 / 0.80 / 1.16 in F (1.65 A / 2.7 A / inner) on a half-shadowed lysozyme
sweep, with R as good as XDS's within every shell and 0.02 worse over all, and
no shell-wise statistic can see it. Each observation now counts once.

REFMAC R-work / R-free on the depositor's flags (d_min of the baseline), model
as deposited, 31 sets chosen by the largest removed ramp and clamped cells,
largest R vs per-shell R gap, plus controls: 5 better by more than 0.0015, 23
the same, 3 worse.
  half-shadowed lysozyme   0.3266 / 0.3306 -> 0.2996 / 0.3069 (XDS on the same
                           model 0.3036 / 0.3083)
  soft X-ray insulin       0.2720 / 0.2841 -> 0.2423 / 0.2570
  4p8r                     0.2143 / 0.2518 -> 0.1873 / 0.2216 (depositor's data
                           0.2216)
  6i3j -0.0029; insulin reference -0.0017; 2wnz +0.0025, 6jgh +0.0016, 5j23 +0.0015.
The ramp removed on the first falls from 186x to 2x. The merge's own hash test
records change as they must.

A mean over the ring still lowers its lit part by the shadow it averages in
(the new test reads 0.72 where 0.57 before and 1 would be exact); a median of the
counts did not do better on the half-shadowed sweep (R-free 0.3090), the only one
it was measured on.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQs733XDnbQQyUbV1vs

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -4519,9 +4519,13 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
     // that each shell still holds tens of thousands of observations.
     constexpr int NSHELL = 20;
     const int nshell = n_terms >= 100 * NSHELL ? NSHELL : 0;
-    // Weight of cell c in shell s: the statistical weight its observations carry there. Computed ONCE
-    // from the incoming sigmas rather than per round, so the projection is a fixed linear operator and
-    // the alternating fit below converges against it.
+    // Weight of cell c in shell s: the number of observations it holds there. Computed ONCE rather than
+    // per round, so the projection is a fixed linear operator and the alternating fit below converges
+    // against it. Not their statistical weight 1/sigma^2: sigma is smallest where little light arrives,
+    // so a cell behind a shadow would carry the shell's mean, its large factor going onto the lit cells
+    // as a small one. How much of a ring is shadowed varies with resolution, so that is a wave in the
+    // merged intensities, one no shell-wise statistic can see. Measured on a sweep with a third of the
+    // detector shadowed: 186x of ramp removed, 3x in <I> between shells against the deposited model.
     std::vector<double> shw(static_cast<size_t>(ncell) * nshell, 0.0);
     std::vector<double> shw_shell(nshell, 0.0), shw_cell(nshell > 0 ? ncell : 0, 0.0);
     std::vector<int32_t> group_shell(n_groups, 0);   // the shell each ASU group sits in, for the gate
@@ -4552,7 +4556,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
             group_shell[t.group] = s;
             const double sc = static_cast<double>(t.sigma) * t.corr;
             if (!(sc > 0.0)) continue;
-            shw[static_cast<size_t>(t.cell) * nshell + s] += 1.0 / (sc * sc);
+            shw[static_cast<size_t>(t.cell) * nshell + s] += 1.0;
         }
     }
     for (int c = 0; c < ncell && nshell > 0; ++c)

image_analysis/scale_merge/RotationScaleMergeGPU.cu

@@ -687,7 +687,7 @@ namespace {
             const auto &t = term[order[q]];
             const double sc = __dmul_rn(double(t.sigma), double(t.corr));
             if (!(sc > 0.0)) continue;
-            w = __dadd_rn(w, 1.0 / __dmul_rn(sc, sc));
+            w = __dadd_rn(w, 1.0);
         }
         shw[b] = w;
     }

image_analysis/scale_merge/RotationScaleMergeGPU.h

@@ -308,9 +308,9 @@ public:
                           int block);
 
     // The resolution gauge's inputs for the terms SurfaceBuildTerms left: each group's shell (that of its
-    // last term in fulls order; length n_groups) and each (cell, shell) bin's statistical weight
-    // sum 1/(sigma corr)^2 in fulls order (length ncell * nshell), over nshell equal-occupancy shells in
-    // 1/d^2 - ApplyCellSurface's numbers to the bit.
+    // last term in fulls order; length n_groups) and each (cell, shell) bin's observation count
+    // (length ncell * nshell), over nshell equal-occupancy shells in 1/d^2 - ApplyCellSurface's numbers
+    // to the bit.
     void SurfaceShellWeights(int n_terms, int nshell, int32_t *group_shell, double *shw);
 
     // One subset of the terms (0 = even frames, 1 = odd, 2 = all), cut into the host's n_blocks

tests/RotationScaleMergeSegmentsTest.cpp

@@ -202,8 +202,8 @@ TEST_CASE("RotationScaleMergeSegments_SegmentsMustCoverTheFrames", "[RotationSca
 
 TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[RotationScale]") {
     // The same synthetic sweep merged without segments gives, to the bit, what the engine gave before
-    // segments existed: FNV-1a over the merged I and sigma, recorded on the engine of rc181 e56b19abb -
-    // on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host
+    // segments existed: FNV-1a over the merged I and sigma, recorded on the engine with the surface gauge
+    // counting observations - on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host
     // path to the bit only where every reduction does; this one carries its own, so it is checked against
     // its own record).
     const std::vector<SweepSpec> sweeps{{40, 1.0, 0.0}, {40, 0.7, 3.0}};
@@ -226,9 +226,54 @@ TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[Rotation
     CHECK(r.segments.empty());
     CHECK(r.merged.size() == 4106);
     if (!on_gpu)
-        CHECK(hash == 3376995616949734726ULL);
+        CHECK(hash == 15900277682318835699ULL);
     else if (get_gpu_description().find("RTX 5080") != std::string::npos)
-        CHECK(hash == 6594155396455850264ULL);
+        CHECK(hash == 6640579017825061097ULL);
+}
+
+TEST_CASE("RotationScaleMerge_ShadowDoesNotBendTheResolutionShape", "[RotationScale]") {
+    // Every reflection lands on a ring of radius proportional to 1/d, at a random azimuth. Past a fixed
+    // detector column a shadow passes a fifth of the light and takes the background with it, so those
+    // observations have the small sigmas of a dark region - and it reaches only the rings wide enough to
+    // cross that column. The surface takes the shadow out; what it must not do is leave the lit part of
+    // the shadowed rings dimmer than the rings the shadow never touches.
+    const std::vector<SweepSpec> sweeps{{400, 1.0, 0.0}};
+    auto data = MakeSweeps(sweeps, 21);
+    std::mt19937 rng(22);
+    std::uniform_real_distribution<float> azimuth(0.0f, 2.0f * static_cast<float>(M_PI));
+    std::normal_distribution<double> noise(0.0, 1.0);
+    for (auto &frame: data.frames)
+        for (auto &r: frame.reflections) {
+            const float radius = 300.0f / r.d;
+            const float phi = azimuth(rng);
+            r.predicted_x = 500.0f + radius * std::cos(phi);
+            r.predicted_y = 500.0f + radius * std::sin(phi);
+            if (r.predicted_x > 590.0f) {
+                const double expected = data.I_true.at(FriedelKey(r.h, r.k, r.l)) * 0.2;
+                const double sigma = std::sqrt(expected + 1.0);
+                r.I = static_cast<float>(expected + sigma * noise(rng));
+                r.sigma = static_cast<float>(sigma);
+                r.var_bkg = static_cast<float>(sigma * sigma - std::max(0.0, static_cast<double>(r.I)));
+                r.bkg = 1.0f;
+            }
+        }
+    const auto result = Merge(data, SyntheticExperiment(true), {});
+
+    // <I_merged> / <I_true> in the rings the shadow misses (d > 3.4 A) against those it crosses
+    double in_merged = 0.0, in_true = 0.0, out_merged = 0.0, out_true = 0.0;
+    for (const auto &m: result.merged) {
+        const auto t = data.I_true.find(FriedelKey(m.h, m.k, m.l));
+        REQUIRE(t != data.I_true.end());
+        (m.d > 3.4 ? in_merged : out_merged) += m.I;
+        (m.d > 3.4 ? in_true : out_true) += t->second;
+    }
+    REQUIRE(in_true > 0.0);
+    REQUIRE(out_true > 0.0);
+    const double shape = (out_merged / out_true) / (in_merged / in_true);
+    // A mean over the ring still lowers its lit part by the shadow it averages in: 0.72 with every
+    // observation counted once, 0.57 where the dark ones' small sigmas gave them the weight.
+    CHECK(shape > 0.65);
+    CHECK(shape < 1.25);
 }
 
 // A sweep in which no rocking event is captured whole enough to make a full - as in one exposed 0.2 deg out