Integration: fall back to the fixed radius on a pattern too dense for it

48008e144 widens the signal radius on crystals with wide spots. On one
battery crystal - simultaneously the widest-spot and among the highest in
mosaicity - the wider radius left its neighbours' background rings with too
few clean pixels and cost 28.5% of its observations. Pass 1 now measures how
often that happens and, above a bound, hands pass 2 the settings from before
the pre-scan widened them.

The obvious quantity does not work. On the total rate of reflections dropped
for a starved ring, the losing crystal reads 4.08% and the rule's four
biggest winners read 1.89-2.04% - and a crystal at the shipped radius reads
2.37%, above all of them. Re-running the winners at r1 = 4 shows why: they
read 2.20-2.32% there too, and widening moves them down. That floor is module
gaps, the beam stop and the resolution mask, which are properties of the
detector and do not move with the radius.

So the counter separates the two. A ring is neighbour-starved when it would
have kept more than five pixels but for the pixels a neighbouring
reflection's signal region occupies. That is exact rather than estimated: the
reflection mask marks the disk inside r2 and the ring is everything outside
it, so a masked ring pixel always belongs to some other reflection's core.
The separation goes from a factor of 2 to a factor of 13 - over the twelve
crystals the radius moves, the rate is 0.000 five times, 0.001 three times,
then 0.004, 0.235, 0.315 and 4.082 - and the bound is the log-space midpoint
of that one gap, 0.0113, a factor 3.6 clear of the nearest measurement on
either side.

Predicted reflection spacing does not separate them at all: the losing
crystal is 19th of 38, a winner sits at 21.9 px, and the loosest pattern in
the battery starves 1.93% of its rings.

Battery: the space group is identical on all 38 and the merged .hkl is
byte-identical on 37, so it is inert wherever it does not fire. On the one
crystal it fires on, <I/sigma> is up 18.3%, R_meas down 29.4%, observations
up 6.3%, CC1/2 0.944 to 0.974, and its two empty top shells come back as
numbers. Its indexing rate, refined distance, beam centre and cell are
bit-identical between the two arms, so this is the guard and not the two-pass
gate.

The counters are a shared channel through both engines, summed across
workers and logged once per pass; on the GPU it is one atomic add per dropped
reflection. The profile-fit runaway guard reports on the same channel, which
is the first measurement of its trip rate.

This does not recover that crystal fully. With the adaptive radius on, pass 1
reaches a different lattice and pass 2 indexes 21% fewer frames - which
happens before the measurement this guard reads exists, and is unaffected by
it. At matched indexing rate the guard recovers 96% of the baseline's
observations against 90.5% without it.

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

image_analysis/MXAnalysisWithoutFPGA.cpp

@@ -357,6 +357,10 @@ void MXAnalysisWithoutFPGA::RunROIOnly(DataMessage &output) {
         roi->Run(*preprocessor_buffer, output.roi);
 }
 
+BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const {
+    return bragg_engine ? bragg_engine->Counts() : BraggIntegrationCounts{};
+}
+
 void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) {
     mask_low_res = settings.low_resolution_limit;
     mask_high_res = settings.high_resolution_limit;

image_analysis/MXAnalysisWithoutFPGA.h

@@ -93,6 +93,10 @@ public:
     void RebuildROI();
     void AnalyzeROIOnly(DataMessage &output);
     void RunROIOnly(DataMessage &output);
+
+    // What this worker's Bragg integrator counted (BraggIntegrationCounts). Each worker builds its own
+    // analysis, so a caller that wants the run's totals sums this over the workers it started.
+    [[nodiscard]] BraggIntegrationCounts BraggCounts() const;
 };
 
 

image_analysis/bragg_integration/BraggIntegrationEngine.h

