Stills integration: background ring reach measured on empty apertures, radial curvature from the whole frame

Two changes to the stills integration background, both stills-only (rotation keeps its path bit for bit):

1. The r2..r3 background ring's outer radius is measured per run. The pre-scan lays apertures on a
   grid over its frames, away from every found spot, and for each candidate r3 (r2+2 .. 14) takes the
   mean square difference between the r1 disk mean and the clipped ring mean. With nothing in the
   aperture that is exactly the background error the ring makes (noise + bias); the run keeps the
   minimum. On a photon-counting pixel detector the noise wins and the shipped r3=14 is kept; on a
   binned CCD the background structure wins and the narrowest ring is taken. r1 stays 6: a smaller
   disk (4,6,14) costs 0.0013 dimple R_free on a JUNGFRAU set, so it does not pass the
   photon-counting gate. --adaptive-integration-radius now defaults on for stills too (=off restores
   the fixed 6/8/14). The pre-scan bandwidth estimate is not applied to stills.

2. The radial background-curvature correction on stills takes its curve from the adaptive spot
   finder's peak-excluded per-ring background (whole frame) instead of the predicted reflections'
   own ring pixels, and auto applies it to every stills image (no ice-score gate: a solvent/carrier
   hump curves the background as much as ice). New engine call RadialBackgroundFromRings; CPU, GPU
   and Metal read the handed-over curve. The GPU box sum now sizes the radial bin index separately
   from accumulating, so a handed-over curve is read at the right radius.

Measured (dimple R_free, deposited free set, battery arguments, all frames unless noted):
  pal86  (Rayonix 4x4, 156 um): 0.354 -> 0.218   (ring 10; at 4,6,10 0.213)
  lcls222 (Rayonix 4x4, 177 um): 0.205 -> 0.196  (ring 10)
  sacla236 (MPCCD): 0.2588 -> 0.2558             (ring 11)
  alvra213_int16 (JUNGFRAU): 0.1986 -> 0.1986    (ring 14, unchanged)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu

image_analysis/MXAnalysisWithoutFPGA.cpp

@@ -235,22 +235,37 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
         // soon as it has detected - so this must be called after detection and before integration.
         // Where no such background exists (no adaptive finder), auto leaves the correction off: the
         // plain profile carries the Bragg peaks and cannot support an absolute threshold.
+        //
+        // On stills the correction also takes its CURVE from that background, and auto applies it to
+        // every image that has one. A serial frame has a few dozen predicted reflections, so their own
+        // ring pixels - where the curve comes from otherwise - sample the radial background sparsely,
+        // while the finder's rings hold every pixel of the frame. And the correction is the measured
+        // curvature itself: where the background has none it is zero, so there is nothing for a gate
+        // to decide. The ice score is the wrong gate there: a smooth hump that is not ice - the
+        // solvent or a carrier medium - curves the background as much as an ice ring does, and the
+        // score does not see it.
         const auto decide_radial_background = [this, &spot_finding_settings, &output]() {
-            if (!bragg_engine->IsBackgroundRadialAuto())
-                return;
             // Same condition as the score at the end of this function: the adaptive finder holds its
             // ring background from whenever it last ran, so requiring that it ran for THIS image is
             // what keeps a stale curve out.
-            if (!spot_finding_settings.adaptive_threshold)
-                return;
-            const std::vector<float> &ring_bkg = adaptiveSpotFinder->GetRingBackground();
-            if (ring_bkg.empty())
+            const std::vector<float> no_rings;
+            const std::vector<float> &ring_bkg = spot_finding_settings.adaptive_threshold
+                                                 ? adaptiveSpotFinder->GetRingBackground() : no_rings;
+            const bool stills = !experiment.IsRotationIndexing();
+            if (stills)
+                bragg_engine->RadialBackgroundFromRings(ring_bkg, integration.GetQBinCount(),
+                                                        integration.Settings().GetLowQ_recipA(),
+                                                        integration.Settings().GetQSpacing_recipA());
+            if (!bragg_engine->IsBackgroundRadialAuto() || ring_bkg.empty())
                 return;
             output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore(
                     ring_bkg, integration.GetQBinCount(), integration.Settings(),
                     spot_finding_settings.ice_ring_width_Q_recipA);
-            bragg_engine->BackgroundRadial(*output.ice_ring_score
-                                           >= experiment.GetScalingSettings().GetIceMinScore());
+            if (stills)
+                bragg_engine->BackgroundRadial(!bragg_engine->GetRadialBackground().empty());
+            else
+                bragg_engine->BackgroundRadial(*output.ice_ring_score
+                                               >= experiment.GetScalingSettings().GetIceMinScore());
         };
 
         // A missing min-pix (std::nullopt) means "choose it per image". This applies only to the stills

image_analysis/bragg_integration/BackgroundAperture.cpp

