ImageSpotFinderGPU: the local-box test is the CPU's to the bit (RESULT-CHANGING)

RESULT-CHANGING for the GPU spot finder (classic and adaptive, which shares this kernel for its
local-box passes) - only on pixels where the two used to disagree.

Two differences between pixel_result and the CPU's StrongInWindow:

- SNR test. Both form T = float(var) * strong2 in float. The CPU compares the int64 square
  (in_minus_mean)^2 with (int64) ceil(T); the GPU compared it with T directly, which converts the
  square to float. Below 2^24 the two differ only when the square equals ceil(T) with T not an
  integer (an exact-integer strong2 such as 9 or 16 never gives that); past 2^24, where T is an
  integer, the GPU's rounded square can land on T and fail a pixel the CPU accepts. The GPU now
  uses the CPU's integer-vs-ceil comparison.
- Count threshold. At photon_count_threshold = 0 the CPU still requires val > 0; the GPU took 0 as
  "no threshold". It matters only for a pixel <= 0 above a negative local mean (a
  pedestal-subtracted detector); the adaptive finder always runs this pass with threshold 0. The GPU
  now applies val > threshold unconditionally, with the threshold carried as int64 like the setting.

Measured: a diagnostic build evaluating both semantics in every pixel_result found zero
disagreements in ~2e10 first-pass and ~2.9e7 candidate tests over lyso_x06da_ref, cytc_x10sa,
myob_x10sa, yag_x10sa_20keV, 4nwv and thau_bl1a_3p8keV (JUNGFRAU); the battery GPU runs before and
after give bit-identical p.hkl on all six sets. The new knife-edge test constructs both cases (a
checkerboard window whose square is var * 25 exactly and rounds onto the float threshold; a 0 above
a negative checkerboard) and failed on the GPU before this change.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4

image_analysis/spot_finding/ImageSpotFinderGPU.cu

@@ -11,7 +11,7 @@ struct spot_parameters {
     int32_t width;
     int32_t height;
     float strong_pixel_threshold2;
-    int32_t count_threshold;
+    int64_t count_threshold;
 };
 
 // input      X x Y pixels array
@@ -66,6 +66,8 @@ __device__ __forceinline__ uint8_t pixel_result(const spot_parameters& params, c
     const int64_t var = count * sum2 - (sum * sum);   // This should be divided by ((2*NBX+1) * (2*NBY+1)-1)*((2*NBX+1) * (2*NBY+1))
     const int64_t in_minus_mean = val * count - sum;  // Should be divided by ((2*NBX+1) * (2*NBY+1));
 
+    // The CPU's test (StrongInWindow in ImageSpotFinderCPU.cpp), to the bit: the square stays an
+    // integer and meets the ceiling of the float threshold. Compared as float, it rounded past 2^24.
     const int64_t tmp1 = in_minus_mean * in_minus_mean;
     const float tmp2 = var * params.strong_pixel_threshold2;
     bool snr_criterion;
@@ -73,9 +75,11 @@ __device__ __forceinline__ uint8_t pixel_result(const spot_parameters& params, c
     if (params.strong_pixel_threshold2 == 0.0f)
         snr_criterion = true;
     else
-        snr_criterion = (count > ImageSpotFinder::MIN_VALID_PIXELS) && (in_minus_mean > 0) && (tmp1 > tmp2);
+        snr_criterion = (count > ImageSpotFinder::MIN_VALID_PIXELS) && (in_minus_mean > 0)
+                        && (tmp1 > static_cast<int64_t>(ceilf(tmp2)));
 
-    bool count_criterion = (params.count_threshold == 0.0f) || (val > params.count_threshold);
+    // As on the CPU, a threshold of 0 still asks for a positive value.
+    bool count_criterion = val > params.count_threshold;
 
     bool strong_pixel = snr_criterion && count_criterion;
     if (val == INT32_MAX)

tests/ImageSpotFinderGPUTest.cpp

@@ -280,4 +280,62 @@ TEST_CASE("ImageSpotFinder_CPU_GPU_Parity_TallImage", "[ImageSpotFinder]") {
     }
 }
 
+
+// The local-box SNR test on the knife edge, where the GPU used to answer differently from the CPU.
+// Both compare (in_minus_mean)^2 against var * strong2 computed in float; the CPU compares the
+// integer square with the ceiling of that float, the GPU compared the square converted to float -
+// which, for squares past 2^24, can round onto the threshold and fail a pixel the CPU accepts.
+//
+// The centre of a checkerboard of 0 and 57 holding 171 is such a pixel at 5 sigma: its 960
+// neighbours are 480 zeros and 480 of 57, so var = 748569600 and in_minus_mean = 136800, whose
+// square 18714240000 is var * 25 exactly. var * 25 rounds down in float to 18714238976, 1024 below
+// the square - half a float spacing - so the square rounds onto the threshold too. The CPU takes
+// the pixel (18714240000 > 18714238976), the GPU used not to.
+//
+// The second image is the count threshold: at photon_count_threshold = 0 the CPU still asks for a
+// positive value, the GPU took 0 to mean no threshold. A 0 above a negative background (the
+// checkerboard of -3 and -1, as a pedestal-subtracted detector gives) is 2 sigma above its mean.
+TEST_CASE("ImageSpotFinder_CPU_GPU_Parity_KnifeEdge", "[ImageSpotFinder]") {
+    if (get_gpu_count() == 0)
+        SKIP("No CUDA GPU present");
+
+    const size_t width = 64, height = 64;
+    const std::vector<bool> res_mask(width * height, false);
+
+    auto run_both = [&](int32_t even, int32_t odd, int32_t centre, float snr, size_t expected) {
+        ImagePreprocessorBufferGPU gpu_buffer(width * height);
+        ImagePreprocessorBuffer cpu_buffer(width * height);
+        for (size_t row = 0; row < height; row++)
+            for (size_t col = 0; col < width; col++) {
+                const int32_t v = ((row + col) % 2 == 0) ? even : odd;
+                gpu_buffer[row * width + col] = v;
+                cpu_buffer[row * width + col] = v;
+            }
+        gpu_buffer[32 * width + 32] = centre;
+        cpu_buffer[32 * width + 32] = centre;
+
+        SpotFindingSettings settings{
+            .signal_to_noise_threshold = snr,
+            .photon_count_threshold = 0,
+            .min_pix_per_spot = 1,
+            .max_pix_per_spot = 1000,
+            .high_resolution_limit = 0.5,
+            .low_resolution_limit = 3.0,
+        };
+        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);
+        CHECK(cpu_spots.size() == expected);
+        CHECK(gpu_spots.size() == cpu_spots.size());
+    };
+
+    SECTION("SNR test with the square past 2^24") {
+        run_both(0, 57, 171, 5.0f, 1);
+    }
+    SECTION("Zero count threshold on a negative background") {
+        run_both(-3, -1, 0, 1.5f, 0);
+    }
+}
+
 #endif