Merge branch 'integ-crowded-6jgj' into rc173 (integration: a ring starved only by neighbours is taken whole)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

# Conflicts:
#	image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

image_analysis/bragg_integration/BraggIntegrationEngine.h

@@ -104,16 +104,17 @@ constexpr double PROFILE_SUMMATION_MAX_NSIGMA = 10.0;
 constexpr double MINPK_MAX_MISSING_PEAK = 0.9;
 }   // namespace bragg_engine
 
-// How often the integrator silently dropped a reflection, or silently declined its own fit, over
+// How often the integrator lost a reflection's clean background ring, or silently declined its own fit, over
 // every image an engine has run. Both engines keep the same two counts, so a caller sees the same
 // numbers whichever one it got. Nothing inside the engine reads them: they exist so a caller that
 // CHOSE the stencil can find out what that choice cost, which no quantity available before
 // integration measures. rugnux reads them out of its first pass (see Rugnux::RunAllPasses).
 struct BraggIntegrationCounts {
     uint64_t predicted = 0;          // reflections offered to the engine
-    uint64_t bkg_starved = 0;        // dropped whole: the r2..r3 ring kept 5 or fewer clean pixels
+    uint64_t bkg_starved = 0;        // the r2..r3 ring kept 5 or fewer clean pixels
     // Of those, the ones the NEIGHBOURS starved: their ring would have kept more than five pixels but
-    // for the ones a neighbouring reflection's signal region occupies. This is the count that answers
+    // for the ones a neighbouring reflection's signal region occupies. These are integrated against
+    // the whole ring rather than dropped; the rest of bkg_starved is dropped whole. This is the count that answers
     // "is the aperture too wide for this pattern", because the rest of bkg_starved is module gaps, the
     // beam stop and the resolution mask - a property of the detector that a wider r1 does not change.
     // Measured over the battery, that floor reaches 2.3% of all reflections while the widening that

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -155,6 +155,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     std::vector<Rough> rough(npredicted);
     double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
     std::vector<int32_t> bkg_vals;   // reused per reflection for the trimmed-mean background (idea 1)
+    std::vector<int32_t> bkg_vals_neighbour;   // the ring pixels a neighbour's region holds, likewise
 
     // 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,
@@ -211,7 +212,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         // someone else's - never this reflection's own core. That makes this an exact count of what
         // the pattern's density took, separable from what the detector took.
         int n_bkg_neighbour = 0;
+        // The readable pixels among those, kept apart so a ring the neighbours starved can fall back
+        // on them (below).
+        double bkg_sum_neighbour = 0.0;
+        int n_bkg_neighbour_valid = 0;
         bkg_vals.clear();
+        bkg_vals_neighbour.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);
@@ -237,7 +243,15 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     y_sum += y;
                     ++n_inner_valid;
                 } else if (d.inner >= r2_sq && d.outer < r3_sq) {
-                    if (refl_mask.Get(x, y)) { ++n_bkg_neighbour; continue; }
+                    if (refl_mask.Get(x, y)) {
+                        ++n_bkg_neighbour;
+                        if (valid(px)) {
+                            bkg_sum_neighbour += static_cast<double>(px);
+                            ++n_bkg_neighbour_valid;
+                            if (bkg_trim_frac > 0.0) bkg_vals_neighbour.push_back(px);
+                        }
+                        continue;
+                    }
                     if (!valid(px)) continue;
                     bkg_sum += static_cast<double>(px);
                     if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
@@ -254,9 +268,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         // survived to constrain the amplitude (XDS's MINPK, dials' valid_foreground_threshold). A box
         // sum has no profile to renormalise with, so there it stays all or nothing.
         const bool full = n_inner_valid == n_inner;
-        // A ring left with five or fewer clean pixels cannot estimate a background, so the reflection
-        // is dropped whole - the one thing the stencil geometry does to the DATA rather than to a
-        // measurement. Counted here because it is the only direct evidence of a radius that has
+        // A ring left with five or fewer readable pixels cannot estimate a background, so the
+        // reflection is dropped whole - the one thing the stencil geometry does to the DATA rather
+        // than to a measurement. Counted here because it is the only direct evidence of a radius that has
         // outgrown the pattern it is integrating (BraggIntegrationCounts).
         const bool keep_partial = full || mode != IntegratorMode::BoxSum;
         if (keep_partial && n_bkg <= 5) {
@@ -265,6 +279,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             // detector, that took it.
             if (n_bkg + n_bkg_neighbour > 5) ++counts.bkg_starved_by_neighbour;
         }
+        // A ring the NEIGHBOURS starved still holds readable pixels - only the mask set them aside.
+        // Most of those neighbours are predictions in the tail of their rocking curve, which put a
+        // few percent of their flux on this frame: on a finely sliced, dense pattern they fill the
+        // whole ring while the frame shows no spot there. Dropping the reflection loses a measurement
+        // to protect a background from flux that is mostly not there, so the ring is taken whole
+        // instead, and the high-side clip below removes the neighbour cores that are.
+        const bool ring_unmasked = keep_partial && n_bkg <= 5 && n_bkg + n_bkg_neighbour_valid > 5;
+        if (ring_unmasked) {
+            bkg_sum += bkg_sum_neighbour;
+            n_bkg += n_bkg_neighbour_valid;
+            bkg_vals.insert(bkg_vals.end(), bkg_vals_neighbour.begin(), bkg_vals_neighbour.end());
+        }
         if (keep_partial && n_bkg > 5) {
             out.bkg = bkg_sum / n_bkg;
             if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
@@ -291,7 +317,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     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.Get(x, y)) continue;
+                        if (refl_mask.Get(x, y) && !ring_unmasked) continue;
                         const int32_t px = img[y * W + x];
                         if (!valid(px)) continue;
                         if (static_cast<double>(px) <= thr) {

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -166,6 +166,9 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     __shared__ unsigned long long s_x, s_y;   // positions behind s_Isum, to take the background out of the centroid
     __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown;
     __shared__ int s_nbkg_nb;   // ring pixels a NEIGHBOUR's signal region holds; see the CPU engine
+    __shared__ int s_nbkg_nbv;  // the readable ones among them, and their sum
+    __shared__ unsigned long long s_bkgsum_nb;
+    __shared__ int s_unmasked;  // the ring was starved by neighbours and is taken whole; see the CPU engine
     __shared__ unsigned long long s_bkgsum;
     __shared__ int s_accept, s_full;
     __shared__ double s_bkg, s_thr;
@@ -190,7 +193,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
         s_r0 = sqrtf(rx * rx + ry * ry);
         s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0;
         s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0;
-        s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0;
+        s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0; s_nbkg_nbv = 0; s_bkgsum_nb = 0;
     }
     __syncthreads();
 
@@ -207,8 +210,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
 
     long long l_Isum = 0, l_Ix = 0, l_Iy = 0;
     long long l_x = 0, l_y = 0;   // positions behind l_Isum, for the background-free centroid
-    int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0;
-    long long l_bkg = 0;
+    int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0, l_nbnbv = 0;
+    long long l_bkg = 0, l_bkg_nb = 0;
     for (int t = threadIdx.x; t < area; t += blockDim.x) {
         const int x = x0 + t % bw, y = y0 + t / bw;
         const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
@@ -227,8 +230,12 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px;
                              l_x += x; l_y += y; ++l_niv; }
         } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) {
-            if (mask[y * p.W + x]) { ++l_nbnb; continue; }
             const int32_t px = img[y * p.W + x];
+            if (mask[y * p.W + x]) {
+                ++l_nbnb;
+                if (valid(px)) { l_bkg_nb += px; ++l_nbnbv; }
+                continue;
+            }
             if (!valid(px)) continue;
             l_bkg += px; ++l_nb;
         }
