mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases 67 Wiki Activity

Adaptive spot finder: sum the rings across blocks in double

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-30 11:05:54 +02:00
co-authored by Claude Opus 5
parent 8dc7254624
commit 04450eb618
2 changed files with 37 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
+31 -20
View File
@@ -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<double>(s_sum[i]));
atomicAdd(&sum2[i], static_cast<double>(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<double>(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<float>(m);
sigma[b] = static_cast<float>(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<double>(host_sum[b]) / host_count[b];
const double var = std::max(0.0, static_cast<double>(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<float>(m),
static_cast<float>(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));
image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
+6 -4
View File
@@ -51,8 +51,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
CudaDevicePtr<float> gpu_corrections;
// Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold.
CudaDevicePtr<float> gpu_sum;
CudaDevicePtr<float> 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<double> gpu_sum;
CudaDevicePtr<double> gpu_sum2;
CudaDevicePtr<uint32_t> gpu_count;
CudaDevicePtr<float> gpu_mean; // per-ring raw mean (clip predicate)
CudaDevicePtr<float> gpu_sigma; // per-ring raw sigma (clip predicate)
@@ -66,8 +68,8 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
CudaDevicePtr<uint32_t> gpu_strong;
// Host mirrors of the small per-ring transfers.
std::vector<float> host_sum; // clipped raw sum } input to the host threshold computation
std::vector<float> host_sum2; // clipped raw sum^2 }
std::vector<double> host_sum; // clipped raw sum } input to the host threshold computation
std::vector<double> host_sum2; // clipped raw sum^2 }
std::vector<uint32_t> host_count; // clipped raw count }
std::vector<float> host_thr; // per-ring threshold (empty -> frame had no valid pixels)
std::vector<float> prof_sum; // plain corrected sum } azimuthal-integration profile