bragg integration: trimmed-mean background, on by default for rotation

The local Bragg background is the mean over the r2..r3 ring. That mean reads
high because the contaminants that survive the signal-disk mask - neighbour-
spot wings, tails, zingers - are one-sided (positive), so it over-subtracts.
Since a weak intensity is a small difference of large numbers (I = S - nS*b),
a per-pixel background bias is fractionally largest at the resolution edge,
exactly where it hurts most.

Replace the ring mean with a symmetric trimmed mean (sort the ring, drop the
lowest and highest fraction f, average the rest), controlled by a new
BraggIntegrationSettings field and the rugnux `--background-trim <f>` option
(default f=0.10; 0 restores the plain mean). Default on for monochromatic
(rotation) data; broadband (stills) keep their tuned high-side sigma-clip, so
the base engine forces the trim to 0 there. Implemented in both the CPU engine
and the GPU kernel (shared-memory bitonic sort per block, flat-mean fallback
above BKG_TRIM_MAX ring pixels); the two agree.

25-crystal rotation battery (fixed SG/cell): <I/sigma> improved on every
crystal (median +50%), ISa on 20/22, resolution-edge R_meas fell several-fold
(e.g. lyso_ref 1.0 A 108%->43%). Last-shell CC1/2 is rescued where the plain
mean had collapsed to noise (Thau_9 at ~2.0 A 3.8%->64%, ~0.5 A of resolution
regained; cytC_10 0.2%->10%) at a small cost (1-3%) in already-clean shells -
it flattens the CC1/2 fall-off rather than shifting it. Stills unchanged.
Documented in CPU_DATA_ANALYSIS.md section 9.2.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

common/BraggIntegrationSettings.cpp

@@ -105,3 +105,15 @@ BraggIntegrationSettings &BraggIntegrationSettings::StillPartiality(bool input)
 bool BraggIntegrationSettings::GetStillPartiality() const {
     return still_partiality;
 }
+
+BraggIntegrationSettings &BraggIntegrationSettings::BackgroundTrimFraction(float input) {
+    check_finite("Background trim fraction", input);
+    check_min("Background trim fraction", input, 0.0);
+    check_max("Background trim fraction", input, 0.49);   // must leave a central majority after trimming
+    bkg_trim_fraction = input;
+    return *this;
+}
+
+float BraggIntegrationSettings::GetBackgroundTrimFraction() const {
+    return bkg_trim_fraction;
+}

common/BraggIntegrationSettings.h

@@ -21,6 +21,12 @@ class BraggIntegrationSettings {
     std::optional<float> fixed_profile_radius;
     float minimum_sigma_in_regards_to_i = 0.02;
     bool still_partiality = false;   // experimental stills excitation-error partiality (rugnux --still-partiality)
+    // Symmetric trimmed-mean fraction for the r2..r3 background ring: drop the lowest and highest this
+    // fraction of ring pixels before averaging. Resists the high-side contamination (neighbour-spot
+    // wings, tails, zingers) that biases the plain ring mean up and makes it over-subtract weak
+    // high-angle reflections. Applied to monochromatic (rotation) data; the integration engine keeps
+    // the tuned high-side sigma-clip for stills instead. 0 = plain ring mean (rugnux --background-trim).
+    float bkg_trim_fraction = 0.10f;
 
 public:
     BraggIntegrationSettings& R1(float input);
@@ -30,6 +36,7 @@ public:
     BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional<float> input);
     BraggIntegrationSettings& Integrator(IntegratorMode input);
     BraggIntegrationSettings& StillPartiality(bool input);
+    BraggIntegrationSettings& BackgroundTrimFraction(float input);
 
 
     [[nodiscard]] IntegratorMode GetIntegrator() const;
@@ -41,4 +48,5 @@ public:
 
     [[nodiscard]] float GetMinimumSigmaInRegardsToI() const;
     [[nodiscard]] bool GetStillPartiality() const;
+    [[nodiscard]] float GetBackgroundTrimFraction() const;
 };

docs/CHANGELOG.md

