spot_finding: fix the GPU finder's main loop, which the tests could not reach

Two bugs in analyze_pixel, both confined to the middle stage of the wave.

The kernel walks each wave's rows in three stages. The priming and drain
loops read prev_out and substitute INT32_MAX for a pixel the previous pass
found strong, exactly as the CPU finder's value_at() does on every read. The
main loop did not - it read the image raw. So in the second pass the pixels
the first pass found strong stayed in the background statistics, inflating the
local mean and variance, and the halo of every broad spot failed the
signal-to-noise test. The two engines therefore did not agree, despite
1a1e05ad1 having set out to make them.

Separately, shared_sum2 is an int64 accumulator but val*val and old*old were
computed in int32 at three sites. That wraps above 46340 while the detector
overloads around 1e6, so any window containing a bright pixel got a corrupted
variance. pixel_result already did the same arithmetic in 64 bits.

The reason this survived is worth recording: numberOfWaves is fixed at 32, so
a wave owns ceil(height/32) rows, and the main loop only runs while
front < rmax with front starting NBX+1 = 16 rows ahead. At the existing test's
100 rows a wave owns 4 rows and THE MAIN LOOP NEVER EXECUTES - every row goes
through priming or drain, and the parity test passed on the broken kernel. On
a 4362-row detector frame a wave owns 137 rows and the main loop carries about
121 of them, so the bug covered roughly 88% of a real image.

The new tall-image test is sized to the partitioning rather than to
convenience: 1024 rows, spots placed inside the main-loop region, one core at
100000 to exercise the overflow. Against the unfixed kernel it reports GPU 25
pixels where the CPU finds 49 and fails six assertions.

The default rugnux path is unaffected because it uses the adaptive finder, and
the 37-crystal battery is identical crystal for crystal. The classic finder is
what SpotFindingSettings defaults to, so this is the online receiver's path;
measured there with --no-adaptive-spots, R-meas improves 11.2% -> 10.9% and
CC1/2 98.8% -> 98.9% on one crystal, with multiplicity up on both tried.

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

image_analysis/spot_finding/ImageSpotFinderGPU.cu

@@ -153,7 +153,7 @@ __global__ void analyze_pixel(const int32_t *in, uint32_t *prev_out, uint32_t *o
             shared_val[(front % window) * blockDim.x + threadIdx.x] = val;
             if (val != INT32_MAX && val != INT32_MIN) {
                 shared_sum[threadIdx.x] += val;
-                shared_sum2[threadIdx.x] += val * val;
+                shared_sum2[threadIdx.x] += static_cast<int64_t>(val) * val;   // see the main loop
                 shared_count[threadIdx.x] += 1;
             }
         }
@@ -183,15 +183,26 @@ __global__ void analyze_pixel(const int32_t *in, uint32_t *prev_out, uint32_t *o
                 cnt = 1;    // bring count to normal
             }
             int32_t val = INT32_MIN;    // past the last row nothing enters the window
-            if (front < static_cast<int32_t>(params.height))
-                val = in[front * params.width + col];
+            if (front < static_cast<int32_t>(params.height)) {
+                // A pixel the previous pass found strong reads as INT32_MAX, exactly as the priming
+                // and drain loops above and below do, and as the CPU finder's value_at() does on
+                // every read. Without it the second pass counts pass-1 spot pixels as background
+                // over the whole middle of the image, which is where nearly all rows are.
+                const int32_t npixel = front * params.width + col;
+                const bool sat = ((prev_out[npixel / 32] & (1U << (npixel % 32))) != 0);
+                val = sat ? INT32_MAX : in[npixel];
+            }
             shared_val[(front % window) * blockDim.x + threadIdx.x] = val;
             if (val == INT32_MAX || val == INT32_MIN) {
                 val = 0;    // no effect value
                 cnt -= 1;   // count diff from normal
             }
             shared_sum[threadIdx.x] += val - old;
