integration: the background is fitted over the ring that survives, not averaged over it

The signal disk and the background ring are concentric, which is the whole reason a linear
background cancels between them. Within the outer ring radius of the edge of the sensor array
that concentricity is gone: the ring loses its outer part while the disk barely loses anything,
so what is left of the ring sits further into the detector, where the radial background is
higher, and the reflection reads low.

Measured at signal-free positions four pixels from a border: the ring reads 162.35 counts per
pixel against a true background over the disk of 160.57, which over a hundred disk pixels is 182
counts of deficit, against 216 to 239 observed. <I/sigma> runs -1.83, -2.34 and -1.18 at nought
to three, three to six and six to nine pixels from the border, and recovers exactly at the outer
ring radius. The same reflection measured at a border reads 179 counts lower than in the
interior over seven thousand matched pairs.

A masked module gap does the same thing but signed by the direction of the displacement, which is
why nothing has caught this: at a gap the two populations cancel in the mean, while at the sensor
border the truncation is always inward, so the bias is always negative.

The background is now the intercept of a straight line in radial offset over whatever ring pixels
survive, read at the reflection's centre. Three extra sums per ring pixel and no extra reads; the
radial distance was already computed there. It is exact under any truncation and reduces to the
mean when the ring is whole, so it is unconditional rather than a mode: a badly truncated ring
pays in sigma, through the variance the fit honestly reports, rather than in a rejection.

On the crystal where this surfaced the outermost shell's correlation with a deposited model goes
from -0.234 to +0.004, and the shell above it from -0.091 to +0.179. Correcting beats discarding:
dropping every observation within fifteen pixels of a border reached only -0.019 and +0.127,
because the corrected observations still carry signal. Interior reflections do not move.

The cost is one geometry: where a neighbour mask has already truncated the ring almost everywhere,
the fit roughly doubles the variance of the background estimate while finding no gradient worth
removing, and a crowded small detector loses one to two points of CC1/2 in its finest shells.

Also: the MINPK denominator counted only profile mass that lands on the detector, so a reflection
whose peak is off the sensor scored a perfect one and no guard could fire.

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

docs/CHANGELOG.md

@@ -3,6 +3,8 @@
 
 ### 1.0.0-rc.167
 