@@ -23,6 +23,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * rugnux: Write anomalous data as a standard CCP4 anomalous MTZ - one row per reflection with `IMEAN`, `I(+)`/`I(-)` (and the matching `F`/`F(+)`/`F(-)` amplitudes) - instead of two `IMEAN` rows per Bijvoet pair, so `aimless`/`mtz2sca`/ANODE read it directly. The non-anomalous MTZ output is unchanged.
 * rugnux: Always write the merged reflections as both `<prefix>.mtz` and `<prefix>.cif` (each has its uses downstream - MTZ for the CCP4/phenix tools, mmCIF for deposition). The `--scaling-output` format selector and the plain-text `.hkl` output are removed.
 * rugnux: Add a detector-plane **modulation** (flat-field) correction surface to rotation scaling - a smooth multiplicative factor over the detector position where a reflection lands, fitted against the merged reference (symmetry-equivalents land at different positions as the crystal rotates) and cross-validated. It corrects detector-response/geometric systematics that inflate R-meas; on the rotation battery it lowers R-meas by several to tens of percent (e.g. lysoC 23->16%, lyso_2 47->28%, EcwtAL500 53->28%) with CC1/2 held or improved and the anomalous signal preserved. On by default with decay/absorption (`--no-scaling-corrections` disables all). The correction-surface cross-validation now scores a sigma-independent R-meas-like agreement instead of a studentized chi^2, so a surface can no longer pass by reshaping the sigmas (which also hardens the absorption surface against a mis-indexed-data regression).
+* Bragg integration: The local background under each Bragg spot is now a **symmetric trimmed mean** of the background ring rather than the plain mean, on by default for monochromatic/rotation data (`--background-trim <f>`, default 0.10; 0 restores the plain mean). The plain ring mean reads high because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided contaminants, so it over-subtracts - and since a weak intensity is a small difference of large numbers, that bias is fractionally largest at the resolution edge. Dropping the lowest and highest fraction of ring pixels removes it while leaving a clean Poisson ring essentially unchanged. On the rotation battery `<I/sigma>` improved on every crystal (median +50%) and resolution-edge R-meas fell several-fold (e.g. lyso_ref 1.0 A 108->43%); the high-resolution CC1/2 is rescued where the plain mean had collapsed it (e.g. a thaumatin case regains ~0.5 A) at a small CC1/2 cost in already-clean shells. Implemented identically in the CPU and GPU integrators; stills keep their existing high-side sigma-clip (§9.2 of CPU_DATA_ANALYSIS.md). The dead, never-read `Reflection.completeness` field was removed.
 
 ### 1.0.0-rc.159
 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use.

docs/CPU_DATA_ANALYSIS.md

@@ -446,7 +446,9 @@ $
 $
 with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below.
 
-For the **profile-fit path on broadband (still) data**, the background mean is additionally computed with a single high-outlier reject (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, then recompute): a bandwidth-streaked high-resolution spot or a close neighbour can leak into the ring and bias the mean high, over-subtracting and driving weak high-resolution intensities negative. A clean Poisson background is essentially unchanged by the cut. The reject is **not** applied to plain box summation (`--integrator boxsum`) or to monochromatic/rotation data.
+**Trimmed-mean background (rotation, default on).** On monochromatic/rotation data the ring mean $\hat{b}$ is by default replaced with a **symmetric trimmed mean**: the ring pixels are sorted, the lowest and highest fraction $f$ are dropped, and the central $(1-2f)$ are averaged ($f=0.10$ by default, `--background-trim`; $f=0$ restores the plain mean). 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. The plain mean reads high there because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants; dropping the extreme ring pixels removes that bias while a clean Poisson ring is essentially unchanged. In practice this lowers the resolution-edge $R_\text{meas}$ several-fold, raises $\langle I/\sigma\rangle$, and rescues the high-resolution CC$_{1/2}$ on data where the plain mean had collapsed it (at a small CC$_{1/2}$ cost in already-clean shells, from the slightly higher variance of the trimmed estimate). It is applied to the shared background used by both the box sum and the profile fit.
+
+For the **profile-fit path on broadband (still) data**, the trimmed mean is *not* used; instead the background mean is computed with a single high-outlier reject (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, then recompute): a bandwidth-streaked high-resolution spot or a close neighbour can leak into the ring and bias the mean high, over-subtracting and driving weak high-resolution intensities negative. A clean Poisson background is essentially unchanged by the cut. Neither robustification is applied to plain box summation (`--integrator boxsum`).
 
 ### 9.3 Profile-fitted extraction (default)
 

