Stills: the radial curvature correction decided per run on the same empty apertures

The first stills battery of e8145d8 lost on one Rayonix set (lcls242_s, dimple R_free 0.2302 ->
0.2323): there the whole-frame radial correction hurts. On 40,000 of its frames: battery default
0.2462, ring 12 with the correction 0.2444, ring 12 without it 0.2336, ring 14 with it 0.2467. Its
beam is strongly polarised and the azimuthally averaged ring background does not describe the local
background.

So the empty apertures now also carry the corrected difference: the frame's ring background
(RingBackground, factored out of the engine) averaged over the ring pixels minus over the disk pixels
- what the integrator's kernel takes off a reflection - is subtracted from each disk-minus-ring
difference. The ring is still chosen on the plain differences, so its reach does not rest on the
radial model; where the correction was left to auto it is then kept only if it makes that ring predict
the disk better. Decisions on 15-60 pre-scan frames:
  pal86 ring 10, correction on (240.6 -> 221.1)   lcls222 ring 10, on (6.34e4 -> 3.41e4)
  lcls242_s ring 12, off (277.2 -> 285.8)         alvra213_int16 ring 14, on (0.01424 -> 0.01422)
Choosing ring and correction jointly instead picked ring 14 with the correction on lcls222, where dimple
reads 0.204 against 0.196 at ring 10; a median instead of the mean square picked ring 14 with the
correction on lcls242_s, the worst arm there. The sequential rule is the one that agrees with dimple on
every set measured.

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

image_analysis/MXAnalysisWithoutFPGA.cpp

@@ -237,7 +237,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
         // 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
+        // every image that has one (rugnux's pre-scan turns it off for a run whose empty apertures say
+        // the radial model does not hold there). 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

image_analysis/bragg_integration/BackgroundAperture.cpp

@@ -6,6 +6,8 @@
 #include <algorithm>
 #include <cmath>
 
+#include "../../common/JFJochMath.h"
+
 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)
@@ -16,10 +18,12 @@ std::vector<float> background_aperture::CandidateR3(float r2, float max_r3) {
 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) {
+                                       const RingBackground &rings, const DiffractionGeometry &geometry,
+                                       const DetectorModules &modules, Tally &tally) {
     if (r3.empty())
         return;
     tally.sum_d2.resize(r3.size(), 0.0);
+    tally.sum_d2_radial.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;
@@ -48,45 +52,57 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width
         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;
+    struct Pixel {
+        int32_t v;
+        float d2;     // squared distance from the aperture centre
+        double b;     // the ring background at this pixel's resolution
+    };
+    std::vector<Pixel> disk, ring;
     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;
+        const int64_t panel = modules.ModuleOfImagePixel(cx[a], cy[a]);
         bool readable = true;
-        ring_px.clear();
-        ring_d2.clear();
+        disk.clear();
+        ring.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)
+                if (d2 >= r_max * r_max || (d2 >= r1 * r1 && d2 < r2 * r2))
                     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);
+                double b = 0.0;
+                if (!rings.empty()) {
+                    const auto [px, py] = modules.ImageToPlane(panel, cx[a] + dx, cy[a] + dy);
+                    b = rings.At(2.0 * PI / geometry.PxlToRes(static_cast<float>(px), static_cast<float>(py)));
                 }
+                (d2 < r1 * r1 ? disk : ring).push_back({v, d2, b});
             }
-        if (!readable || disk_n == 0)
+        if (!readable || disk.empty())
             continue;