+* The background under a reflection is now fitted as a straight line in radius over the background ring and read off at the reflection's centre, instead of being taken as the ring mean, so a reflection whose ring is cut by the edge of the array or by a mask is no longer measured against the background somewhere else.
+* A reflection whose expected profile falls largely outside the detector is now rejected by `--overlap-minpk`, which previously measured that fraction only over the part of the profile that was on the array and so could never fire.
 * `rugnux --model` reports CC(model, data) - the correlation of the merged intensities with the placed, scaled model - by resolution shell, on the same shells as CC1/2, with the reflection count and a significance for each.
 * `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale.
 * The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off.

docs/CPU_DATA_ANALYSIS_INTEGRATION.md

@@ -119,7 +119,15 @@ $
 \hat{b} = \frac{B}{n_B},\qquad
 \hat{I} = S - n_S \hat{b},
 $
-with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 count (pixel values are photon counts) and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates the **variance** by $1+n_S/n_B$ — 1.11 with the shipped circular stencil ($n_S = 45$, $n_B = 408$) — and so understates $\sigma$ by up to $\sqrt{1+n_S/n_B} \approx 1.05$, a bound attained on background-limited (weak) reflections and falling towards 1 on strong ones, where $S$ dominates; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run. The same term is carried into the profile fit (§9.3), where it adds $\big(\sum P/v \,\big/ \sum P^2/v\big)^2\,\mathrm{var}(\hat b)$ — the square of $\partial I/\partial\hat b$ for that fit; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below, where $S$ and $n_S$ then count only the pixels that were actually read.
+with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 count (pixel values are photon counts) and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b\,(1/n_B + \bar\rho^2/S_{\rho\rho})$ (the fitted background just below; $\hat b/n_B$ whenever the ring is whole), and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates the **variance** by $1+n_S/n_B$ — 1.11 with the shipped circular stencil ($n_S = 45$, $n_B = 408$) — and so understates $\sigma$ by up to $\sqrt{1+n_S/n_B} \approx 1.05$, a bound attained on background-limited (weak) reflections and falling towards 1 on strong ones, where $S$ dominates; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run. The same term is carried into the profile fit (§9.3), where it adds $\big(\sum P/v \,\big/ \sum P^2/v\big)^2\,\mathrm{var}(\hat b)$ — the square of $\partial I/\partial\hat b$ for that fit; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below, where $S$ and $n_S$ then count only the pixels that were actually read.
+
+**Fitted, not averaged (always on).** The ring mean is the background at the ring's *centroid*, and that is the value under the signal disk only while the ring is whole — which is exactly the concentricity the paragraph on the radial correction below rests on. Within $r_3$ of the edge of the array (or of a large mask) the ring loses its outer part while the much smaller $r_1$ disk keeps nearly all of its own, so the surviving ring sits off-centre along the radius, where a radial background has a different value; the mean then carries that difference into every one of the $n_S$ disk pixels it is subtracted from. The bias scales with the background *level*, so it is fractions of a count on ordinary data and large where the background is not: measured at a sensor border on a dataset running 164 counts/px, the ring read $+1.8$ counts/px high over a 100-pixel disk, i.e. $-180$ counts on every affected reflection — and always with the same sign, because a sensor boundary truncates the ring inward and nothing averages it away.
+
+So $\hat b$ is the **intercept of a straight line in the radial offset**, fitted over whatever ring pixels survive and evaluated at the reflection's own centre:
+$
+\hat b \;=\; \bar B \;-\; \frac{S_{\rho B}}{S_{\rho\rho}}\,\bar\rho ,
+$
+with $\rho$ each ring pixel's projection on the beam→reflection direction (already computed for the stencil, so this costs three more sums and no extra pixel reads). This is exact for a radial background under *any* truncation, needs no threshold and rejects nothing; a whole ring has $\bar\rho\approx0$ and the estimate reduces to the mean. The price is paid where it is due: the intercept's variance is the mean's inflated by $1+n_B\bar\rho^2/S_{\rho\rho}$, so a badly truncated ring is extrapolated back to the centre and reports the extra uncertainty rather than being discarded. The same fit removes the corresponding bias at module gaps and masked regions, which the boundary case only makes visible because there the truncation always has one sign. A *tangential* background gradient at a truncation is not corrected — the fit is deliberately one-dimensional, because the background is a function of radius and a second free direction would only take in noise.
 
 **High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), lowered by `rugnux` to $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. That is only a default — the flag sets $n$ whatever the bandwidth is. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px.
 
@@ -127,7 +135,7 @@ The clip cuts only the high tail, which matters: the **symmetric** trimmed mean
 
 Both estimators are computed in the shared background pass, but only the trim reaches plain box summation: the high-side clip is skipped for `--integrator boxsum`, which therefore uses the plain ring mean unless `--background-trim` is given.
 
-**Radial background correction (opt-in).** A ring mean estimates the background *under* the signal disk correctly only if the background is flat there. The signal disk and the ring are concentric, so for a background that is **linear** in position $\langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk}$ identically — a plane or gradient fit buys exactly nothing. The leading error is the **curvature** of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. That error is a kernel over radial offset,
+**Radial background correction (opt-in).** The signal disk and an **untruncated** ring are concentric, so for a background that is **linear** in position $\langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk}$ identically — and where the ring *is* truncated, the linear fit above has already taken that term out. What neither removes is the **curvature** of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. That error is a kernel over radial offset,
 $
 \delta \hat b \;=\; \textstyle\sum_k \kappa_k\, \bar B(r_0+k),
 $
@@ -153,7 +161,7 @@ v = \max\!\left(B + I\,P,\ \tfrac{1}{2}B\right),
 $
 where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). The plug-in $I$ enters **as it is**: half-wave rectifying it, $v=B+\max(I,0)P$, lets $v$ — and with it the reported $1/\sum P^2/v$ — respond only to *upward* fluctuations of a noisy estimate, which adds $\approx0.4\,\sigma\sum P^3/(\sum P^2)^2$ to every $\sigma$ whatever the count rate. That offset is invisible on strong reflections and a large fractional inflation on weak ones; the $\tfrac12 B$ clamp keeps $v$ positive without reintroducing it. As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the background term is floored at $0.01$ ct/px — enough to keep $P^2/v$ finite when the ring mean reads exactly zero, which a ring of $n_B$ pixels cannot distinguish from any background below $\approx1/n_B$. The rotation/excitation partiality is carried exactly as in the box-sum path.
 
-**Pixels the fit cannot use (MINPK).** A profile fit is the amplitude of a *normalised* profile, so a pixel left out of the sum renormalises the estimator by construction: it costs information — $\sum P^2/v$ shrinks and $\sigma$ grows — but biases nothing. That is what keeps a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload: those pixels are simply not read, and the fit is taken over the rest, exactly as the shared pixels of a crowded reflection are (`--overlap exclude`). The reflection is kept only while enough of the expected profile survives — at least `--overlap-minpk` of the profile mass that falls on the detector at all, default 0.75, which is XDS's `MINPK` and dials' `valid_foreground_threshold`. The complete reflections alone teach the profile, its resolution shells and their widths. `--integrator boxsum` has no profile to renormalise with and keeps the all-or-nothing rule of §9.2.
+**Pixels the fit cannot use (MINPK).** A profile fit is the amplitude of a *normalised* profile, so a pixel left out of the sum renormalises the estimator by construction: it costs information — $\sum P^2/v$ shrinks and $\sigma$ grows — but biases nothing. That is what keeps a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload: those pixels are simply not read, and the fit is taken over the rest, exactly as the shared pixels of a crowded reflection are (`--overlap exclude`). The reflection is kept only while enough of the expected profile survives — at least `--overlap-minpk` of the profile's **whole** mass, default 0.75, which is XDS's `MINPK` and dials' `valid_foreground_threshold`. Mass the array does not reach counts as lost, like any other: measured over the on-detector part alone the ratio is 1 by construction for a reflection whose peak is off the sensor, and no threshold on it can fire. The complete reflections alone teach the profile, its resolution shells and their widths. `--integrator boxsum` has no profile to renormalise with and keeps the all-or-nothing rule of §9.2.
 
 "Biases nothing" holds only while the profile *model* is exact. Lose the peak and the amplitude is set by the wings alone, so the result stops being a measurement of the reflection and becomes a measurement of how well the fitted shape describes it. The worst case is a pixel invalidated *by the flux it saw* — a detector's per-frame overload marker: that pixel goes missing **because** the reflection was bright, so the loss concentrates on the strong low-resolution reflections that are the largest terms of $R_\mathrm{meas}$, where the fit reads $-50\%$ against the symmetry mates. MINPK cannot catch it, because it cuts on profile *mass* and the peak of a broad spot is a few percent of the mass. So a second condition applies alongside it, on any unreadable pixel whatever made it unreadable: **no unreadable pixel may carry more than 0.9 of the profile's own peak value**. As a fraction of the peak rather than a radius in pixels, that scales with the spot — for a Gaussian it is a cut at $\sqrt{-2\ln f}\,\sigma = 0.46\sigma$, the peak pixel alone where $\sigma$ is 0.8 px and the crest of the ridge where the profile is a bandwidth streak — and it costs well under 0.1 % of the recovered observations.
 

docs/RUGNUX_ADVANCED.md

@@ -352,7 +352,7 @@ Integration:
 | `--integration-stencil <k>` | Push the `r2..r3` background ring out by `k` times the beam's radial streak `bandwidth·R_px`, per reflection (default `0` = the 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 — and the growth is capped at `2·r3`. The neighbour exclusion grows with it, so on a crowded pattern a few reflections can be left with too little background and dropped. Needs `--bandwidth`: on a monochromatic beam the streak is zero and this does nothing |
 | `--background-clip <n>` | Monochromatic (rotation + still): high-side clip of the background ring at `mean + n·√mean` (default 4; 0 = off). The default background estimator — it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3σ; ignored by `--integrator boxsum` |
 | `--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. Rings holding more than 512 pixels fall back to the plain mean (the GPU sorts the ring in shared memory and the CPU now matches it), which the default radii never reach but wide ones do |
