rugnux: set the integration radius from the crystal's own spot width

On rotation data the signal radius is now r1 = clamp(round(2*r80), 4, 6),
where r80 is the 80% encircled-flux radius of the crystal's own spots. The
background ring keeps its area (r3 = sqrt(r2^2 + 133)), so r1 = 4 is the
shipped default bit for bit and 26 of the 38 battery crystals come out
byte-identical.

The width had to be measured somewhere new. rugnux already has one -
shell_sigma2[].tan - but it is a second moment taken inside the r1 disk it
would be setting, and it saturates at r1/2, so feeding it back measures the
cap and not the crystal. SpotWidth instead measures encircled flux over an
aperture fixed for the whole file (14 px, normalised at 8), in the pre-scan,
from spots the finder already produces on the frames the beam-stop projection
already reads. It touches no integrator output and runs before the first
integration pass, so there is no loop, it costs no extra frame reads, and
both passes - including the space-group search, which runs in pass 1 - see
the same radius. A default run pays a median 1.9 s.

k = 2 is not fitted. For a Gaussian r80 = 1.794 sigma, so r1 = 2*r80 is
3.59 sigma, where the truncated second moment recovers 0.990 of sigma^2. The
new test checks the estimator returns 1.794 sigma on a known Gaussian.

Battery, 38 crystals, both arms run twice: the space group is identical on
all 38 and 35 agree with the reference in both arms. Per shell on the 12
crystals the rule moves, 6 win and 4 tie, with mean per-shell <I/sigma> up
30.6, 24.6, 15.8, 9.1, 8.0 and 5.1 per cent and R_meas down as much as 23.8.
Runtime is neutral - 19m13s against 22m00s warm.

One crystal is a real cost and is named in docs/RUGNUX.md with its
workaround: an I222 case that is simultaneously the widest-spot and among the
highest-mosaicity in the set loses 28.5% of its observations at unchanged
completeness, because at r1 = 6 its predicted reflection density leaves the
background ring too few clean pixels. No cheap guard separates it - its
predicted spacing is mid-table, larger than five crystals that survive r1 =
12 - and the guard that would, on the measured drop rate out of pass 1, needs
a diagnostic channel out of both integration engines and is not yet
validated.

This depends on 3ea120677: at the previous twin-law bound of 1.70 the
wider radius costs one crystal its 422, refused on an H ratio of 1.73 even
though every operator correlation confirms the point group.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P

docs/CHANGELOG.md

@@ -9,6 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * rugnux: a lattice centring the data could not test - the crystal was integrated on the primitive sub-cell, so the reflections the centring extinguishes were never measured - is marked `UNTESTED` in the space-group candidate table and, where it is adopted, is warned about as coming from the lattice metric rather than from the intensities.
 * rugnux: the space-group search prints the twin-law disagreement H for every operator it tested and the H ratio of the point group it adopted, on every run, instead of only when that ratio refuses a promotion.
 * rugnux: `--mode scale` works on a `_process.h5` whose space group came from re-seating the lattice; a file written before this stops with a message naming the two cells instead of failing inside the merge.
+* rugnux: on rotation data the integration signal radius is set from how wide the crystal's own spots are, measured in the pre-scan, instead of the fixed 4 px; `--adaptive-integration-radius=off` restores the fixed radius, and an explicit `--integration-radius` still overrides both.
 * rugnux reports diffraction anisotropy: the anisotropic deltaB and the diffraction limit along each
 * rugnux: `--mode scale` reports the Wilson B-factor estimate instead of `WILSON_B= nan`.
 * rugnux: `--export-unmerged` writes the integrated observations as `<prefix>_unmerged.mtz`, an unmerged MTZ readable by aimless, pointless, careless and `iotbx.merging_statistics`, in `--mode mx` and `--mode scale` alike; each rotation reflection's partials are summed into one full, and `--export-unmerged-partials` writes one row per image instead. Lattice-centring absences are not written; screw and glide absences are.

docs/RUGNUX.md

@@ -534,6 +534,7 @@ Integration:
 | --- | --- |
 | `--integrator <txt>` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) |
 | `--integration-radius <r>` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` |
+| `--adaptive-integration-radius[=on\|off]` | Set the signal radius `r1` from how wide this crystal's spots actually are (default: **on** for rotation, off for stills). `r1` is not the integration domain — that is the profile-fit grid — but it *is* the aperture the profile **width** is learned over, and a second moment over a disk of radius `a` saturates at `a²/4`, so at the shipped `r1 = 4` the learned σ can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. The pre-scan reads `r80`, the radius holding 80 % of a spot's flux, off isolated strong spots over a fixed 14 px aperture that owes nothing to `r1`, fits it against `1/d` and evaluates it at 5 Å; then `r1 = clamp(round(2·r80), 4, 6)`, `r2 = r1 + 2`, and `r3` is taken so the `r2..r3` background ring keeps the area it has at the default `4,6,13`. The ceiling of 6 is pattern density: `r2` also drives the neighbour-ownership radius and the ring's inner edge, and past it a dense pattern starts losing reflections whose ring falls below six clean pixels. Ignored when `--integration-radius` is given. Over the rotation regression battery it moves 12 of 38 crystals and leaves the merged intensities of the other 26 unchanged; where it moves them, per-shell ⟨I/σ⟩ improves by up to 31 % and R_meas by up to 24 %. One crystal with unusually broad spots *and* a large mosaicity is a known cost: its background rings are starved by the wider `r2` and it loses a quarter of its observations, so if multiplicity collapses on a very mosaic crystal, switch this off |
 | `--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 |

