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