diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 64a864868..9a28c8775 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -148,6 +148,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, std::vector rough(npredicted); double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; std::vector bkg_vals; // reused per reflection for the trimmed-mean background (idea 1) + // The background ring's usable pixels and their radial offsets, in the order pass A reads them - what + // the sigma-clip pass walks, instead of the stencil a second time. + struct RingPixel { int32_t px; float rad; }; + std::vector ring_pixels; // Radial background curve, accumulated from the annulus pixels this pass already reads. A pixel's // radius is the reflection's radius plus the pixel's projection on the beam->reflection direction, @@ -188,6 +192,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // the pattern's density took, separable from what the detector took. int n_bkg_neighbour = 0; bkg_vals.clear(); + ring_pixels.clear(); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); @@ -217,6 +222,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (!valid(px)) continue; bkg_sum += static_cast(px); if (bkg_trim_frac > 0.0) bkg_vals.push_back(px); + if (do_clip) ring_pixels.push_back({px, d.rad}); ++n_bkg; } } @@ -263,23 +269,16 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(std::max(out.bkg, 1.0)); double s = 0.0; int n = 0; - for (int y = y0; y <= y1; ++y) - for (int x = x0; x <= x1; ++x) { - const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); - if (!(d.inner >= r2_sq && d.outer < r3_sq)) continue; - if (refl_mask[y * W + x]) continue; - const int32_t px = img[y * W + x]; - if (!valid(px)) continue; - if (static_cast(px) <= thr) { - s += px; - ++n; - if (bkg_radial) { - // The radial curve is binned on the TRUE detector radius, so the - // offset here is the unshrunk radial projection. - const int b = std::clamp(static_cast(std::lround(st.r0 + d.rad)), 0, n_rad - 1); - rad_sum[b] += static_cast(px); - ++rad_cnt[b]; - } + for (const auto &[px, rad] : ring_pixels) + if (static_cast(px) <= thr) { + s += px; + ++n; + if (bkg_radial) { + // The radial curve is binned on the TRUE detector radius, so the + // offset here is the unshrunk radial projection. + const int b = std::clamp(static_cast(std::lround(st.r0 + rad)), 0, n_rad - 1); + rad_sum[b] += static_cast(px); + ++rad_cnt[b]; } } if (n > 5) { out.bkg = s / n; n_bkg_used = n; }