-        // Every candidate ring as the integrator reads it: its mean, then one high-side clip pass.
+        double disk_v = 0.0, disk_b = 0.0;
+        for (const auto &p : disk) {
+            disk_v += p.v;
+            disk_b += p.b;
+        }
+        disk_v /= static_cast<double>(disk.size());
+        disk_b /= static_cast<double>(disk.size());
+
+        // Every candidate ring as the integrator reads it: the mean, then one high-side clip pass. The
+        // radial correction is taken over the ring's geometry, as the integrator's kernel is.
+        std::vector<double> d_plain(r3.size()), d_radial(r3.size());
         bool measured = true;
         for (size_t k = 0; k < r3.size() && measured; ++k) {
-            double sum = 0.0;
+            const float r3_sq = r3[k] * r3[k];
+            double sum = 0.0, sum_b = 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];
+            for (const auto &p : ring)
+                if (p.d2 < r3_sq) {
+                    sum += p.v;
+                    sum_b += p.b;
                     ++n;
                 }
             if (n == 0) {
@@ -95,28 +111,55 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width
             }
             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];
+            double kept = 0.0;
+            int n_kept = 0;
+            for (const auto &p : ring)
+                if (p.d2 < r3_sq && (clip_nsigma <= 0.0f || p.v <= thr)) {
+                    kept += p.v;
+                    ++n_kept;
                 }
-            measured = ring_n[k] > 0;
+            measured = n_kept > 0;
+            if (measured) {
+                d_plain[k] = disk_v - kept / n_kept;
+                d_radial[k] = d_plain[k] - (disk_b - sum_b / n);
+            }
         }
         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.sum_d2[k] += d_plain[k] * d_plain[k];
+            if (!rings.empty())
+                tally.sum_d2_radial[k] += d_radial[k] * d_radial[k];
         }
         ++tally.apertures;
+        if (!rings.empty())
+            ++tally.apertures_radial;
     }
 }
 
