Integration: the spot footprint includes how far the spots sit from their predicted reflections

A crystal of slightly misaligned domains (a ferroelastic domain twin below a phase transition, a
split crystal) records each reflection as two or more compact spots around the averaged lattice's
prediction, moving apart with resolution. The pre-scan footprint is measured about each spot, so it
saw compact spots; the r1 disk held the gap between them and the background ring sat on them.

The geometry pre-pass now compares every indexed spot with the predicted position of its own
reflection on the same frame and adds the mean square offset (radial and tangential, by distance
from the beam) to the pre-scan widths; the canonical pass integrates with that table. Where spots sit
on their predictions this moves the widths by the prediction error alone (lysozyme: 0.2-1.0 px, no
reflection outgrows r1); on a 100 K KDP domain twin the offsets reach 8-16 px at the edge.

KDP (kdp_x10sa_20keV), battery SHELXL recipe on the COD model: R1 0.272 -> 0.048, wR2 0.685 ->
0.129, EXTI 26.9 -> 0.025, GooF 3.2 -> 1.29 (XDS: 0.112 / 0.333 / 0.054 / 1.25); R_meas 21.6% -> 5.3%.

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

image_analysis/bragg_integration/SpotFootprint.cpp

@@ -5,6 +5,8 @@
 
 #include <algorithm>
 #include <cmath>
+#include <map>
+#include <tuple>
 
 namespace {
 
@@ -124,3 +126,50 @@ SpotFootprint FootprintFromSpots(const std::vector<FootprintSpot> &spots, float
     }
     return fp;
 }
+
+void MeasureFootprintOffsets(const std::vector<SpotToSave> &spots, const std::vector<Reflection> &reflections,
+                             float beam_x, float beam_y, std::vector<FootprintOffset> &out) {
+    std::map<std::tuple<int, int, int>, const Reflection *> predicted;
+    for (const auto &r : reflections)
+        predicted[{r.h, r.k, r.l}] = &r;
+    for (const auto &s : spots) {
+        if (!s.indexed || s.lattice != 0) continue;
+        const auto it = predicted.find({static_cast<int>(s.h), static_cast<int>(s.k), static_cast<int>(s.l)});
+        if (it == predicted.end()) continue;
+        const float px = it->second->predicted_x, py = it->second->predicted_y;
+        const float rx = px - beam_x, ry = py - beam_y;
+        const float r = std::sqrt(rx * rx + ry * ry);
+        if (!(r > 1.0f)) continue;
+        const float ux = rx / r, uy = ry / r;
+        const float dx = s.x - px, dy = s.y - py;
+        out.push_back({r, dx * ux + dy * uy, -dx * uy + dy * ux});
+    }
+}
+
+SpotFootprint FootprintWithOffsets(const SpotFootprint &widths, std::vector<FootprintOffset> offsets) {
+    if (widths.empty()) return widths;
+    // Sorted, so the sums below do not depend on the order the frames were measured in.
+    std::sort(offsets.begin(), offsets.end());
+    const int n = static_cast<int>(widths.sigma_rad.size());
+    std::vector<double> s2r(n, 0.0), s2t(n, 0.0);
+    std::vector<int> cnt(n, 0);
+    for (const auto &o : offsets) {
+        const int b = std::clamp(static_cast<int>(o.r_px / widths.bin_px), 0, n - 1);
+        s2r[b] += static_cast<double>(o.off_rad) * o.off_rad;
+        s2t[b] += static_cast<double>(o.off_tan) * o.off_tan;
+        ++cnt[b];
+    }
+    std::vector<int> filled;
+    for (int b = 0; b < n; ++b)
+        if (cnt[b] >= FOOTPRINT_MIN_SPOTS_PER_BIN) filled.push_back(b);
+    if (filled.empty()) return widths;
+    SpotFootprint fp = widths;
+    for (int b = 0; b < n; ++b) {
+        int best = filled.front();
+        for (int f : filled)
+            if (std::abs(f - b) < std::abs(best - b)) best = f;
+        fp.sigma_rad[b] = static_cast<float>(std::sqrt(widths.sigma_rad[b] * widths.sigma_rad[b] + s2r[best] / cnt[best]));
+        fp.sigma_tan[b] = static_cast<float>(std::sqrt(widths.sigma_tan[b] * widths.sigma_tan[b] + s2t[best] / cnt[best]));
+    }
+    return fp;
+}

image_analysis/bragg_integration/SpotFootprint.h

@@ -23,12 +23,25 @@
 // to the integrator through BraggIntegrationSettings; where it says a spot outgrows the r1 disk, the
 // background ring is moved clear of the spot and the profile is fitted at the measured width
 // (BraggStencil.h). Compact spots leave the integration exactly as it was.