@@ -245,6 +252,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     WARP_ATOMIC_ADD(s_ndisk, l_nd);
     WARP_ATOMIC_ADD(s_nown, l_no);
     WARP_ATOMIC_ADD(s_nbkg_nb, l_nbnb);
+    WARP_ATOMIC_ADD(s_nbkg_nbv, l_nbnbv);
+    WARP_ATOMIC_ADD(s_bkgsum_nb, (unsigned long long) l_bkg_nb);
     __syncthreads();
 
     for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; }
@@ -254,14 +263,21 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
         // renormalises to the pixels it can read and Pass B cuts on how much of the profile survived
         // (XDS's MINPK). A box sum has no profile to renormalise with. See the CPU engine.
         s_full = (s_ninner_valid == s_ninner) ? 1 : 0;
-        s_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 1 : 0;
+        s_unmasked = 0;
         // A reflection the disk would have kept but the ring cannot support: dropped whole for want of
         // a background. One atomic per DROPPED reflection, so nothing is paid where none is dropped.
         // See BraggIntegrationCounts; the CPU engine counts the same condition.
         if ((s_full || p.partial_ok) && s_nbkg <= 5) {
             atomicAdd(counts + COUNT_BKG_STARVED, 1ull);
             if (s_nbkg + s_nbkg_nb > 5) atomicAdd(counts + COUNT_BKG_STARVED_NEIGHBOUR, 1ull);
+            // A ring the neighbours starved is taken whole instead of dropping the reflection.
+            if (s_nbkg + s_nbkg_nbv > 5) {
+                s_unmasked = 1;
+                s_nbkg += s_nbkg_nbv;
+                s_bkgsum += s_bkgsum_nb;
+            }
         }