-int background_aperture::Best(const Tally &tally) {
-    if (tally.apertures == 0 || tally.sum_d2.empty())
+void background_aperture::Tally::Add(const Tally &other) {
+    sum_d2.resize(std::max(sum_d2.size(), other.sum_d2.size()), 0.0);
+    sum_d2_radial.resize(std::max(sum_d2_radial.size(), other.sum_d2_radial.size()), 0.0);
+    for (size_t k = 0; k < other.sum_d2.size(); ++k)
+        sum_d2[k] += other.sum_d2[k];
+    for (size_t k = 0; k < other.sum_d2_radial.size(); ++k)
+        sum_d2_radial[k] += other.sum_d2_radial[k];
+    apertures += other.apertures;
+    apertures_radial += other.apertures_radial;
+}
+
+std::vector<double> background_aperture::Tally::MeanSquare(bool radial) const {
+    const size_t n = radial ? apertures_radial : apertures;
+    std::vector<double> out;
+    if (n == 0)
+        return out;
+    for (const double v : radial ? sum_d2_radial : sum_d2)
+        out.push_back(v / static_cast<double>(n));
+    return out;
+}
+
+int background_aperture::Best(const std::vector<double> &mean_square) {
+    if (mean_square.empty())
         return -1;
-    return static_cast<int>(std::min_element(tally.sum_d2.begin(), tally.sum_d2.end()) - tally.sum_d2.begin());
+    return static_cast<int>(std::min_element(mean_square.begin(), mean_square.end()) - mean_square.begin());
 }

image_analysis/bragg_integration/BackgroundAperture.h

@@ -7,8 +7,11 @@
 #include <cstdint>
 #include <vector>
 
+#include "../../common/DetectorModules.h"
+#include "../../common/DiffractionGeometry.h"
 #include "../../common/DiffractionSpot.h"
 #include "../image_preprocessing/ImagePreprocessorBuffer.h"
+#include "RingBackground.h"
 
 // How far out the background ring should reach, measured on apertures with nothing in them.
 //
@@ -27,12 +30,24 @@
 // 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)).
+//
+// The same apertures also say whether the radial background-curvature correction helps: each difference
+// is taken a second time with the correction applied - the frame's own ring background (RingBackground)
+// averaged over the ring's pixels minus over the disk's, which is what the integrator's kernel takes off a
+// reflection. Where the background really is a function of the radius the corrected ring predicts the
+// disk better; where it is not - a strongly polarised beam, panels of different response - it predicts it
+// worse, and the run measures which. The ring is chosen on the plain differences first, so its reach does
+// not rest on the radial model; the correction is then kept where it improves that ring.
 namespace background_aperture {
 
-// The squared disk-minus-ring differences, summed per candidate r3, over every aperture measured.
+// The squared disk-minus-ring differences, summed per candidate r3 over every aperture measured: plain,
+// and with the radial correction (counted separately: a frame without a ring background adds none).
 struct Tally {
-    std::vector<double> sum_d2;
-    size_t apertures = 0;
+    std::vector<double> sum_d2, sum_d2_radial;
+    size_t apertures = 0, apertures_radial = 0;
+    void Add(const Tally &other);
+    // The mean square difference of each candidate, plain or corrected; empty where nothing was measured.
+    [[nodiscard]] std::vector<double> MeanSquare(bool radial) const;
 };
 
 // The candidate outer radii: from r2 + 2 up to max_r3, one pixel apart.
@@ -40,13 +55,16 @@ struct Tally {
 
 // 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.
+// that radius plus SPOT_CLEARANCE_PX of its centre. `rings` is this frame's ring background, or empty;
+// `geometry` and `modules` place the pixels on the plane to read it at their resolution.
 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);
+                  const std::vector<float> &r3, float clip_nsigma, float gain,
+                  const RingBackground &rings, const DiffractionGeometry &geometry,
+                  const DetectorModules &modules, Tally &tally);
 
-// The candidate with the smallest mean square difference, or -1 where no aperture was measured.
-[[nodiscard]] int Best(const Tally &tally);
+// The candidate with the smallest mean square difference, or -1 where nothing was measured.
+[[nodiscard]] int Best(const std::vector<double> &mean_square);
 
 constexpr float SPOT_CLEARANCE_PX = 4.0f;
 // At most this many apertures a frame: a 16 Mpx detector would otherwise lay twenty thousand.

image_analysis/bragg_integration/BraggIntegrationEngine.cpp

@@ -13,6 +13,7 @@
 
 #include "../../common/JFJochMath.h"   // PI (M_PI is not standard, and MSVC does not define it)
 #include "../SensorAbsorption.h"
+#include "RingBackground.h"
 
 namespace {
 
@@ -221,45 +222,16 @@ void BraggIntegrationEngine::BuildRadialKernel(float grow) {
 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))
+    const RingBackground rings(ring_bkg, q_bins, low_q, dq);
+    if (rings.empty())
         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.
+    // corner.
     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);
-    }
+    for (int r = 0; r < n; ++r)
+        radial_curve[r] = static_cast<float>(rings.At(two_k * std::sin(0.5 * std::atan(r / F_px))));
 }
 
 std::vector<BraggPanel> BraggIntegrationEngine::Panels(const std::vector<Reflection> &predicted,

image_analysis/bragg_integration/CMakeLists.txt

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

image_analysis/bragg_integration/RingBackground.cpp

@@ -0,0 +1,33 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RingBackground.h"
+
+#include <algorithm>
+#include <cmath>
+
+RingBackground::RingBackground(const std::vector<float> &ring_bkg, int32_t q_bins, float low_q, float dq) {
+    if (q_bins <= 0 || !(dq > 0.0f))
+        return;
+    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];
+        }
+    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]);
+        }
+}
+
+double RingBackground::At(double q_recipA) const {
+    if (q_recipA <= q.front())
+        return b.front();
+    if (q_recipA >= q.back())
+        return b.back();
+    const size_t j = std::upper_bound(q.begin(), q.end(), q_recipA) - q.begin() - 1;
+    return b[j] + (b[j + 1] - b[j]) * (q_recipA - q[j]) / (q[j + 1] - q[j]);
+}

image_analysis/bragg_integration/RingBackground.h