-| `--background-radial[=on\|off\|auto]` | Correct the background ring for the **curvature** of the radial background (default **off**). Disk and 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, since the model is a function of radius alone: on ice made of discrete crystallite spots there is no smooth ring and the correction makes the bias worse. Ignored by `--integrator boxsum` (no clip pass to take the curve from) |
+| `--background-radial[=on\|off\|auto]` | Correct the background ring for the **curvature** of the radial background (default **off**). Disk and ring are concentric, and the linear part of the background is fitted out of the ring in any case, so 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, since the model is a function of radius alone: on ice made of discrete crystallite spots there is no smooth ring and the correction makes the bias worse. Ignored by `--integrator boxsum` (no clip pass to take the curve from) |
 | `--integration-high-resolution <num>` | High-resolution limit for prediction and integration. Omitted (or 0) means integration extends as far as the detector reaches — which is what the predictor can place on the detector anyway, since it rejects reflections that miss it. Set a value to integrate less than the detector offers |
 | `--max-hkl <n>` | Predict reflections with \|h\|,\|k\|,\|l\| ≤ `n` (max 511). By default this is derived per crystal from the refined cell as `ceil(max(a,b,c)/d_min) + 1`, which is the exact bound: the predictor keeps only \|q\| ≤ 1/d_min and `h = a·q`, so no reflection can lie outside it and no candidate inside it is wasted on a shorter axis. Set it only to override that |
 | `--bandwidth <num>` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) |

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -168,6 +168,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         int n_disk = 0, n_own = 0;   // pixels in the signal disk, and how many are this reflection's
         double bkg_sum = 0.0;
         int n_bkg = 0;
+        // Sums for the radial-linear background fit below, over n_fit ring pixels: sum(rad),
+        // sum(rad^2), sum(rad*px), with rad the pixel's offset along the beam->reflection direction
+        // (already computed for the stencil). The clip pass below replaces them with its own.
+        double rad_s1 = 0.0, rad_s2 = 0.0, rad_sp = 0.0;
         // Ring pixels a NEIGHBOUR's signal region occupies. refl_mask marks d.inner < r2_sq and the
         // ring is d.inner >= r2_sq, so the two are complementary and a masked ring pixel is always
         // someone else's - never this reflection's own core. That makes this an exact count of what
@@ -202,12 +206,16 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; }
                     if (!valid(px)) continue;
                     bkg_sum += static_cast<double>(px);
+                    rad_s1 += d.rad;
+                    rad_s2 += static_cast<double>(d.rad) * d.rad;
+                    rad_sp += static_cast<double>(d.rad) * px;
                     if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
                     ++n_bkg;
                 }
             }
 
         int n_bkg_used = n_bkg;   // pixels behind the FINAL background value (trim/clip shrink it)
+        int n_fit = n_bkg;        // pixels behind rad_s1/rad_s2/rad_sp (the clip pass shrinks them)
         // A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the
         // reflection outright. A profile fit does not need it to: the fit is the amplitude of a
         // NORMALISED profile, so leaving pixels out renormalises the estimator by construction and
@@ -249,6 +257,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                 const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(std::max(out.bkg, 1.0));
                 double s = 0.0;
                 int n = 0;