rugnux/CMakeLists.txt

@@ -14,6 +14,8 @@ ADD_LIBRARY(Rugnux STATIC
         ModelValidation.h
         ResultReport.cpp
         ResultReport.h
+        SpotWidth.cpp
+        SpotWidth.h
 )
 
 TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi)

rugnux/Rugnux.cpp

@@ -3,6 +3,7 @@
 
 #include "Rugnux.h"
 #include "ModelValidation.h"
+#include "SpotWidth.h"
 
 #include <algorithm>
 #include <atomic>
@@ -339,7 +340,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
                                          [](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; });
     const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow;
     const bool want_beam_center = config_.estimate_beam_center && !beam_center_placed_;
-    if (!want_shadow && !want_beam_center)
+    // The spot width is a property of the crystal, not of the geometry, so the second pass reuses the
+    // radius the first one measured rather than re-reading the frames.
+    const bool want_width = config_.adaptive_integration_radius && !spot_width_measured_;
+    if (!want_shadow && !want_beam_center && !want_width)
         return;
 
     // The two consumers want different frames and each gets its own set. The shadow is built from
@@ -381,15 +385,16 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     // tens of seconds once the beam centre asks for its own frames on top of the projection.
     if (observer)
         observer->OnPhase(want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop")
-                                      : "Beam centre");
+                                      : (want_beam_center ? "Beam centre" : "Spot width"));
 
     ShadowFinder finder(experiment_, pixel_mask_);
 
-    // The spot symmetry is read in the same pass as the projection rather than in a pre-pass of its
-    // own. The image is preprocessed a second time for it: the shadow is accumulated from raw counts
-    // and the spot finder works on the converted image.
+    // The spot symmetry and the spot width are read in the same pass as the projection rather than in
+    // a pre-pass of their own. The image is preprocessed a second time for them: the shadow is
+    // accumulated from raw counts and the spot finder works on the converted image.
+    const bool want_spots = want_spot_symmetry || want_width;
     std::unique_ptr<AzimuthalIntegrationMapping> prescan_mapping;
-    if (want_spot_symmetry)
+    if (want_spots)
         prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(experiment_, pixel_mask_);
 
     // Everything reading one frame needs that cannot be shared between threads: its own scratch for
@@ -403,7 +408,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     };
     const auto make_worker = [&] {
         PreScanWorker w;
-        if (want_spot_symmetry) {
+        if (want_spots) {
             w.preprocessor = std::make_unique<ImagePreprocessorCPU>(experiment_, pixel_mask_);
             w.preprocessed = std::make_unique<ImagePreprocessorBuffer>(experiment_.GetPixelsNum());
             if (config_.spot_finding.adaptive_threshold)
@@ -415,11 +420,13 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         return w;
     };
 
-    // The spots of one image. Taken as a step of its own because the sparse-sweep pass below reads
-    // images the shadow does not. The frame they belong to is stamped by the caller, which is what
-    // keeps the pool independent of the order the workers finish in.
+    // The spots of one image, and what the two consumers take off them. Taken as a step of its own
+    // because the sparse-sweep pass below reads images the shadow does not. The frame they belong to
+    // is stamped by the caller, which is what keeps the pool independent of the order the workers
+    // finish in.
     const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx,
-                                std::vector<BeamCenterSpot> &out) {
+                                bool for_beam_center, std::vector<BeamCenterSpot> &out,
+                                bool for_width, std::vector<spot_width::FluxCurve> &curves) {
         try {
             w.preprocessor->Analyze(*w.preprocessed,
                                     image.GetUncompressedPtr(w.decompression_buffer), image.GetMode());
@@ -428,6 +435,13 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             return;
         }
         auto spots = w.spot_finder->Run(*w.preprocessed, config_.spot_finding);
+        // The width goes first: it wants the whole spot list, to know which spots are isolated.
+        if (for_width)
+            MeasureSpotFluxCurves(*w.preprocessed, static_cast<int>(experiment_.GetXPixelsNumConv()),
+                                  static_cast<int>(experiment_.GetYPixelsNumConv()),
+                                  experiment_.GetDiffractionGeometry(), spots, curves);
+        if (!for_beam_center)
+            return;
         // The strongest of a crowded frame: the symmetry is over-determined either way, and the
         // matching is quadratic in the spots of one frame.
         if (spots.size() > BEAM_CENTER_SPOTS_PER_IMAGE) {
@@ -444,6 +458,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
 
     std::vector<float> frame_angle_deg;
     std::vector<BeamCenterSpot> beam_center_spots;
+    std::vector<spot_width::FluxCurve> width_curves;
 
     // Read the sample on several workers. The reader serialises on the HDF5 lock, but the
     // decompression, the projection and the spot finding - which is all of the cost on a large
@@ -453,6 +468,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     {
         const std::vector<int> ordinals(sample.begin(), sample.end());
         std::vector<std::vector<BeamCenterSpot>> spots_of(ordinals.size());
+        std::vector<std::vector<spot_width::FluxCurve>> curves_of(ordinals.size());
         // A frame joins the pool only if it could be read, as it did when this was a serial loop.
         std::vector<char> spot_read(ordinals.size(), 0);
         const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
@@ -484,14 +500,22 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
                     if (shadow_set.contains(ordinal))
                         finder.AddImage(msg, w.shadow_buffer, t);
 
-                    if (spot_set.contains(ordinal)) {
-                        spot_read[i] = 1;
-                        find_spots(w, msg.image, image_idx, spots_of[i]);
-                    }
+                    // The width is measured on the projection's own even spread over the sweep, so
+                    // asking for a beam centre cannot move it either.
+                    const bool for_beam_center = spot_set.contains(ordinal);
+                    const bool for_width = want_width && shadow_set.contains(ordinal);
+                    if (for_beam_center) spot_read[i] = 1;
+                    if (for_beam_center || for_width)
+                        find_spots(w, msg.image, image_idx, for_beam_center, spots_of[i],
+                                   for_width, curves_of[i]);
                 }
             }));
         for (auto &f : futures) f.get();
 