+//
+// A spot can also be in the wrong PLACE for the disk: a crystal of slightly misaligned domains - the
+// ferroelastic domains of a crystal below a phase transition, a cracked or split crystal - records
+// each reflection as two or more spots around the position the averaged lattice predicts, and they
+// move apart with resolution. Each of them is compact, so the widths above, measured about each spot,
+// say nothing; the r1 disk holds the gap between them and the background ring lands on them. What
+// does see it is where the spots sit against the prediction: once the lattice is known, every
+// indexed spot is compared with the predicted position of its own reflection on the same frame, and
+// the mean square of that offset is added to the width of the spot, radially and tangentially apart.
+// The table then describes the reflection rather than the spot - and on a crystal whose spots sit on
+// their predictions it moves the widths by the prediction error alone, a fraction of a pixel.
 // =============================================================================
 
 #include <cstdint>
 #include <vector>
 
 #include "../../common/BraggIntegrationSettings.h"
+#include "../../common/Reflection.h"
+#include "../../common/SpotToSave.h"
 
 // One spot's measured widths, and its distance from the beam centre [px].
 struct FootprintSpot {
@@ -48,6 +61,27 @@ void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam
 // holding too few spots takes the nearest bin that has enough; with no such bin the table is empty.
 SpotFootprint FootprintFromSpots(const std::vector<FootprintSpot> &spots, float r_max);
 
+// One indexed spot's offset from the predicted position of its reflection [px], along and across the
+// radius, and the prediction's distance from the beam centre [px].
+struct FootprintOffset {
+    float r_px = 0.0f;
+    float off_rad = 0.0f;
+    float off_tan = 0.0f;
+    bool operator<(const FootprintOffset &o) const {
+        return r_px != o.r_px ? r_px < o.r_px : off_rad != o.off_rad ? off_rad < o.off_rad : off_tan < o.off_tan;
+    }
+};
+
+// The offsets of one frame: every spot indexed on the main lattice against the reflection of the same
+// hkl predicted on that frame. Spots whose reflection is not among the predictions are left out.
+void MeasureFootprintOffsets(const std::vector<SpotToSave> &spots, const std::vector<Reflection> &reflections,
+                             float beam_x, float beam_y, std::vector<FootprintOffset> &out);
+
+// The width table with the offsets added in quadrature: in each distance bin of `widths`, the mean
+// square offset along and across the radius. A bin holding too few offsets takes the nearest bin
+// that has enough; with no such bin, or no widths, the widths come back unchanged.
+SpotFootprint FootprintWithOffsets(const SpotFootprint &widths, std::vector<FootprintOffset> offsets);
+
 // Spots per distance bin a table needs before it is believed.
 constexpr int FOOTPRINT_MIN_SPOTS_PER_BIN = 20;
 constexpr int FOOTPRINT_BINS = 12;

rugnux/Rugnux.cpp

@@ -1311,6 +1311,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
             const float H = static_cast<float>(experiment_.GetYPixelsNum());
             const SpotFootprint fp = FootprintFromSpots(footprint_spots,
                                                         std::hypot(std::max(bx, W - bx), std::max(by, H - by)));
+            if (!fp.empty())
+                prescan_footprint_ = fp;
             // It acts only where a spot outgrows the r1 disk (BraggStencil.h), so a pattern of compact
             // spots is left on exactly the settings - and the passes - it had without it.
             const float r1_now = experiment_.GetBraggIntegrationSettings().GetR1();
@@ -3719,8 +3721,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
         indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A();
         indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
         indexer = std::make_unique<IndexAndRefine>(experiment_, indexer_pool.get());
-        // With no per-image file to write, nothing reads the message's reflection list.
-        if (!writer_queue)
+        // With no per-image file to write, nothing reads the message's reflection list - except the
+        // geometry pre-pass, which measures the spots against it (SpotFootprint.h).
+        if (!writer_queue && !(geometry_prepass && prescan_footprint_))
             indexer->KeepReflectionsInMessage(false);
         // The reference that breaks the per-image indexing ambiguity: a reference MTZ where there is
         // one, otherwise intensities computed from the model, which carry the same information. Serial
@@ -6117,6 +6120,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     std::mutex harmonic_m;
     HarmonicEvidence harmonic_evidence;
 
+    // Where the indexed spots sit against their predicted reflections, measured on the geometry
+    // pre-pass for the canonical pass's footprint (SpotFootprint.h).
+    const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value();
+    std::mutex footprint_m;
+    std::vector<FootprintOffset> footprint_offsets;
+
     auto azint_worker = [&]() {
         std::vector<uint8_t> decompression_buffer;
         std::shared_ptr<JFJochReaderRawImage> img;
@@ -6192,6 +6201,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
                                        /*enable_fused_adaptive_gpu=*/true);
         AzimuthalIntegrationProfile profile(mapping);
         HarmonicEvidence worker_harmonic;
+        std::vector<FootprintOffset> worker_offsets;
 
         while (!cancelled_) {
             const int ordinal = next_ordinal.fetch_add(1);
@@ -6240,6 +6250,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             }
 
             AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic);
+            if (measure_offsets)
+                MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(),
+                                        experiment_.GetBeamY_pxl(), worker_offsets);
 
             plots.Add(msg, profile);
             if (writer_queue) writer_queue->Post(msg, img);
@@ -6253,6 +6266,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             std::lock_guard lock(harmonic_m);
             harmonic_evidence.Add(worker_harmonic);
         }
