spot_finding: run the same two passes on the CPU as on the GPU

ImageSpotFinderGPU::Detect launches its kernel twice, feeding the first
pass's strong-pixel bitmap back in so the second recomputes each local
background with those pixels excluded and keeps them strong. The CPU
finder ran a single pass, so the two returned different spot lists for the
same frame and a dataset processed without a GPU did not match one
processed with it.

It matters for any spot wide enough to reach into its own 31x31 background
box: the spot inflates the mean and variance it is then tested against, so
its outer pixels fail the SNR test. On the test image added here - a 5x5
core at 300 counts with a one-pixel ring at 25 - a single pass returns the
25-pixel core and 7500 counts where two passes return the full 49 pixels
and 8100.

pxl_val also becomes int64_t, matching the GPU's pixel_result signature.
It was int32_t, so pxl_val * pxl_val overflowed above 46341 counts even
though the surrounding sums were already 64-bit.

The new parity test compares PixelCount and Count, not just the centroid,
which does not move for a symmetric spot whether or not the ring was
picked up; it was confirmed to fail against the old single-pass CPU.

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

image_analysis/spot_finding/ImageSpotFinderCPU.cpp

@@ -7,12 +7,35 @@
 #include "StrongPixelSet.h"
 
 ImageSpotFinderCPU::ImageSpotFinderCPU(int32_t in_width, int32_t in_height)
-    : ImageSpotFinder(in_width, in_height) {}
+    : ImageSpotFinder(in_width, in_height), first_pass_buffer(OutputSize(), 0) {}
 
 void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                                 const SpotFindingSettings &settings) {
+    // Two passes, as ImageSpotFinderGPU::Detect does. The second recomputes every local background
+    // with the pixels the first found strong taken out of it, and keeps those pixels strong. It
+    // matters because a spot wide enough to reach into its own background window inflates the mean
+    // and variance it is then tested against, so its outer pixels fail the SNR test on a single
+    // pass. The GPU has always done this; running one pass here made the two finders return
+    // different spot lists for the same frame.
+    DetectPass(image, settings, nullptr, first_pass_buffer);
+    DetectPass(image, settings, first_pass_buffer.data(), output_buffer);
+}
+
+void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image,
+                                    const SpotFindingSettings &settings,
+                                    const uint32_t *prev_strong,
+                                    std::vector<uint32_t> &out_buffer) {
     for (int i = 0; i < OutputSize(); i++)
-        output_buffer[i] = 0;
+        out_buffer[i] = 0;
+
+    // A pixel found strong by the previous pass reads as INT32_MAX, which the accumulation below
+    // already skips and the acceptance test below already takes as strong - the same substitution
+    // the GPU kernel makes when it reads prev_out.
+    auto value_at = [&](int32_t pxl) -> int32_t {
+        if (prev_strong && (prev_strong[pxl / 32] & (1U << (pxl % 32))))
+            return INT32_MAX;
+        return image[pxl];
+    };
 
     std::bitset<32> out = 0;
 
@@ -20,12 +43,12 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
         if (settings.photon_count_threshold > 0) {
             for (int pxl = 0; pxl < height * width; pxl++) {
                 int32_t bit = pxl % 32;
-                int32_t pxl_val = image[pxl];
+                int32_t pxl_val = value_at(pxl);
                 if (pxl_val == INT32_MAX || (pxl_val > settings.photon_count_threshold && pxl_val != INT32_MIN))
                     out.set(bit);
 
                 if (bit == 31) {
-                    output_buffer[pxl / 32] = out.to_ulong();
+                    out_buffer[pxl / 32] = out.to_ulong();
                     out.reset();
                 }
             }
@@ -44,8 +67,8 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
             for (int col = 0; col < width; col++) {
                 auto pxl = line * width + col;
 
-                if (image[pxl] != INT32_MAX && image[pxl] != INT32_MIN) {
-                    int64_t tmp = image[pxl];
+                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;
@@ -58,8 +81,8 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                 if (line < height - NBX) {
                     auto pxl = (line + NBX) * width + col;
 
-                    if (image[pxl] != INT32_MAX && image[pxl] != INT32_MIN) {
-                        int64_t tmp = image[pxl];
+                    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;
@@ -68,8 +91,8 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
 
                 if (line >= NBX + 1) {
                     auto pxl = (line - (NBX + 1)) * width + col;
-                    if (image[pxl] != INT32_MAX && image[pxl] != INT32_MIN) {
-                        int64_t tmp = image[pxl];
+                    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;
@@ -101,7 +124,7 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                 }
                 const int32_t pxl = line * width + col;
 
-                int32_t pxl_val = image[pxl];
+                const int64_t pxl_val = value_at(pxl);
                 int64_t sum_local = sum - pxl_val;
                 int64_t sum2_local = sum2 - pxl_val * pxl_val;
                 int64_t valid_local = valid - 1;
@@ -111,7 +134,7 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
 
                 const int32_t bit = pxl % 32;
 
-                if ((pxl_val == INT32_MAX) // saturated pixel is accepted always
+                if ((pxl_val == INT32_MAX) // saturated pixel, or strong in the previous pass, is accepted always
                     || ((pxl_val != INT32_MIN && // pixel is not bad pixel
                          valid_local > MIN_VALID_PIXELS && // too many bad pixels around will give poor statistics
                          (pxl_val > settings.photon_count_threshold) && // pixel is above count threshold
@@ -121,7 +144,7 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                     out.set(bit);
 
                 if (bit == 31) {
-                    output_buffer[pxl / 32] = out.to_ulong();
+                    out_buffer[pxl / 32] = out.to_ulong();
                     out.reset() ;
                 }
             }
@@ -129,5 +152,5 @@ void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
     }
 
     if (height * width % 32 != 0)
-        output_buffer[OutputSize() - 1] = out.to_ulong();
+        out_buffer[OutputSize() - 1] = out.to_ulong();
 }