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.