Merge branch 'cpu-loops-tail' into rc173 (CPU: vectorised spot-finder loops, prediction prefilter, compact scaling partials, tiled integration masks)

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

common/PixelMask.cpp

@@ -59,6 +59,9 @@ void PixelMask::UpdateBinaryMask() {
     for (size_t i = 0; i < mask.size(); i++)
         binary_mask[i] = (mask[i] != 0);
     binary_mask_checksum = TableChecksum(binary_mask.data(), binary_mask.size());
+    packed_mask.assign((mask.size() + 31) / 32, 0);
+    for (size_t i = 0; i < mask.size(); i++)
+        packed_mask[i / 32] |= static_cast<uint32_t>(binary_mask[i]) << (i % 32);
 }
 
 void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) {
@@ -127,6 +130,10 @@ const std::vector<uint8_t> &PixelMask::GetBinaryMask() const {
     return binary_mask;
 }
 
+const std::vector<uint32_t> &PixelMask::GetPackedMask() const {
+    return packed_mask;
+}
+
 uint64_t PixelMask::GetBinaryMaskChecksum() const {
     return binary_mask_checksum;
 }

common/PixelMask.h

@@ -26,6 +26,9 @@ class PixelMask {
     // derived here once instead of in every one of those engines.
     std::vector<uint8_t> binary_mask;
     uint64_t binary_mask_checksum = 0;
+    // The same flag packed 32 pixels to a word (pixel i is bit i % 32 of word i / 32) - the form the CPU
+    // image preprocessor reads, derived here for the same reason.
+    std::vector<uint32_t> packed_mask;
 
     uint32_t LoadMask(const std::vector<uint32_t>& mask, uint8_t bit);
     // Everything that follows from the mask, recomputed wherever the mask changes.
@@ -68,6 +71,7 @@ public:
     [[nodiscard]] const std::vector<uint32_t> &GetMask(const DiffractionExperiment& experiment) const;
     [[nodiscard]] const std::vector<uint32_t> &GetMask() const;
     [[nodiscard]] const std::vector<uint8_t> &GetBinaryMask() const;
+    [[nodiscard]] const std::vector<uint32_t> &GetPackedMask() const;
     [[nodiscard]] uint64_t GetBinaryMaskChecksum() const;
 
     [[nodiscard]] std::vector<uint32_t> GetUserMask(const DiffractionExperiment& experiment) const;

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -44,7 +44,9 @@ struct ImageSampler {
 }   // namespace
 
 BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment)
-    : BraggIntegrationEngine(experiment) {}
+    : BraggIntegrationEngine(experiment),
+      refl_mask(static_cast<int>(xpixel), static_cast<int>(ypixel), 0),
+      owner(static_cast<int>(xpixel), static_cast<int>(ypixel), BRAGG_OWNER_NONE) {}
 
 template <class Sampler>
 std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
@@ -71,17 +73,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     //     signal is excluded from this reflection's r2..r3 background ring. The region is the INNER
     //     stencil ellipse, taken in each neighbour's OWN frame - an elongated neighbour whose streak
     //     is still marked as a disk would leak its tails into this reflection's ring. ---
-    // Back to empty what the previous call wrote (see the header).
-    for (const Rect &b : refl_mask_dirty)
-        for (int y = b.y0; y <= b.y1; ++y)
-            std::fill(refl_mask.begin() + y * W + b.x0, refl_mask.begin() + y * W + b.x1 + 1, 0);
-    refl_mask_dirty.clear();
-    for (const Rect &b : owner_dirty)
-        for (int y = b.y0; y <= b.y1; ++y)
-            std::fill(owner.begin() + y * W + b.x0, owner.begin() + y * W + b.x1 + 1, BRAGG_OWNER_NONE);
-    owner_dirty.clear();
-
-    refl_mask.resize(npixel, 0);
+    refl_mask.Clear();
+    owner.Clear();
     for (size_t i = 0; i < npredicted; ++i) {
         const auto &r = predicted[i];
         const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
@@ -89,12 +82,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + st.ex_in + 1.0f)));
         const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_in - 1.0f)));
         const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_in + 1.0f)));