+                double c_s1 = 0.0, c_s2 = 0.0, c_sp = 0.0;
                 for (int y = y0; y <= y1; ++y)
                     for (int x = x0; x <= x1; ++x) {
                         const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
@@ -259,6 +268,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                         if (static_cast<double>(px) <= thr) {
                             s += px;
                             ++n;
+                            c_s1 += d.rad;
+                            c_s2 += static_cast<double>(d.rad) * d.rad;
+                            c_sp += static_cast<double>(d.rad) * px;
                             if (bkg_radial) {
                                 // The radial curve is binned on the TRUE detector radius, so the
                                 // offset here is the unshrunk radial projection.
@@ -268,7 +280,44 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                             }
                         }
                     }
-                if (n > 5) { out.bkg = s / n; n_bkg_used = n; }
+                if (n > 5) {
+                    out.bkg = s / n;
+                    n_bkg_used = n;
+                    bkg_sum = s;
+                    rad_s1 = c_s1; rad_s2 = c_s2; rad_sp = c_sp; n_fit = n;
+                }
+            }
+            // --- Fit the background rather than average it. The ring mean is the background AT THE
+            //     RING'S CENTROID, and that is the value under the disk only while the ring is
+            //     whole: the two are concentric, which is the whole reason a background linear in
+            //     position cancels between them. Within r3 of the array edge (or of a large mask)
+            //     the ring loses its outer part while the much smaller disk keeps nearly all of its
+            //     own, so the surviving ring sits off-centre along the radius, where the background
+            //     is different - and the mean carries that difference into every one of the ~n_inner
+            //     pixels it is subtracted from. Measured at a sensor border on a 164 counts/px
+            //     background: +1.8 counts/px over a 100 px disk, -180 counts on the intensity, and
+            //     always the same sign, because the truncation is always inward.
+            //     A straight line in the radial offset, evaluated at rad = 0 (the reflection centre),
+            //     is exact for a linear background under ANY truncation. It reduces to the mean when
+            //     the ring is whole (the centroid is then the centre), needs no threshold and rejects
+            //     nothing. The curvature that survives even an untruncated ring is a separate,
+            //     smaller term and is what --background-radial corrects; a linear fit is in that
+            //     kernel's null space, so the two do not overlap.
+            const double rad_mean = rad_s1 / n_fit;
+            const double Sxx = rad_s2 - rad_s1 * rad_mean;
+            const double Sxy = rad_sp - rad_mean * bkg_sum;
+            // The intercept is estimated from the ring pixels' spread in radius, so its variance is
+            // the mean's inflated by how far the ring's centroid had to be extrapolated back to the
+            // centre. A whole ring has rad_mean ~ 0 and gets bkg/n_bkg_used, as before; a badly
+            // truncated one pays for the extrapolation in sigma instead of in a rejection.
+            // The Poisson scale of that variance stays the ring MEAN, which is an average of counts
+            // and cannot be negative. The fitted intercept can be, on a background of well under one
+            // count per pixel, and a negative variance is not a variance - it reaches sigma as a NaN.
+            const double bkg_mean = out.bkg;
+            double bkg_var_factor = 1.0 / n_bkg_used;
+            if (Sxx > 0.0) {
+                out.bkg -= (Sxy / Sxx) * rad_mean;
+                bkg_var_factor += rad_mean * rad_mean / Sxx;
             }
             // The sum is over the pixels actually READ, so that is the count the background is
             // subtracted with; with nothing missing it is the whole disk, exactly as before.
@@ -277,8 +326,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
             // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
             // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
-            out.bkg_var = out.bkg / n_bkg_used;
-            out.var_bkg = static_cast<double>(n_inner_valid) * out.bkg
+            out.bkg_var = bkg_mean * bkg_var_factor;
+            out.var_bkg = static_cast<double>(n_inner_valid) * bkg_mean
                           + static_cast<double>(n_inner_valid) * n_inner_valid * out.bkg_var;
             out.I_sum = I_sum;
             out.n_inner = static_cast<int>(n_inner_valid);
