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>

docs/CPU_DATA_ANALYSIS.md

@@ -181,6 +181,8 @@ A pixel is considered strong if:
 
 This is equivalent to a local z-score criterion but implemented in integer arithmetic to be robust and fast.
 
+The test is applied in **two passes** over the image. The first is as described above. The second repeats it with every pixel found strong by the first excluded from the local background — it is treated exactly like a saturated pixel, so it contributes to no window it falls into and stays strong itself. This matters for any spot wide enough to reach into its own background box: on a single pass such a spot inflates the mean and variance it is then tested against, and its outer pixels fail the criterion. Excluding the core recovers them, so the spot is reported with its true extent rather than its brightest few pixels. Both the CPU and GPU implementations run these two passes and return the same spot list for the same frame.
+
 Special cases:
 - saturated pixels can be forced to “strong” (useful for detecting overloaded Bragg spots),
 - invalid pixels are never strong.

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();
 }

image_analysis/spot_finding/ImageSpotFinderCPU.h

@@ -18,6 +18,15 @@
 // and it aims for 100 ms execution
 
 class ImageSpotFinderCPU : public ImageSpotFinder {
+    // Output of the first pass, and the exclusion mask for the second.
+    std::vector<uint32_t> first_pass_buffer;
+
+    // One detection pass. prev_strong is the previous pass's bitmap, or nullptr for the first pass;
+    // a pixel set in it is treated as INT32_MAX, i.e. kept out of every local background it would
+    // fall into, and reported strong.
+    void DetectPass(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings,
+                    const uint32_t *prev_strong, std::vector<uint32_t> &out);
+
 public:
     ImageSpotFinderCPU(int32_t width, int32_t height);
     void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override;

tests/ImageSpotFinderGPUTest.cpp

@@ -7,6 +7,7 @@
 #ifdef JFJOCH_USE_CUDA
 
 #include "../image_analysis/spot_finding/ImageSpotFinderGPU.h"
+#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
 #include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h"
 
 static void fill_test_image(ImagePreprocessorBuffer& buffer, size_t width, size_t height) {
@@ -133,4 +134,78 @@ TEST_CASE("ImageSpotFinderGPU_20M") {
     REQUIRE(spots[1].RawCoord().y == 50);
 }
 
+// The two finders must return the same spots for the same frame - a dataset processed on a machine
+// without a GPU has to give the same answer as one processed with it.
+//
+// The spots here are deliberately broad. Both finders measure a pixel against a 31x31 local
+// background, so a blob several pixels across sits inside its own background window and inflates
+// the mean and variance it is tested against. That is what the second pass exists to undo: it
+// recomputes the background with the pixels found strong by the first pass excluded. A single-pass
+// finder loses the outer pixels of every broad spot, so this comparison fails unless both sides run
+// the same two passes.
+TEST_CASE("ImageSpotFinder_CPU_GPU_Parity", "[ImageSpotFinder]") {
+    if (get_gpu_count() == 0)
+        SKIP("No CUDA GPU present");
+
+    const size_t width = 100, height = 100;
+
+    ImagePreprocessorBufferGPU gpu_buffer(width * height);
+    ImagePreprocessorBuffer cpu_buffer(width * height);
+
+    // Background alternating 5/10 (mean 7.5, sd 2.5), plus two spots shaped to make the second pass
+    // matter: a bright 5x5 core (300) with a thin one-pixel ring around it (25). The ring is well
+    // above the clean background, but the core sitting inside the ring's own 31x31 window drags that
+    // window's mean to ~16 and its sd to ~46, so on a single pass the ring fails the SNR test and the
+    // spot comes out as the 25-pixel core. The second pass takes the core out of the background and
+    // the ring passes, giving 49 pixels. The ring is kept thin on purpose: a wide halo would swamp
+    // its own background and stay undetectable either way.
+    auto fill = [&](ImagePreprocessorBuffer &b) {
+        for (size_t i = 0; i < width * height; i++)
+            b[i] = (i % 2) * 5 + 5;
+
+        const struct { int cx, cy; } spots[] = {{50, 50}, {22, 74}};
+        for (const auto &s : spots) {
+            for (int dy = -3; dy <= 3; dy++) {
+                for (int dx = -3; dx <= 3; dx++) {
+                    const bool ring = std::abs(dx) == 3 || std::abs(dy) == 3;
+                    b[(s.cy + dy) * width + (s.cx + dx)] = ring ? 25 : 300;
+                }
+            }
+        }
+    };
+    fill(gpu_buffer);
+    fill(cpu_buffer);
+
+    SpotFindingSettings settings{
+        .signal_to_noise_threshold = 3.0,
+        .photon_count_threshold = 0,
+        .min_pix_per_spot = 1,
+        .max_pix_per_spot = 1000,
+        .high_resolution_limit = 0.5,
+        .low_resolution_limit = 3.0,
+    };
+    const std::vector<bool> res_mask(width * height, false);
+
+    ImageSpotFinderCPU cpu(static_cast<int32_t>(width), static_cast<int32_t>(height));
+    const auto cpu_spots = cpu.Run(cpu_buffer, settings, res_mask);
+
+    const auto gpu_spots = run_gpu_and_collect_spots(gpu_buffer, width, height, settings, res_mask);
+
+    REQUIRE(cpu_spots.size() == gpu_spots.size());
+    REQUIRE(cpu_spots.size() == 2);   // guard against both finding nothing and "agreeing"
+
+    for (size_t i = 0; i < cpu_spots.size(); i++) {
+        CHECK(cpu_spots[i].RawCoord().x == Catch::Approx(gpu_spots[i].RawCoord().x));
+        CHECK(cpu_spots[i].RawCoord().y == Catch::Approx(gpu_spots[i].RawCoord().y));
+        // The pixel set is what the second pass changes, so compare it rather than the centroid,
+        // which stays put for a symmetric spot whether or not the halo was picked up.
+        CHECK(cpu_spots[i].PixelCount() == gpu_spots[i].PixelCount());
+        CHECK(cpu_spots[i].Count() == gpu_spots[i].Count());
+    }
+
+    // Both finders must reach past the bright core into the ring: 49 pixels, not the 25 of the core
+    // alone. Without this the comparison above would still pass if both ran a single pass.
+    CHECK(cpu_spots[0].PixelCount() == 49);
+}
+
 #endif