CPU adaptive spot finder: sigma clips from a histogram, first pass only where read

Two whole-image passes per frame out of the CPU spot finder (the pre-scan's
finder on every build, and every image on the CPU-only build):

- AccumulateRings ran three passes over the frame - the plain ring statistics
  and two sigma clips. The plain pass now also counts each ring's valid
  values in a histogram (0..1023, the rest in a short list), and the clip
  passes sum over the distinct values: each meets the same float test its
  pixels would, and the sums are integers, so the totals are the same.
- The local test's first pass is read by DetectAt only inside a candidate's
  window. It now marks the row/32-column blocks those windows reach and
  keeps its sliding sums everywhere but skips the per-pixel test elsewhere;
  the bits it leaves unset are never read.

md5-identical output on four sets (GPU) and three (CPU-only). CPU-only
16M: 203 -> 177 s and 303 -> 263 s; GPU 16M 23.8 -> 22.8 s (the pre-scan's
finder).

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

image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp

@@ -21,6 +21,7 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &
     ring_thr.assign(nbins, 0.0f);
     ring_bkg.assign(nbins, NAN);
     ring_bits.assign(OutputSize(), 0);
+    ring_hist.assign(nbins * HIST_VALUES, 0);
 }
 
 // Per-ring background statistics with iterated peak exclusion following peakfinder8:
@@ -37,19 +38,44 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image
     std::fill(ring_sum2.begin(), ring_sum2.end(), 0);
     std::fill(ring_cnt.begin(), ring_cnt.end(), 0);
 
-    for (size_t pxl = 0; pxl < npix; ++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 (clip_k > 0.0f) {
+    if (clip_k <= 0.0f) {
+        std::fill(ring_hist.begin(), ring_hist.end(), 0);
+        ring_overflow.clear();
+        for (size_t pxl = 0; pxl < npix; ++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)
+            ring_sum[b] += v;
+            ring_sum2[b] += static_cast<uint64_t>(static_cast<int64_t>(v) * v);
+            ring_cnt[b] += 1;
+            if (v >= 0 && v < HIST_VALUES)
+                ring_hist[b * HIST_VALUES + v] += 1;
+            else
+                ring_overflow.emplace_back(b, v);
+        }
+    } else {
+        // A sigma-clip pass over the plain pass's values: each distinct value of a ring meets the same
+        // test the pixels holding it would, and its pixels are added as a count.
+        const auto keep = [&](uint16_t b, int32_t v) {
             const float lo = ring_mean[b] - clip_k * ring_sigma[b];
             const float hi = ring_mean[b] + clip_k * ring_sigma[b];
-            if (v < lo || v > hi) continue;                // exclude peaks / outliers
+            return !(v < lo || v > hi);                        // exclude peaks / outliers
+        };
+        for (size_t b = 0; b < nbins; ++b)
+            for (int32_t v = 0; v < HIST_VALUES; ++v) {
+                const uint32_t n = ring_hist[b * HIST_VALUES + v];
+                if (n == 0 || !keep(static_cast<uint16_t>(b), v)) continue;
+                ring_sum[b] += static_cast<int64_t>(n) * v;
+                ring_sum2[b] += static_cast<uint64_t>(n) * static_cast<uint64_t>(static_cast<int64_t>(v) * v);
+                ring_cnt[b] += n;
+            }
+        for (const auto &[b, v] : ring_overflow) {
+            if (!keep(b, v)) continue;
+            ring_sum[b] += v;
+            ring_sum2[b] += static_cast<uint64_t>(static_cast<int64_t>(v) * v);
+            ring_cnt[b] += 1;
         }
-        ring_sum[b] += v;
-        ring_sum2[b] += static_cast<uint64_t>(static_cast<int64_t>(v) * v);
-        ring_cnt[b] += 1;
     }
 
     for (size_t b = 0; b < nbins; ++b) {

image_analysis/spot_finding/AdaptiveSpotFinderCPU.h

@@ -63,6 +63,12 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinderCPU {
     // Pixels at or above their ring's threshold, packed like output_buffer. Intersected with the
     // local-box mask that ImageSpotFinderCPU::Detect leaves in output_buffer.
     std::vector<uint32_t> ring_bits;
+    // The plain pass's valid pixels as a per-ring histogram of their values (HIST_VALUES bins per
+    // ring) plus a list of the values outside it, so the two sigma-clip passes sum over distinct
+    // values instead of over the image again. Integer sums, so the same totals.
+    static constexpr int32_t HIST_VALUES = 1024;
+    std::vector<uint32_t> ring_hist;
+    std::vector<std::pair<uint16_t, int32_t>> ring_overflow;   // (ring, value)
 
     void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k);
     // Fill ring_bits from the thresholds of the current frame.

image_analysis/spot_finding/ImageSpotFinderCPU.cpp

@@ -50,6 +50,21 @@ bool StrongInWindow(int64_t pxl_val, int64_t sum, int64_t sum2, int64_t valid,
 void ImageSpotFinderCPU::DetectAt(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings,
                                   const std::vector<uint32_t> &candidates) {
     candidate_windows.clear();
+    // The first pass's bits are read only inside the window of a candidate (and at the candidate), so
+    // it tests only the pixels of the row/column blocks such a window reaches; the bits it leaves
+    // unset there are never read.
+    const int32_t nblocks = (width + 31) / 32;
+    first_pass_needed.assign(static_cast<size_t>(height) * nblocks, 0);
+    for (size_t w = 0; w < candidates.size(); w++)
+        for (uint32_t bits = candidates[w]; bits; bits &= bits - 1) {
+            const int64_t pxl = static_cast<int64_t>(w) * 32 + std::countr_zero(bits);
+            if (pxl >= static_cast<int64_t>(width) * height)
+                break;
+            const int32_t line = static_cast<int32_t>(pxl / width), col = static_cast<int32_t>(pxl % width);
+            for (int32_t y = std::max(line - NBX, 0); y <= std::min(line + NBX, height - 1); y++)
+                for (int32_t b = std::max(col - NBX, 0) / 32; b <= std::min(col + NBX, width - 1) / 32; b++)
+                    first_pass_needed[static_cast<size_t>(y) * nblocks + b] = 1;
+        }
     DetectPass(image, settings, nullptr, first_pass_buffer, candidates.data());
 
     std::fill(output_buffer.begin(), output_buffer.end(), 0);
@@ -196,7 +211,9 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image,
                 if (candidates && (candidates[pxl / 32] & (1U << bit)))
                     candidate_windows.push_back({pxl, sum, sum2, valid});
 
-                if (StrongInWindow(value_at(pxl), sum, sum2, valid, settings, strong2))
+                const bool needed = !candidates
+                                    || first_pass_needed[static_cast<size_t>(line) * ((width + 31) / 32) + col / 32];
+                if (needed && StrongInWindow(value_at(pxl), sum, sum2, valid, settings, strong2))
                     out.set(bit);
 
                 if (bit == 31) {

image_analysis/spot_finding/ImageSpotFinderCPU.h

@@ -37,6 +37,9 @@ class ImageSpotFinderCPU : public ImageSpotFinder {
         int64_t sum, sum2, valid;
     };
     std::vector<CandidateWindow> candidate_windows;
+    // Where the first pass's own result is read by DetectAt - within NBX of a candidate - per row and
+    // 32-column block. Elsewhere that pass keeps its sums but skips the test.
+    std::vector<uint8_t> first_pass_needed;
 
 protected:
     // Detect, but with the result wanted only at the candidate pixels (packed like output_buffer):