@@ -521,10 +570,14 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             for (int dx = -Rf; dx <= Rf; ++dx) {
                 const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
                 if (Pp <= 0.0) continue;
+                // Profile mass the array does not reach is mass this fit has LOST, so it belongs in
+                // the denominator. Counted after the bounds test it cancels out of p_valid / p_grid,
+                // which is then exactly 1 for a reflection whose peak is off the sensor entirely -
+                // and no threshold on a ratio that is always 1 can fire.
+                p_grid += Pp;
                 const int x = rh.cx + dx, y = rh.cy + dy;
                 if (x < 0 || y < 0 || x >= W || y >= H) continue;
                 const bool in_disk = dx * dx + dy * dy < r1_sq;
-                p_grid += Pp;
                 p_peak = std::max(p_peak, Pp);
                 if (in_disk) m_all += Pp;
                 if (!valid(img[y * W + x])) {

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -185,6 +185,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     unsigned long long *s_radv = s_rad;                       // signed value carried as two's complement
     int   *s_radn = (int *) (s_rad + p.rad_w);
     __shared__ int   s_radbase;
+    __shared__ float s_slot[3][BRAGG_WARPS_MAX];   // per-warp partials, see block_sum
     if (threadIdx.x == 0) {
         const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y;
         s_r0 = sqrtf(rx * rx + ry * ry);
@@ -209,6 +210,9 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     long long l_x = 0, l_y = 0;   // positions behind l_Isum, for the background-free centroid
     int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0;
     long long l_bkg = 0;
+    // Sums for the radial-linear background fit: sum(rad), sum(rad^2), sum(rad*px) over the ring.
+    // See the CPU engine for what they are for.
+    float l_rs1 = 0.0f, l_rs2 = 0.0f, l_rsp = 0.0f;
     for (int t = threadIdx.x; t < area; t += blockDim.x) {
         const int x = x0 + t % bw, y = y0 + t / bw;
         const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
@@ -231,6 +235,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             const int32_t px = img[y * p.W + x];
             if (!valid(px)) continue;
             l_bkg += px; ++l_nb;
+            l_rs1 += d.rad; l_rs2 += d.rad * d.rad; l_rsp += d.rad * (float) px;
         }
     }
     WARP_ATOMIC_ADD(s_Isum, (unsigned long long) l_Isum);
@@ -245,6 +250,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     WARP_ATOMIC_ADD(s_ndisk, l_nd);
     WARP_ATOMIC_ADD(s_nown, l_no);
     WARP_ATOMIC_ADD(s_nbkg_nb, l_nbnb);
+    // Float sums, so a fixed-order block reduction rather than an atomic; see block_sum.
+    const float ring_s1 = block_sum(l_rs1, s_slot[0]);
+    const float ring_s2 = block_sum(l_rs2, s_slot[1]);
+    const float ring_sp = block_sum(l_rsp, s_slot[2]);
     __syncthreads();
 
     for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; }
@@ -321,8 +330,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     }
 
     // Second ring pass for the high-side sigma-clip (re-reads the annulus; avoids storing bkg values).
+    float clip_s1 = 0.0f, clip_s2 = 0.0f, clip_sp = 0.0f;
     if (s_accept && p.bkg_clip_nsigma > 0.0f && !do_trim) {
         long long c_l = 0; int cn_l = 0;
+        float c_s1 = 0.0f, c_s2 = 0.0f, c_sp = 0.0f;
         for (int t = threadIdx.x; t < area; t += blockDim.x) {
             const int x = x0 + t % bw, y = y0 + t / bw;
             const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
@@ -332,6 +343,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             if (!valid(px)) continue;
             if ((long long) px <= s_thr_i) {
                 c_l += px; ++cn_l;
+                c_s1 += d.rad; c_s2 += d.rad * d.rad; c_sp += d.rad * (float) px;
                 if (n_rad > 0) {
                     // The radial curve is binned on the TRUE detector radius, so this is the
                     // unshrunk radial projection - no per-pixel sqrt.
@@ -344,6 +356,9 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
             }
         }
         WARP_ATOMIC_ADD(s_clipsum, (unsigned long long) c_l); WARP_ATOMIC_ADD(s_clipn, cn_l);
+        clip_s1 = block_sum(c_s1, s_slot[0]);
+        clip_s2 = block_sum(c_s2, s_slot[1]);
+        clip_sp = block_sum(c_sp, s_slot[2]);
     }
     __syncthreads();
 
@@ -367,19 +382,39 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
 
     double bkg = s_bkg;
     int n_bkg_used = s_nbkg;   // pixels behind the FINAL background value (trim/clip shrink it)
+    // The fit sums and the pixel sum they go with: the whole ring, or the clip's survivors.
+    double fit_s1 = ring_s1, fit_s2 = ring_s2, fit_sp = ring_sp;
+    double fit_sum = (double) (long long) s_bkgsum;
+    int n_fit = s_nbkg;
     if (do_trim) {
         const int nb = min(s_bn, BKG_TRIM_MAX);
         const int lo = (int) (nb * p.bkg_trim), hi = nb - lo;
         if (hi > lo) n_bkg_used = hi - lo;
     }
-    if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) { bkg = (double) (long long) s_clipsum / (double) s_clipn; n_bkg_used = s_clipn; }
+    if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) {
+        bkg = (double) (long long) s_clipsum / (double) s_clipn; n_bkg_used = s_clipn;
+        fit_s1 = clip_s1; fit_s2 = clip_s2; fit_sp = clip_sp;
+        fit_sum = (double) (long long) s_clipsum; n_fit = s_clipn;
+    }
+    // Fit the background rather than average it: the ring mean is the background at the RING'S
+    // centroid, which is the value under the disk only while the ring is whole. See the CPU engine.
+    const double rad_mean = fit_s1 / (double) n_fit;
+    const double Sxx = fit_s2 - fit_s1 * rad_mean;
+    const double Sxy = fit_sp - rad_mean * fit_sum;
+    // The variance's Poisson scale stays the ring MEAN, which cannot be negative; see the CPU engine.
+    const double bkg_mean = bkg;
+    double bkg_var_factor = 1.0 / (double) n_bkg_used;
+    if (Sxx > 0.0) {
+        bkg -= (Sxy / Sxx) * rad_mean;
+        bkg_var_factor += rad_mean * rad_mean / Sxx;
+    }
     const long long Isum = (long long) s_Isum;
     // The sum is over the pixels actually READ, so that is the count the background is subtracted
     // with; with nothing missing it is the whole disk, exactly as before.
     const double I = (double) Isum - (double) s_ninner_valid * bkg;
     // See the CPU engine: bkg is estimated from n_bkg_used pixels and subtracted n_inner times, so
     // var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree.
