Metal Bragg integration engine on macOS, output unchanged (bit-identical to the CPU engine)

BraggIntegrationEngineMetal is BraggIntegrationEngineCPU with its integer pixel work run as Metal
kernels, one thread per reflection: the neighbour mask, the owner map, pass A's walk over every box
and the background clip. Everything in double, the profile learning and pass B stay on the CPU, in
the CPU engine's own code - Apple GPUs have no double, and no GPU exp matches macOS expf bit for bit.
The sums are integers, so the GPU's order does not matter; the stencil arithmetic is the same
BraggStencil.h, compiled at run time (MSL 3.0, no offline Metal toolchain needed) with fast math off
(sqrt and division correctly rounded, measured) and fp contraction off, as the host already builds
(-ffp-contract=off, now extended to Objective-C++). The image is read in place (unified memory).

The CPU engine is split so both engines share the finishing code: pass A counts into BraggBoxSums
(integer; the double ring sums were exact integers already), the clip threshold is its integer floor,
FinishPassA and ProfilePasses are shared. Verified byte-identical p.hkl and MTZ data before/after the
split and with the Metal engine on lysozyme, cytochrome C, myoglobin, 4nwv and a small-molecule set.

Measured on an Apple M6 (wall, CPU engine -> Metal): 4nwv 180.5 -> 151.2 s, cytochrome C 72.7 ->
59.8 s, lysozyme 29.2 -> 22.5 s, myoglobin 43.3 -> 41.2 s; the sparse small-molecule set +1.5 s
(few reflections per frame leave the GPU idle).

CI: a macOS step runs tests/BraggIntegrationEngineMetalTest.cpp, which requires every float of every
reflection and every count to match the CPU engine over 19 setting combinations on random crowded
frames (negative pixels, exact-pixel ties, hot, masked and overloaded pixels), three frames per
engine; it fails without a Metal device and if the GPU path fell back to the CPU. Ten planted kernel
bugs were each caught. rugnux --no-gpu hides the Metal device too (gpus_are_hidden).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -36,13 +36,6 @@ inline void prefetch_row(const void *first, const void *last) {
     prefetch(last);
 }
 
-// Identity sampler over the preprocessed int32 buffer (already in that convention).
-struct BufferSampler {
-    const int32_t *p;
-    int32_t operator[](size_t i) const { return p[i]; }
-    const void *address(size_t i) const { return p + i; }
-};
-
 // Sampler over a raw detector image of pixel type T: masked pixels carry the type minimum, saturated
 // the type maximum (the FPGA image has no lossy-codec +/-1 band). Only the pixels actually read - the
 // reflection disks - are converted, so there is no whole-image pass.
@@ -65,9 +58,9 @@ struct ImageSampler {
 BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment,
                                                      const PixelMask &mask)
     : BraggIntegrationEngine(experiment),
-      static_mask(mask.GetPackedMask()),
       refl_mask(static_cast<int>(xpixel), static_cast<int>(ypixel), 0),
-      owner(static_cast<int>(xpixel), static_cast<int>(ypixel), BRAGG_OWNER_NONE) {}
+      owner(static_cast<int>(xpixel), static_cast<int>(ypixel), BRAGG_OWNER_NONE),
+      static_mask(mask.GetPackedMask()) {}
 
 void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector<Reflection> &predicted, size_t npredicted) {
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
@@ -133,23 +126,13 @@ template <class Sampler>
 std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                                                            const std::vector<Reflection> &predicted,
                                                            size_t npredicted, int64_t image_number) {
-    std::vector<BraggFitResult> results(npredicted);
     if (npredicted == 0)
-        return Finalize(predicted, npredicted, results, image_number);
+        return Finalize(predicted, npredicted, std::vector<BraggFitResult>(npredicted), 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;
 
-    // Symmetric trimmed-mean background fraction (BraggIntegrationSettings, rugnux --background-trim):
-    // replaces the r2..r3 ring MEAN with an f-trimmed mean, robust to the high-side contamination
-    // (neighbour wings, tails) that biases the plain mean up and makes it over-subtract weak high-angle
-    // reflections. The base ctor has already forced it to 0 whenever the clip below is in force, which
-    // is the default - the two are alternatives.
-    const double bkg_trim_frac = bkg_trim;
-
-    auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); };
-
     MarkNeighbours(predicted, npredicted);
     owner.Clear();
 