image_analysis/bragg_integration/BraggIntegrationEngine.cpp

@@ -68,6 +68,10 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe
     beam_y = geom.GetBeamY_pxl();
     use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0);
 
+    // Trimmed-mean background applies to monochromatic (rotation) data; broadband (stills) keep the tuned
+    // high-side sigma-clip, so the trim is forced off there. The fraction itself comes from the settings.
+    bkg_trim = broadband ? 0.0f : settings.GetBackgroundTrimFraction();
+
     polarization = experiment.GetPolarizationFactor();
 }
 

image_analysis/bragg_integration/BraggIntegrationEngine.h

@@ -97,6 +97,10 @@ protected:
     double c_radial;             // radial variance coefficient of tan^2(2theta): parallax + capture
     double F_px;                 // detector distance expressed in pixels
     float  beam_x, beam_y;
+    // Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings): the configured
+    // fraction for monochromatic (rotation) data, forced to 0 for broadband (stills, which keep their
+    // high-side sigma-clip). 0 = plain ring mean. Read by both the CPU and GPU engines.
+    float  bkg_trim = 0.0f;
 
     DiffractionGeometry geom;    // kept for the per-reflection polarization correction
     std::optional<float> polarization;

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -7,6 +7,7 @@
 #include <cmath>
 #include <cstdint>
 #include <limits>
+#include <vector>
 
 #include "../../common/CompressedImage.h"
 #include "../../common/JFJochException.h"
@@ -56,6 +57,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
     const bool do_clip = apply_bkg_clip && mode != IntegratorMode::BoxSum;
 
+    // Symmetric trimmed-mean background fraction (BraggIntegrationSettings, rugnux --background-trim):
+    // replaces the r2..r3 ring MEAN with an f-trimmed mean, robust to the high-side contamination
+    // (neighbour wings, tails) that biases the plain mean up and makes it over-subtract weak high-angle
+    // reflections. Already forced to 0 for broadband/stills by the base ctor (they keep the clip below).
+    const double bkg_trim_frac = bkg_trim;
+
     auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); };
 
     // --- Reflection mask: mark the r2 signal disk of every predicted reflection so a neighbour's
@@ -82,6 +89,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     };
     std::vector<Rough> rough(npredicted);
     double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