-    const double bkg_var = bkg / (double) n_bkg_used;
+    const double bkg_var = bkg_mean * bkg_var_factor;
     const double var_bkg_term = (double) s_ninner_valid * (double) s_ninner_valid * bkg_var;
     double sigma = 1.0;
     uint8_t hasobs = 0; double ox = 0.0, oy = 0.0;
@@ -404,7 +439,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     cy_o[i] = (int) lroundf(cy);
     I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg;
     bkgvar_o[i] = (float) bkg_var;
-    varbkg_o[i] = (float) ((double) s_ninner_valid * bkg + var_bkg_term);
+    varbkg_o[i] = (float) ((double) s_ninner_valid * bkg_mean + var_bkg_term);
     isum_o[i] = (float) Isum;
     ninner_o[i] = s_ninner_valid;
     rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0);
@@ -684,11 +719,13 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x
     for (int k = threadIdx.x; k < GfGf; k += blockDim.x) {
         const float Pp = Pbuf[k];
         if (Pp <= 0.0f) continue;
+        // Profile mass the array does not reach is mass this fit has LOST, so it belongs in the
+        // denominator. See the CPU engine.
+        l_pgrid += Pp;
         const int dx = k % Gf - Rf, dy = k / Gf - Rf;
         const int x = cx + dx, y = cy + dy;
         if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue;
         const bool in_disk = (float) (dx * dx + dy * dy) < p.r1_sq;
-        l_pgrid += Pp;
         l_ppeak = fmaxf(l_ppeak, Pp);
         if (in_disk) l_mall += Pp;
         if (!valid(img[y * p.W + x])) {
@@ -826,7 +863,7 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment
                               "BraggIntegrationEngineGPU: profile grid exceeds shared memory (r2 too large)");
     // boxsum's dynamic window sits alongside its static shared arrays (the trimmed-mean buffer
     // dominates them), so leave room for those rather than budgeting the whole block.
-    if (boxsum_shared_bytes + sizeof(int) * BKG_TRIM_MAX + 256 > prop.sharedMemPerBlock)
+    if (boxsum_shared_bytes + sizeof(int) * BKG_TRIM_MAX + 512 > prop.sharedMemPerBlock)
         throw JFJochException(JFJochExceptionCategory::GPUCUDAError,
                               "BraggIntegrationEngineGPU: background ring exceeds shared memory");
 

rugnux/rugnux_cli.cpp

@@ -190,10 +190,10 @@ 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, or 3 when --bandwidth is set; 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 off). 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 off). The signal disk and the background ring are concentric, and the linear part is fitted out of the ring in any case, so 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-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. Applies whatever --bandwidth is set to" << 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;
+    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, overloaded, or off the array altogether - 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;
     std::cout << "   --integrator <txt>              Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
     std::cout << "   --simple-stills                 stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl;
     std::cout << "   -q, --azim-q-spacing <num>      Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;

tests/BraggBackgroundTest.cpp