@@ -0,0 +1,122 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "BackgroundAperture.h"
+
+#include <algorithm>
+#include <cmath>
+
+std::vector<float> background_aperture::CandidateR3(float r2, float max_r3) {
+    std::vector<float> out;
+    for (float r3 = r2 + 2.0f; r3 <= max_r3 + 1e-3f; r3 += 1.0f)
+        out.push_back(r3);
+    return out;
+}
+
+void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width, int height,
+                                       const std::vector<DiffractionSpot> &spots, float r1, float r2,
+                                       const std::vector<float> &r3, float clip_nsigma, float gain,
+                                       Tally &tally) {
+    if (r3.empty())
+        return;
+    tally.sum_d2.resize(r3.size(), 0.0);
+    const float r_max = r3.back();
+    const int half = static_cast<int>(std::ceil(r_max));
+    const int side = 2 * half + 1;
+
+    // A grid of centres, at least one aperture apart, kept clear of every found spot.
+    const double area = static_cast<double>(width) * height;
+    const int step = std::max(side, static_cast<int>(std::ceil(std::sqrt(area / MAX_APERTURES_PER_FRAME))));
+    const float clear = r_max + SPOT_CLEARANCE_PX;
+    std::vector<int32_t> cx, cy;
+    for (int y = half; y + half < height; y += step)
+        for (int x = half; x + half < width; x += step) {
+            bool near_spot = false;
+            for (const auto &s : spots) {
+                const Coord c = s.RawCoord();
+                if ((c.x - x) * (c.x - x) + (c.y - y) * (c.y - y) < clear * clear) {
+                    near_spot = true;
+                    break;
+                }
+            }
+            if (!near_spot) {
+                cx.push_back(x);
+                cy.push_back(y);
+            }
+        }
+    if (cx.empty())
+        return;
+    const int32_t *tiles = image.GatherTiles(cx, cy, width, half);
+
+    std::vector<double> ring_sum(r3.size());
+    std::vector<int> ring_n(r3.size());
+    std::vector<int32_t> ring_px;
+    std::vector<float> ring_d2;
+    for (size_t a = 0; a < cx.size(); ++a) {
+        const int32_t *t = tiles + a * static_cast<size_t>(side) * side;
+        double disk_sum = 0.0;
+        int disk_n = 0;
+        bool readable = true;
+        ring_px.clear();
+        ring_d2.clear();
+        for (int dy = -half; dy <= half && readable; ++dy)
+            for (int dx = -half; dx <= half; ++dx) {
+                const float d2 = static_cast<float>(dx * dx + dy * dy);
+                if (d2 >= r_max * r_max)
+                    continue;
+                const int32_t v = t[(dy + half) * side + dx + half];
+                if (v == INT32_MIN || v == INT32_MAX) {
+                    readable = false;
+                    break;
+                }
+                if (d2 < r1 * r1) {
+                    disk_sum += v;
+                    ++disk_n;
+                } else if (d2 >= r2 * r2) {
+                    ring_px.push_back(v);
+                    ring_d2.push_back(d2);
+                }
+            }
+        if (!readable || disk_n == 0)
+            continue;
+        // Every candidate ring as the integrator reads it: its mean, then one high-side clip pass.
+        bool measured = true;
+        for (size_t k = 0; k < r3.size() && measured; ++k) {
+            double sum = 0.0;
+            int n = 0;
+            for (size_t p = 0; p < ring_px.size(); ++p)
+                if (ring_d2[p] < r3[k] * r3[k]) {
+                    sum += ring_px[p];
+                    ++n;
+                }
+            if (n == 0) {
+                measured = false;
+                break;
+            }
+            const double mean = sum / n;
+            const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, 1.0));
+            ring_sum[k] = 0.0;
+            ring_n[k] = 0;
+            for (size_t p = 0; p < ring_px.size(); ++p)
+                if (ring_d2[p] < r3[k] * r3[k] && (clip_nsigma <= 0.0f || ring_px[p] <= thr)) {
+                    ring_sum[k] += ring_px[p];
+                    ++ring_n[k];
+                }
+            measured = ring_n[k] > 0;
+        }
+        if (!measured)
+            continue;
+        const double disk = disk_sum / disk_n;
+        for (size_t k = 0; k < r3.size(); ++k) {
+            const double d = disk - ring_sum[k] / ring_n[k];
+            tally.sum_d2[k] += d * d;
+        }
+        ++tally.apertures;
+    }
+}
+
+int background_aperture::Best(const Tally &tally) {
+    if (tally.apertures == 0 || tally.sum_d2.empty())
+        return -1;
+    return static_cast<int>(std::min_element(tally.sum_d2.begin(), tally.sum_d2.end()) - tally.sum_d2.begin());
+}

image_analysis/bragg_integration/BackgroundAperture.h

@@ -0,0 +1,55 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+#include "../../common/DiffractionSpot.h"
+#include "../image_preprocessing/ImagePreprocessorBuffer.h"
+
+// How far out the background ring should reach, measured on apertures with nothing in them.
+//
+// The ring mean stands in for the background under the signal disk. A wider ring averages more pixels,
+// so its own noise falls; but it also reaches further from the reflection, so whatever the background
+// does between the disk and the ring - a curved solvent ring, the tile structure of a fibre-optic CCD,
+// the tails of reflections nobody predicted - enters the intensity as a bias. Which of the two wins is
+// a property of the detector and the sample, not a constant: on a pixel detector counting a fraction of
+// a photon per pixel the noise wins and the widest ring is best, on a binned CCD reading hundreds of
+// ADU a pixel the structure wins and the narrowest is.
+//
+// So it is measured. Apertures are laid on a grid over pre-scan frames, away from every spot found; on
+// each, the mean of the r1 disk and the mean of the r2..r3 ring are taken for every candidate r3, and
+// the run keeps the r3 whose ring predicts the disk best - the smallest mean square of the difference.
+// With nothing in the aperture the difference is exactly the error the ring makes in the background of
+// a reflection: its noise plus its bias. The disk's own noise is the same for every candidate, so it
+// does not move the choice. The ring is clipped as the integrator clips it (high side, at
+// mean + nsigma * sqrt(gain * mean)).
+namespace background_aperture {
+
+// The squared disk-minus-ring differences, summed per candidate r3, over every aperture measured.
+struct Tally {
+    std::vector<double> sum_d2;
+    size_t apertures = 0;
+};
+
+// The candidate outer radii: from r2 + 2 up to max_r3, one pixel apart.
+[[nodiscard]] std::vector<float> CandidateR3(float r2, float max_r3);
+
+// The apertures of one frame, appended to `tally`. Centres lie on a grid over the frame, and an aperture
+// is used only where every pixel out to the largest candidate is readable and no found spot is within
+// that radius plus SPOT_CLEARANCE_PX of its centre.
+void MeasureFrame(ImagePreprocessorBuffer &image, int width, int height,
+                  const std::vector<DiffractionSpot> &spots, float r1, float r2,
+                  const std::vector<float> &r3, float clip_nsigma, float gain, Tally &tally);
+
+// The candidate with the smallest mean square difference, or -1 where no aperture was measured.
+[[nodiscard]] int Best(const Tally &tally);
+
+constexpr float SPOT_CLEARANCE_PX = 4.0f;
+// At most this many apertures a frame: a 16 Mpx detector would otherwise lay twenty thousand.
+constexpr size_t MAX_APERTURES_PER_FRAME = 2000;
+
+}   // namespace background_aperture

image_analysis/bragg_integration/BraggIntegrationEngine.cpp

