From 0e5a15cb3c73b650333754e49d891d92dfdab261 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 7 Nov 2026 11:17:21 +0100 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_01HVeDQs733XDnbQQyUbV1vs --- .../BraggIntegrationEngineCPU.cpp | 8 +++-- .../BraggIntegrationEngineGPU.cu | 5 +-- tests/BraggIntegrationEngineCPUTest.cpp | 33 +++++++++++++++++++ 3 files changed, 42 insertions(+), 4 deletions(-) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 684d2d608..e10a3a6aa 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/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(s.I_sum) - s.n_inner_valid * out.bkg; - if (net > 0.0) { + if (net > out.sigma) { out.obs_x = (static_cast(s.I_sum_x) - out.bkg * s.x_sum) / net; out.obs_y = (static_cast(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(s.I_sum_x) / static_cast(s.I_sum); out.obs_y = static_cast(s.I_sum_y) / static_cast(s.I_sum); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index ec33b726f..41f8d3cef 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/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 { diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 845f73072..81ea24326 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/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{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.