@@ -0,0 +1,167 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <cstdint>
+#include <vector>
+
+#include "../common/BraggIntegrationSettings.h"
+#include "../common/DetectorSetup.h"
+#include "../common/DiffractionExperiment.h"
+#include "../common/Reflection.h"
+#include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h"
+#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
+
+// The background under a reflection is estimated from a ring that is CONCENTRIC with the signal disk,
+// which is the whole reason a background varying across the reflection cancels between them. Within r3
+// of the edge of the array the ring loses its outer part and stops being concentric: what is left sits
+// off-centre along the radius, where a radial background has a different value, and the ring mean
+// carries that difference into every disk pixel it is subtracted from. These tests put a known radial
+// ramp under reflections at a range of distances from the border and ask for the background at the
+// reflection's own centre back.
+
+namespace {
+
+constexpr double BKG_LEVEL = 200.0;      // counts/px at the beam centre
+constexpr double BKG_GRADIENT = 2.0;     // counts/px per pixel of radius
+
+Reflection MakeReflection(float x, float y, int hkl) {
+    Reflection r{};
+    r.h = hkl; r.k = hkl; r.l = hkl;
+    r.predicted_x = x;
+    r.predicted_y = y;
+    r.d = 2.0f;
+    r.prescaling_corr = 1.0f;
+    r.partiality = 1.0f;
+    return r;
+}
+
+double TrueBackground(double x, double y, double beam_x, double beam_y) {
+    return BKG_LEVEL + BKG_GRADIENT * std::hypot(x - beam_x, y - beam_y);
+}
+
+DiffractionExperiment MakeExperiment(IntegratorMode mode, float beam_x, float beam_y) {
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV)
+              .BeamX_pxl(beam_x).BeamY_pxl(beam_y);
+    BraggIntegrationSettings settings;
+    settings.Integrator(mode);
+    experiment.ImportBraggIntegrationSettings(settings);
+    return experiment;
+}
+
+}   // namespace
+
+// The ramp is radial and there is no signal anywhere, so the answer is known exactly: the background
+// the engine reports has to be the ramp's value at the reflection's own centre, at the border as much
+// as in the middle of the array. Averaging the surviving ring instead reports its centroid's value,
+// which at the border is several pixels of radius away.
+TEST_CASE("BraggBackground_TruncatedRingIsNotBiased", "[Integration]") {
+    const float beam_x = 400.0f, beam_y = 400.0f;
+    const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::BoxSum, beam_x, beam_y);
+    const int W = static_cast<int>(experiment.GetXPixelsNum());
+    const int H = static_cast<int>(experiment.GetYPixelsNum());
+
+    ImagePreprocessorBuffer image(experiment.GetPixelsNum());
+    for (int y = 0; y < H; ++y)
+        for (int x = 0; x < W; ++x)
+            image[static_cast<size_t>(y) * W + x] =
+                    static_cast<int32_t>(std::lround(TrueBackground(x, y, beam_x, beam_y)));
+
+    // Reflections marching in towards the array from the bottom edge, at three azimuths so the border
+    // cuts the ring at a different angle to the radius each time, plus an interior control.
+    std::vector<Reflection> predicted;
+    int hkl = 1;
+    for (int x : {400, 700, 1000})
+        for (int inset : {2, 4, 6, 9, 12, 16, 40})
+            predicted.push_back(MakeReflection(static_cast<float>(x) + 0.3f,
+                                               static_cast<float>(H - 1 - inset) - 0.2f, hkl++));
+
+    BraggIntegrationEngineCPU engine(experiment);
+    const auto out = engine.Run(image, predicted, predicted.size(), 1);
+    REQUIRE(out.size() == predicted.size());
+
+    for (const auto &r : out) {
+        const double expected = TrueBackground(r.predicted_x, r.predicted_y, beam_x, beam_y);
+        INFO("reflection at " << r.predicted_x << "," << r.predicted_y
+             << " expected bkg " << expected << " got " << r.bkg);
+        CHECK(r.bkg == Catch::Approx(expected).margin(0.5));
+        // A flat background under the disk means the box sum has nothing above it to report.
+        CHECK(std::abs(r.I) < 6.0f * r.sigma);
+    }
+}
+
+// The same reflections read against a FLAT background: the fit must not invent a correction where
+// there is no gradient to correct, at the border or anywhere else.
+TEST_CASE("BraggBackground_FlatBackgroundIsUnchangedAtTheBorder", "[Integration]") {
+    const float beam_x = 400.0f, beam_y = 400.0f;
+    const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::BoxSum, beam_x, beam_y);
+    const int H = static_cast<int>(experiment.GetYPixelsNum());
+
+    ImagePreprocessorBuffer image(experiment.GetPixelsNum());
+    for (size_t i = 0; i < experiment.GetPixelsNum(); ++i)
+        image[i] = static_cast<int32_t>(BKG_LEVEL);
+
+    std::vector<Reflection> predicted;
+    int hkl = 1;
+    for (int inset : {2, 4, 6, 9, 12, 16, 40}) {
+        predicted.push_back(MakeReflection(700.3f, static_cast<float>(H - 1 - inset) - 0.2f, hkl++));
+        predicted.push_back(MakeReflection(static_cast<float>(inset) + 0.3f, 700.2f, hkl++));
+    }
+
+    BraggIntegrationEngineCPU engine(experiment);
+    const auto out = engine.Run(image, predicted, predicted.size(), 1);
+    REQUIRE(out.size() == predicted.size());
+    for (const auto &r : out) {
+        INFO("reflection at " << r.predicted_x << "," << r.predicted_y);
+        CHECK(r.bkg == Catch::Approx(BKG_LEVEL).margin(1e-3));
+    }
+}
+
+// MINPK is a fraction of the profile the fit can see, and the denominator has to be the WHOLE profile:
+// counted over the grid cells that land on the array it is 1 by construction for a reflection whose
+// peak is off the sensor, and then no threshold on it can fire.
+TEST_CASE("BraggBackground_ProfileMassOffTheArrayIsRejected", "[Integration]") {
+    const float beam_x = 400.0f, beam_y = 400.0f;
+    const DiffractionExperiment experiment =
+            MakeExperiment(IntegratorMode::ProfileGaussian, beam_x, beam_y);
+    const int W = static_cast<int>(experiment.GetXPixelsNum());
+    const int H = static_cast<int>(experiment.GetYPixelsNum());
+
+    ImagePreprocessorBuffer image(experiment.GetPixelsNum());
+    for (size_t i = 0; i < experiment.GetPixelsNum(); ++i)
+        image[i] = static_cast<int32_t>(BKG_LEVEL);
+    // Strong, well-formed spots so the profile is learned and the fit has something to work on.
+    auto add_spot = [&](float cx, float cy) {
+        for (int dy = -6; dy <= 6; ++dy)
+            for (int dx = -6; dx <= 6; ++dx) {
+                const int x = static_cast<int>(std::lround(cx)) + dx;
+                const int y = static_cast<int>(std::lround(cy)) + dy;
+                if (x < 0 || y < 0 || x >= W || y >= H) continue;
+                const double ex = x - cx, ey = y - cy;
+                image[static_cast<size_t>(y) * W + x] +=
+                        static_cast<int32_t>(std::lround(4000.0 * std::exp(-(ex * ex + ey * ey) / (2.0 * 1.3 * 1.3))));
+            }
+    };
+
+    std::vector<Reflection> predicted;
+    int hkl = 1;
+    for (int gy = 0; gy < 8; ++gy)
+        for (int gx = 0; gx < 8; ++gx) {
+            const float cx = 100.0f + 60.0f * gx, cy = 100.0f + 60.0f * gy;
+            add_spot(cx, cy);
+            predicted.push_back(MakeReflection(cx, cy, hkl++));
+        }
+    // A reflection whose predicted centre sits just outside the array: nearly all of its profile,
+    // its peak included, is off the sensor.
+    const size_t off_array = predicted.size();
+    predicted.push_back(MakeReflection(700.0f, static_cast<float>(H) + 2.0f, hkl++));
+
+    BraggIntegrationEngineCPU engine(experiment);
+    const auto out = engine.Run(image, predicted, predicted.size(), 1);
+    for (const auto &r : out)
+        CHECK(r.h != static_cast<int>(off_array) + 1);
+    CHECK(out.size() == off_array);
+}

