diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index 8ccdcf06..ca63d69a 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -13,8 +13,8 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping & static_cast(in_mapping.GetHeight())), mapping(in_mapping) { const size_t nbins = mapping.GetBinNumber(); - ring_sum.assign(nbins, 0.0); - ring_sum2.assign(nbins, 0.0); + ring_sum.assign(nbins, 0); + ring_sum2.assign(nbins, 0); ring_cnt.assign(nbins, 0); ring_mean.assign(nbins, 0.0f); ring_sigma.assign(nbins, 0.0f); @@ -29,8 +29,8 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image const size_t nbins = ring_sum.size(); const size_t npix = static_cast(width) * height; - std::fill(ring_sum.begin(), ring_sum.end(), 0.0); - std::fill(ring_sum2.begin(), ring_sum2.end(), 0.0); + std::fill(ring_sum.begin(), ring_sum.end(), 0); + std::fill(ring_sum2.begin(), ring_sum2.end(), 0); std::fill(ring_cnt.begin(), ring_cnt.end(), 0); for (size_t pxl = 0; pxl < npix; ++pxl) { @@ -44,14 +44,14 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image if (v < lo || v > hi) continue; // exclude peaks / outliers } ring_sum[b] += v; - ring_sum2[b] += static_cast(v) * v; + ring_sum2[b] += static_cast(static_cast(v) * v); ring_cnt[b] += 1; } for (size_t b = 0; b < nbins; ++b) { if (ring_cnt[b] > 0) { - const double m = ring_sum[b] / ring_cnt[b]; - const double var = std::max(0.0, ring_sum2[b] / ring_cnt[b] - m * m); + const double m = static_cast(ring_sum[b]) / ring_cnt[b]; + const double var = std::max(0.0, static_cast(ring_sum2[b]) / ring_cnt[b] - m * m); ring_mean[b] = static_cast(m); ring_sigma[b] = static_cast(std::sqrt(var)); } @@ -74,8 +74,8 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, double g_sum = 0.0, g_sum2 = 0.0; for (size_t b = 0; b < nbins; ++b) { n_total += ring_cnt[b]; - g_sum += ring_sum[b]; - g_sum2 += ring_sum2[b]; + g_sum += static_cast(ring_sum[b]); + g_sum2 += static_cast(ring_sum2[b]); } if (n_total == 0) { // Nothing valid to threshold against: leave no strong pixels for ExtractSpots to build on. diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index bb4aa220..fa67ebfc 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -28,8 +28,11 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { const AzimuthalIntegrationMapping &mapping; // per-ring scratch, sized to the mapping's bin count - std::vector ring_sum; - std::vector ring_sum2; + // Exact integers: a preprocessed pixel is an int32 and the sentinels are skipped, so v and v*v + // are exact in 64 bits. That is what lets the GPU engine reproduce these bit for bit - integer + // addition is associative, so its block atomics can arrive in any order. + std::vector ring_sum; + std::vector ring_sum2; std::vector ring_cnt; std::vector ring_mean; std::vector ring_sigma; diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 7dd43306..aa403eb0 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -24,24 +24,29 @@ __global__ void reduce_rings_shared( const float *__restrict__ sigma, float clip_k, bool accumulate_corrected, - double *__restrict__ sum, double *__restrict__ sum2, uint32_t *__restrict__ count, + unsigned long long *__restrict__ sum, unsigned long long *__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]; + // The raw accumulators are INTEGERS, not floats. A preprocessed pixel is an exact int32 (the + // masked and saturated sentinels are skipped below), so v and v*v are exact in 64 bits and + // integer addition is associative - which makes the ring mean and sigma, and therefore the + // detection threshold, independent of the order the atomics happen to arrive in. With float + // accumulators the threshold moved in its last bits between runs, and because detection compares + // an integer pixel value against it, a threshold crossing an integer flipped every pixel of that + // value in the ring at once. + // The CORRECTED sums stay float: they are a pixel value times a float correction, so there is no + // exact integer form. They feed the reported azimuthal profile, not the detection decision. + extern __shared__ unsigned long long sh[]; + unsigned long long *s_sum = sh; // signed value carried as two's complement + unsigned long long *s_sum2 = &s_sum[nbins]; uint32_t *s_count = reinterpret_cast(&s_sum2[nbins]); float *s_sum_corr = reinterpret_cast(&s_count[nbins]); float *s_sum2_corr = &s_sum_corr[nbins]; for (int i = threadIdx.x; i < nbins; i += blockDim.x) { - s_sum[i] = 0.0f; - s_sum2[i] = 0.0f; + s_sum[i] = 0; + s_sum2[i] = 0; s_count[i] = 0; if (accumulate_corrected) { s_sum_corr[i] = 0.0f; @@ -61,8 +66,8 @@ __global__ void reduce_rings_shared( const float hi = mean[b] + clip_k * sigma[b]; if (fv < lo || fv > hi) continue; } - atomicAdd(&s_sum[b], fv); - atomicAdd(&s_sum2[b], fv * fv); + atomicAdd(&s_sum[b], static_cast(static_cast(v))); + atomicAdd(&s_sum2[b], static_cast(static_cast(v) * v)); atomicAdd(&s_count[b], 1u); if (accumulate_corrected) { const float cv = fv * corrections[idx]; @@ -73,8 +78,8 @@ __global__ void reduce_rings_shared( __syncthreads(); for (int i = threadIdx.x; i < nbins; i += blockDim.x) { - atomicAdd(&sum[i], static_cast(s_sum[i])); - atomicAdd(&sum2[i], static_cast(s_sum2[i])); + atomicAdd(&sum[i], s_sum[i]); + atomicAdd(&sum2[i], s_sum2[i]); atomicAdd(&count[i], s_count[i]); if (accumulate_corrected) { atomicAdd(&sum_corr[i], s_sum_corr[i]); @@ -93,7 +98,7 @@ __global__ void reduce_rings_global( const float *__restrict__ sigma, float clip_k, bool accumulate_corrected, - double *__restrict__ sum, double *__restrict__ sum2, uint32_t *__restrict__ count, + unsigned long long *__restrict__ sum, unsigned long long *__restrict__ sum2, uint32_t *__restrict__ count, float *__restrict__ sum_corr, float *__restrict__ sum2_corr, size_t npix, int nbins) { @@ -108,9 +113,8 @@ __global__ void reduce_rings_global( const float hi = mean[b] + clip_k * sigma[b]; if (fv < lo || fv > hi) continue; } - const double dv = static_cast(v); - atomicAdd(&sum[b], dv); - atomicAdd(&sum2[b], dv * dv); + atomicAdd(&sum[b], static_cast(static_cast(v))); + atomicAdd(&sum2[b], static_cast(static_cast(v) * v)); atomicAdd(&count[b], 1u); if (accumulate_corrected) { const float cv = fv * corrections[idx]; @@ -122,15 +126,16 @@ __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 double *__restrict__ sum, const double *__restrict__ sum2, +__global__ void finalize_rings(const unsigned long long *__restrict__ sum, + const unsigned long long *__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) { // 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); + const double m = static_cast(static_cast(sum[b])) / count[b]; + const double var = fmax(0.0, static_cast(sum2[b]) / count[b] - m * m); mean[b] = static_cast(m); sigma[b] = static_cast(sqrt(var)); } @@ -206,8 +211,8 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & // measured no better (181 vs 175 us/launch). flag_blocks = 4 * prop.multiProcessorCount; - shared_plain = static_cast(nbins) * (4 * sizeof(float) + sizeof(uint32_t)); - shared_clip = static_cast(nbins) * (2 * sizeof(float) + sizeof(uint32_t)); + shared_plain = static_cast(nbins) * (2 * sizeof(unsigned long long) + 2 * sizeof(float) + sizeof(uint32_t)); + shared_clip = static_cast(nbins) * (2 * sizeof(unsigned long long) + sizeof(uint32_t)); use_shared = (shared_plain < prop.sharedMemPerBlock); // Both tables are functions of the detector geometry alone, so they are uploaded once per GPU and @@ -246,8 +251,8 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting double g_sum = 0.0, g_sum2 = 0.0; for (int b = 0; b < nbins; ++b) { n_total += host_count[b]; - g_sum += host_sum[b]; - g_sum2 += host_sum2[b]; + g_sum += static_cast(static_cast(host_sum[b])); + g_sum2 += static_cast(host_sum2[b]); } if (n_total == 0) { host_thr.clear(); @@ -269,8 +274,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 = host_sum[b] / host_count[b]; - const double var = std::max(0.0, host_sum2[b] / host_count[b] - m * m); + const double m = static_cast(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); host_thr[b] = adaptive_threshold::RingThreshold(static_cast(m), static_cast(std::sqrt(var)), p, z); } @@ -284,8 +289,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(double) * nbins, *stream)); - cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(double) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * 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)); @@ -302,16 +307,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(double) * nbins, *stream)); - cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(double) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * 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(double) * nbins, cudaMemcpyDeviceToHost, *stream)); - cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(double) * nbins, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(unsigned long long) * nbins, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(unsigned long long) * 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 7b5d2700..fb7dc749 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -56,8 +56,8 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { // Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold. // 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_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) @@ -71,8 +71,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 } (exact integers - see the kernel) 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 diff --git a/tests/AdaptiveSpotFinderGPUTest.cpp b/tests/AdaptiveSpotFinderGPUTest.cpp index 61b8653d..ce38ede7 100644 --- a/tests/AdaptiveSpotFinderGPUTest.cpp +++ b/tests/AdaptiveSpotFinderGPUTest.cpp @@ -187,15 +187,51 @@ TEST_CASE("AdaptiveSpotFinderGPU_RunToRunReproducible", "[AdaptiveSpotFinderGPU] auto stream = std::make_shared(); AdaptiveSpotFinderGPU gpu(mapping, stream); + // The ring accumulators are exact integers, so the threshold does not depend on the order the + // block atomics arrive in and the spot list has to be bit-identical every time - not merely + // close. Repeat enough times to give a scheduling-dependent threshold a chance to show itself: + // the effect it used to have was ~1 changed observation in a million, so a handful of repeats on + // a quiet background would not have caught it. const auto first = gpu.Run(buffer, settings); REQUIRE(first.size() > 0); - for (int repeat = 0; repeat < 4; repeat++) { + for (int repeat = 0; repeat < 50; repeat++) { const auto again = gpu.Run(buffer, settings); REQUIRE(again.size() == first.size()); REQUIRE(SortedCoords(again) == SortedCoords(first)); } } +// The threshold is computed from sums of int32 pixel values, so the two engines can agree EXACTLY +// rather than approximately - and that is the property worth locking, because it is what makes the +// GPU path's spot list independent of how the reduction happened to be scheduled. +TEST_CASE("AdaptiveSpotFinderGPU_RingStatsMatchCPUExactly", "[AdaptiveSpotFinderGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_RingStatsMatchCPUExactly"); + return; + } + + DiffractionExperiment x = MakeExperiment(); + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + FillTestImage(buffer, x); + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), + x.GetPixelsNum() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); + REQUIRE(cudaDeviceSynchronize() == cudaSuccess); + + const SpotFindingSettings settings = AdaptiveSettings(); + + auto stream = std::make_shared(); + AdaptiveSpotFinderGPU gpu(mapping, stream); + AdaptiveSpotFinderCPU cpu(mapping); + + const auto gpu_spots = gpu.Run(buffer, settings); + const auto cpu_spots = cpu.Run(buffer, settings); + REQUIRE(gpu_spots.size() == cpu_spots.size()); + REQUIRE(SortedCoords(gpu_spots) == SortedCoords(cpu_spots)); +} + TEST_CASE("AdaptiveSpotFinderGPU_Speed", "[AdaptiveSpotFinderGPU][.benchmark]") { if (get_gpu_count() == 0) { WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_Speed");