@@ -218,6 +218,50 @@ void BraggIntegrationEngine::BuildRadialKernel(float grow) {
         k_diff.push_back(static_cast<float>(hist_ann[k] / sa - hist_disk[k] / sd));
 }
 
+void BraggIntegrationEngine::RadialBackgroundFromRings(const std::vector<float> &ring_bkg, int32_t q_bins,
+                                                       float low_q, float dq) {
+    radial_curve.clear();
+    if (q_bins <= 0 || !(dq > 0.0f))
+        return;
+    // The rings averaged over their azimuthal bins, where the profile has any.
+    std::vector<double> sum(q_bins, 0.0);
+    std::vector<int> cnt(q_bins, 0);
+    for (size_t i = 0; i < ring_bkg.size(); ++i)
+        if (std::isfinite(ring_bkg[i])) {
+            sum[i % q_bins] += ring_bkg[i];
+            ++cnt[i % q_bins];
+        }
+    std::vector<double> q, b;
+    for (int i = 0; i < q_bins; ++i)
+        if (cnt[i] > 0) {
+            q.push_back(low_q + (i + 0.5) * dq);
+            b.push_back(sum[i] / cnt[i]);
+        }
+    if (q.size() < 2)
+        return;
+
+    // The same bins as the engines' own curve: one per pixel of distance from the beam, out to the far
+    // corner. Each is read off the rings by linear interpolation in q; beyond the first and the last
+    // ring the end value is held, which has no curvature to correct.
+    const int n = static_cast<int>(std::ceil(r_max)) + 2;
+    const double two_k = 4.0 * PI / geom.GetWavelength_A();
+    radial_curve.resize(n);
+    size_t j = 0;
+    for (int r = 0; r < n; ++r) {
+        const double qr = two_k * std::sin(0.5 * std::atan(r / F_px));
+        while (j + 1 < q.size() && q[j + 1] < qr)
+            ++j;
+        double v;
+        if (qr <= q.front())
+            v = b.front();
+        else if (j + 1 >= q.size())
+            v = b.back();
+        else
+            v = b[j] + (b[j + 1] - b[j]) * (qr - q[j]) / (q[j + 1] - q[j]);
+        radial_curve[r] = static_cast<float>(v);
+    }
+}
+
 std::vector<BraggPanel> BraggIntegrationEngine::Panels(const std::vector<Reflection> &predicted,
                                                        size_t npredicted) const {
     std::vector<BraggPanel> out(npredicted);

image_analysis/bragg_integration/BraggIntegrationEngine.h

@@ -229,6 +229,10 @@ protected:
     int   k_len = 0;                 // entries per kernel
     int   n_kern = 1;                // kernels in the table (1 = circular stencil)
     std::vector<float> k_diff;       // n_kern * k_len, annulus-minus-disk weight per radial offset
+    // The radial background B the kernels are applied to, one value per pixel of distance from the
+    // beam, where the caller handed one over (RadialBackgroundFromRings). Empty: the curve is
+    // accumulated from the clipped ring pixels of the predicted reflections themselves.
+    std::vector<float> radial_curve;
 
     // One radial-offset kernel for a ring grown by `grow` px, appended to k_diff. Kept out of line
     // so the circular and elongated cases cannot drift apart. The signal disk does not change with
@@ -277,6 +281,16 @@ public:
     void BackgroundRadial(bool on) { bkg_radial = on && bkg_radial_built; }
     [[nodiscard]] bool IsBackgroundRadialAuto() const { return bkg_radial_auto; }
 
+    // The radial background curve for the images that follow, from a peak-excluded background per
+    // resolution ring of the whole detector (the adaptive spot finder's GetRingBackground): q bin i
+    // holds the ring at low_q + (i + 0.5) * dq, entries past the first q_bins are further azimuthal
+    // bins and are averaged in, and a ring with no value (NAN) is bridged by its neighbours. On a
+    // serial frame the predicted reflections' own ring pixels are a sparse sample of the same curve -
+    // a few dozen reflections, against every pixel of the frame here. An empty profile goes back to
+    // accumulating the curve from the reflections.
+    void RadialBackgroundFromRings(const std::vector<float> &ring_bkg, int32_t q_bins, float low_q, float dq);
+    [[nodiscard]] const std::vector<float> &GetRadialBackground() const { return radial_curve; }
+
     // 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.

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -311,7 +311,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                 if (px <= thr) {
                     clip.sum += px;
                     ++clip.n;
-                    if (bkg_radial) {
+                    if (bkg_radial && radial_curve.empty()) {
                         // 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);
@@ -335,7 +335,7 @@ BraggIntegrationEngineCPU::PassA BraggIntegrationEngineCPU::StartPassA(size_t np
     const float by0 = static_cast<float>(plane_bounds.y0), by1 = static_cast<float>(plane_bounds.y1);
     a.n_rad = static_cast<int>(std::ceil(std::hypot(std::max<double>(beam_x - bx0, bx1 - beam_x),
                                                     std::max<double>(beam_y - by0, by1 - beam_y)))) + 2;
-    if (bkg_radial) {
+    if (bkg_radial && radial_curve.empty()) {
         a.rad_sum.assign(a.n_rad, 0);
         a.rad_cnt.assign(a.n_rad, 0);
     }
@@ -612,7 +612,11 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
             double corr = 0.0;
             for (int k = 0; k < k_len; ++k) {
                 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;
+                double v;
+                if (!radial_curve.empty())
+                    v = radial_curve[b];
+                else
+                    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

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -168,7 +168,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
                        float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o,
                        uint8_t *hasobs_o, uint8_t *overloaded_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, int accumulate_rad,
                        unsigned long long *counts) {
     const int i = blockIdx.x;
     if (i >= n) return;
@@ -360,7 +360,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
             if (!valid(px)) continue;
             if ((long long) px <= s_thr_i) {
                 c_l += px; ++cn_l;
-                if (n_rad > 0) {
+                if (accumulate_rad) {
                     // The radial curve is binned on the TRUE detector radius, so this is the
                     // unshrunk radial projection - no per-pixel sqrt.
                     const int idx = (int) lroundf(s_r0 + d.rad) - s_radbase;
@@ -375,7 +375,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
     }
     __syncthreads();
 
-    if (n_rad > 0) {
+    if (accumulate_rad) {
         for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) {
             if (s_radn[t] == 0) continue;
             const int b = min(max(s_radbase + t, 0), n_rad - 1);
@@ -602,7 +602,7 @@ __global__ void build_profiles(const unsigned long long *shell_grid, const unsig
 //     correction because the kernel weights sum to zero. The kernel depends on how far this
 //     reflection's ring was elongated, so k_diff is a table and kbin_a picks the row.
 //     Mirrors BraggIntegrationEngineCPU. ---
-__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, int n_rad,
+__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,
@@ -615,7 +615,8 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad
     for (int k = 0; k < k_len; ++k) {
         int b = rbin_a[i] + k - k_off;
         b = min(max(b, 0), n_rad - 1);
-        const float v = rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg;
+        const float v = curve ? curve[b]
+                              : rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg;
         corr += kern[k] * v;
     }
     const float bkg_new = bkg - corr;
@@ -906,6 +907,7 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment
         n_rad = static_cast<int>(std::ceil(r_max)) + 2;
         d_rad_sum = CudaDevicePtr<unsigned long long>(n_rad);
         d_rad_cnt = CudaDevicePtr<int>(n_rad);
+        d_rad_curve = CudaDevicePtr<float>(n_rad);
         d_k_diff = CudaDevicePtr<float>(k_diff.size());
         cuda_err(cudaMemcpyAsync(d_k_diff, k_diff.data(), sizeof(float) * k_diff.size(),
                                  cudaMemcpyHostToDevice, *this->stream));
@@ -1015,7 +1017,11 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
 
     // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under
     // the auto mode they are allocated for every image and used for the ones the ice score selects.
-    const int rad_n = bkg_radial ? n_rad : 0;
+    // Where the caller handed over a curve (RadialBackgroundFromRings) it is read instead, and the
+    // box sums accumulate nothing.
+    const bool radial_on = bkg_radial && n_rad > 0;
+    const bool radial_given = radial_on && !radial_curve.empty();
+    const int rad_n = radial_on ? n_rad : 0;
 
     // The mask and the owner map can be taken back either by clearing the whole frame or by revisiting
     // the boxes that were marked, and which is cheaper depends on the detector: on 18 Mpx the marks are
@@ -1034,7 +1040,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
         if (p.overlap)
             cuda_err(cudaMemsetAsync(d_owner, 0xff, sizeof(uint32_t) * npixel, *stream));   // BRAGG_OWNER_NONE
     }
-    if (rad_n > 0) {
+    if (rad_n > 0 && !radial_given) {
         cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(unsigned long long) * n_rad, *stream));
         cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream));
     }
@@ -1050,14 +1056,17 @@ 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_overloaded, d_invd2,
                                        d_isum, d_ninner, d_rbin, d_kbin,
-                                       d_rad_sum, d_rad_cnt, rad_n, d_counts);
+                                       d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given, d_counts);
     cuda_err(cudaGetLastError());
 
     // Correct the flat annulus background for the curvature of the radial background before anything
     // downstream (profile fit, variance) reads it.
-    if (rad_n > 0) {
+    if (radial_given)
+        cuda_err(cudaMemcpyAsync(d_rad_curve, radial_curve.data(), sizeof(float) * n_rad,
+                                 cudaMemcpyHostToDevice, *stream));
+    if (radial_on) {
         radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>(
-            d_rad_sum, d_rad_cnt, rad_n, d_k_diff, k_len, k_off,
+            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);
         cuda_err(cudaGetLastError());
     }

image_analysis/bragg_integration/BraggIntegrationEngineGPU.h

@@ -54,6 +54,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
     CudaDevicePtr<unsigned long long> d_rad_sum;   // integer pixel sums, see boxsum
     CudaDevicePtr<float>   d_k_diff;
     CudaDevicePtr<int>     d_rad_cnt;
+    CudaDevicePtr<float>   d_rad_curve;              // the curve handed over (RadialBackgroundFromRings)
 
     // --- fixed-size device arrays ---
     // The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the

image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm

@@ -220,7 +220,9 @@ std::vector<Reflection> BraggIntegrationEngineMetal::Run(const ImagePreprocessor
         f.inv_claim = inv_claim;
         f.use_owner = use_owner;
         f.exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum;
-        f.n_rad = bkg_radial ? a.n_rad : 0;
+        // The curve is added up here only where no ring profile was handed over (RadialBackgroundFromRings).
+        const bool accumulate_radial = bkg_radial && radial_curve.empty();
+        f.n_rad = accumulate_radial ? a.n_rad : 0;
         f.static_mask_words = static_cast<uint32_t>(static_mask.size());
         f.n = static_cast<uint32_t>(npredicted);
 
@@ -285,12 +287,12 @@ std::vector<Reflection> BraggIntegrationEngineMetal::Run(const ImagePreprocessor
             }
         }
         if (any_clip) {
-            if (bkg_radial && dev->rad_sum == nil) {
+            if (accumulate_radial && dev->rad_sum == nil) {
                 dev->rad_sum = buffer(a.n_rad * 2 * sizeof(uint32_t));
                 dev->rad_cnt = buffer(a.n_rad * sizeof(int32_t));
             }
             cb = [dev->queue commandBuffer];
-            if (bkg_radial) {
+            if (accumulate_radial) {
                 blit = [cb blitCommandEncoder];
                 [blit fillBuffer:dev->rad_sum range:NSMakeRange(0, a.n_rad * 2 * sizeof(uint32_t)) value:0];
                 [blit fillBuffer:dev->rad_cnt range:NSMakeRange(0, a.n_rad * sizeof(int32_t)) value:0];
@@ -306,13 +308,13 @@ std::vector<Reflection> BraggIntegrationEngineMetal::Run(const ImagePreprocessor
             [enc setBuffer:dev->clips offset:0 atIndex:5];
             [enc setBuffer:dev->thr offset:0 atIndex:6];
             [enc setBuffer:dev->clip_out offset:0 atIndex:7];
-            [enc setBuffer:(bkg_radial ? dev->rad_sum : dev->clips) offset:0 atIndex:8];   // not written without the curve
-            [enc setBuffer:(bkg_radial ? dev->rad_cnt : dev->clips) offset:0 atIndex:9];
+            [enc setBuffer:(accumulate_radial ? dev->rad_sum : dev->clips) offset:0 atIndex:8];   // not written without the curve
+            [enc setBuffer:(accumulate_radial ? dev->rad_cnt : dev->clips) offset:0 atIndex:9];
             [enc setBuffer:dev->panels offset:0 atIndex:10];
             [enc dispatchThreads:grid threadsPerThreadgroup:group];
             [enc endEncoding];
             Wait(cb);
-            if (bkg_radial) {
+            if (accumulate_radial) {
                 const auto *w = static_cast<const uint32_t *>(dev->rad_sum.contents);
                 const auto *c = static_cast<const int32_t *>(dev->rad_cnt.contents);
                 for (int b = 0; b < a.n_rad; ++b) {

image_analysis/bragg_integration/CMakeLists.txt

@@ -18,6 +18,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC
         SpotWidth.h
         BackgroundGain.cpp
         BackgroundGain.h
+        BackgroundAperture.cpp
+        BackgroundAperture.h
 )
 
 TARGET_LINK_LIBRARIES(JFJochBraggIntegration JFJochImagePreprocessing JFJochCommon)

rugnux/Rugnux.h

@@ -45,6 +45,7 @@
 #include "../image_analysis/WriteReflections.h"   // OutputProvenance
 #include "../image_analysis/bragg_integration/SpotWidth.h"       // FluxCurve
 #include "../image_analysis/bragg_integration/SpotFootprint.h"   // FootprintSpot
+#include "../image_analysis/bragg_integration/BackgroundAperture.h"   // background_aperture::Tally
 #include "DiagnosticOutput.h"                                     // PerImagePlot
 
 class JFJochReader;
@@ -189,7 +190,9 @@ struct ProcessConfig {
     // aperture that owes nothing to the integrator (SpotWidth.h) - and sets r1 from it. r1 is the
     // aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius
     // a saturates at a^2/4, so at the shipped r1 = 4 a spot wider than sigma = 2 px is fitted with a
-    // profile the model cannot represent. Ignored where the radii were set by hand.
+    // profile the model cannot represent. On stills the signal disk is kept and the pre-scan measures
+    // the background ring's outer radius instead, on empty apertures (BackgroundAperture.h). Ignored
+    // where the radii were set by hand.
     bool adaptive_integration_radius = false;
 
     // The input's pixel values are ADU, not counted photons (SMV and marCCD: CCDs). The pre-scan then
@@ -1331,6 +1334,7 @@ class Rugnux {
             std::vector<FootprintSpot> footprint;
             std::vector<float> spot_q;
             size_t width_images = 0;
+            background_aperture::Tally apertures;
         };
         struct Defects {
             size_t pixels = 0;

rugnux/RugnuxPreScan.cpp

@@ -7,6 +7,7 @@
 #include "../image_analysis/indexing/SpindleCuspLoss.h"
 #include "../image_analysis/bragg_integration/SpotWidth.h"
 #include "../image_analysis/bragg_integration/BackgroundGain.h"
+#include "../image_analysis/bragg_integration/BackgroundAperture.h"
 #include "../image_analysis/bragg_integration/SpotFootprint.h"
 #include "../image_analysis/bragg_integration/BraggStencil.h"   // BRAGG_FOOTPRINT_NSIGMA
 #include "../image_analysis/hot_pixels/HotPixels.h"
@@ -366,10 +367,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         return w;
     };
 
+    // Stills: how far the background ring should reach, measured on empty apertures of the same frames
+    // (BackgroundAperture.h) at the signal disk and gap the run integrates with. The candidates run up
+    // to the ring the run starts from, so where the noise wins that ring is what is kept. The clip gain
+    // is read once the gain below has been measured.
+    const bool stills = !experiment_.IsRotationIndexing();
+    const BraggIntegrationSettings aperture_bis = prescan_x.GetBraggIntegrationSettings();
+    const std::vector<float> aperture_r3 = stills
+            ? background_aperture::CandidateR3(aperture_bis.GetR2(), aperture_bis.GetR3()) : std::vector<float>();
+    float aperture_gain = 1.0f;
+
     // The spots of one image, and what is taken off them.
     const auto find_spots = [&](PreScanWorker &w, const CompressedImage &image, int image_idx,
                                 bool for_width, std::vector<spot_width::FluxCurve> &curves,
-                                std::vector<float> &spot_q, std::vector<FootprintSpot> &footprint) {
+                                std::vector<float> &spot_q, std::vector<FootprintSpot> &footprint,
+                                background_aperture::Tally &apertures) {
         try {
             // As in the image loops: a frame the device cannot decode goes to the host decoder.
             ImageStatistics stats;
@@ -413,6 +425,11 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             MeasureSpotsOnTiles(*w.preprocessed, static_cast<int>(prescan_x.GetXPixelsNumConv()),
                                 static_cast<int>(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(),
                                 prescan_x.GetModules(), spots, curves, footprint);
+        if (for_width && !aperture_r3.empty())
+            background_aperture::MeasureFrame(*w.preprocessed, static_cast<int>(prescan_x.GetXPixelsNumConv()),
+                                              static_cast<int>(prescan_x.GetYPixelsNumConv()), spots,
+                                              aperture_bis.GetR1(), aperture_bis.GetR2(), aperture_r3,
+                                              aperture_bis.GetBackgroundClipNSigma(), aperture_gain, apertures);
     };
 
     std::vector<float> prescan_spot_q;
@@ -422,6 +439,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     std::vector<spot_width::FluxCurve> shape_curves;
     // Every spot's measured widths, for the spot footprint (SpotFootprint.h).
     std::vector<FootprintSpot> footprint_spots;
+    // The empty apertures of the frames the width is measured on (stills).
+    background_aperture::Tally aperture_tally;
     // Images the width ended up being measured on, for the log: the tiers below stop as soon as the
     // answer has settled, so this is a property of the crystal and worth reporting.
     size_t width_images = 0;
@@ -595,6 +614,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         std::vector<std::vector<spot_width::FluxCurve>> curves_of(ordinals.size());
         std::vector<std::vector<float>> spot_q_of(ordinals.size());
         std::vector<std::vector<FootprintSpot>> footprint_of(ordinals.size());
+        std::vector<background_aperture::Tally> apertures_of(ordinals.size());
         // Frames read for the bandwidth alone, past the tier the width settled on. Their spots stay
         // out of the powder measurement, which is taken on the frames it always was.
         std::vector<char> shape_only(ordinals.size(), 0);
@@ -641,7 +661,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
                         const size_t i = visit[v];
                         const int image_idx = start_image + ordinals[i] * config_.stride;
                         if (const CompressedImage *image = frame(i, w.raw_image))
-                            find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i]);
+                            find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i],
+                                       apertures_of[i]);
                         {
                             std::lock_guard lock(m);
                             done_in_tier[tier_of(v)]++;
@@ -708,11 +729,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         for (size_t i = 0; i < ordinals.size(); i++)
             if (!shape_only[i])
                 prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end());
+        // In sample order, over the frames the width was taken on.
+        aperture_tally.sum_d2.assign(aperture_r3.size(), 0.0);
+        for (size_t i = 0; i < ordinals.size(); i++)
+            if (position[i] < width_images && !apertures_of[i].sum_d2.empty()) {
+                for (size_t k = 0; k < aperture_r3.size(); k++)
+                    aperture_tally.sum_d2[k] += apertures_of[i].sum_d2[k];
+                aperture_tally.apertures += apertures_of[i].apertures;
+            }
         // In sample order, so the table does not depend on how the workers interleaved.
         for (const auto &f : footprint_of)
             footprint_spots.insert(footprint_spots.end(), f.begin(), f.end());
     };
     std::future<void> spots_measured;
+    // The gain is applied above, before the spots are measured, so the apertures clip as the run will.
+    aperture_gain = experiment_.GetBraggIntegrationSettings().GetBackgroundPixelGain();
     if (want_spots)
         spots_measured = std::async(std::launch::async, spot_measurement);
 
@@ -725,7 +756,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             if (!measured)
                 measured.emplace();
             measured->spots = PreScanMemo::Spots{width_curves, shape_curves, footprint_spots, prescan_spot_q,
-                                                 width_images};
+                                                 width_images, aperture_tally};
         } else if (want_width) {
             logger.Info("Spot width: the spots kept from an earlier pre-scan of this run on the same frames, "
                         "geometry, mask and spot finding");
@@ -735,6 +766,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             footprint_spots = k.footprint;
             prescan_spot_q = k.spot_q;
             width_images = k.width_images;
+            aperture_tally = k.apertures;
         }
 
         // The integration radius from the recorded spot width. r1 is the aperture the integrator learns
@@ -743,7 +775,33 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         // the model cannot represent. r2 follows r1, and r3 is taken so that the r2..r3 background ring
         // keeps the area it has at the shipped 4/6/13 - a ring that shrinks with the disk is what makes a
         // bare --integration-radius worse than the default it replaces.
-        if (want_width) {
+        //
+        // Stills keep their signal disk and take the ring's outer radius from the empty apertures
+        // instead (BackgroundAperture.h): the candidate whose ring predicted the background under the
+        // disk best. Where the noise wins that is the widest, which is the ring the run started from.
+        if (want_width && stills) {
+            spot_width_measured_ = true;
+            const int best = background_aperture::Best(aperture_tally);
+            if (best < 0) {
+                logger.Info("Background ring: no empty aperture measured on {} images, keeping r3={:.1f}",
+                            width_images, experiment_.GetBraggIntegrationSettings().GetR3());
+            } else {
+                std::string table;
+                for (size_t k = 0; k < aperture_r3.size(); k++)
+                    table += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k],
+                                         aperture_tally.sum_d2[k] / static_cast<double>(aperture_tally.apertures));
+                BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
+                bis.R3(aperture_r3[best]);
+                experiment_.ImportBraggIntegrationSettings(bis);
+                logger.Info("Background ring: on {} empty apertures from {} images the r{:.0f}..r3 ring predicts "
+                            "the r{:.0f} disk with a mean square error of (r3:value){} => integration radii "
+                            "r1={:.1f} r2={:.1f} r3={:.1f}", aperture_tally.apertures, width_images,
+                            bis.GetR2(), bis.GetR1(), table, bis.GetR1(), bis.GetR2(), bis.GetR3());
+                Note(fmt::format("background ring {:.0f}..{:.0f} px, measured on empty apertures",
+                                 bis.GetR2(), bis.GetR3()));
+            }
+        }
+        if (want_width && !stills) {
             spot_width_measured_ = true;
             const auto r80 = spot_width::R80AtReference(width_curves);
             if (!r80) {
@@ -824,8 +882,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         // The X-ray bandwidth, read off the shapes of the same spots (spot_width::EstimateBandwidth). A
         // value the file states, or --bandwidth, wins; otherwise a significant estimate is the bandwidth
         // of the run - prediction, partiality, the profile's radial width and the background stencil all
-        // take it - and anything short of significant leaves the beam monochromatic.
-        if (want_width) {
+        // take it - and anything short of significant leaves the beam monochromatic. Not on stills, where the
+        // estimate has not been shown to measure the beam (on the serial sets it reads not significant).
+        if (want_width && !stills) {
             const auto &det = experiment_.GetDetectorSetup();
             const double lambda = experiment_.GetWavelength_A();
             const auto estimate = spot_width::EstimateBandwidth(

rugnux/rugnux_cli.cpp

@@ -249,12 +249,12 @@ void print_usage() {
     std::cout << "  Integration" << std::endl;
     std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM). Default: the file's value, else measured in the pre-scan from the radial elongation of strong spots and used where significant, else 0; --bandwidth 0 forces a monochromatic beam" << std::endl;
     std::cout << "   --integration-radius <r>        Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
-    std::cout << "   --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default: on for rotation, off for stills). r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl;
+    std::cout << "   --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default on). On stills it measures instead how far the background ring should reach: on empty apertures of the pre-scan frames, the r3 (up to 14) whose ring predicts the background under the signal disk best. r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl;
     std::cout << "   --integration-stencil <k>       Push the r2..r3 background ring out by k times the beam's radial streak (bandwidth*Rpx), per reflection (default 0 = a fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially - the r1 signal box stays a circle. On a monochromatic beam the streak is zero and this does nothing" << std::endl;
     std::cout << "   --integration-high-resolution <num> High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl;
     std::cout << "   --max-hkl <n>                   Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl;
     std::cout << "   --background-clip <n>           High-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl;
-    std::cout << "   --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl;
+    std::cout << "   --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. On stills =auto applies it to every image, with the curve taken from the spot finder's peak-excluded background of the whole frame: a smooth solvent or carrier hump curves the background as much as ice does. Costs one short dot product per reflection and no extra pixel reads" << std::endl;
     std::cout << "   --background-trim <f>           Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl;
     std::cout << "   --overlap <txt>                 What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl;
     std::cout << "   --overlap-minpk <f>             Least fraction of a reflection's expected profile that must be usable for the reflection to be kept (default 0.75, XDS MINPK). Governs both ways part of a profile is lost: the fraction that must be READABLE - not masked, untrusted, in a detector gap or overloaded - in every profile mode, and, under --overlap reject, the fraction that must be cleanly the reflection's own. --integrator boxsum has no profile to renormalise with, so there a signal disk with any unreadable pixel is still discarded outright and the reject fraction is by disk AREA, which cuts harder" << std::endl;
@@ -3308,15 +3308,15 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     } else if (!rotation_indexing) {
         // An explicit --integration-radius always wins over the stills default.
         ApplyRugnuxIntegrationRadii(experiment, rotation_indexing);
-        logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=14.0 (override with --integration-radius)");
+        logger.Info("Stills integration radii start at r1=6.0 r2=8.0 r3=14.0, the ring's reach measured in the pre-scan (override with --integration-radius)");
     }
 
-    // Adaptive signal radius: measured in the pre-scan and applied there (Rugnux::PreScan). Radii the
-    // user set by hand always win, and stills already default to a wide box, so the rule is on for
-    // rotation only.
+    // Adaptive radii: measured in the pre-scan and applied there (Rugnux::PreScan) - on rotation the
+    // signal radius from the spot width, on stills the background ring's reach from empty apertures.
+    // Radii the user set by hand always win.
     const bool adaptive_integration_radius =
             !integration_radius_arg
-            && adaptive_integration_radius_arg.value_or(rotation_indexing);
+            && adaptive_integration_radius_arg.value_or(true);
 
     if (forced_prediction_mosaicity_arg) {
         BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
@@ -3427,7 +3427,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         const auto radial = experiment.GetBraggIntegrationSettings().GetBackgroundRadialCorrection();
         logger.Info("Background ring: radial curvature correction {}",
                     radial.has_value() ? (*radial ? "on" : "off")
-                                       : "auto (per image, where the ice score sees a smooth ring)");
+                                       : rotation_indexing
+                                             ? "auto (per image, where the ice score sees a smooth ring)"
+                                             : "auto (every image, from the radial background of the whole frame)");
     }
 
     // Detection defaults differ by workflow (RugnuxDefaultSpotFinding); each flag always wins.

tests/BackgroundApertureTest.cpp

@@ -0,0 +1,72 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <random>
+
+#include "../image_analysis/bragg_integration/BackgroundAperture.h"
+#include "../common/JFJochMath.h"
+
+namespace {
+
+constexpr int W = 600, H = 600;
+
+// A background of `level` counts a pixel with Poisson noise, modulated by `ripple` (relative) with a
+// period of 20 px along x - the structure a ring further from the disk reads wrong.
+ImagePreprocessorBuffer Frame(double level, double ripple, uint32_t seed) {
+    std::mt19937 rng(seed);
+    ImagePreprocessorBuffer img(static_cast<size_t>(W) * H);
+    for (int y = 0; y < H; y++)
+        for (int x = 0; x < W; x++) {
+            std::poisson_distribution<int> p(level * (1.0 + ripple * std::sin(2.0 * PI * x / 20.0)));
+            img[static_cast<size_t>(y) * W + x] = p(rng);
+        }
+    return img;
+}
+
+int BestR3(double level, double ripple) {
+    const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f);
+    background_aperture::Tally tally;
+    for (uint32_t seed = 1; seed <= 4; seed++) {
+        auto img = Frame(level, ripple, seed);
+        background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, tally);
+    }
+    REQUIRE(tally.apertures > 100);
+    return static_cast<int>(r3[background_aperture::Best(tally)]);
+}
+
+}   // namespace
+
+TEST_CASE("BackgroundAperture_CandidateR3", "[Integration]") {
+    const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f);
+    REQUIRE(r3.size() == 5);
+    CHECK(r3.front() == 10.0f);
+    CHECK(r3.back() == 14.0f);
+}
+
+// A few counts a pixel on a flat background: the noise wins, and the widest ring is kept.
+TEST_CASE("BackgroundAperture_FlatNoisyKeepsWidestRing", "[Integration]") {
+    CHECK(BestR3(0.5, 0.0) == 14);
+}
+
+// Thousands of counts a pixel on a background that varies over the stencil's own scale: the structure
+// wins, and the narrowest ring is kept.
+TEST_CASE("BackgroundAperture_StructuredBrightTakesNarrowestRing", "[Integration]") {
+    CHECK(BestR3(2000.0, 0.05) == 10);
+}
+
+// An aperture over an unreadable pixel is not used.
+TEST_CASE("BackgroundAperture_SkipsUnreadablePixels", "[Integration]") {
+    const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f);
+    auto img = Frame(1.0, 0.0, 7);
+    background_aperture::Tally all;
+    background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, all);
+    for (size_t i = 0; i < static_cast<size_t>(W) * H; i += 97)
+        img[i] = INT32_MIN;
+    background_aperture::Tally masked;
+    background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, masked);
+    CHECK(masked.apertures == 0);
+    CHECK(all.apertures > 0);
+}

tests/BraggIntegrationEngineCPUTest.cpp

@@ -10,6 +10,7 @@
 #include "../common/BraggIntegrationSettings.h"
 #include "../common/DetectorSetup.h"
 #include "../common/DiffractionExperiment.h"
+#include "../common/JFJochMath.h"
 #include "../common/Reflection.h"
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h"
 #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
@@ -148,3 +149,57 @@ TEST_CASE("BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped", "[Integr
     CHECK(saturated[0].overloaded);
     CHECK(saturated[0].clipped);
 }
+
+// An empty aperture on a background that is a curved function of the radius alone: the flat ring mean
+// reads it wrong by the curvature, and the radial correction with the curve taken from a per-ring
+// profile of the whole frame (as the spot finder hands it over on stills) takes that back out.
+TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[Integration][portable]") {
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f);
+    const float bx = 400.0f, by = 400.0f;
+    const auto background = [](double r) { return 50.0 + 0.05 * (r - 300.0) * (r - 300.0); };
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+    ImagePreprocessorBuffer image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        image[i] = static_cast<int32_t>(std::lround(
+                background(std::hypot(static_cast<double>(i % width) - bx, static_cast<double>(i / width) - by))));
+
+    // The profile per q bin, as the spot finder's rings would hold it.
+    const float low_q = 0.05f, dq = 0.01f;
+    const double F_px = 100.0 / 0.075, lambda = 1.0;
+    std::vector<float> rings(500);
+    for (size_t i = 0; i < rings.size(); ++i) {
+        const double q = low_q + (i + 0.5) * dq;
+        rings[i] = static_cast<float>(background(F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI)))));
+    }
+    rings[100] = NAN;   // a ring with no value is bridged
+
+    Reflection r{};
+    r.h = 1; r.k = 0; r.l = 0;
+    r.predicted_x = bx + 300.0f;
+    r.predicted_y = by;
+    r.d = 2.0f;
+    r.prescaling_corr = 1.0f;
+    r.partiality = 1.0f;
+    const std::vector<Reflection> predicted{r};
+
+    auto integrate = [&](bool correct) {
+        experiment.ImportBraggIntegrationSettings(
+                BraggIntegrationSettings().Integrator(IntegratorMode::BoxSum).BackgroundRadialCorrection(correct));
+        BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment));
+        engine.RadialBackgroundFromRings(rings, static_cast<int32_t>(rings.size()), low_q, dq);
+        REQUIRE(!engine.GetRadialBackground().empty());
+        const auto out = engine.Run(image, predicted, predicted.size(), 0);
+        REQUIRE(out.size() == 1);
+        return out[0].I;
+    };
+    // The curvature is along the radius only, so the error is n_inner x 0.05 x (<u^2>_ring - <u^2>_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);
+
+    // An empty profile hands the curve back to the reflections' own rings.
+    BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment));
+    engine.RadialBackgroundFromRings({}, 500, low_q, dq);
+    CHECK(engine.GetRadialBackground().empty());
+}