+    std::vector<int32_t> bkg_vals;   // reused per reflection for the trimmed-mean background (idea 1)
 
     for (size_t i = 0; i < npredicted; ++i) {
         const auto &r = predicted[i];
@@ -94,6 +102,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0;
         double bkg_sum = 0.0;
         int n_bkg = 0;
+        bkg_vals.clear();
         for (int y = y0; y <= y1; ++y)
             for (int x = x0; x <= x1; ++x) {
                 const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
@@ -109,15 +118,28 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     if (refl_mask[y * W + x]) continue;
                     if (!valid(px)) continue;
                     bkg_sum += static_cast<double>(px);
+                    if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
                     ++n_bkg;
                 }
             }
 
         if (n_inner_valid == n_inner && n_bkg > 5) {
             out.bkg = bkg_sum / n_bkg;
-            // One high-outlier sigma-clip pass on the background ring (stills-only): reject pixels
-            // above mean + 3*sqrt(mean) to strip a bandwidth-streaked neighbour that biases the mean.
-            if (do_clip) {
+            if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5) {
+                // Symmetric trimmed mean over the background ring (idea 1): drop the lowest and highest
+                // bkg_trim_frac of the pixels, average the rest. Robust to the high-side contamination
+                // that biases the plain ring mean and makes it over-subtract at high resolution.
+                std::sort(bkg_vals.begin(), bkg_vals.end());
+                const size_t lo = static_cast<size_t>(bkg_vals.size() * bkg_trim_frac);
+                const size_t hi = bkg_vals.size() - lo;
+                if (hi > lo) {
+                    double s = 0.0;
+                    for (size_t t = lo; t < hi; ++t) s += bkg_vals[t];
+                    out.bkg = s / static_cast<double>(hi - lo);
+                }
+            } else if (do_clip) {
+                // One high-outlier sigma-clip pass on the background ring (stills-only): reject pixels
+                // above mean + 3*sqrt(mean) to strip a bandwidth-streaked neighbour that biases the mean.
                 const double thr = out.bkg + 3.0 * std::sqrt(std::max(out.bkg, 1.0));
                 double s = 0.0;
                 int n = 0;

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -26,10 +26,14 @@ struct BraggGpuParams {
     float c_radial;
     float F_px;
     float beam_x, beam_y;
+    float bkg_trim;       // idea 1: symmetric trimmed-mean background fraction (0 = plain ring mean)
 };
 
 __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; }
 
+// idea 1: max annulus pixels held in shared memory for the per-block trimmed-mean sort (else flat mean).
+constexpr int BKG_TRIM_MAX = 512;
+
 // --- Mark the r2 signal disk of every predicted reflection (race-free: all writes are 1). ---
 __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, BraggGpuParams p, int n) {
     const int i = blockIdx.x;
@@ -111,8 +115,57 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
     }
     __syncthreads();
 
+    // Trimmed-mean background (idea 1): collect the annulus into shared memory, bitonic-sort the block, and
+    // average the middle (1 - 2*bkg_trim) fraction - robust to the high-side contamination that biases the
+    // plain ring mean. Falls back to the flat mean when the ring exceeds the shared buffer.
+    __shared__ int s_bvals[BKG_TRIM_MAX];
+    __shared__ int s_bn;
+    if (threadIdx.x == 0) s_bn = 0;
+    __syncthreads();
+    const bool do_trim = s_accept && p.bkg_trim > 0.0f && s_nbkg > 5 && s_nbkg <= BKG_TRIM_MAX;
+    if (do_trim) {
+        for (int t = threadIdx.x; t < area; t += blockDim.x) {
+            const int x = x0 + t % bw, y = y0 + t / bw;
+            const float ddx = (float) x - cx, ddy = (float) y - cy;
+            const float d2 = ddx * ddx + ddy * ddy;
+            if (!(d2 >= p.r2_sq && d2 < p.r3_sq)) continue;
+            if (mask[y * p.W + x]) continue;
+            const int32_t px = img[y * p.W + x];
+            if (!valid(px)) continue;
+            const int slot = atomicAdd(&s_bn, 1);
+            if (slot < BKG_TRIM_MAX) s_bvals[slot] = px;
+        }
+        __syncthreads();
+        const int nb = min(s_bn, BKG_TRIM_MAX);
+        int n2 = 1;
+        while (n2 < nb) n2 <<= 1;
+        for (int t = threadIdx.x + nb; t < n2; t += blockDim.x) s_bvals[t] = INT32_MAX;   // pad to power of 2
+        __syncthreads();
+        for (int k = 2; k <= n2; k <<= 1)                       // ascending bitonic sort of s_bvals[0..n2)
+            for (int j = k >> 1; j > 0; j >>= 1) {
+                for (int idx = threadIdx.x; idx < n2; idx += blockDim.x) {
+                    const int ixj = idx ^ j;
+                    if (ixj > idx) {
+                        const bool up = ((idx & k) == 0);
+                        const int a = s_bvals[idx], b = s_bvals[ixj];
+                        if ((up && a > b) || (!up && a < b)) { s_bvals[idx] = b; s_bvals[ixj] = a; }
+                    }
+                }
+                __syncthreads();
+            }
+        if (threadIdx.x == 0) {
+            const int lo = (int) (nb * p.bkg_trim), hi = nb - lo;
+            if (hi > lo) {
+                double s = 0.0;
+                for (int t = lo; t < hi; ++t) s += (double) s_bvals[t];
+                s_bkg = s / (double) (hi - lo);
+            }
+        }
+        __syncthreads();
+    }
+
     // Second ring pass for the stills sigma-clip (re-reads the annulus; avoids storing bkg values).