@@ -184,20 +167,6 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     };
 
     // --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). ---
-    struct Rough {
-        double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0;
-        double bkg_var = 0.0;   // variance of the background ESTIMATE itself, bkg / n_bkg
-        double var_bkg = 0.0;   // total non-signal variance carried to the merge
-        int64_t I_sum = 0;      // kept so I can be rebuilt after the radial background correction
-        int n_inner = 0;
-        int n_disk = 0, n_own = 0;   // signal-disk pixels, and how many of them are this reflection's
-        int r_bin = 0;          // rounded distance from the beam centre, indexes the radial curve
-        int k_bin = 0;          // which radial-background kernel this reflection's stencil needs
-        int cx = 0, cy = 0, shell = -1;
-        bool ok = false, strong = false, has_obs = false;
-        bool full = false;      // every pixel of the signal disk was readable
-        bool overloaded = false;   // a pixel of the signal disk was saturated
-    };
     // A saturated pixel, or one unreadable on this frame alone - the overload marker of a detector
     // that writes its error value for a pixel it could not count (EIGER), which the preprocessor turns
     // into a masked pixel like a gap's. The run's mask holds the gaps, so whatever is unreadable beyond
@@ -205,26 +174,13 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     auto static_masked = [&](size_t idx) {
         return idx / 32 < static_mask.size() && ((static_mask[idx / 32] >> (idx % 32)) & 1U);
     };
-    std::vector<Rough> rough(npredicted);
-    double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
+    PassA a = StartPassA(npredicted);
     std::vector<int32_t> bkg_vals;   // reused per reflection for the trimmed-mean background (idea 1)
     // The ring pixels the background clip below reads, with their radial offsets, gathered while pass
-    // A walks the box anyway - in the order the clip sums them, so it sums the same values.
+    // A walks the box anyway.
     struct ClipPixel { int32_t px; float rad; };
     std::vector<ClipPixel> clip_ring;
 
-    // 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,
-    // so no per-pixel sqrt is needed.
-    const int n_rad = static_cast<int>(std::ceil(std::hypot(std::max<double>(beam_x, W - beam_x),
-                                                            std::max<double>(beam_y, H - beam_y)))) + 2;
-    std::vector<double> rad_sum;
-    std::vector<int> rad_cnt;
-    if (bkg_radial) {
-        rad_sum.assign(n_rad, 0.0);
-        rad_cnt.assign(n_rad, 0);
-    }
-
     // The box pass A reads: the OUTER ellipse, tightly - taking the largest semi-axis in both
     // directions instead would read up to 60% more pixels for nothing.
     struct Box { int x0, x1, y0, y1; };
@@ -247,28 +203,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                              img.address(static_cast<size_t>(y) * W + b.x1));
         }
         const auto &r = predicted[i];
-        Rough out;
         // This reflection's stencil: the r1 signal disk, and the r2..r3 ring elongated radially by
         // the analytic smear.
         const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
         const Box box = box_of(r, st);
         const int x0 = box.x0, x1 = box.x1, y0 = box.y0, y1 = box.y1;
 
-        // Both from the stencil's own radius, so nothing downstream is derived from a second one.
-        out.r_bin = std::clamp(static_cast<int>(std::lround(st.r0)), 0, n_rad - 1);
-        out.k_bin = BraggStencilKernelIndex(st, n_kern);
-
-        int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0;
-        int64_t x_sum = 0, y_sum = 0;   // positions of the pixels behind I_sum, for the centroid
-        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;
-        // Readable ring pixels, and those of them no predicted neighbour's signal region covers.
-        // 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
+        // n_ring / n_clean: readable ring pixels, and those of them no predicted neighbour's signal
+        // region covers. 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 the difference an exact count of what the pattern's density took,
         // separable from what the detector took.