tests/BraggIntegrationEngineGPUTest.cpp

@@ -15,6 +15,7 @@
 #include "../common/BraggIntegrationSettings.h"
 #include "../common/DetectorSetup.h"
 #include "../common/DiffractionExperiment.h"
+#include "../common/JFJochMath.h"
 #include "../common/Reflection.h"
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h"
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineGPU.h"
@@ -413,6 +414,60 @@ TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") {
     }
 }
 
+// The radial correction with its curve handed over from a per-ring profile (RadialBackgroundFromRings), as
+// on stills: the GPU reads the same curve and corrects the same empty aperture as the CPU.
+TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integration]") {
+    if (get_gpu_count() == 0) {
+        WARN("No CUDA GPU present. Skipping");
+        return;
+    }
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f);
+    const float bx = 400.0f, by = 400.0f;
+    const auto background = [](double r) { return 50.0 + 0.05 * (r - 300.0) * (r - 300.0); };
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+    ImagePreprocessorBuffer cpu_image(npixel);
+    ImagePreprocessorBufferGPU gpu_image(npixel);
+    for (size_t i = 0; i < npixel; ++i) {
+        cpu_image[i] = static_cast<int32_t>(std::lround(
+                background(std::hypot(static_cast<double>(i % width) - bx, static_cast<double>(i / width) - by))));
+        gpu_image[i] = cpu_image[i];
+    }
+    auto stream = std::make_shared<CudaStream>();
+    REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(),
+                            npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess);
+
+    const float low_q = 0.05f, dq = 0.01f;
+    const double F_px = 100.0 / 0.075, lambda = 1.0;
+    std::vector<float> rings(500);
+    for (size_t i = 0; i < rings.size(); ++i) {
+        const double q = low_q + (i + 0.5) * dq;
+        rings[i] = static_cast<float>(background(F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI)))));
+    }
+
+    std::vector<Reflection> predicted;
+    for (int i = 0; i < 8; ++i)
+        predicted.push_back(MakeReflection(bx + 300.0f * std::cos(0.7f * i), by + 300.0f * std::sin(0.7f * i), 2.0f, i + 1));
+
+    for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
+        experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings().Integrator(mode).BackgroundRadialCorrection(true));
+        BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment));
+        BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment));
+        cpu.RadialBackgroundFromRings(rings, static_cast<int32_t>(rings.size()), low_q, dq);
+        gpu.RadialBackgroundFromRings(rings, static_cast<int32_t>(rings.size()), low_q, dq);
+        const auto out_cpu = cpu.Run(cpu_image, predicted, predicted.size(), 0);
+        const auto out_gpu = gpu.Run(gpu_image, predicted, predicted.size(), 0);
+        REQUIRE(out_cpu.size() == predicted.size());
+        REQUIRE(out_gpu.size() == out_cpu.size());
+        for (size_t i = 0; i < out_cpu.size(); ++i) {
+            INFO("mode " << static_cast<int>(mode) << " reflection " << i);
+            CHECK(out_gpu[i].bkg == Catch::Approx(out_cpu[i].bkg).epsilon(0.001).margin(0.05));
+            CHECK(std::fabs(out_gpu[i].I) < 10.0f);
+            CHECK(std::fabs(out_cpu[i].I) < 10.0f);
+        }
+    }
+}
+
 // Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU
 // for ProfileIntegrate2D). Run explicitly with:  ./jfjoch_test "[bragg_bench]"
 TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") {

tests/CMakeLists.txt

@@ -92,6 +92,7 @@ ADD_EXECUTABLE(jfjoch_test
         BraggIntegrationEngineCompressedImageTest.cpp
         BraggIntegrationEngineCPUTest.cpp
         BackgroundGainTest.cpp
+        BackgroundApertureTest.cpp
         BraggStencilTest.cpp
         SpotFootprintTest.cpp
         SpotTilesTest.cpp