+        // Sample order, so the pooled width does not depend on how the workers interleaved.
+        for (auto &c : curves_of)
+            width_curves.insert(width_curves.end(), c.begin(), c.end());
+
         // Frame numbering follows the sample order, exactly as the serial read did.
         for (size_t i = 0; i < ordinals.size(); i++) {
             if (!spot_read[i]) continue;
@@ -507,6 +531,32 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         }
     }
 
+    // The integration radius from the recorded spot width. r1 is the aperture the integrator learns
+    // the profile width over and a second moment over a disk of radius a saturates at a^2/4, so at the
+    // shipped r1 = 4 the learned sigma cannot exceed 2 px and a wider spot is fitted with a profile
+    // the model cannot represent. r2 follows r1, and r3 is taken so that the r2..r3 background ring
+    // keeps the area it has at the shipped 4/6/13 - a ring that shrinks with the disk is what makes a
+    // bare --integration-radius worse than the default it replaces.
+    if (want_width) {
+        spot_width_measured_ = true;
+        const auto r80 = spot_width::R80AtReference(width_curves);
+        if (!r80) {
+            logger.Info("Spot width: not measurable on {} spots, keeping the integration radius at "
+                        "r1={:.1f}", width_curves.size(),
+                        experiment_.GetBraggIntegrationSettings().GetR1());
+        } else {
+            const float r1 = spot_width::R1ForWidth(*r80);
+            const float r2 = r1 + 2.0f;
+            const float r3 = std::sqrt(r2 * r2 + 133.0f);
+            BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
+            bis.R1(r1).R2(r2).R3(r3);
+            experiment_.ImportBraggIntegrationSettings(bis);
+            logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots) => integration radii "
+                        "r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A,
+                        width_curves.size(), r1, r2, r3);
+        }
+    }
+
     if (finder.GetFrameCount() == 0) {
         logger.Warning("Pre-scan: no image could be read. Skipping.");
         return;
@@ -591,7 +641,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
                         }
                         if (img) {
                             extra_read[i] = 1;
-                            find_spots(w, img->image, image_idx, extra_spots[i]);
+                            // The extra frames go to the beam centre alone; the width has already
+                            // been measured on the projection's own sample.
+                            find_spots(w, img->image, image_idx, true, extra_spots[i],
+                                       false, width_curves);
                         }
                     }
                 }));
@@ -1140,7 +1193,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
 
     // After any geometry refinement, so the shadow is found about the beam centre actually used,
     // and before the azimuthal mapping and the output mask below, which both read pixel_mask_.
-    if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center)
+    if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center
+        || config_.adaptive_integration_radius)
         PreScan(start_image, images_to_process,
                 config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES), observer);
 

rugnux/Rugnux.h