-        if (x0 <= x1 && y0 <= y1)
-            refl_mask_dirty.push_back({x0, x1, y0, y1});
         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);
-                if (d.inner < r2_sq) refl_mask[y * W + x] = 1;
+                if (d.inner < r2_sq) refl_mask.At(x, y) = 1;
             }
     }
 
@@ -106,7 +97,6 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     const bool exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum;
     const bool use_owner = overlap != OverlapMode::Off;
     if (use_owner) {
-        owner.resize(npixel, BRAGG_OWNER_NONE);
         const float claim_sq = claim * claim;
         for (size_t i = 0; i < npredicted; ++i) {
             const auto &r = predicted[i];
@@ -114,21 +104,19 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + claim)));
             const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - claim)));
             const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + claim)));
-            if (x0 <= x1 && y0 <= y1)
-                owner_dirty.push_back({x0, x1, y0, y1});
             for (int y = y0; y <= y1; ++y)
                 for (int x = x0; x <= x1; ++x) {
                     const float dx = x - r.predicted_x, dy = y - r.predicted_y;
                     const float d2 = dx * dx + dy * dy;
                     if (d2 >= claim_sq) continue;
-                    uint32_t &o = owner[y * W + x];
+                    uint32_t &o = owner.At(x, y);
                     o = std::min(o, BraggOwnerKey(std::sqrt(d2), inv_claim, static_cast<int>(i)));
                 }
         }
     }
     auto clean = [&](int x, int y, size_t i) {
         return !use_owner
-               || BraggOwnedBy(owner[static_cast<size_t>(y) * W + x], static_cast<int>(i));
+               || BraggOwnedBy(owner.Get(x, y), static_cast<int>(i));
     };
 
     // --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). ---
@@ -213,7 +201,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     y_sum += y;
                     ++n_inner_valid;
                 } else if (d.inner >= r2_sq && d.outer < r3_sq) {
-                    if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; }
+                    if (refl_mask.Get(x, y)) { ++n_bkg_neighbour; continue; }
                     if (!valid(px)) continue;
                     bkg_sum += static_cast<double>(px);
                     if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