-        int n_ring = 0, n_clean = 0;
+        BraggBoxSums s{};
         bkg_vals.clear();
         clip_ring.clear();
         for (int y = y0; y <= y1; ++y)
@@ -282,154 +228,223 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     // mode it leaves the disk entirely - the sum, the pixel count the background is
                     // subtracted with, and the all-or-nothing validity gate alike. The box sum only
                     // counts how many were lost, which is all it can act on.
-                    ++n_disk;
+                    ++s.n_disk;
                     if (overlap != OverlapMode::Off) {
-                        if (clean(x, y, i)) ++n_own;
+                        if (clean(x, y, i)) ++s.n_own;
                         else if (exclude) continue;
                     }
-                    ++n_inner;
+                    ++s.n_inner;
                     if (px == INT32_MAX
                         || (px == INT32_MIN && !static_masked(static_cast<size_t>(y) * W + x)))
-                        out.overloaded = true;
+                        s.overloaded = 1;
                     if (!valid(px)) continue;
-                    I_sum += px;
-                    I_sum_x += static_cast<int64_t>(x) * px;
-                    I_sum_y += static_cast<int64_t>(y) * px;
+                    s.I_sum += px;
+                    s.I_sum_x += static_cast<int64_t>(x) * px;
+                    s.I_sum_y += static_cast<int64_t>(y) * px;
                     // Position sums over the very same pixels, so the background can be taken out of
                     // the centroid below. It is not known yet - the ring is read in this same loop.
-                    x_sum += x;
-                    y_sum += y;
-                    ++n_inner_valid;
+                    s.x_sum += x;
+                    s.y_sum += y;
+                    ++s.n_inner_valid;
                 } else if (d.inner >= r2_sq && d.outer < r3_sq) {
                     if (!valid(px)) continue;
-                    ++n_ring;
+                    ++s.n_ring;
                     const uint8_t m = refl_mask.Get(x, y);
-                    if (m == 0) ++n_clean;
+                    if (m == 0) ++s.n_clean;
                     if (m == MASK_FLUX) continue;
-                    bkg_sum += static_cast<double>(px);
-                    if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
-                    ++n_bkg;
+                    s.bkg_sum += px;
+                    if (bkg_trim > 0.0f) bkg_vals.push_back(px);
+                    ++s.n_bkg;
                 }
             }
 