-    if (s_accept && p.do_clip) {
+    if (s_accept && p.do_clip && !do_trim) {
         double c_l = 0.0; int cn_l = 0;
         for (int t = threadIdx.x; t < area; t += blockDim.x) {
             const int x = x0 + t % bw, y = y0 + t / bw;
@@ -437,6 +490,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
         .bw_sigma = static_cast<float>(bw_sigma), .c_radial = static_cast<float>(c_radial),
         .F_px = static_cast<float>(F_px),
         .beam_x = beam_x, .beam_y = beam_y,
+        .bkg_trim = bkg_trim,   // effective trim fraction (0 for stills), set by the base ctor from settings
     };
 
     // Pass A: reset accumulators, mask, then box-sum.

rugnux/rugnux_cli.cpp

@@ -111,6 +111,7 @@ void print_usage() {
     std::cout << "  Integration" << std::endl;
     std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
     std::cout << "   --integration-radius <r>        Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
+    std::cout << "   --background-trim <f>           Rotation: symmetric trimmed-mean fraction for the background ring (0<=f<0.5, default 0.10; 0 = plain mean). Removes the high-side bias that over-subtracts weak high-angle spots" << std::endl;
     std::cout << "   --integrator <txt>              Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
     std::cout << "   --still-partiality              Experimental: weight stills reflections by a Gaussian excitation-error partiality exp(-dist_ewald^2/2sigma^2) instead of treating each as a full" << std::endl;
     std::cout << "   -q, --azim-q-spacing <num>      Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;
@@ -151,6 +152,7 @@ enum {
     OPT_REFINE_GEOMETRY,
     OPT_BANDWIDTH,
     OPT_INTEGRATION_RADIUS,
+    OPT_BACKGROUND_TRIM,
     OPT_REJECT_OUTLIERS,
     OPT_REJECT_DELTA_CCHALF,
     OPT_REFERENCE_COLUMN,
@@ -252,6 +254,7 @@ static option long_options[] = {
     {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS},
     {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
     {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
+    {"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
     {"integrator", required_argument, nullptr, OPT_INTEGRATOR},
     {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY},
     {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
@@ -523,6 +526,7 @@ int main(int argc, char **argv) {
     std::optional<double> resolution_cc_target;                      // --resolution-cc-target
     std::optional<int> report_shell_count;                           // --resolution-shells
     std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
+    std::optional<double> background_trim_arg;         // --background-trim: background-ring trimmed-mean fraction
     std::optional<IntegratorMode> integrator_mode;     // --integrator boxsum|gaussian|empirical
     bool still_partiality_flag = false;                // --still-partiality (experimental stills partiality)
     std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
@@ -786,6 +790,9 @@ int main(int argc, char **argv) {
         case OPT_INTEGRATION_RADIUS:
             integration_radius_arg = optarg;
             break;
+        case OPT_BACKGROUND_TRIM:
+            background_trim_arg = parse_double_arg(optarg, "--background-trim", logger);
+            break;
         case OPT_INTEGRATOR:
             if (strcmp(optarg, "boxsum") == 0) integrator_mode = IntegratorMode::BoxSum;
             else if (strcmp(optarg, "gaussian") == 0) integrator_mode = IntegratorMode::ProfileGaussian;
@@ -1398,6 +1405,14 @@ int main(int argc, char **argv) {
                     : "profile (empirical)");
     }
 
+    if (background_trim_arg) {
+        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
+        bis.BackgroundTrimFraction(static_cast<float>(*background_trim_arg));
+        experiment.ImportBraggIntegrationSettings(bis);
+        logger.Info("Background-ring trimmed-mean fraction set to {:.2f} (rotation; stills keep the sigma-clip)",
+                    *background_trim_arg);
+    }
+
     if (still_partiality_flag) {
         BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
         bis.StillPartiality(true);