@@ -267,7 +255,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     for (int x = x0; x <= x1; ++x) {
                         const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
                         if (!(d.inner >= r2_sq && d.outer < r3_sq)) continue;
-                        if (refl_mask[y * W + x]) continue;
+                        if (refl_mask.Get(x, y)) continue;
                         const int32_t px = img[y * W + x];
                         if (!valid(px)) continue;
                         if (static_cast<double>(px) <= thr) {

image_analysis/bragg_integration/BraggIntegrationEngineCPU.h

@@ -14,15 +14,44 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine {
     // Core integrator, templated on a pixel sampler so it reads either the preprocessed int32 buffer
     // or a raw CompressedImage of any pixel type - both presented per-pixel in the INT32_MIN(masked)/
     // INT32_MAX(saturated) convention - without ever materialising a second full-image copy.
-    // Full-frame scratch of RunImpl: the reflection mask and the signal-region owner map. Kept from
-    // call to call rather than allocated per image - at 18 Mpixel that was 90 MB, a fresh mmap and a
-    // page fault per 4 kB every call - and put back to their empty state one rectangle at a time: the
-    // rectangles a call writes are recorded, and the next call clears exactly those before it starts.
-    struct Rect { int x0, x1, y0, y1; };
-    std::vector<uint8_t> refl_mask;
-    std::vector<uint32_t> owner;
-    std::vector<Rect> refl_mask_dirty;
-    std::vector<Rect> owner_dirty;
+    // Full-frame scratch of RunImpl: the reflection mask and the signal-region owner map. A call writes
+    // only around its reflections, so these keep only the 16x16-pixel tiles written since the last
+    // Clear(): a frame-sized array was mostly never read, yet over a sweep every page of it got
+    // touched, in every worker. Reading a tile nothing wrote gives `empty`.
+    template <class T>
+    class TiledFrame {
+        static constexpr int TILE = 16;
+        int tiles_x;
+        T empty;
+        std::vector<int32_t> tile_start;   // per tile: where it starts in `pixels`, -1 = not written
+        std::vector<int32_t> written;      // the tiles written, for Clear()
+        std::vector<T> pixels;
+    public:
+        TiledFrame(int width, int height, T empty)
+            : tiles_x((width + TILE - 1) / TILE), empty(empty),
+              tile_start(static_cast<size_t>(tiles_x) * ((height + TILE - 1) / TILE), -1) {}
+        T Get(int x, int y) const {
+            const int32_t start = tile_start[(y / TILE) * tiles_x + x / TILE];
+            return start < 0 ? empty : pixels[start + (y % TILE) * TILE + x % TILE];
+        }
+        T &At(int x, int y) {
+            const int t = (y / TILE) * tiles_x + x / TILE;
+            if (tile_start[t] < 0) {
+                tile_start[t] = static_cast<int32_t>(pixels.size());
+                pixels.resize(pixels.size() + TILE * TILE, empty);
+                written.push_back(t);
+            }
+            return pixels[tile_start[t] + (y % TILE) * TILE + x % TILE];
+        }
+        void Clear() {
+            for (int t : written)
+                tile_start[t] = -1;
+            written.clear();
+            pixels.clear();
+        }
+    };
+    TiledFrame<uint8_t> refl_mask;
+    TiledFrame<uint32_t> owner;
 
     template <class Sampler>
     std::vector<Reflection> RunImpl(const Sampler &img, const std::vector<Reflection> &predicted,

image_analysis/bragg_prediction/BraggPredictionRot.cpp

@@ -77,6 +77,28 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
         experiment.GetBraggIntegrationSettings().GetFlightPath(), geom.GetDetectorDistance_mm(),
         geom.GetWavelength_A());
 
+    // The phi test after the atan2 below, taken before it: a solution with |phi| above phi_limit fails
+    // that test whatever its own bandwidth term, since mos_eff_rad is largest at the resolution limit.
+    // |phi| > phi_limit is cos(phi) < cos(phi_limit), which the (cos, sin) pair answers without the
+    // angle. phi_limit carries a 1e-3 rad margin, far above float rounding, so every solution rejected
+    // here is rejected by that test too. -1 = nothing is rejected here.
+    double cos_phi_limit = -1.0;
+    if (settings.min_zeta > 0.0f) {
+        double mos_max_rad = mos_angle_rad;
+        if (bandwidth_sigma > 0.0f) {
+            const double sin_theta = half_wavelength_A * static_cast<double>(one_over_dmax);
+            const double dphi_bw = sin_theta < 1.0
+                                   ? bandwidth_sigma * sin_theta / std::sqrt(1.0 - sin_theta * sin_theta)
+                                   : INFINITY;
+            mos_max_rad = std::sqrt(mos_max_rad * mos_max_rad + dphi_bw * dphi_bw);
+        }
+        const double phi_limit = half_wedge_angle_rad
+                                 + (settings.mosaicity_multiplier * mos_max_rad + 1e-5) / settings.min_zeta
+                                 + 1e-3;
+        if (phi_limit < PI)
+            cos_phi_limit = std::cos(phi_limit);
+    }
+
     for (int h = -settings.max_h; h <= settings.max_h; h++) {
         // Precompute A* h contribution
 
@@ -127,6 +149,9 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
 
                     const float cosphi = (p_m1 * p0_m1 + p_m3 * p0_m3) / rho_sq;
                     const float sinphi = (p_m1 * p0_m3 - p_m3 * p0_m1) / rho_sq;
+                    const double c = cosphi, s = sinphi;
+                    if (c < cos_phi_limit * std::sqrt(c * c + s * s))
+                        continue;
                     float phi = -1.0f * std::atan2(sinphi, cosphi);
 
                     // Most solutions are far from this frame. zeta >= min_zeta, so the rocking-curve test

image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp

@@ -5,10 +5,7 @@
 
 ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &experiment, const PixelMask &mask)
     : ImagePreprocessor(experiment),
-      mask_1bit(npixels, false) {
-    for (int i = 0; i < npixels; i++)
-        mask_1bit[i] = (mask.GetMask().at(i) != 0);
-}
+      mask_bits(mask.GetPackedMask()) {}
 
 ImageStatistics ImagePreprocessorCPU::Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *image_ptr, CompressedImageMode image_mode) {
     if (processed_image.size() != npixels)
@@ -49,7 +46,7 @@ void ImagePreprocessorCPU::AnalyzeBlock(ImagePreprocessorBuffer &processed_image
         sat_pixel_val = static_cast<T>(saturation_limit);
 
     for (size_t i = 0; i < n; i++) {
-        if (mask_1bit[first + i] != 0) {
+        if ((mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U) {
             out[i] = INT32_MIN;
             ++ret.masked_pixel_count;
         } else if (image[i] == err_pixel_val) {

image_analysis/image_preprocessing/ImagePreprocessorCPU.h

@@ -7,7 +7,9 @@
 
 class ImagePreprocessorCPU : public ImagePreprocessor {
 
-    std::vector<bool> mask_1bit;
+    // The pixel mask one bit per pixel, 32 to a word - PixelMask derives it once, so each engine only
+    // copies it.
+    std::vector<uint32_t> mask_bits;
 
     template <class T> void AnalyzeBlock(ImagePreprocessorBuffer &processed_image, size_t first, const uint8_t *input,
                                          size_t n, T err_value, T sat_value, ImageStatistics &stats);

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -1336,12 +1336,13 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
     return n_groups;
 }
 
-double RotationScaleMerge::ReferenceWeight(const Obs &o, double min_partiality, float &I_corr) {
-    if (o.group < 0) return 0.0;
-    if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) return 0.0;
-    if (o.partiality < min_partiality) return 0.0;
-    const float Ic = o.I * o.corr;
-    const float sigma_corr = o.sigma * o.corr;
+double RotationScaleMerge::ReferenceWeight(int32_t group, float corr, float partiality, float I, float sigma,
+                                           double min_partiality, float &I_corr) {
+    if (group < 0) return 0.0;
+    if (!(corr > 0.0f) || !std::isfinite(corr)) return 0.0;
+    if (partiality < min_partiality) return 0.0;
+    const float Ic = I * corr;
+    const float sigma_corr = sigma * corr;
     if (!std::isfinite(Ic) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) return 0.0;
     I_corr = Ic;
     return 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
@@ -1383,7 +1384,30 @@ void RotationScaleMerge::ReducePartialGroupMeans(int n_groups, std::vector<doubl
     });
 }
 
-void RotationScaleMerge::FitPerFrameG(std::vector<Obs> &obs, const std::vector<int32_t> &fstart,
+void RotationScaleMerge::ReduceScalingGroupMeans(int n_groups, const std::vector<ScalingRef> &refs,
+                                                 const std::vector<float> &corr,
+                                                 std::vector<double> &out_mean) const {
+    out_mean.assign(n_groups, NAN);
+    ParallelChunks(n_groups, ThreadsForWork(corr.size(), nthreads), [&](int glo, int ghi) {
+        for (int g = glo; g < ghi; ++g) {
+            double sw = 0.0, swI = 0.0;
+            const int lo = partial_group_start[g], hi = lo + partial_group_count[g];
+            for (int p = lo; p < hi; ++p) {
+                const ScalingRef &r = refs[p];
+                float I_corr = 0.0f;
+                const double w = ReferenceWeight(g, corr[partial_group_perm[p]], r.partiality, r.I, r.sigma,
+                                                 min_partiality, I_corr);
+                if (w <= 0.0) continue;
+                sw += w;
+                swI += w * I_corr;
+            }
+            if (sw > 0.0) out_mean[g] = swI / sw;
+        }
+    });
+}
+
+template <class O>
+void RotationScaleMerge::FitPerFrameG(const std::vector<O> &obs, const std::vector<int32_t> &fstart,
                                       const std::vector<int32_t> &fcount,
                                       const std::vector<double> &group_mean_in,
                                       bool unity, std::vector<double> &g) {
@@ -5560,6 +5584,30 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
             logger.Info("Per-frame scaling (partials): the same grouping and settings as an earlier merge, "
                         "so its {} iterations stand", partial_loop.iterations);
         } else {
+            // The loop works on its own compact copy of the partials (ScalingObs in the header) and hands
+            // corr back to them when it is done.
+            const int n_partials = static_cast<int>(partials.size());
+            const size_t partial_threads = ThreadsForWork(partials.size(), nthreads);
+            std::vector<ScalingObs> scaling_obs(n_partials);
+            std::vector<float> scaling_corr(n_partials);
+            ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
+                for (int i = lo; i < hi; ++i) {
+                    const Obs &o = partials[i];
+                    scaling_obs[i] = {o.frame, o.group, o.I, o.sigma, o.partiality, o.prescaling_corr, o.zeta,
+                                      o.on_ice};
+                    scaling_corr[i] = o.corr;
+                }
+            });
+            std::vector<ScalingRef> scaling_refs(partial_group_start.empty()
+                                                 ? 0 : partial_group_start.back() + partial_group_count.back());
+            ParallelChunks(n_groups, partial_threads, [&](int glo, int ghi) {
+                for (int g = glo; g < ghi; ++g)
+                    for (int p = partial_group_start[g]; p < partial_group_start[g] + partial_group_count[g]; ++p) {
+                        const ScalingObs &o = scaling_obs[partial_group_perm[p]];
+                        scaling_refs[p] = {o.I, o.sigma, o.partiality};
+                    }
+            });
+
             // RunScalingLoop rescales after every iteration before anything reads corr, so the iteration's
             // UpdateCorr is folded into that pass over the partials: the G and the fitted frames as the
             // iteration left them, then the rescale on top - the same two roundings, one pass instead of two.
@@ -5567,34 +5615,32 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
             std::vector<uint8_t> frame_fitted;
             partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window,
                                           /*release_restraint=*/false, [&] {
-                ReducePartialGroupMeans(n_groups, partial_mean);
-                FitPerFrameG(partials, frame_start, frame_count, partial_mean, /*unity=*/false, g_partial);
+                ReduceScalingGroupMeans(n_groups, scaling_refs, scaling_corr, partial_mean);
+                FitPerFrameG(scaling_obs, frame_start, frame_count, partial_mean, /*unity=*/false, g_partial);
                 g_fitted = g_partial;
                 frame_fitted = frame_scaled_scratch;
             }, [&](const std::vector<uint8_t> &apply, const std::vector<double> &ratio) {
-                ParallelChunks(static_cast<int>(partials.size()), ThreadsForWork(partials.size(), nthreads),
-                               [&](int lo, int hi) {
+                ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
                     for (int i = lo; i < hi; ++i) {
-                        Obs &o = partials[i];
+                        const ScalingObs &o = scaling_obs[i];
+                        float &corr = scaling_corr[i];
                         if (frame_fitted[o.frame]) {   // UpdateCorr
                             const double denom = static_cast<double>(o.partiality) * g_fitted[o.frame];
                             if (std::isfinite(o.prescaling_corr) && std::isfinite(denom) && denom > 0.0)
-                                o.corr = static_cast<float>(o.prescaling_corr / denom);
+                                corr = static_cast<float>(o.prescaling_corr / denom);
                             else
-                                o.corr = NAN;
+                                corr = NAN;
                         }
-                        if (apply[o.frame] && std::isfinite(o.corr))
-                            o.corr = static_cast<float>(o.corr * ratio[o.frame]);
+                        if (apply[o.frame] && std::isfinite(corr))
+                            corr = static_cast<float>(corr * ratio[o.frame]);
                     }
                 });
             });
 
-            PartialLoopMemo m{key, partial_loop, frame_scaled_scratch, g_partial, {}};
-            m.corr.resize(partials.size());
-            ParallelChunks(static_cast<int>(partials.size()), ThreadsForWork(partials.size(), nthreads),
-                           [&](int lo, int hi) {
-                for (int i = lo; i < hi; ++i) m.corr[i] = partials[i].corr;
+            ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
+                for (int i = lo; i < hi; ++i) partials[i].corr = scaling_corr[i];
             });
+            PartialLoopMemo m{key, partial_loop, frame_scaled_scratch, g_partial, std::move(scaling_corr)};
             if (partial_loop_memos.size() == MAX_PARTIAL_LOOP_MEMOS)
                 partial_loop_memos.erase(partial_loop_memos.begin());
             partial_loop_memos.push_back(std::move(m));

image_analysis/scale_merge/RotationScaleMerge.h

@@ -446,17 +446,40 @@ private:
     std::vector<int32_t> partial_group_start, partial_group_count;
     std::unique_ptr<int32_t[]> partial_group_perm;
 
+    // The CPU partial scaling loop's working copy of the partials. It passes over them several times
+    // per iteration and reads a handful of fields each time, so it works on these rather than on the
+    // eighty-byte Obs: ScalingObs holds, in the partials' own order, what FitPerFrameG and the corr
+    // update read; ScalingRef holds, in the CSR's group order, what the group means read besides corr,
+    // so that only corr is gathered. The same values, so the same arithmetic in the same order.
+    struct ScalingObs {
+        int32_t frame, group;
+        float I, sigma, partiality, prescaling_corr, zeta;
+        uint8_t on_ice;
+    };
+    struct ScalingRef {
+        float I, sigma, partiality;
+    };
+
     // The observation's inverse-variance weight in the scaling reference, and its scaled intensity;
     // 0 when the observation is not in the reference at all.
-    static double ReferenceWeight(const Obs &o, double min_partiality, float &I_corr);
+    static double ReferenceWeight(int32_t group, float corr, float partiality, float I, float sigma,
+                                  double min_partiality, float &I_corr);
+    static double ReferenceWeight(const Obs &o, double min_partiality, float &I_corr) {
+        return ReferenceWeight(o.group, o.corr, o.partiality, o.I, o.sigma, min_partiality, I_corr);
+    }
     // Inverse-variance per-group mean of I*corr over `obs` (the merge reference).
     void ReduceGroupMeans(const std::vector<Obs> &obs, int n_groups, std::vector<double> &out_mean) const;
     // The same over the partials, using the CSR above - the same sums in the same order.
     void ReducePartialGroupMeans(int n_groups, std::vector<double> &out_mean) const;
+    // The same again over the scaling loop's copy: refs in group order, corr in the partials' order.
+    void ReduceScalingGroupMeans(int n_groups, const std::vector<ScalingRef> &refs,
+                                 const std::vector<float> &corr, std::vector<double> &out_mean) const;
 
     // Robust per-frame G fit (IRLS, Cauchy k=3), unity=false uses the rotation partiality, unity=true the
     // scale-fulls (partiality already folded in). Reads out_mean[group] as the reference intensity.
-    void FitPerFrameG(std::vector<Obs> &obs, const std::vector<int32_t> &fstart,
+    // O is Obs or ScalingObs.
+    template <class O>
+    void FitPerFrameG(const std::vector<O> &obs, const std::vector<int32_t> &fstart,
                       const std::vector<int32_t> &fcount, const std::vector<double> &group_mean_in,
                       bool unity, std::vector<double> &g);
 

image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp

@@ -60,26 +60,89 @@ void AdaptiveSpotFinderCPU::AccumulateRingsBlock(const ImagePreprocessorBuffer &
     const size_t nbins = ring_sum.size();
     const float *corrections = mapping.Corrections().data();
 
+    // Consecutive pixels mostly share a ring, so a ring's sums are held in locals while they do and
+    // written back when the ring changes: the same additions in the same order, without a store and a
+    // reload of the same address on every pixel. The azimuthal-integration sums get a loop of their own
+    // over the block, so that neither loop runs out of registers for its sums.
+    if (fuse_azint) {
+        size_t cur = nbins;   // the ring held in the locals below; nbins = none
+        float az_sum = 0.0f, az_sum2 = 0.0f;
+        uint32_t az_cnt = 0;
+        for (size_t pxl = first; pxl < first + n; ++pxl) {
+            const int32_t v = image[pxl];
+            if (v == INT32_MIN || v == INT32_MAX) continue;   // bad / saturated
+            const uint16_t b = pixel_to_bin[pxl];
+            if (b >= nbins) continue;                          // masked / out of range (UINT16_MAX)
+            if (b != cur) {
+                if (cur != nbins) {
+                    azint_sum[cur] = az_sum;
+                    azint_sum2[cur] = az_sum2;
+                    azint_count[cur] = az_cnt;
+                }
+                cur = b;
+                az_sum = azint_sum[b];
+                az_sum2 = azint_sum2[b];
+                az_cnt = azint_count[b];
+            }
+            const float val = static_cast<float>(v) * corrections[pxl];
+            const float val_sq = val * val;
+            az_sum += val;
+            az_sum2 += val_sq;
+            ++az_cnt;
+        }
+        if (cur != nbins) {
+            azint_sum[cur] = az_sum;
+            azint_sum2[cur] = az_sum2;
+            azint_count[cur] = az_cnt;
+        }
+    }
+
+    // Values outside the histogram are listed by a second loop, run only when the block has any: a
+    // call in this loop would leave the sums in memory again.
+    uint32_t *hist = ring_hist.data();
+    bool overflow = false;
+    size_t cur = nbins;
+    int64_t sum = 0, cnt = 0;
+    uint64_t sum2 = 0;
     for (size_t pxl = first; pxl < first + n; ++pxl) {
         const int32_t v = image[pxl];
         if (v == INT32_MIN || v == INT32_MAX) continue;   // bad / saturated
         const uint16_t b = pixel_to_bin[pxl];
         if (b >= nbins) continue;                          // masked / out of range (UINT16_MAX)
-        if (fuse_azint) {
-            const float val = static_cast<float>(v) * corrections[pxl];
-            const float val_sq = val * val;
-            azint_sum[b] += val;
-            azint_sum2[b] += val_sq;
-            ++azint_count[b];
+        if (b != cur) {
+            if (cur != nbins) {
+                ring_sum[cur] = sum;
+                ring_sum2[cur] = sum2;
+                ring_cnt[cur] = cnt;
+            }
+            cur = b;
+            sum = ring_sum[b];
+            sum2 = ring_sum2[b];
+            cnt = ring_cnt[b];
         }
-        ring_sum[b] += v;
-        ring_sum2[b] += static_cast<uint64_t>(static_cast<int64_t>(v) * v);
-        ring_cnt[b] += 1;
+        sum += v;
+        sum2 += static_cast<uint64_t>(static_cast<int64_t>(v) * v);
+        cnt += 1;
         if (v >= 0 && v < HIST_VALUES)
-            ring_hist[b * HIST_VALUES + v] += 1;
+            hist[b * HIST_VALUES + v] += 1;
         else
-            ring_overflow.emplace_back(b, v);
+            overflow = true;
     }
+    if (cur != nbins) {
+        ring_sum[cur] = sum;
+        ring_sum2[cur] = sum2;
+        ring_cnt[cur] = cnt;
+    }
+
+    if (overflow)
+        for (size_t pxl = first; pxl < first + n; ++pxl) {
+            const int32_t v = image[pxl];
+            if (v == INT32_MIN || v == INT32_MAX) continue;
+            const uint16_t b = pixel_to_bin[pxl];
+            if (b >= nbins) continue;
+            if (v < 0 || v >= HIST_VALUES)
+                ring_overflow.emplace_back(b, v);
+        }
 }
 
 // A sigma-clip pass over the plain pass's values: each distinct value of a ring meets the same test the

image_analysis/spot_finding/ImageSpotFinderCPU.cpp

@@ -45,6 +45,24 @@ bool StrongInWindow(int64_t pxl_val, int64_t sum, int64_t sum2, int64_t valid,
                 // pixel is above SNR threshold
 }
 
+// Adds row `in` to the vertical sums and takes row `out` out of them; nullptr is no row. A bad or
+// saturated pixel adds 0. Written without branches, and as a function of its own so that the loop
+// vectorises whatever it is inlined into; integer sums, so the same totals as pixel by pixel.
+void SlideVertical(const int32_t *in, const int32_t *out, int32_t width,
+                   int64_t *sum, int64_t *sum2, uint16_t *valid) {
+    for (int32_t col = 0; col < width; col++) {
+        const int64_t a = in ? in[col] : INT32_MIN;
+        const int64_t r = out ? out[col] : INT32_MIN;
+        const bool a_ok = a != INT32_MAX && a != INT32_MIN;
+        const bool r_ok = r != INT32_MAX && r != INT32_MIN;
+        const int64_t a_v = a_ok ? a : 0;
+        const int64_t r_v = r_ok ? r : 0;
+        sum[col] += a_v - r_v;
+        sum2[col] += a_v * a_v - r_v * r_v;
+        valid[col] += static_cast<uint16_t>(a_ok - r_ok);
+    }
+}
+
 } // namespace
 
 void ImageSpotFinderCPU::DetectAt(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings,
@@ -163,45 +181,55 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image,
         std::vector<int64_t> sum2_vert(width, 0);
         std::vector<uint16_t> valid_vert(width, 0);
 
+        // Add the row entering the window (line_in) and take out the one leaving it (line_out); -1 is
+        // no row. A pixel strong in the previous pass is not counted either; SlideVertical adds it with
+        // the rest of its row and it is taken out again below - there are few, and the sums are
+        // integers, so the totals are the same as never adding it.
+        const int32_t *img = image.data();
+        const int32_t w = width;
+        int64_t *sv = sum_vert.data();
+        int64_t *sv2 = sum2_vert.data();
+        uint16_t *vv = valid_vert.data();
+        auto slide_vert = [&](int line_in, int line_out) {
+            SlideVertical(line_in >= 0 ? img + static_cast<size_t>(line_in) * w : nullptr,
+                          line_out >= 0 ? img + static_cast<size_t>(line_out) * w : nullptr, w, sv, sv2, vv);
+            if (!prev_strong)
+                return;
+            // sign = +1 takes out what line_in added, -1 puts back what line_out took out.
+            const auto uncount_strong = [&](int line, int64_t sign) {
+                const int32_t first = line * w, last = first + w - 1;
+                for (int32_t word = first / 32; word <= last / 32; word++) {
+                    uint32_t bits = prev_strong[word];
+                    if (word == first / 32)
+                        bits &= UINT32_MAX << (first % 32);
+                    if (word == last / 32 && last % 32 != 31)
+                        bits &= (1U << (last % 32 + 1)) - 1;
+                    for (; bits; bits &= bits - 1) {
+                        const int32_t pxl = word * 32 + std::countr_zero(bits);
+                        const int64_t v = img[pxl];
+                        if (v == INT32_MAX || v == INT32_MIN)
+                            continue;
+                        sv[pxl - first] -= sign * v;
+                        sv2[pxl - first] -= sign * v * v;
+                        vv[pxl - first] -= static_cast<uint16_t>(sign);
+                    }
+                }
+            };
+            if (line_in >= 0)
+                uncount_strong(line_in, 1);
+            if (line_out >= 0)
+                uncount_strong(line_out, -1);
+        };
+
         for (int line = 0; line < NBX; line++) {
             new_row(line);
-            for (int col = 0; col < width; col++) {
-                auto pxl = line * width + col;
-
-                int64_t tmp = value_at(pxl);
-                if (tmp != INT32_MAX && tmp != INT32_MIN) {
-                    sum_vert[col] += tmp;
-                    sum2_vert[col] += tmp * tmp;
-                    valid_vert[col] += 1;
-                }
-            }
+            slide_vert(line, -1);
         }
 
         for (int line = 0; line < height; line++) {
             if (line < height - NBX)
                 new_row(line + NBX);
-            for (int col = 0; col < width; col++) {
-                if (line < height - NBX) {
-                    auto pxl = (line + NBX) * width + col;
-
-                    int64_t tmp = value_at(pxl);
-                    if (tmp != INT32_MAX && tmp != INT32_MIN) {
-                        sum_vert[col] += tmp;
-                        sum2_vert[col] += tmp * tmp;
-                        valid_vert[col] += 1;
-                    }
-                }
-
-                if (line >= NBX + 1) {
-                    auto pxl = (line - (NBX + 1)) * width + col;
-                    int64_t tmp = value_at(pxl);
-                    if (tmp != INT32_MAX && tmp != INT32_MIN) {
-                        sum_vert[col] -= tmp;
-                        sum2_vert[col] -= tmp * tmp;
-                        valid_vert[col] -= 1;
-                    }
-                }
-            }
+            slide_vert(line < height - NBX ? line + NBX : -1, line >= NBX + 1 ? line - (NBX + 1) : -1);
 
             if (candidates) {
                 // Only the runs of 32-column blocks whose result is read (first_pass_needed; every