@@ -92,6 +92,14 @@ struct ProcessConfig {
     bool estimate_beam_center = false;
     bool fit_spindle = false;
 
+    // Adaptive integration radius (--adaptive-integration-radius). When set, the pre-scan measures how
+    // wide the recorded spots are - r80, the radius holding 80 % of a spot's flux, over a fixed
+    // aperture that owes nothing to the integrator (SpotWidth.h) - and sets r1 from it. r1 is the
+    // aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius
+    // a saturates at a^2/4, so at the shipped r1 = 4 a spot wider than sigma = 2 px is fitted with a
+    // profile the model cannot represent. Ignored where the radii were set by hand.
+    bool adaptive_integration_radius = false;
+
     // Rotation two-pass geometry post-refinement (FullAnalysis, rotation only; on by default in the rugnux
     // CLI, --rotation-no-postrefine disables it). When set, a first pass integrates and post-refines the
     // detector distance + beam (from the observed spot positions) and the cell scale + rotation axis (from the
@@ -298,11 +306,16 @@ class Rugnux {
     // a starting point for indexing, so it must not overwrite either of those.
     bool beam_center_placed_ = false;
 
-    // Pre-scan: read a spread sample of frames and take two things off them - the shadow of the beam
-    // stop and its holder, added to the pixel mask (config_.detect_beam_stop), and the beam centre
+    // The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the
+    // two-pass rotation run measures it once and both passes integrate at the same radius.
+    bool spot_width_measured_ = false;
+
+    // Pre-scan: read a spread sample of frames and take three things off them - the shadow of the beam
+    // stop and its holder, added to the pixel mask (config_.detect_beam_stop), the beam centre
     // (config_.estimate_beam_center), from the symmetry of the spots where the sweep reaches half a
-    // turn and from the isotropy of the scattered background where it does not. Either may be asked
-    // for without the other.
+    // turn and from the isotropy of the scattered background where it does not, and the recorded spot
+    // width that sets the integration radius (config_.adaptive_integration_radius). Any of the three
+    // may be asked for without the others.
     void PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer);
 
     // Stills global geometry-refinement first pass (config_.refine_geometry): index a spread sample of

rugnux/SpotWidth.cpp

@@ -0,0 +1,242 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "SpotWidth.h"
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <utility>
+
+using namespace spot_width;
+
+namespace {
+
+// The engine reads pixels in the INT32_MIN(masked)/INT32_MAX(saturated) convention.
+inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; }
+
+// Nothing inside this radius of the beam centre: the beam stop and its halo are not spots.
+constexpr float MIN_BEAM_DISTANCE_PX = 60.0f;
+// A neighbour this close puts its own flux inside the aperture, which would read as extra width.
+constexpr float ISOLATION_PX = 28.0f;
+// Spots taken per resolution band per image, strongest first.
+constexpr int PER_BAND_PER_IMAGE = 40;
+// The r <= 4 px sum must be this many sigma above the background before the tail is believed.
+constexpr double SNR_MIN = 15.0;
+constexpr int R_CENTROID = 4;
+constexpr double MAX_CENTROID_OFFSET_PX = 2.0;
+// Spots needed before a band, and the crystal, are characterised at all.
+constexpr size_t MIN_SPOTS_PER_BAND = 15;
+constexpr size_t MIN_SPOTS_TOTAL = 20;
+
+// Resolution bands, A. The quota is per band, so a crystal is characterised over its whole range
+// and not wherever its strongest spots happen to sit.
+constexpr int N_BAND = 5;
+constexpr std::array<std::pair<float, float>, N_BAND> BANDS = {{
+        {2.0f, 3.0f}, {3.0f, 4.5f}, {4.5f, 7.0f}, {7.0f, 12.0f}, {12.0f, 30.0f}}};
+
+int band_of(float d_A) {
+    for (int b = 0; b < N_BAND; b++)
+        if (d_A >= BANDS[b].first && d_A < BANDS[b].second) return b;
+    return -1;
+}
+
+// The radius at which the curve reaches `frac`, linearly interpolated. prof[i] is the flux inside
+// radius i+1.
+float interpolate_radius(double frac, const std::array<float, R_MAX> &prof) {
+    if (prof[0] >= frac)
+        return prof[0] > 0.0f ? static_cast<float>(frac / prof[0]) : 1.0f;
+    for (int i = 1; i < R_MAX; i++)
+        if (prof[i] >= frac)
+            return static_cast<float>(i + (frac - prof[i - 1]) / (prof[i] - prof[i - 1]));
+    return static_cast<float>(R_MAX);
+}
+
+double median_of(std::vector<double> &v) {
+    if (v.empty()) return 0.0;
+    const size_t mid = v.size() / 2;
+    std::nth_element(v.begin(), v.begin() + mid, v.end());
+    const double hi = v[mid];
+    if (v.size() % 2 == 1) return hi;
+    return 0.5 * (hi + *std::max_element(v.begin(), v.begin() + mid));
+}
+
+}   // namespace
+
+void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height,
+                           const DiffractionGeometry &geometry,
+                           const std::vector<DiffractionSpot> &spots,
+                           std::vector<FluxCurve> &out) {
+    if (spots.empty()) return;
+
+    const float beam_x = geometry.GetBeamX_pxl(), beam_y = geometry.GetBeamY_pxl();
+
+    // Where every spot of this image sits, so isolation can be tested against all of them and not
+    // only against the ones that survive the gates below.
+    std::vector<Coord> centre(spots.size());
+    for (size_t i = 0; i < spots.size(); i++)
+        centre[i] = spots[i].RawCoord();
+
+    // Isolation on a grid of ISOLATION_PX cells: a neighbour within that distance is in this cell or
+    // one of the eight around it.
+    const int gw = static_cast<int>(width / ISOLATION_PX) + 1;
+    const int gh = static_cast<int>(height / ISOLATION_PX) + 1;
+    std::vector<std::vector<uint32_t>> cell(static_cast<size_t>(gw) * gh);
+    const auto cell_of = [&](const Coord &c) {
+        const int gx = std::clamp(static_cast<int>(c.x / ISOLATION_PX), 0, gw - 1);
+        const int gy = std::clamp(static_cast<int>(c.y / ISOLATION_PX), 0, gh - 1);
+        return std::pair<int, int>(gx, gy);
+    };
+    for (size_t i = 0; i < spots.size(); i++) {
+        const auto [gx, gy] = cell_of(centre[i]);
+        cell[static_cast<size_t>(gy) * gw + gx].push_back(static_cast<uint32_t>(i));
+    }
+    const auto isolated = [&](size_t i) {
+        const auto [gx, gy] = cell_of(centre[i]);
+        for (int y = std::max(0, gy - 1); y <= std::min(gh - 1, gy + 1); y++)
+            for (int x = std::max(0, gx - 1); x <= std::min(gw - 1, gx + 1); x++)
+                for (uint32_t j : cell[static_cast<size_t>(y) * gw + x]) {
+                    if (j == i) continue;
+                    if (std::hypot(centre[j].x - centre[i].x, centre[j].y - centre[i].y) < ISOLATION_PX)
+                        return false;
+                }
+        return true;
+    };
+
+    // Candidates that pass the geometric gates, by band, strongest first.
+    struct Candidate { size_t index; int64_t count; float d_A; };
+    std::array<std::vector<Candidate>, N_BAND> candidates;
+    for (size_t i = 0; i < spots.size(); i++) {
+        const Coord &c = centre[i];
+        const int cx = static_cast<int>(std::lround(c.x)), cy = static_cast<int>(std::lround(c.y));
+        if (cx < R_BKG_OUT || cy < R_BKG_OUT || cx >= width - R_BKG_OUT || cy >= height - R_BKG_OUT)
+            continue;
+        if (std::hypot(c.x - beam_x, c.y - beam_y) < MIN_BEAM_DISTANCE_PX) continue;
+        const float d_A = geometry.PxlToRes(c.x, c.y);
+        const int band = band_of(d_A);
+        if (band < 0) continue;
+        if (!isolated(i)) continue;
+        candidates[band].push_back({i, spots[i].Count(), d_A});
+    }
+
+    std::vector<double> ring;
+    for (int band = 0; band < N_BAND; band++) {
+        auto &cand = candidates[band];
+        const size_t take = std::min<size_t>(cand.size(), PER_BAND_PER_IMAGE);
+        std::partial_sort(cand.begin(), cand.begin() + take, cand.end(),
+                          [](const Candidate &a, const Candidate &b) { return a.count > b.count; });
+        for (size_t k = 0; k < take; k++) {
+            const Coord &c = centre[cand[k].index];
+            const int cx = static_cast<int>(std::lround(c.x)), cy = static_cast<int>(std::lround(c.y));
+
+            // The background under the spot, and a check that the whole aperture is readable: a hole
+            // in it removes flux from one radius and not another, which is exactly the shape this
+            // measures.
+            ring.clear();
+            bool readable = true;
+            for (int dy = -R_BKG_OUT; dy <= R_BKG_OUT && readable; dy++)
+                for (int dx = -R_BKG_OUT; dx <= R_BKG_OUT; dx++) {
+                    const int d2 = dx * dx + dy * dy;
+                    if (d2 > R_BKG_OUT * R_BKG_OUT) continue;
+                    const int32_t px = image[static_cast<size_t>(cy + dy) * width + (cx + dx)];
+                    if (!valid(px)) { readable = false; break; }
+                    if (d2 >= R_BKG_IN * R_BKG_IN) ring.push_back(px);
+                }
+            if (!readable || ring.size() < 20) continue;
+            const size_t n_ring = ring.size();
+            const double bkg = median_of(ring);
+
+            // Flux and centroid over the r <= 4 px core, then the signal-to-noise gate. A weak spot's
+            // tail is background, and an encircled-flux curve built on it measures the background.
+            double core = 0.0, mx = 0.0, my = 0.0;
+            int n_core = 0;
+            for (int dy = -R_CENTROID; dy <= R_CENTROID; dy++)
+                for (int dx = -R_CENTROID; dx <= R_CENTROID; dx++) {
+                    if (dx * dx + dy * dy > R_CENTROID * R_CENTROID) continue;
+                    const double v = image[static_cast<size_t>(cy + dy) * width + (cx + dx)] - bkg;
+                    core += v;
+                    mx += v * dx;
+                    my += v * dy;
+                    ++n_core;
+                }
+            if (core <= 0.0) continue;
+            const double noise = std::sqrt(core + n_core * std::max(bkg, 0.05)
+                                                  * (1.0 + static_cast<double>(n_core) / n_ring));
+            if (core / noise < SNR_MIN) continue;
+            mx /= core;
+            my /= core;
+            if (std::abs(mx) > MAX_CENTROID_OFFSET_PX || std::abs(my) > MAX_CENTROID_OFFSET_PX)
+                continue;
+
+            // The encircled flux about that centroid, out to the fixed aperture.
+            FluxCurve curve;
+            curve.d_A = cand[k].d_A;
+            for (int dy = -R_MAX; dy <= R_MAX; dy++)
+                for (int dx = -R_MAX; dx <= R_MAX; dx++) {
+                    const double rc = std::hypot(dx - mx, dy - my);
+                    if (rc > R_MAX) continue;
+                    const double v = image[static_cast<size_t>(cy + dy) * width + (cx + dx)] - bkg;
+                    for (int t = std::max(1, static_cast<int>(std::ceil(rc))); t <= R_MAX; t++)
+                        curve.c[t - 1] += static_cast<float>(v);
+                }
+            if (!(curve.c[R_NORM - 1] > 0.0f) || !(curve.c[R_MAX - 1] > 0.0f)) continue;
+            const float norm = curve.c[R_NORM - 1];
+            for (float &v : curve.c) v /= norm;
+            out.push_back(curve);
+        }
+    }
+}
+
+std::optional<float> spot_width::R80AtReference(const std::vector<FluxCurve> &curves) {
+    if (curves.size() < MIN_SPOTS_TOTAL) return std::nullopt;
+
+    // One point per band: the median curve of the band, the radius it holds 80 % of its flux at, and
+    // the median resolution it was measured at.
+    struct Point { double inv_d; double r80; double weight; };
+    std::vector<Point> points;
+    std::vector<double> values, band_d;
+    for (int b = 0; b < N_BAND; b++) {
+        band_d.clear();
+        for (const auto &c : curves)
+            if (c.d_A >= BANDS[b].first && c.d_A < BANDS[b].second) band_d.push_back(c.d_A);
+        if (band_d.size() < MIN_SPOTS_PER_BAND) continue;
+        std::array<float, R_MAX> profile{};
+        for (int t = 0; t < R_MAX; t++) {
+            values.clear();
+            for (const auto &c : curves)
+                if (c.d_A >= BANDS[b].first && c.d_A < BANDS[b].second) values.push_back(c.c[t]);
+            profile[t] = static_cast<float>(median_of(values));
+        }
+        const double d_med = median_of(band_d);
+        if (d_med <= 0.0) continue;
+        points.push_back({1.0 / d_med, interpolate_radius(0.8, profile),
+                          static_cast<double>(band_d.size())});
+    }
+    if (points.empty()) return std::nullopt;
+    if (points.size() == 1) return static_cast<float>(points[0].r80);
+
+    // The mosaic contribution to the detector footprint grows as 1/d, so r80 is linear in 1/d.
+    double sw = 0.0, sx = 0.0, sxx = 0.0, sy = 0.0, sxy = 0.0;
+    for (const auto &p : points) {
+        sw += p.weight;
+        sx += p.weight * p.inv_d;
+        sxx += p.weight * p.inv_d * p.inv_d;
+        sy += p.weight * p.r80;
+        sxy += p.weight * p.inv_d * p.r80;
+    }
+    const double det = sw * sxx - sx * sx;
+    double value = sy / sw;
+    if (std::abs(det) > 1e-12) {
+        const double c1 = (sw * sxy - sx * sy) / det;
+        value = (sy - c1 * sx) / sw + c1 / D_REF_A;
+    }
+    // Never extrapolate outside what the bands actually measured.
+    double lo = std::numeric_limits<double>::max(), hi = 0.0;
+    for (const auto &p : points) { lo = std::min(lo, p.r80); hi = std::max(hi, p.r80); }
+    return static_cast<float>(std::clamp(value, 0.8 * lo, 1.25 * hi));
+}
+
+float spot_width::R1ForWidth(float r80) {
+    return std::clamp(std::round(2.0f * r80), 4.0f, 6.0f);
+}