-        if (use_owner && n_own < n_disk) {
-            ++counts.overlapped;
-            if (n_own < overlap_min_peak * n_disk) ++counts.overlapped_major;
-        }
-
-        int n_bkg_used = n_bkg;   // pixels behind the FINAL background value (trim/clip shrink it)
-        // A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the
-        // reflection outright. A profile fit does not need it to: the fit is the amplitude of a
-        // NORMALISED profile, so leaving pixels out renormalises the estimator by construction and
-        // widens sigma by the information they carried - the same argument Exclude makes for a
-        // neighbour's pixels. Pass B keeps the reflection only while enough of the expected profile
-        // 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 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 || out.overloaded;
-        if (keep_partial && n_bkg <= 5) ++counts.bkg_starved;
-        // Would every predicted neighbour, the tails included, leave the ring starved that the
-        // detector alone would not? That is the pattern's density, what the guard on a widened
-        // radius reads (BraggIntegrationCounts).
-        if (keep_partial && n_clean <= 5 && n_ring > 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)) {
-                // Symmetric trimmed mean over the background ring (idea 1): drop the lowest and highest
-                // bkg_trim_frac of the pixels, average the rest. Robust to the high-side contamination
-                // that biases the plain ring mean and makes it over-subtract at high resolution.
-                std::sort(bkg_vals.begin(), bkg_vals.end());
-                const size_t lo = static_cast<size_t>(bkg_vals.size() * bkg_trim_frac);
-                const size_t hi = bkg_vals.size() - lo;
-                if (hi > lo) {
-                    double s = 0.0;
-                    for (size_t t = lo; t < hi; ++t) s += bkg_vals[t];
-                    out.bkg = s / static_cast<double>(hi - lo);
-                    n_bkg_used = static_cast<int>(hi - lo);
-                }
-            } else if (do_clip) {
-                // One high-outlier sigma-clip pass on the background ring: reject pixels above
-                // mean + n*sqrt(gain*mean) to strip a neighbour core or a zinger that biases the mean.
-                // The gain is 1 for counted photons, so this is then exactly the Poisson clip.
-                const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(static_cast<double>(bkg_pixel_gain)
-                                                                         * std::max(out.bkg, 1.0));
-                double s = 0.0;
-                int n = 0;
-                for (const auto &[px, rad] : clip_ring)
-                    if (static_cast<double>(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<int>(std::lround(st.r0 + rad)), 0, n_rad - 1);
-                            rad_sum[b] += static_cast<double>(px);
-                            ++rad_cnt[b];
-                        }
+        // One high-outlier sigma-clip pass on the background ring (see ClipThreshold), and the radial
+        // background curve from the pixels it keeps.
+        BraggClipSums clip{};
+        if (ClipsRing(s)) {
+            const int32_t thr = ClipThreshold(s);
+            for (const auto &[px, rad] : clip_ring)
+                if (px <= thr) {
+                    clip.sum += px;
+                    ++clip.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<int>(std::lround(st.r0 + rad)), 0, a.n_rad - 1);
+                        a.rad_sum[b] += px;
+                        ++a.rad_cnt[b];
                     }
-                if (n > 5) { out.bkg = s / n; n_bkg_used = n; }
-            }
-            // The sum is over the pixels actually READ, so that is the count the background is
-            // subtracted with; with nothing missing it is the whole disk, exactly as before.
-            out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner_valid) * out.bkg;
-            // I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its
-            // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
-            // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
-            // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
-            out.bkg_var = out.bkg / n_bkg_used;
-            out.var_bkg = static_cast<double>(n_inner_valid) * out.bkg
-                          + static_cast<double>(n_inner_valid) * n_inner_valid * out.bkg_var;
-            out.I_sum = I_sum;
-            out.n_inner = static_cast<int>(n_inner_valid);
-            out.n_disk = n_disk;
-            out.n_own = n_own;
-            const double var_bkg_term = static_cast<double>(n_inner_valid) * n_inner_valid * out.bkg_var;
-            out.sigma = 1.0;
-            if (I_sum > 0) {
-                out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(I_sum) + var_bkg_term));
-                // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by
-                // signal plus background, and the background is flat over a disk centred on the
-                // PREDICTION - so its own centroid is the prediction exactly, and it pulls the answer
-                // there. The measured displacement comes out shrunk by I/(I + n*bkg), which is worst
-                // where the background dominates, i.e. at high resolution. Post-refinement fits the
-                // 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.
-                const double net = static_cast<double>(I_sum) - n_inner_valid * out.bkg;
-                if (net > 0.0) {
-                    out.obs_x = (static_cast<double>(I_sum_x) - out.bkg * x_sum) / net;
-                    out.obs_y = (static_cast<double>(I_sum_y) - out.bkg * y_sum) / net;
-                } else {
-                    // Nothing above background to take a centroid of; the raw one is all there is.
-                    out.obs_x = static_cast<double>(I_sum_x) / static_cast<double>(I_sum);
-                    out.obs_y = static_cast<double>(I_sum_y) / static_cast<double>(I_sum);
                 }
-                // A disk with a hole in it gives a centroid pulled away from the hole, and the hole
-                // sits at a fixed place on the detector - post-refinement would read that as geometry.
-                out.has_obs = full;
-            }
-            out.cx = static_cast<int>(std::lround(r.predicted_x));
-            out.cy = static_cast<int>(std::lround(r.predicted_y));
-            out.ok = true;
-            out.full = full;
-            // The profile, its resolution shells and their widths are learned from COMPLETE
-            // reflections: a partial one has no intensity to normalise its grid by and a hole to
-            // stack into it. So the learning sees exactly what it saw before this rescue existed.
-            out.strong = full && out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA;
-            if (full && r.d > 0.0f) {
-                const double inv_d2 = 1.0 / (static_cast<double>(r.d) * r.d);
-                inv_d2_min = std::min(inv_d2_min, inv_d2);
-                inv_d2_max = std::max(inv_d2_max, inv_d2);
-            }
         }
