Merge sf-ceil into rc176: the GPU spot finder's local-box test, count threshold and sigma-clip bounds are the CPU's to the bit

RESULT-CHANGING in principle, unchanged in practice: diagnostic builds counted 0 disagreements of
the SNR and count tests over ~1.9e10 first-pass and 2.9e7 candidate tests on lyso_x06da_ref,
cytc_x10sa, myob_x10sa, yag_x10sa_20keV, 4nwv and thau_bl1a_3p8keV, and 2 clip-bound flips each on
lyso, yag and 4nwv; p.hkl bit-identical before and after on all six. New knife-edge parity test
fails on the old GPU code in both directions.

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

image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu

@@ -44,12 +44,14 @@ __device__ __forceinline__ void flush_ring(unsigned long long *s_sum, unsigned l
     }
 }
 
-// The sigma-clip test, with both bounds as single fused multiply-adds. Written out rather than left to
-// the compiler, which fuses `mean - clip_k * sigma` in one kernel and may not in another - and the
-// histogram clip below must draw exactly the line the image pass draws.
+// The sigma-clip test, with both bounds rounded as AdaptiveSpotFinderCPU::ClipRings rounds them: the
+// product, then the sum, never fused. Written out rather than left to the compiler, so that every
+// kernel - the histogram clip below and the image passes - draws exactly the CPU's line. A fused
+// bound differs in the last bit, and where that crosses an integer every pixel of that value in the
+// ring moves in or out of the clip.
 __device__ __forceinline__ bool clip_keep(float fv, float mean, float sigma, float clip_k) {
-    const float lo = __fmaf_rn(-clip_k, sigma, mean);
-    const float hi = __fmaf_rn(clip_k, sigma, mean);
+    const float lo = __fsub_rn(mean, __fmul_rn(clip_k, sigma));
+    const float hi = __fadd_rn(mean, __fmul_rn(clip_k, sigma));
     return !(fv < lo || fv > hi);
 }
 

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