diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index 4219bc96..27f3f82a 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/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(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(params.height)) - val = in[front * params.width + col]; + if (front < static_cast(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(val) * val + - static_cast(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(val) * val + - static_cast(old) * old; // see the main loop shared_count[threadIdx.x] += cnt; } front++; diff --git a/tests/ImageSpotFinderGPUTest.cpp b/tests/ImageSpotFinderGPUTest.cpp index 84226cca..280af1c2 100644 --- a/tests/ImageSpotFinderGPUTest.cpp +++ b/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 res_mask(width * height, false); + + 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); + + 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