-        rough[i] = out;
+        FinishPassA(a, i, r, st, s, clip, bkg_vals);
     }
 
+    return ProfilePasses(img, owner, predicted, npredicted, a, image_number);
+}
+
+BraggIntegrationEngineCPU::PassA BraggIntegrationEngineCPU::StartPassA(size_t npredicted) const {
+    const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
+    PassA a;
+    a.rough.resize(npredicted);
+    // A pixel's radius is the reflection's radius plus the pixel's projection on the beam->reflection
+    // direction, so no per-pixel sqrt is needed.
+    a.n_rad = static_cast<int>(std::ceil(std::hypot(std::max<double>(beam_x, W - beam_x),
+                                                    std::max<double>(beam_y, H - beam_y)))) + 2;
+    if (bkg_radial) {
+        a.rad_sum.assign(a.n_rad, 0);
+        a.rad_cnt.assign(a.n_rad, 0);
+    }
+    a.inv_d2_min = std::numeric_limits<double>::max();
+    a.inv_d2_max = 0.0;
+    return a;
+}
+
+bool BraggIntegrationEngineCPU::ClipsRing(const BraggBoxSums &s) const {
+    // The clip is the alternative to the trimmed mean, so where it is on the trim is off, and the
+    // profile modes keep every reflection with a background (FinishPassA).
+    return bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum && s.n_bkg > 5;
+}
+
+int32_t BraggIntegrationEngineCPU::ClipThreshold(const BraggBoxSums &s) const {
+    // Reject pixels above mean + n*sqrt(gain*mean) to strip a neighbour core or a zinger that biases
+    // the mean. The gain is 1 for counted photons, so this is then exactly the Poisson clip. A pixel
+    // value is an integer, so it is at or below the threshold exactly when it is at or below its floor.
+    const double bkg = static_cast<double>(s.bkg_sum) / s.n_bkg;
+    const double thr = bkg + bkg_clip_nsigma * std::sqrt(static_cast<double>(bkg_pixel_gain) * std::max(bkg, 1.0));
+    return static_cast<int32_t>(std::clamp(std::floor(thr), static_cast<double>(INT32_MIN),
+                                           static_cast<double>(INT32_MAX)));
+}
+
+void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection &r, const BraggStencil &st,
+                                            const BraggBoxSums &s, const BraggClipSums &clip,
+                                            std::vector<int32_t> &bkg_vals) {
+    Rough out;
+    // Both from the stencil's own radius, so nothing downstream is derived from a second one.
+    out.r_bin = std::clamp(static_cast<int>(std::lround(st.r0)), 0, a.n_rad - 1);
+    out.k_bin = BraggStencilKernelIndex(st, n_kern);
+    out.overloaded = s.overloaded != 0;
+
+    if (overlap != OverlapMode::Off && s.n_own < s.n_disk) {
+        ++counts.overlapped;
+        if (s.n_own < overlap_min_peak * s.n_disk) ++counts.overlapped_major;
+    }
+
+    int n_bkg_used = s.n_bkg;   // pixels behind the FINAL background value (trim/clip shrink it)
+    // A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the
+    // reflection outright. A profile fit does not need it to: the fit is the amplitude of a
+    // NORMALISED profile, so leaving pixels out renormalises the estimator by construction and
+    // widens sigma by the information they carried - the same argument Exclude makes for a
+    // neighbour's pixels. Pass B keeps the reflection only while enough of the expected profile
+    // 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 = s.n_inner_valid == s.n_inner;
+    // 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 || out.overloaded;
+    if (keep_partial && s.n_bkg <= 5) ++counts.bkg_starved;
+    // Would every predicted neighbour, the tails included, leave the ring starved that the
+    // detector alone would not? That is the pattern's density, what the guard on a widened
+    // radius reads (BraggIntegrationCounts).
+    if (keep_partial && s.n_clean <= 5 && s.n_ring > 5) ++counts.bkg_starved_by_neighbour;
+    if (keep_partial && s.n_bkg > 5) {
+        out.bkg = static_cast<double>(s.bkg_sum) / s.n_bkg;
+        // Symmetric trimmed-mean background fraction (BraggIntegrationSettings, rugnux --background-trim):
+        // replaces the r2..r3 ring MEAN with an f-trimmed mean, robust to the high-side contamination
+        // (neighbour wings, tails) that biases the plain mean up and makes it over-subtract weak
+        // high-angle reflections. The base ctor has already forced it to 0 whenever the clip is in
+        // force, which is the default - the two are alternatives.
+        const double bkg_trim_frac = bkg_trim;
+        if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
+            && bkg_vals.size() <= static_cast<size_t>(bragg_engine::BKG_TRIM_MAX)) {
+            // Symmetric trimmed mean over the background ring (idea 1): drop the lowest and highest
+            // bkg_trim_frac of the pixels, average the rest. Robust to the high-side contamination
+            // that biases the plain ring mean and makes it over-subtract at high resolution.
+            std::sort(bkg_vals.begin(), bkg_vals.end());
+            const size_t lo = static_cast<size_t>(bkg_vals.size() * bkg_trim_frac);
+            const size_t hi = bkg_vals.size() - lo;
+            if (hi > lo) {
+                double sum = 0.0;
+                for (size_t t = lo; t < hi; ++t) sum += bkg_vals[t];
+                out.bkg = sum / static_cast<double>(hi - lo);
+                n_bkg_used = static_cast<int>(hi - lo);
+            }
+        } else if (ClipsRing(s)) {
+            if (clip.n > 5) { out.bkg = static_cast<double>(clip.sum) / clip.n; n_bkg_used = clip.n; }
+        }
+        // The sum is over the pixels actually READ, so that is the count the background is
+        // subtracted with; with nothing missing it is the whole disk, exactly as before.
+        out.I = static_cast<double>(s.I_sum) - static_cast<double>(s.n_inner_valid) * out.bkg;
+        // I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its
+        // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
+        // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
+        // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
+        out.bkg_var = out.bkg / n_bkg_used;
+        out.var_bkg = static_cast<double>(s.n_inner_valid) * out.bkg
+                      + static_cast<double>(s.n_inner_valid) * s.n_inner_valid * out.bkg_var;
+        out.I_sum = s.I_sum;
+        out.n_inner = static_cast<int>(s.n_inner_valid);
+        out.n_disk = s.n_disk;
+        out.n_own = s.n_own;
+        const double var_bkg_term = static_cast<double>(s.n_inner_valid) * s.n_inner_valid * out.bkg_var;
+        out.sigma = 1.0;
+        if (s.I_sum > 0) {
+            out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(s.I_sum) + var_bkg_term));
+            // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by
+            // signal plus background, and the background is flat over a disk centred on the
+            // PREDICTION - so its own centroid is the prediction exactly, and it pulls the answer
+            // there. The measured displacement comes out shrunk by I/(I + n*bkg), which is worst
+            // where the background dominates, i.e. at high resolution. Post-refinement fits the
+            // 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.
+            const double net = static_cast<double>(s.I_sum) - s.n_inner_valid * out.bkg;
+            if (net > 0.0) {
+                out.obs_x = (static_cast<double>(s.I_sum_x) - out.bkg * s.x_sum) / net;
+                out.obs_y = (static_cast<double>(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.
+                out.obs_x = static_cast<double>(s.I_sum_x) / static_cast<double>(s.I_sum);
+                out.obs_y = static_cast<double>(s.I_sum_y) / static_cast<double>(s.I_sum);
+            }
+            // A disk with a hole in it gives a centroid pulled away from the hole, and the hole
+            // sits at a fixed place on the detector - post-refinement would read that as geometry.
+            out.has_obs = full;
+        }
+        out.cx = static_cast<int>(std::lround(r.predicted_x));
+        out.cy = static_cast<int>(std::lround(r.predicted_y));
+        out.ok = true;
+        out.full = full;
+        // The profile, its resolution shells and their widths are learned from COMPLETE
+        // reflections: a partial one has no intensity to normalise its grid by and a hole to
+        // stack into it. So the learning sees exactly what it saw before this rescue existed.
+        out.strong = full && out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA;
+        if (full && r.d > 0.0f) {
+            const double inv_d2 = 1.0 / (static_cast<double>(r.d) * r.d);
+            a.inv_d2_min = std::min(a.inv_d2_min, inv_d2);
+            a.inv_d2_max = std::max(a.inv_d2_max, inv_d2);
+        }
+    }
+    a.rough[i] = out;
+}
+
+template <class Sampler, class OwnerMap>
+std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &img, const OwnerMap &owner_map,
+                                                                 const std::vector<Reflection> &predicted,
+                                                                 size_t npredicted, PassA &a,
+                                                                 int64_t image_number) {
+    const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
+    std::vector<BraggFitResult> results(npredicted);
+    std::vector<Rough> &rough = a.rough;
+    const bool exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum;
+    const bool use_owner = overlap != OverlapMode::Off;
+    auto clean = [&](int x, int y, size_t i) {
+        return !use_owner
+               || BraggOwnedBy(owner_map.Get(x, y), static_cast<int>(i));
+    };
+    auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); };
+
     // --- Radial background curvature correction. The annulus mean is blind to the curvature of the
     //     radial background (a linear background cancels between the concentric disk and annulus), so
     //     correct it by mean_annulus(B) - mean_disk(B) taken from the curve just accumulated. Reads no
