From 04450eb618d2a95fddcf381c08b2074360c2d5ef Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 30 Jul 2026 11:05:54 +0200 Subject: [PATCH] Adaptive spot finder: sum the rings across blocks in double The ring sigma is the cancelling difference sum2/n - m^2, and both sums were float accumulated by atomics whose order is arbitrary. Two costs: the cancellation left only ~4 digits in the variance, and the ordering moved the resulting threshold by ~0.05 counts between runs - enough to flip a pixel sitting on the hard "value >= threshold" test, and with it a connected component's size. So the GPU engine did not reproduce the CPU one and did not reproduce itself. Only the accumulators that span blocks are widened. The per-block staging stays float, because a block contributes a few dozen similar-magnitude pixels to a ring and there is nothing to lose there - that also keeps the shared-memory footprint of the hot loop, and hence its occupancy, exactly as it was: measured on a 4.5 MP frame, 0.960 vs 0.966 ms/frame (40.9x over the CPU path, unchanged). finalize_rings now does the cancellation in double and rounds to float last, which is what AdaptiveSpotFinderCPU::AccumulateRings does. The device properties are also read from the current device rather than device 0; callers round-robin engines across GPUs. Co-Authored-By: Claude Opus 5 (1M context) --- .../spot_finding/AdaptiveSpotFinderGPU.cu | 51 +++++++++++-------- .../spot_finding/AdaptiveSpotFinderGPU.h | 10 ++-- 2 files changed, 37 insertions(+), 24 deletions(-) diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 3a0b4d5d..6fd29494 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -24,10 +24,14 @@ __global__ void reduce_rings_shared( const float *__restrict__ sigma, float clip_k, bool accumulate_corrected, - float *__restrict__ sum, float *__restrict__ sum2, uint32_t *__restrict__ count, + double *__restrict__ sum, double *__restrict__ sum2, uint32_t *__restrict__ count, float *__restrict__ sum_corr, float *__restrict__ sum2_corr, size_t npix, int nbins) { + // The per-block staging stays float: a block contributes only a few dozen pixels to a given ring, + // all of similar magnitude, so there is nothing to lose there - and float keeps the shared footprint + // (and hence the occupancy) of the hot loop unchanged. The precision that matters is in the sum over + // ALL blocks and in the cancelling difference sum2/n - m^2 that follows it, so those are double. extern __shared__ float sh[]; float *s_sum = sh; float *s_sum2 = &s_sum[nbins]; @@ -69,8 +73,8 @@ __global__ void reduce_rings_shared( __syncthreads(); for (int i = threadIdx.x; i < nbins; i += blockDim.x) { - atomicAdd(&sum[i], s_sum[i]); - atomicAdd(&sum2[i], s_sum2[i]); + atomicAdd(&sum[i], static_cast(s_sum[i])); + atomicAdd(&sum2[i], static_cast(s_sum2[i])); atomicAdd(&count[i], s_count[i]); if (accumulate_corrected) { atomicAdd(&sum_corr[i], s_sum_corr[i]); @@ -89,7 +93,7 @@ __global__ void reduce_rings_global( const float *__restrict__ sigma, float clip_k, bool accumulate_corrected, - float *__restrict__ sum, float *__restrict__ sum2, uint32_t *__restrict__ count, + double *__restrict__ sum, double *__restrict__ sum2, uint32_t *__restrict__ count, float *__restrict__ sum_corr, float *__restrict__ sum2_corr, size_t npix, int nbins) { @@ -104,8 +108,9 @@ __global__ void reduce_rings_global( const float hi = mean[b] + clip_k * sigma[b]; if (fv < lo || fv > hi) continue; } - atomicAdd(&sum[b], fv); - atomicAdd(&sum2[b], fv * fv); + const double dv = static_cast(v); + atomicAdd(&sum[b], dv); + atomicAdd(&sum2[b], dv * dv); atomicAdd(&count[b], 1u); if (accumulate_corrected) { const float cv = fv * corrections[idx]; @@ -117,15 +122,17 @@ __global__ void reduce_rings_global( // Per-ring mean/sigma from the current raw accumulators. Rings with no pixels this pass keep their // previous value (matches the CPU, which leaves ring_mean/ring_sigma untouched when the count is 0). -__global__ void finalize_rings(const float *__restrict__ sum, const float *__restrict__ sum2, +__global__ void finalize_rings(const double *__restrict__ sum, const double *__restrict__ sum2, const uint32_t *__restrict__ count, float *__restrict__ mean, float *__restrict__ sigma, int nbins) { for (int b = blockIdx.x * blockDim.x + threadIdx.x; b < nbins; b += blockDim.x * gridDim.x) { if (count[b] > 0) { - const float m = sum[b] / count[b]; - const float var = fmaxf(0.0f, sum2[b] / count[b] - m * m); - mean[b] = m; - sigma[b] = sqrtf(var); + // In double, then rounded to float for the clip predicate - the same two steps, in the same + // order and the same types, as AdaptiveSpotFinderCPU::AccumulateRings. + const double m = sum[b] / count[b]; + const double var = fmax(0.0, sum2[b] / count[b] - m * m); + mean[b] = static_cast(m); + sigma[b] = static_cast(sqrt(var)); } } } @@ -183,8 +190,12 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & output_buffer_reg(output_buffer), last_profile(in_mapping) { + // The current device, not device 0: callers round-robin engines across GPUs, so device 0's shared + // memory and SM count can belong to a different card than the one these kernels launch on. + int device = 0; + cuda_err(cudaGetDevice(&device)); cudaDeviceProp prop{}; - cuda_err(cudaGetDeviceProperties(&prop, 0)); + cuda_err(cudaGetDeviceProperties(&prop, device)); reduce_blocks = 4 * prop.multiProcessorCount; flag_blocks = 4 * prop.multiProcessorCount; @@ -249,8 +260,8 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting if (host_count[b] < adaptive_threshold::MIN_RING_PIXELS) { host_thr[b] = g_thr; } else { - const double m = static_cast(host_sum[b]) / host_count[b]; - const double var = std::max(0.0, static_cast(host_sum2[b]) / host_count[b] - m * m); + const double m = host_sum[b] / host_count[b]; + const double var = std::max(0.0, host_sum2[b] / host_count[b] - m * m); host_thr[b] = adaptive_threshold::RingThreshold(static_cast(m), static_cast(std::sqrt(var)), p, z); } @@ -264,8 +275,8 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, "AdaptiveSpotFinderGPU::Detect: mismatch in pixel size"); // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) --- - cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(float) * nbins, *stream)); - cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(float) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(double) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(double) * nbins, *stream)); cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * nbins, *stream)); cuda_err(cudaMemsetAsync(gpu_mean, 0, sizeof(float) * nbins, *stream)); cuda_err(cudaMemsetAsync(gpu_sigma, 0, sizeof(float) * nbins, *stream)); @@ -282,16 +293,16 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, cuda_err(cudaMemcpyAsync(prof_count.data(), gpu_count, sizeof(uint32_t) * nbins, cudaMemcpyDeviceToHost, *stream)); for (int pass = 0; pass < 2; ++pass) { - cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(float) * nbins, *stream)); - cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(float) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(double) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(double) * nbins, *stream)); cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * nbins, *stream)); ReducePass(image, 3.0f, false); FinalizeStats(); } // Snapshot the clipped raw stats that drive the threshold. - cuda_err(cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(float) * nbins, cudaMemcpyDeviceToHost, *stream)); - cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(float) * nbins, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(double) * nbins, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(double) * nbins, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(host_count.data(), gpu_count, sizeof(uint32_t) * nbins, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaStreamSynchronize(*stream)); diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index 78740f2a..9481de42 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -51,8 +51,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaDevicePtr gpu_corrections; // Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold. - CudaDevicePtr gpu_sum; - CudaDevicePtr gpu_sum2; + // double, like the CPU engine's ring accumulators: the ring sigma is the cancelling difference + // sum2/n - m^2, and the block atomics that fill these arrive in an arbitrary order. + CudaDevicePtr gpu_sum; + CudaDevicePtr gpu_sum2; CudaDevicePtr gpu_count; CudaDevicePtr gpu_mean; // per-ring raw mean (clip predicate) CudaDevicePtr gpu_sigma; // per-ring raw sigma (clip predicate) @@ -66,8 +68,8 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaDevicePtr gpu_strong; // Host mirrors of the small per-ring transfers. - std::vector host_sum; // clipped raw sum } input to the host threshold computation - std::vector host_sum2; // clipped raw sum^2 } + std::vector host_sum; // clipped raw sum } input to the host threshold computation + std::vector host_sum2; // clipped raw sum^2 } std::vector host_count; // clipped raw count } std::vector host_thr; // per-ring threshold (empty -> frame had no valid pixels) std::vector prof_sum; // plain corrected sum } azimuthal-integration profile