From 9c4e22465a0e0beb1b39c0794fbb6a653b923fc4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 17 Oct 2026 13:18:19 +0200 Subject: [PATCH 1/2] 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 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- .../spot_finding/ImageSpotFinderGPU.cu | 10 +++- tests/ImageSpotFinderGPUTest.cpp | 58 +++++++++++++++++++ 2 files changed, 65 insertions(+), 3 deletions(-) diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index 72abcfe2d..4d8ece36c 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/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(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) diff --git a/tests/ImageSpotFinderGPUTest.cpp b/tests/ImageSpotFinderGPUTest.cpp index 280af1c2c..9b4a409ff 100644 --- a/tests/ImageSpotFinderGPUTest.cpp +++ b/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 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(width), static_cast(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 From 7539f216dd95c805a4f1d3a70fcabf315ec5be15 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 17 Oct 2026 16:06:58 +0200 Subject: [PATCH 2/2] AdaptiveSpotFinderGPU: sigma-clip bounds rounded as the CPU rounds them (RESULT-CHANGING) RESULT-CHANGING for the GPU adaptive spot finder, only where a clip bound sat within its last bit of an integer. clip_keep formed mean -/+ clip_k * sigma as one fused multiply-add each; AdaptiveSpotFinderCPU::ClipRings computes the product and then the sum (contraction is off project-wide). The two bounds differ in the last bit, and where that crosses an integer pixel value every pixel of that value in the ring moves in or out of the clip, so the ring mean/sigma - and the ring's threshold - differ. clip_keep now rounds product then sum (__fmul_rn, __fsub_rn / __fadd_rn), still in one place for every kernel. Measured: a diagnostic build evaluating both forms in every clip_keep found 2 differing tests each on lyso_x06da_ref (of 4.7e7), yag_x10sa_20keV (of 2.2e8) and 4nwv (of 7.0e7), none on cytc_x10sa, myob_x10sa and thau_bl1a_3p8keV. The battery GPU runs before and after give bit-identical p.hkl on all six sets. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 1c9123f23..7dc5e36b0 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/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); }