Integration: no signal-centroid from a net below sigma

Pass A's background-free centroid divided the disk's first moment by
net = I_sum - n*bkg whenever net > 0. An integer sum against a rational
background reaches zero to rounding (net ~ 1e-14), and any small net gives a
position off the disk: the first moment's noise is radius x sigma, so the
centroid's is radius x sigma/net. Below sigma the raw centroid is taken, as it
already was for net <= 0. Both engines; Metal shares the CPU code.

XDET/YDET up to 1e16 in the unmerged MTZ was the visible part: 67 of 68
rotation sets scanned had rows off the detector. Their I is right (profile fit
and box sum are centred on the prediction; z against symmetry mates has rms
0.92 vs 1.04 for the control).

The real harm was in post-refinement. Such a centroid can reach the strongest
20000 positions, where the held-out cost is a plain mean with a noise from the
fourth moment: one outlier inflates the noise and the cross-validation refuses
a fit that the paired change supports at -17 sd. Measured, --no-gpu, this commit
against its parent, 14 sets: 5 change their verdicts and so their geometry;
the 9 others have merged MTZs bitwise identical. On the half-image in-house
lysozyme: joint noise 6.8e-8 -> 2.4e-8, reject -> COMMIT, cell a 78.44 -> 78.35
(XDS 77.74), R_meas 56.0 -> 46.8 %, ISa 4.3 -> 4.7. Open-arm 8c3e: 1.74 -> 1.67 A.
Small-molecule L-alanine single sweep commits too, ISa 57.3 -> 55.6.

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

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -461,12 +461,16 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection
             // beam and the distance on these centroids, so it under-corrects by that factor and
             // does so resolution-dependently. Subtracting a flat pedestal from a first moment is
             // exact: sum(x*(px-bkg)) = I_sum_x - bkg*x_sum over the same pixels.
+            //
+            // The first moment's noise is the disk radius times sigma, so the centroid's is the radius
+            // times sigma/net: a net below sigma places it off the disk, and a net that is zero to
+            // rounding, which an integer sum against a rational background reaches, at 1e15 pixels.
             const double net = static_cast<double>(s.I_sum) - s.n_inner_valid * out.bkg;
-            if (net > 0.0) {
+            if (net > out.sigma) {
                 out.obs_x = (static_cast<double>(s.I_sum_x) - out.bkg * s.x_sum) / net;
                 out.obs_y = (static_cast<double>(s.I_sum_y) - out.bkg * s.y_sum) / net;
             } else {
-                // Nothing above background to take a centroid of; the raw one is all there is.
+                // No signal to take a centroid of; the raw one is all there is.
                 out.obs_x = static_cast<double>(s.I_sum_x) / static_cast<double>(s.I_sum);
                 out.obs_y = static_cast<double>(s.I_sum_y) / static_cast<double>(s.I_sum);
             }

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -448,9 +448,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
         // The centroid of the SIGNAL. Weighting by raw counts weights by signal plus background,
         // and the background is flat over a disk centred on the PREDICTION, so its centroid is the
         // prediction and it pulls the answer there - the displacement comes out shrunk by
-        // I/(I + n*bkg), worst where the background dominates. See the CPU engine.
+        // I/(I + n*bkg), worst where the background dominates. A net below sigma has no centroid to
+        // give. See the CPU engine.
         const double net = (double) Isum - (double) s_ninner_valid * bkg;
-        if (net > 0.0) {
+        if (net > sigma) {
             ox = ((double) (long long) s_Ix - bkg * (double) (long long) s_x) / net;
             oy = ((double) (long long) s_Iy - bkg * (double) (long long) s_y) / net;
         } else {

tests/BraggIntegrationEngineCPUTest.cpp

@@ -155,6 +155,39 @@ TEST_CASE("BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped", "[Integr
     CHECK(saturated[0].clipped);
 }
 
+// A reflection with next to nothing above the background: a first moment over a net of a few counts
+// is no position, and the centroid of the disk has to stay on the disk.
+TEST_CASE("BraggIntegrationEngineCPU_NoSignalCentroidStaysOnTheDisk", "[Integration][portable]") {
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f)
+            .Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
+    experiment.ImportIndexingSettings(IndexingSettings().RotationIndexing(true));
+    experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings());
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+
+    const int cx = 600, cy = 300;
+    ImagePreprocessorBuffer image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        image[i] = 150;
+    // A noise dipole with a net of one count.
+    image[cy * width + cx + 2] += 100;
+    image[cy * width + cx - 2] -= 99;
+
+    Reflection r{};
+    r.h = 1; r.k = 2; r.l = 3;
+    r.predicted_x = cx; r.predicted_y = cy;
+    r.d = 2.0f;
+    r.prescaling_corr = 1.0f;
+    r.partiality = 1.0f;
+
+    BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment));
+    const auto out = engine.Run(image, std::vector<Reflection>{r}, 1, 0);
+    REQUIRE(out.size() == 1);
+    REQUIRE(std::isfinite(out[0].observed_x));
+    CHECK(std::fabs(out[0].observed_x - cx) < 2.0f);
+    CHECK(std::fabs(out[0].observed_y - cy) < 2.0f);
+}
+
 // An empty aperture on a background that is a curved function of the radius alone: the flat ring mean
 // reads it wrong by the curvature, and the radial correction with the curve taken from a per-ring
 // profile of the whole frame (as the spot finder hands it over on stills) takes that back out.