@@ -104,6 +104,24 @@ 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
+// 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
+    // 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
+    // "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
+    // costs data adds 4.1 percentage points on top of a floor of 0.02%.
+    uint64_t bkg_starved_by_neighbour = 0;
+    uint64_t profile_fallback = 0;   // the profile fit disagreed with the box sum; the box sum was kept
+};
+
 // One reflection's extracted intensity, produced by the derived engine and turned into a
 // Reflection by Finalize() (which owns the polarization correction and scale bookkeeping).
 struct BraggFitResult {
@@ -199,6 +217,10 @@ protected:
     DiffractionGeometry geom;    // kept for the per-reflection polarization correction
     std::optional<float> polarization;
 
+    // Accumulated over every image this engine has run. An engine belongs to one worker thread, so
+    // these are plain counters and the caller sums over the engines it made.
+    BraggIntegrationCounts counts;
+
     // Assemble output reflections from the per-reflection fit results (polarization + scale corr).
     std::vector<Reflection> Finalize(const std::vector<Reflection> &predicted, size_t npredicted,
                                      const std::vector<BraggFitResult> &results,
@@ -222,4 +244,9 @@ public:
     // radial buffers were never allocated, so the CPU would correct and the GPU would not.
     void BackgroundRadial(bool on) { bkg_radial = on && bkg_radial_built; }
     [[nodiscard]] bool IsBackgroundRadialAuto() const { return bkg_radial_auto; }
+
+    // What this engine has counted since it was built (BraggIntegrationCounts). The GPU engine keeps
+    // its counts on the device and brings them back here, so this is a synchronising call - ask for it
+    // once a pass, not once an image.
+    [[nodiscard]] virtual BraggIntegrationCounts Counts() const { return counts; }
 };

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -53,6 +53,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     std::vector<BraggFitResult> results(npredicted);
     if (npredicted == 0)
         return Finalize(predicted, npredicted, results, image_number);
+    counts.predicted += npredicted;
 
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
     const bool do_clip = bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum;
@@ -166,6 +167,11 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         int n_disk = 0, n_own = 0;   // pixels in the signal disk, and how many are this reflection's
         double bkg_sum = 0.0;
         int n_bkg = 0;
+        // Ring pixels a NEIGHBOUR's signal region occupies. refl_mask marks d.inner < r2_sq and the
+        // ring is d.inner >= r2_sq, so the two are complementary and a masked ring pixel is always
+        // 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;
         bkg_vals.clear();
         for (int y = y0; y <= y1; ++y)
             for (int x = x0; x <= x1; ++x) {
@@ -188,7 +194,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     I_sum_y += static_cast<int64_t>(y) * px;
                     ++n_inner_valid;
                 } else if (d.inner >= r2_sq && d.outer < r3_sq) {
-                    if (refl_mask[y * W + x]) continue;
+                    if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; }
                     if (!valid(px)) continue;
                     bkg_sum += static_cast<double>(px);
                     if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
@@ -205,7 +211,18 @@ 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;
-        if ((full || mode != IntegratorMode::BoxSum) && n_bkg > 5) {
+        // 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
+        // outgrown the pattern it is integrating (BraggIntegrationCounts).
+        const bool keep_partial = full || mode != IntegratorMode::BoxSum;
+        if (keep_partial && n_bkg <= 5) {
+            ++counts.bkg_starved;
+            // Would the ring have been enough without the neighbours? Then it is the pattern, not the
+            // detector, that took it.
+            if (n_bkg + n_bkg_neighbour > 5) ++counts.bkg_starved_by_neighbour;
+        }
+        if (keep_partial && n_bkg > 5) {
             out.bkg = bkg_sum / n_bkg;
             if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
                 && bkg_vals.size() <= static_cast<size_t>(bragg_engine::BKG_TRIM_MAX)) {
@@ -555,6 +572,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             I = rh.I;
             sigma = rh.sigma;
             var_bkg = rh.var_bkg;
+            ++counts.profile_fallback;
         }
         // Carry the Pass-A box-sum intensity-weighted centroid (observed spot position) through the
         // profile path too - post-refinement uses it as the observed position (beam-centre / distance).

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -42,6 +42,14 @@ __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MA
 // Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL.
 constexpr int MOM_STRIDE = 3;
 
+// Slots of the device counter array behind BraggIntegrationCounts. They stay on the device for the
+// engine's lifetime and come back only when Counts() is asked for, so an image costs nothing beyond
+// the atomics themselves - and both are on paths taken by a small minority of reflections.
+constexpr int COUNT_BKG_STARVED = 0;
+constexpr int COUNT_BKG_STARVED_NEIGHBOUR = 1;
+constexpr int COUNT_PROFILE_FALLBACK = 2;
+constexpr int COUNT_SLOTS = 3;
+
 // --- Mark the r2 signal region of every predicted reflection (race-free: all writes are 1). The
 //     region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine.
 //     The same sweep fills the owner map when an overlap treatment is on: a pixel two signal regions
@@ -149,12 +157,14 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
                        float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o,
                        uint8_t *hasobs_o, unsigned long long *invd2mm,
                        float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o,
-                       unsigned long long *rad_sum, int *rad_cnt, int n_rad) {
+                       unsigned long long *rad_sum, int *rad_cnt, int n_rad,
+                       unsigned long long *counts) {
     const int i = blockIdx.x;
     if (i >= n) return;
 
     __shared__ unsigned long long s_Isum, s_Ix, s_Iy;
     __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__ unsigned long long s_bkgsum;
     __shared__ int s_accept, s_full;
     __shared__ double s_bkg, s_thr;
@@ -179,7 +189,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_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0;
-        s_ndisk = 0; s_nown = 0;
+        s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0;
     }
     __syncthreads();
 
@@ -195,7 +205,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     const int area = (bw > 0 && bh > 0) ? bw * bh : 0;
 
     long long l_Isum = 0, l_Ix = 0, l_Iy = 0;
-    int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0;
+    int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0;
     long long l_bkg = 0;
     for (int t = threadIdx.x; t < area; t += blockDim.x) {
         const int x = x0 + t % bw, y = y0 + t / bw;
@@ -214,7 +224,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             const int32_t px = img[y * p.W + x];
             if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; ++l_niv; }
         } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) {
-            if (mask[y * p.W + x]) continue;
+            if (mask[y * p.W + x]) { ++l_nbnb; continue; }
             const int32_t px = img[y * p.W + x];
             if (!valid(px)) continue;
             l_bkg += px; ++l_nb;
@@ -229,6 +239,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     WARP_ATOMIC_ADD(s_bkgsum, (unsigned long long) l_bkg);
     WARP_ATOMIC_ADD(s_ndisk, l_nd);
     WARP_ATOMIC_ADD(s_nown, l_no);
+    WARP_ATOMIC_ADD(s_nbkg_nb, l_nbnb);
     __syncthreads();
 
     for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; }