-            shared_sum2[threadIdx.x] += val * val - old * old;
+            // 64-bit squares: the accumulator is int64, but val*val in int32 wraps above 46340 and
+            // the detector saturates far higher (overload ~1e6), which corrupted the variance for
+            // every window containing a bright pixel.
+            shared_sum2[threadIdx.x] += static_cast<int64_t>(val) * val
+                                        - static_cast<int64_t>(old) * old;
             shared_count[threadIdx.x] += cnt;
         }
         front++;
@@ -230,7 +241,8 @@ __global__ void analyze_pixel(const int32_t *in, uint32_t *prev_out, uint32_t *o
                     cnt += 1;   // count diff from normal
             }
             shared_sum[threadIdx.x] += val - old;
-            shared_sum2[threadIdx.x] += val * val - old * old;
+            shared_sum2[threadIdx.x] += static_cast<int64_t>(val) * val
+                                        - static_cast<int64_t>(old) * old;   // see the main loop
             shared_count[threadIdx.x] += cnt;
         }
         front++;

tests/ImageSpotFinderGPUTest.cpp

@@ -210,4 +210,74 @@ TEST_CASE("ImageSpotFinder_CPU_GPU_Parity", "[ImageSpotFinder]") {
     CHECK(cpu_spots[0].PixelCount() == 49);
 }
 
+// The same comparison on a TALL image, which is what makes it bite.
+//
+// The GPU kernel splits the image into ImageSpotFinderGPU::numberOfWaves (32) horizontal waves and
+// each wave walks its rows in three stages: a priming loop, a main loop, and a drain loop. The main
+// loop only runs while front < rmax, and front starts NBX+1 rows ahead of the wave's first row - so
+// for an image of 100 rows, where a wave owns ceil(100/32) = 4 rows, THE MAIN LOOP NEVER EXECUTES.
+// Every row goes through priming/drain. On a real detector frame (4362 rows -> 137 rows per wave)
+// the main loop instead carries ~121 of every 137 rows, i.e. almost the whole image.
+//
+// That is why a bug living only in the main loop survived: it was unreachable at the test's size.
+// This case is 1024 rows (32 per wave, ~16 through the main loop) so the path is covered, and the
+// spots are placed deep inside a wave rather than at its edges. Values above 46340 are included
+// deliberately: the window keeps a sum of squares, and squaring in 32 bits overflows past that,
+// while a real detector saturates around 1e6.
+TEST_CASE("ImageSpotFinder_CPU_GPU_Parity_TallImage", "[ImageSpotFinder]") {
+    if (get_gpu_count() == 0)
+        SKIP("No CUDA GPU present");
+
+    const size_t width = 128, height = 1024;
+
+    ImagePreprocessorBufferGPU gpu_buffer(width * height);
+    ImagePreprocessorBuffer cpu_buffer(width * height);
+
+    // Spots at rows a wave reaches through its MAIN loop, not its priming or drain rows.
+    const struct { int cx, cy; int32_t core; } spots[] = {
+        {64,  80, 300},        // ordinary bright spot
+        {64, 400, 100000},     // above 46340: 32-bit squaring wraps, 64-bit does not
+        {64, 720, 300},
+    };
+    auto fill = [&](ImagePreprocessorBuffer &b) {
+        for (size_t i = 0; i < width * height; i++)
+            b[i] = (i % 2) * 5 + 5;
+        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 : s.core;
+                }
+    };
+    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));
+    cpu.SetResolutionMask(res_mask);
+    const auto cpu_spots = cpu.Run(cpu_buffer, settings);
+    const auto gpu_spots = run_gpu_and_collect_spots(gpu_buffer, width, height, settings, res_mask);
+
+    REQUIRE(cpu_spots.size() == 3);   // all three found, so "agreeing on nothing" cannot pass
+    REQUIRE(gpu_spots.size() == cpu_spots.size());
+
+    for (size_t i = 0; i < cpu_spots.size(); i++) {
+        INFO("spot " << 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));
+        CHECK(cpu_spots[i].PixelCount() == gpu_spots[i].PixelCount());
+        CHECK(cpu_spots[i].Count() == gpu_spots[i].Count());
+        CHECK(cpu_spots[i].PixelCount() == 49);   // core + ring, i.e. the second pass did its job
+    }
+}
+
 #endif