+        s_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 1 : 0;
         s_bkg = s_accept ? ((double) (long long) s_bkgsum / (double) s_nbkg) : 0.0;
         s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0));
         // The pixel is an integer, so comparing it against the floor of the threshold accepts
@@ -285,7 +301,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             const int x = x0 + t % bw, y = y0 + t / bw;
             const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
             if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue;
-            if (mask[y * p.W + x]) continue;
+            if (mask[y * p.W + x] && !s_unmasked) continue;
             const int32_t px = img[y * p.W + x];
             if (!valid(px)) continue;
             const int slot = atomicAdd(&s_bn, 1);
@@ -327,7 +343,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             const int x = x0 + t % bw, y = y0 + t / bw;
             const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
             if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue;
-            if (mask[y * p.W + x]) continue;
+            if (mask[y * p.W + x] && !s_unmasked) continue;
             const int32_t px = img[y * p.W + x];
             if (!valid(px)) continue;
             if ((long long) px <= s_thr_i) {

tests/BraggIntegrationEngineGPUTest.cpp

@@ -140,7 +140,8 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional<float> b
     return experiment;
 }
 
-void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
+// Returns the fraction of the predicted reflections the engines kept.
+double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
                      float clip_nsigma = 4.0f, bool radial = false, int spacing = 60,
                      float stencil_k = 0.0f,
                      float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f,
@@ -200,6 +201,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
         CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.03).margin(2.0));
         CHECK(out_gpu[i].sigma == Catch::Approx(out_cpu[i].sigma).epsilon(0.03).margin(0.5));
     }
+    return static_cast<double>(out_cpu.size()) / static_cast<double>(scene.predicted.size());
 }
 
 } // namespace
@@ -269,6 +271,12 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
         CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f,
                         0.0f, 0.0f, 0.0f, OverlapMode::Exclude);
     }
+    // Spots 8 px apart: the neighbours' r2 regions cover every background ring, so only the edge of
+    // the grid keeps a clean ring pixel. Both engines have to fall back to the whole ring and keep the
+    // reflections rather than drop them for want of a background.
+    SECTION("ProfileGaussian neighbour-starved rings") {
+        CHECK(CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8) > 0.95);
+    }
     SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); }
     // Unreadable pixels inside the signal disks themselves: the MINPK rescue keeps the reflection and
     // fits it over what is left, and the peak-loss rule throws back the ones that lost the profile's