rugnux/SpotWidth.h

@@ -0,0 +1,73 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+// =============================================================================
+// SpotWidth - how wide the recorded spots actually are, in detector pixels
+// =============================================================================
+//
+// The integrator's signal radius r1 is not the integration domain (the profile fit runs over a
+// (2*ceil(r2)+1)^2 grid), but it IS the aperture the profile WIDTH is learned over, and a second
+// moment taken over a disk of radius a saturates at a^2/4 - so at the shipped r1 = 4 the learned
+// Gaussian sigma can never exceed 2 px whatever the spot is. A crystal whose spots are wider than
+// that is integrated with a profile the model cannot represent, and it is exactly those crystals
+// that gain from a larger r1. Measured over the rotation battery, one threshold on the recorded
+// width separates the two populations cleanly.
+//
+// This measures that width WITHOUT asking the integrator, so the answer can be used to set r1:
+// r80, the radius holding 80 % of a spot's flux, is read off the encircled-flux curve of isolated
+// strong spots over a FIXED 14 px aperture normalised at 8 px. Nothing here depends on r1, r2 or
+// r3, or on any quantity the integrator computes; the aperture is a constant of this file. The
+// curve is built about the spot's own flux-weighted centroid, no shape is assumed, and the
+// candidates are stratified by resolution - the strongest spots of a weak crystal sit at high
+// angle and of a strong one at low angle, and the mosaic footprint grows as 1/d, so an unstratified
+// sample compares crystals at different resolutions.
+// =============================================================================
+
+#include <array>
+#include <optional>
+#include <vector>
+
+#include "../common/DiffractionGeometry.h"
+#include "../common/DiffractionSpot.h"
+#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
+
+namespace spot_width {
+// The encircled-flux curve runs to R_MAX px and is normalised at R_NORM, clear of the noisy tail.
+// R_NORM = 8 also fixes what this estimator can see: a spot needs r80 <= 8 px to be measured at
+// all, which is five times the width the decision below is taken at.
+constexpr int R_MAX = 14;
+constexpr int R_NORM = 8;
+// The background under a spot is the median of the R_MAX+2 .. R_MAX+6 px ring around it.
+constexpr int R_BKG_IN = R_MAX + 2;
+constexpr int R_BKG_OUT = R_MAX + 6;
+// Resolution the widths are compared at. Inside the measured range of essentially every crystal.
+constexpr float D_REF_A = 5.0f;
+
+// One spot's encircled flux C(1), C(2), ... C(R_MAX), divided by C(R_NORM), with the resolution it
+// was recorded at.
+struct FluxCurve {
+    float d_A = 0.0f;
+    std::array<float, R_MAX> c{};
+};
+
+// r80 at D_REF_A from a pool of curves, or nothing where too few spots were measurable. Fitted as
+// r80(d) = c0 + c1/d over resolution bands, weighted by the spots in each, so crystals whose strong
+// spots sit at different resolutions are compared like with like.
+[[nodiscard]] std::optional<float> R80AtReference(const std::vector<FluxCurve> &curves);
+
+// The signal radius the measured width asks for: r1 = clamp(round(2*r80), 4, 6). The floor is the
+// shipped default and the ceiling is set by pattern density - r2 = r1 + 2 also drives the
+// neighbour-ownership radius and the inner edge of the background ring, and past 6 px a dense
+// pattern starts losing reflections whose ring falls below six clean pixels.
+[[nodiscard]] float R1ForWidth(float r80);
+}   // namespace spot_width
+
+// Every isolated, strong, fully readable spot of one image, as a normalised encircled-flux curve,
+// appended to `out`. Called on one worker's own vector, so the pool is independent of the order the
+// workers finish in.
+void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height,
+                           const DiffractionGeometry &geometry,
+                           const std::vector<DiffractionSpot> &spots,
+                           std::vector<spot_width::FluxCurve> &out);