@@ -239,6 +250,13 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
         // (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;
+        // 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);
+        }
         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
@@ -567,7 +585,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x
                     const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a,
                     const float *shell_P, const float *global_P,
                     const float *sigma2_r, const float *sigma2_t, const int *shell_n,
-                    float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n) {
+                    float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n,
+                    unsigned long long *counts) {
     const int i = blockIdx.x;
     if (i >= n) return;
     extern __shared__ float Pbuf[];
@@ -737,6 +756,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x
                 I = I_seed[i];
                 sigma = sigma_seed[i];
                 var_bkg = varbkg_seed[i];
+                atomicAdd(counts + COUNT_PROFILE_FALLBACK, 1ull);
             }
             I_o[i] = I; sigma_o[i] = sigma; varbkg_o[i] = var_bkg; ok_o[i] = 1;
         } else ok_o[i] = 0;
@@ -759,8 +779,10 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment
       d_sigma2_t(bragg_engine::N_SHELL + 1),
       d_shell_n(bragg_engine::N_SHELL),
       d_global_n(1),
-      d_invd2(2) {
+      d_invd2(2),
+      d_counts(COUNT_SLOTS) {
     threads = 128;
+    cuda_err(cudaMemset(d_counts, 0, sizeof(unsigned long long) * COUNT_SLOTS));
 
     // Fit profile grid: R for empirical / box, up to 3R (radially elongated) for the Gaussian.
     const int max_Rf = empirical ? R : 3 * R;
@@ -856,6 +878,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
     std::vector<BraggFitResult> results(npredicted);
     if (image.size() != npixel || npredicted == 0)
         return Finalize(predicted, npredicted, results, image_number);
+    counts.predicted += npredicted;
 
     const int32_t *img = image.getGPUBuffer();
     if (img == nullptr)
@@ -925,7 +948,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
                                        d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y,
                                        d_ok, d_strong, d_has_obs, d_invd2,
                                        d_isum, d_ninner, d_rbin, d_kbin,
-                                       d_rad_sum, d_rad_cnt, rad_n);
+                                       d_rad_sum, d_rad_cnt, rad_n, d_counts);
 
     // Correct the flat annulus background for the curvature of the radial background before anything
     // downstream (profile fit, variance) reads it.
@@ -945,7 +968,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
         fit<<<n, threads, fit_shared_bytes, *stream>>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2,
                                                        d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P,
                                                        d_sigma2_r, d_sigma2_t, d_shell_n,
-                                                       d_I, d_sigma, d_var_bkg, d_ok, p, n);
+                                                       d_I, d_sigma, d_var_bkg, d_ok, p, n, d_counts);
     }
 
     // Nothing below reads the mask or the owner map, so put them back now, over the boxes that were
@@ -983,3 +1006,14 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
     }
     return Finalize(predicted, npredicted, results, image_number);
 }
+
+BraggIntegrationCounts BraggIntegrationEngineGPU::Counts() const {
+    BraggIntegrationCounts out = counts;   // `predicted` is counted on the host, as on the CPU engine
+    unsigned long long host[COUNT_SLOTS] = {0, 0, 0};
+    cuda_err(cudaMemcpyAsync(host, d_counts, sizeof(host), cudaMemcpyDeviceToHost, *stream));
+    cuda_err(cudaStreamSynchronize(*stream));
+    out.bkg_starved = host[COUNT_BKG_STARVED];
+    out.bkg_starved_by_neighbour = host[COUNT_BKG_STARVED_NEIGHBOUR];
+    out.profile_fallback = host[COUNT_PROFILE_FALLBACK];
+    return out;
+}

image_analysis/bragg_integration/BraggIntegrationEngineGPU.h

@@ -64,6 +64,9 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
     CudaDevicePtr<float>   d_sigma2_r, d_sigma2_t;        // radial/tangential widths, N_SHELL + global
     CudaDevicePtr<int>     d_shell_n, d_global_n;
     CudaDevicePtr<unsigned long long> d_invd2;           // [min,max] inv-d^2 as monotonic bit patterns
+    // BraggIntegrationCounts, accumulated on the device so an image costs no transfer; brought back
+    // only when Counts() is asked for.
+    CudaDevicePtr<unsigned long long> d_counts;
 
     // --- host staging (copied back once per frame) ---
     // Pinned, like every other engine's staging: a copy out of pageable memory does not return until the
@@ -80,4 +83,8 @@ public:
     std::vector<Reflection> Run(const ImagePreprocessorBuffer &image,
                                 const std::vector<Reflection> &predicted, size_t npredicted,
                                 int64_t image_number) override;
+
+    // Brings the two device counters back before answering. Synchronises the stream, so ask once a
+    // pass rather than once an image.
+    [[nodiscard]] BraggIntegrationCounts Counts() const override;
 };