diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 22bc416c7..205621eb8 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -630,6 +630,11 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & } rh.bkg -= corr; // annulus mean -> mean over the signal disk rh.I = static_cast(rh.I_sum) - static_cast(rh.n_inner) * rh.bkg; + // The background under the disk is now the corrected one, and so is its Poisson share of + // the non-signal variance. The estimate's own variance (bkg_var) stays: it is the noise of + // the ring pixels, which the correction - a smooth curve - does not change; nor does sigma, + // which is taken from the counts in the disk. + rh.var_bkg -= pixel_gain * static_cast(rh.n_inner) * corr; } } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 412562901..78a85d108 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -615,8 +615,8 @@ __global__ void build_profiles(const unsigned long long *shell_grid, const unsig __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve, int n_rad, const float *k_diff, int k_len, int k_off, const float *isum_a, const int *ninner_a, const int *rbin_a, - const int *kbin_a, - const uint8_t *ok_a, float *bkg_o, float *I_o, int n) { + const int *kbin_a, float pixel_gain, + const uint8_t *ok_a, float *bkg_o, float *I_o, float *varbkg_o, int n) { const int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= n || !ok_a[i]) return; const float bkg = bkg_o[i]; @@ -632,6 +632,7 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad const float bkg_new = bkg - corr; bkg_o[i] = bkg_new; I_o[i] = isum_a[i] - (float) ninner_a[i] * bkg_new; + varbkg_o[i] -= pixel_gain * (float) ninner_a[i] * corr; // the disk's background share, as the CPU } // --- Pass B Kabsch profile fit: I = sum P(c-B)/v over sum P^2/v, v = B + max(I,0)P (iterate). The @@ -1081,7 +1082,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (radial_on) { radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>( d_rad_sum, d_rad_cnt, radial_given ? d_rad_curve.get() : nullptr, n_rad, d_k_diff, k_len, k_off, - d_isum, d_ninner, d_rbin, d_kbin, d_ok, d_bkg, d_I, n); + d_isum, d_ninner, d_rbin, d_kbin, static_cast(pixel_gain), d_ok, d_bkg, d_I, d_var_bkg, n); cuda_err(cudaGetLastError()); } diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 7b68a0014..a8b3e9ac6 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -194,12 +194,16 @@ TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[Integration][porta REQUIRE(!engine.GetRadialBackground().empty()); const auto out = engine.Run(image, predicted, predicted.size(), 0); REQUIRE(out.size() == 1); - return out[0].I; + return out[0]; }; // The curvature is along the radius only, so the error is n_inner x 0.05 x (_ring - _disk) // with u the radial offset: 49 px x 0.05 x ((6^2 + 13^2) / 4 - 4^2 / 4) = 116 counts read too low. - CHECK(integrate(false) < -80.0f); - CHECK(std::fabs(integrate(true)) < 5.0f); + const Reflection plain = integrate(false), corrected = integrate(true); + CHECK(plain.I < -80.0f); + CHECK(std::fabs(corrected.I) < 5.0f); + // The non-signal variance follows the corrected background under the signal disk (45 pixels at + // r1 = 4); the ring's own estimate noise is the same either way. + CHECK(corrected.var_bkg - plain.var_bkg == Catch::Approx(45.0 * (corrected.bkg - plain.bkg)).epsilon(1e-4)); // An empty profile hands the curve back to the reflections' own rings. BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index de125599f..889360085 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -530,6 +530,7 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat INFO("mode " << static_cast(mode) << " reflection " << i); CHECK(out_gpu[i].bkg == Catch::Approx(out_cpu[i].bkg).epsilon(0.001).margin(0.05)); CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.01).margin(2.0)); + CHECK(out_gpu[i].var_bkg == Catch::Approx(out_cpu[i].var_bkg).epsilon(0.001)); } } }