rugnux/rugnux_cli.cpp

@@ -171,6 +171,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 << "   --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default: on for rotation, off for stills). r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl;
     std::cout << "   --integration-stencil <k>       Push the r2..r3 background ring out by k times the beam's radial streak (bandwidth*Rpx), per reflection (default 0 = a fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially - the r1 signal box stays a circle. Needs --bandwidth: on a monochromatic beam the streak is zero and this does nothing" << std::endl;
     std::cout << "   --integration-high-resolution <num> High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl;
     std::cout << "   --max-hkl <n>                   Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl;
@@ -231,6 +232,7 @@ enum {
     OPT_NO_FIT_SPINDLE,
     OPT_BANDWIDTH,
     OPT_INTEGRATION_RADIUS,
+    OPT_ADAPTIVE_INTEGRATION_RADIUS,
     OPT_INTEGRATION_STENCIL,
     OPT_BACKGROUND_TRIM,
     OPT_OVERLAP,
@@ -361,6 +363,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},
+    {"adaptive-integration-radius", optional_argument, nullptr, OPT_ADAPTIVE_INTEGRATION_RADIUS},
     {"integration-stencil", required_argument, nullptr, OPT_INTEGRATION_STENCIL},
     {"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
     {"overlap", required_argument, nullptr, OPT_OVERLAP},
@@ -667,6 +670,8 @@ static int RunRugnux(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;
+    // --adaptive-integration-radius: unset means the per-workflow default below (on for rotation).
+    std::optional<bool> adaptive_integration_radius_arg;
     std::optional<std::string> integration_stencil_arg; // --integration-stencil: ring elongation, in sigma
     std::optional<double> background_trim_arg;         // --background-trim: background-ring trimmed-mean fraction
     std::optional<OverlapMode> overlap_arg;            // --overlap: treatment of shared signal pixels
@@ -1038,6 +1043,9 @@ static int RunRugnux(int argc, char **argv) {
         case OPT_INTEGRATION_RADIUS:
             integration_radius_arg = optarg;
             break;
+        case OPT_ADAPTIVE_INTEGRATION_RADIUS:
+            adaptive_integration_radius_arg = !optarg || std::string(optarg) != "off";
+            break;
         case OPT_INTEGRATION_STENCIL:
             integration_stencil_arg = optarg;
             break;
@@ -2058,6 +2066,13 @@ static int RunRugnux(int argc, char **argv) {
         logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=12.0 (override with --integration-radius)");
     }
 
+    // Adaptive signal radius: measured in the pre-scan and applied there (Rugnux::PreScan). Radii the
+    // user set by hand always win, and stills already default to a wide box, so the rule is on for
+    // rotation only.
+    const bool adaptive_integration_radius =
+            !integration_radius_arg
+            && adaptive_integration_radius_arg.value_or(rotation_indexing);
+
     if (forced_prediction_mosaicity_arg) {
         BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
         bis.ForcedPredictionMosaicity_deg(static_cast<float>(*forced_prediction_mosaicity_arg));
@@ -2210,6 +2225,7 @@ static int RunRugnux(int argc, char **argv) {
     config.detect_beam_stop = detect_beam_stop;
     config.estimate_beam_center = estimate_beam_center;
     config.fit_spindle = fit_spindle;
+    config.adaptive_integration_radius = adaptive_integration_radius;
     config.rotation_postrefine_geometry = rotation_postrefine_geometry;
     config.rotation_scale = rotation_scale;
     config.rotation_indexing_image_count = rotation_indexing_image_count;

tests/RugnuxTest.cpp

@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include <catch2/catch_all.hpp>
+#include <cmath>
 #include <filesystem>
 
 #include "../common/DiffractionExperiment.h"
@@ -10,6 +11,7 @@
 #include "../reader/JFJochHDF5Reader.h"
 #include "../rugnux/Rugnux.h"
 #include "../rugnux/RugnuxCommandLine.h"
+#include "../rugnux/SpotWidth.h"
 
 namespace {
     // Write a small VDS dataset of `n` flat images and return nothing (prefix_master.h5 +
@@ -183,3 +185,59 @@ TEST_CASE("RugnuxCommandLine_AzInt", "[process]") {
     CHECK(cmd.find("--azim-min-q") != std::string::npos);
     CHECK(cmd.find("in.h5") != std::string::npos);
 }
+
+namespace {
+    // A field of identical round Gaussian spots on three rings, so that the width estimator sees
+    // several resolution bands with the same true width and its 1/d fit has to come back flat.
+    void PaintGaussianSpots(ImagePreprocessorBuffer &image, int w, double sigma, double total_counts,
+                            std::vector<DiffractionSpot> &spots) {
+        constexpr int BKG = 3;
+        for (size_t i = 0; i < image.size(); i++) image[i] = BKG;
+        const double amp = total_counts / (2.0 * M_PI * sigma * sigma);
+        for (int radius : {150, 350, 550})
+            for (int k = 0; k < 20; k++) {
+                const double phi = 2.0 * M_PI * k / 20.0 + 0.1 * radius;
+                const int cx = static_cast<int>(std::lround(600 + radius * std::cos(phi)));
+                const int cy = static_cast<int>(std::lround(600 + radius * std::sin(phi)));
+                for (int dy = -14; dy <= 14; dy++)
+                    for (int dx = -14; dx <= 14; dx++)
+                        image[static_cast<size_t>(cy + dy) * w + (cx + dx)] +=
+                                static_cast<int32_t>(std::lround(
+                                        amp * std::exp(-(dx * dx + dy * dy) / (2.0 * sigma * sigma))));
+                spots.emplace_back(static_cast<uint32_t>(cx), static_cast<uint32_t>(cy),
+                                   static_cast<int64_t>(total_counts));
+            }
+    }
+}
+
+// The width the adaptive integration radius is set from. A round Gaussian of width sigma holds 80 %
+// of its flux inside sqrt(2 ln 5) * sigma = 1.794 * sigma, and that is what the estimator has to
+// return - over an aperture that owes nothing to the integrator's r1, which is the whole point of
+// measuring it here rather than reading the integrator's own second moment.
+TEST_CASE("SpotWidth_Gaussian", "[process]") {
+    constexpr int W = 1200, H = 1200;
+    DiffractionGeometry geometry;
+    geometry.BeamX_pxl(600).BeamY_pxl(600).DetectorDistance_mm(200).PixelSize_mm(0.075)
+            .Wavelength_A(1.0);
+
+    for (double sigma : {1.0, 2.2}) {
+        ImagePreprocessorBuffer image(static_cast<size_t>(W) * H);
+        std::vector<DiffractionSpot> spots;
+        PaintGaussianSpots(image, W, sigma, 20000.0, spots);
+
+        std::vector<spot_width::FluxCurve> curves;
+        MeasureSpotFluxCurves(image, W, H, geometry, spots, curves);
+        REQUIRE(curves.size() >= 45);
+
+        const auto r80 = spot_width::R80AtReference(curves);
+        REQUIRE(r80.has_value());
+        CHECK(*r80 == Catch::Approx(1.794 * sigma).margin(0.3));
+    }
+
+    // The rule the measurement drives: the shipped radius below the line, the capped one above it.
+    CHECK(spot_width::R1ForWidth(1.0f) == 4.0f);
+    CHECK(spot_width::R1ForWidth(1.794f) == 4.0f);
+    CHECK(spot_width::R1ForWidth(2.4f) == 5.0f);
+    CHECK(spot_width::R1ForWidth(3.947f) == 6.0f);
+    CHECK(spot_width::R1ForWidth(9.0f) == 6.0f);
+}