tests/BraggIntegrationEngineGPUTest.cpp

@@ -45,17 +45,26 @@ Reflection MakeReflection(float x, float y, float d, int hkl) {
 // centre of every 5th, a mid-profile pixel of every 7th and a disk-edge pixel of every 11th. That is
 // the MINPK rescue's own case - a reflection kept and fitted over the pixels it has - and with it the
 // peak-loss rule, which has to fire on the same reflections in both engines.
+// border > 0 puts a radial background ramp of that many counts per pixel of radius under the scene
+// and moves the spot grid onto the edges of the array, so the r2..r3 ring of the outermost spots is
+// truncated by the boundary. That is the case the radial-linear background fit exists for, and the
+// two engines have to fit the same line over the pixels each of them kept.
 Scene BuildScene(size_t width, size_t height, int spacing = 60, float companion_dx = 0.0f,
-                 bool clip_spots = false) {
+                 bool clip_spots = false, float border = 0.0f) {
     Scene s;
     s.width = width;
     s.height = height;
     s.image.assign(width * height, 12);   // flat background
+    if (border > 0.0f)
+        for (size_t y = 0; y < height; ++y)
+            for (size_t x = 0; x < width; ++x)
+                s.image[y * width + x] = 200 + static_cast<int32_t>(std::lround(
+                        border * std::hypot(static_cast<double>(x) - 400.0, static_cast<double>(y) - 400.0)));
 
     // A grid of spots, well separated so background rings do not overlap the neighbours' disks.
     // A spread of intensities (some weak, some very strong) and a spread of d (so several resolution
     // shells are populated) exercises the strong-spot selection, shell learning and the fit.
-    const int margin = 45;
+    const int margin = border > 0.0f ? 3 : 45;
     int hkl = 1;
     for (int gy = 0; margin + gy * spacing < static_cast<int>(height) - margin; ++gy) {
         for (int gx = 0; margin + gx * spacing < static_cast<int>(width) - margin; ++gx) {
@@ -145,7 +154,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
                      float stencil_k = 0.0f,
                      float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f,
                      OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f,
-                     bool clip_spots = false) {
+                     bool clip_spots = false, float border = 0.0f) {
     const DiffractionExperiment experiment =
             MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3,
                            overlap);
@@ -154,7 +163,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
     const size_t npixel = experiment.GetPixelsNum();
     REQUIRE(npixel == width * height);
 
-    const Scene scene = BuildScene(width, height, spacing, companion_dx, clip_spots);
+    const Scene scene = BuildScene(width, height, spacing, companion_dx, clip_spots, border);
     REQUIRE(scene.image.size() == npixel);
     REQUIRE(scene.predicted.size() > 60);
 
@@ -183,7 +192,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
     if (clip_spots) {
         // Guard against the coverage going vacuous: the punched pixels have to actually cost some
         // reflections, or the two engines are being compared on a case neither of them meets.
-        const Scene clean_scene = BuildScene(width, height, spacing, companion_dx, false);
+        const Scene clean_scene = BuildScene(width, height, spacing, companion_dx, false, border);
         ImagePreprocessorBuffer clean_image(npixel);
         for (size_t i = 0; i < npixel; ++i)
             clean_image[i] = clean_scene.image[i];
@@ -293,6 +302,22 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
         CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f,
                         0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, true);
     }
+    // Spots on the edges of the array under a radial background ramp: their background rings are
+    // truncated by the boundary, so the radial-linear fit that replaces the ring mean has a different
+    // pixel set - and a different line - for every one of them. The clip and the trim select that
+    // pixel set differently in the two engines, so both estimators need the case.
+    SECTION("ProfileGaussian border gradient") {
+        CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f,
+                        0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 2.0f);
+    }
+    SECTION("BoxSum border gradient") {
+        CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f,
+                        0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 2.0f);
+    }
+    SECTION("ProfileGaussian border gradient trim") {
+        CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f, false, 60, 0.0f,
+                        0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 2.0f);
+    }
     // The radial background curvature correction is computed independently in the two engines
     // (host loop vs radial_correct kernel), so it needs its own parity coverage.
     SECTION("BoxSum radial")             { CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, true); }

tests/CMakeLists.txt

@@ -38,6 +38,7 @@ ADD_EXECUTABLE(jfjoch_test
         ROIIntegrationGPUTest.cpp
         BraggIntegrationEngineGPUTest.cpp
         BraggIntegrationEngineCompressedImageTest.cpp
+        BraggBackgroundTest.cpp
         BraggStencilTest.cpp
         LossyFilterTest.cpp
         ImageBufferTest.cpp