@@ -443,8 +458,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             const float *kern = k_diff.data() + static_cast<size_t>(rh.k_bin) * k_len;
             double corr = 0.0;
             for (int k = 0; k < k_len; ++k) {
-                const int b = std::clamp(rh.r_bin + k - k_off, 0, n_rad - 1);
-                const double v = rad_cnt[b] > 0 ? rad_sum[b] / rad_cnt[b] : rh.bkg;
+                const int b = std::clamp(rh.r_bin + k - k_off, 0, a.n_rad - 1);
+                const double v = a.rad_cnt[b] > 0 ? static_cast<double>(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg;
                 corr += static_cast<double>(kern[k]) * v;
             }
             rh.bkg -= corr;                       // annulus mean -> mean over the signal disk
@@ -473,8 +488,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     }
 
     auto shell_of = [&](float d) {
-        if (!(d > 0.0f) || inv_d2_max <= inv_d2_min) return 0;
-        const double t = (1.0 / (static_cast<double>(d) * d) - inv_d2_min) / (inv_d2_max - inv_d2_min);
+        if (!(d > 0.0f) || a.inv_d2_max <= a.inv_d2_min) return 0;
+        const double t = (1.0 / (static_cast<double>(d) * d) - a.inv_d2_min) / (a.inv_d2_max - a.inv_d2_min);
         return std::clamp(static_cast<int>(t * N_SHELL), 0, N_SHELL - 1);
     };
     for (size_t i = 0; i < npredicted; ++i)
@@ -749,6 +764,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     return Finalize(predicted, npredicted, results, image_number);
 }
 
+template std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(
+        const BraggIntegrationEngineCPU::BufferSampler &, const BraggIntegrationEngineCPU::FlatOwnerMap &,
+        const std::vector<Reflection> &, size_t, BraggIntegrationEngineCPU::PassA &, int64_t);
+
 std::vector<Reflection> BraggIntegrationEngineCPU::Run(const ImagePreprocessorBuffer &image,
                                                        const std::vector<Reflection> &predicted,
                                                        size_t npredicted, int64_t image_number) {