+        {
+            std::lock_guard lock(footprint_m);
+            footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end());
+        }
 
         const auto bragg = analysis.BraggCounts();
         bragg_predicted += bragg.predicted;
@@ -6306,6 +6323,31 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     result.cancelled = cancelled_;
     result.images_processed = finished_count.load();
 
+    // The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the
+    // spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the
+    // reflection outgrows the r1 disk, and like them it is held back for the canonical pass.
+    if (measure_offsets && !cancelled_) {
+        const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets);
+        const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings());
+        bool outgrows = false;
+        for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
+            outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1();
+        std::string table;
+        for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
+            table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]);
+        logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their "
+                    "predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}",
+                    footprint_offsets.size(), table,
+                    outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1()));
+        if (outgrows) {
+            if (!bragg_before_adaptive_)
+                bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings();
+            BraggIntegrationSettings bis = canonical;
+            bis.Footprint(fp);
+            bragg_adaptive_ = bis;
+        }
+    }
+
     if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) {
         const SupercellProbe probe = indexer->GetSupercellProbe();
         if (probe[0][0].n > 0 && probe[0][1].n > 0) {

rugnux/Rugnux.h

@@ -719,6 +719,11 @@ class Rugnux {
     // radius before the pre-pass and applying it after keeps the geometry - and so the lattice - exactly
     // what the fixed radius produced, and still integrates the canonical pass at the measured width.
     std::optional<BraggIntegrationSettings> bragg_adaptive_;
+    // The spot footprint as the pre-scan measured it, about each spot (SpotFootprint.h), whether or
+    // not it was adopted. The geometry pre-pass adds to it how far the spots sit from their predicted
+    // reflections - which takes a lattice the pre-scan does not have - and the sum goes to the
+    // canonical pass in bragg_adaptive_. Kept apart so a repeated pre-pass starts from it again.
+    std::optional<SpotFootprint> prescan_footprint_;
     // Fraction of the predicted reflections the last completed pass dropped because NEIGHBOURING
     // reflections left their background ring below six clean pixels. Measured, not predicted - it is
     // the only thing that tells a radius the pattern can take from one it cannot (see RunAllPasses).

tests/SpotFootprintTest.cpp

@@ -66,3 +66,39 @@ TEST_CASE("SpotFootprint_TableFillsSparseBinsFromNeighbours", "[Integration][por
     REQUIRE(fp.sigma_rad[FOOTPRINT_BINS - 3] == 3.0f);   // nearer the last
     REQUIRE(FootprintFromSpots({}, 1200.0f).empty());
 }
+
+// A reflection recorded as a doublet: two indexed spots on either side of the prediction, along the
+// radius. The offsets are found against the reflection of the same hkl, and their mean square is
+// added to the widths of the bin they fall in.
+TEST_CASE("SpotFootprint_OffsetsFromPredictionWidenTheTable", "[Integration][portable]") {
+    const float bx = 600.0f, by = 600.0f;
+    std::vector<Reflection> refl(1);
+    refl[0].h = 1; refl[0].k = 2; refl[0].l = 3;
+    refl[0].predicted_x = bx + 1000.0f;   // on +x, so radial = x
+    refl[0].predicted_y = by;
+    std::vector<SpotToSave> spots;
+    spots.push_back({.x = bx + 1004.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true});
+    spots.push_back({.x = bx + 996.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true});
+    spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = 0, .h = 3, .k = 2, .l = 1, .indexed = true});   // no such prediction
+    spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = -1, .h = 1, .k = 2, .l = 3, .indexed = false});  // not indexed
+    std::vector<FootprintOffset> offsets;
+    MeasureFootprintOffsets(spots, refl, bx, by, offsets);
+    REQUIRE(offsets.size() == 2);
+    REQUIRE_THAT(offsets[0].off_rad, Catch::Matchers::WithinAbs(4.0, 1e-4));
+    REQUIRE_THAT(offsets[0].off_tan, Catch::Matchers::WithinAbs(0.0, 1e-4));
+
+    std::vector<FootprintOffset> pool;
+    for (int i = 0; i < FOOTPRINT_MIN_SPOTS_PER_BIN; ++i)
+        pool.insert(pool.end(), offsets.begin(), offsets.end());
+    SpotFootprint widths;
+    widths.bin_px = 100.0f;
+    widths.sigma_rad.assign(12, 1.0f);
+    widths.sigma_tan.assign(12, 2.0f);
+    const SpotFootprint fp = FootprintWithOffsets(widths, pool);
+    for (int b = 0; b < 12; ++b) {   // one filled bin (10) serves them all
+        REQUIRE_THAT(fp.sigma_rad[b], Catch::Matchers::WithinAbs(std::sqrt(17.0), 1e-4));
+        REQUIRE_THAT(fp.sigma_tan[b], Catch::Matchers::WithinAbs(2.0, 1e-4));
+    }
+    REQUIRE(FootprintWithOffsets(widths, {}).sigma_rad[3] == 1.0f);
+    REQUIRE(FootprintWithOffsets(SpotFootprint{}, pool).empty());
+}