@@ -0,0 +1,20 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstdint>
+#include <vector>
+
+// A frame's background as a function of q, from a peak-excluded background per resolution ring of the
+// whole detector (the adaptive spot finder's GetRingBackground): ring i holds 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. Between rings the value is interpolated linearly; beyond
+// the first and the last ring the end value is held, which has no curvature.
+struct RingBackground {
+    std::vector<double> q, b;
+
+    RingBackground(const std::vector<float> &ring_bkg, int32_t q_bins, float low_q, float dq);
+    [[nodiscard]] bool empty() const { return q.size() < 2; }
+    [[nodiscard]] double At(double q_recipA) const;
+};

rugnux/RugnuxPreScan.cpp

@@ -8,6 +8,7 @@
 #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/RingBackground.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"
@@ -425,11 +426,21 @@ 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())
+        if (for_width && !aperture_r3.empty()) {
+            // The frame's ring background, as the integration of this frame will read it - from the
+            // adaptive finder, which has it whole-frame and peak-excluded.
+            const std::vector<float> no_rings;
+            const RingBackground rings(config_.spot_finding.adaptive_threshold ? w.spot_finder->GetRingBackground()
+                                                                               : no_rings,
+                                       prescan_mapping->GetQBinCount(),
+                                       prescan_mapping->Settings().GetLowQ_recipA(),
+                                       prescan_mapping->Settings().GetQSpacing_recipA());
             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);
+                                              aperture_bis.GetBackgroundClipNSigma(), aperture_gain, rings,
+                                              prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), apertures);
+        }
     };
 
     std::vector<float> prescan_spot_q;
@@ -730,13 +741,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             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;
-            }
+            if (position[i] < width_images)
+                aperture_tally.Add(apertures_of[i]);
         // 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());
@@ -779,26 +786,44 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         // 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.
+        // The same apertures decide the radial curvature correction, where it was left to auto: on
+        // where the corrected ring predicts the disk better.
         if (want_width && stills) {
             spot_width_measured_ = true;
-            const int best = background_aperture::Best(aperture_tally);
+            BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
+            const auto radial_setting = bis.GetBackgroundRadialCorrection();
+            // The ring is chosen on the plain differences, so its reach does not rest on the radial
+            // model; the correction is then kept, where it was left to auto, if it makes that ring predict
+            // the disk better. Read only where every frame gave a ring background.
+            const std::vector<double> plain = aperture_tally.MeanSquare(false);
+            const std::vector<double> radial = aperture_tally.apertures_radial == aperture_tally.apertures
+                                               ? aperture_tally.MeanSquare(true) : std::vector<double>();
+            const int best = background_aperture::Best(plain);
+            const bool use_radial = radial_setting.has_value()
+                                    ? *radial_setting
+                                    : best >= 0 && !radial.empty() && radial[best] < plain[best];
             if (best < 0) {
                 logger.Info("Background ring: no empty aperture measured on {} images, keeping r3={:.1f}",
-                            width_images, experiment_.GetBraggIntegrationSettings().GetR3());
+                            width_images, bis.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();
+                const auto table = [&](const std::vector<double> &m) {
+                    std::string t;
+                    for (size_t k = 0; k < m.size(); k++)
+                        t += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k], m[k]);
+                    return t;
+                };
                 bis.R3(aperture_r3[best]);
+                if (!radial_setting)
+                    bis.BackgroundRadialCorrection(use_radial ? std::nullopt : std::optional<bool>(false));
                 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()));
+                            "the r{:.0f} disk with a mean square error of (r3:value){}, and with the radial "
+                            "curvature correction{} => integration radii r1={:.1f} r2={:.1f} r3={:.1f}, radial "
+                            "correction {}", aperture_tally.apertures, width_images, bis.GetR2(), bis.GetR1(),
+                            table(plain), radial.empty() ? std::string(" (no ring background)") : table(radial),
+                            bis.GetR1(), bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off");
+                Note(fmt::format("background ring {:.0f}..{:.0f} px, radial correction {}, measured on empty "
+                                 "apertures", bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off"));
             }
         }
         if (want_width && !stills) {

rugnux/rugnux_cli.cpp

@@ -254,7 +254,7 @@ void print_usage() {
     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. 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-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 the curve is taken from the spot finder's peak-excluded background of the whole frame, and =auto applies it to every image of a run where it makes the background ring predict the background of empty apertures better (measured in the pre-scan): 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;
@@ -3429,7 +3429,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                     radial.has_value() ? (*radial ? "on" : "off")
                                        : rotation_indexing
                                              ? "auto (per image, where the ice score sees a smooth ring)"
-                                             : "auto (every image, from the radial background of the whole frame)");
+                                             : "auto (decided in the pre-scan on empty apertures; the curve from the whole frame)");
     }
 
     // Detection defaults differ by workflow (RugnuxDefaultSpotFinding); each flag always wins.

tests/BackgroundApertureTest.cpp

@@ -26,15 +26,25 @@ ImagePreprocessorBuffer Frame(double level, double ripple, uint32_t seed) {
     return img;
 }
 
+const RingBackground no_rings({}, 0, 0.0f, 0.0f);
+
+DiffractionGeometry Geometry() {
+    DiffractionGeometry g;
+    g.BeamX_pxl(300.0f).BeamY_pxl(300.0f).DetectorDistance_mm(60.0f).PixelSize_mm(0.075f).Wavelength_A(1.0f);
+    return g;
+}
+
 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);
+        background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(),
+                                          DetectorModules(), tally);
     }
     REQUIRE(tally.apertures > 100);
-    return static_cast<int>(r3[background_aperture::Best(tally)]);
+    CHECK(tally.apertures_radial == 0);
+    return static_cast<int>(r3[background_aperture::Best(tally.MeanSquare(false))]);
 }
 
 }   // namespace
@@ -62,11 +72,47 @@ 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);
+    background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(),
+                                      DetectorModules(), 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);
+    background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(),
+                                      DetectorModules(), masked);
     CHECK(masked.apertures == 0);
     CHECK(all.apertures > 0);
 }
+
+// A bright background that is a curved function of the radius alone: the radial correction, read off
+// the frame's ring background, makes the ring predict the disk better than it does without it.
+TEST_CASE("BackgroundAperture_RadialCorrectionHelpsOnARadialBackground", "[Integration]") {
+    const DiffractionGeometry g = Geometry();
+    const auto background = [](double r) { return 2000.0 + 0.2 * (r - 150.0) * (r - 150.0); };
+    std::mt19937 rng(3);
+    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(background(std::hypot(x - 300.0, y - 300.0)));
+            img[static_cast<size_t>(y) * W + x] = p(rng);
+        }
+    const float low_q = 0.05f, dq = 0.005f;
+    std::vector<float> profile(1000);
+    for (size_t i = 0; i < profile.size(); i++) {
+        const double q = low_q + (i + 0.5) * dq;
+        const double r = 60.0 / 0.075 * std::tan(2.0 * std::asin(q * 1.0 / (4.0 * PI)));
+        profile[i] = static_cast<float>(background(r));
+    }
+    const RingBackground rings(profile, static_cast<int32_t>(profile.size()), low_q, dq);
+    const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f);
+    background_aperture::Tally tally;
+    // Unclipped: a gradient this steep would put the high side of the ring above a photon clip.
+    background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 0.0f, 1.0f, rings, g, DetectorModules(), tally);
+    REQUIRE(tally.apertures_radial == tally.apertures);
+    REQUIRE(tally.apertures > 100);
+    const auto plain = tally.MeanSquare(false), radial = tally.MeanSquare(true);
+    // The pre-scan's question: at the ring the plain differences choose, does the correction help.
+    const int best = background_aperture::Best(plain);
+    CHECK(radial[best] < 0.8 * plain[best]);
+    // With the curvature taken out, the noise wins again and a wide ring is best.
+    CHECK(r3[background_aperture::Best(radial)] >= 12.0f);
+}