mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

Merge branch 'spotkern' into rc175: GPU adaptive spot finder, same spots in about half the kernel time

Browse Source 
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
This commit is contained in:
leonarski_f
2026-10-11 05:06:16 +02:00
co-authored by Claude Opus 5.5
parent 65fde24fb6 60ca872bc2
commit 56bad5fdba
6 changed files with 304 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
+105 -51
View File
@@ -81,6 +81,34 @@ __device__ __forceinline__ void count_value(uint32_t *hist, int2 *overflow, uint
}
}
// A pixel value as an unsigned key with the same order, the masked sentinel INT32_MIN mapping to 0 -
// so a zeroed word maximum is "every pixel masked" and atomicMax on unsigned keeps the order.
__device__ __forceinline__ uint32_t order_key(int32_t v) {
return static_cast<uint32_t>(v) ^ 0x80000000u;
}
// The maximum of a quad's key over the eight lanes that hold one 32-pixel word (quads 8w .. 8w + 7;
// a warp starts on a multiple of 32 quads, so the eight are lanes 8j .. 8j + 7), written by the first
// of them. Only the warp's last pass can have lanes past the end, and those are a suffix: a partner
// outside the active mask is simply not taken. A word whose eight quads are all here has no other
// writer, so it is stored outright; only the last word of the image, which may also hold the pixels
// past the last quad, is merged with an atomic (an atomic per word measured ~10 % of the kernel).
__device__ __forceinline__ void note_word_max(uint32_t *word_max, size_t q, uint32_t key) {
const unsigned active = __activemask();
const int lane = threadIdx.x % 32;
for (int offset = 1; offset < 8; offset *= 2) {
const uint32_t other = __shfl_xor_sync(active, key, offset);
if (active & (1u << (lane ^ offset)))
key = max(key, other);
}
if (lane % 8 == 0) {
if (((active >> lane) & 0xFFu) == 0xFFu)
word_max[q / 8] = key;
else
atomicMax(&word_max[q / 8], key);
}
}
// One ring reduction, staging per-ring sums in shared memory (fast path). Shared layout:
// [ sum | sum2 | sum_corr | sum2_corr | count(uint32) ] x nbins
// The corrected arrays exist only when accumulate_corrected is true (the plain first pass); on the
@@ -96,7 +124,7 @@ __global__ void reduce_rings_shared(
unsigned long long *__restrict__ sum, unsigned long long *__restrict__ sum2, uint32_t *__restrict__ count,
unsigned long long *__restrict__ sum_corr, unsigned long long *__restrict__ sum2_corr,
uint32_t *__restrict__ hist, int2 *__restrict__ overflow, uint32_t *__restrict__ overflow_n,
uint32_t overflow_cap,
uint32_t overflow_cap, uint32_t *__restrict__ word_max,
size_t npix, int nbins) {
// A sigma-clip pass reads the image only when the plain pass's overflow list did not fit;
@@ -150,6 +178,11 @@ __global__ void reduce_rings_shared(
const uint16_t bq[4] = {b4.x, b4.y, b4.z, b4.w};
const float cq[4] = {c4.x, c4.y, c4.z, c4.w};
// The plain pass also notes the largest value of every 32-pixel word, for flag_strong.
if (accumulate_corrected)
note_word_max(word_max, q, max(max(order_key(v4.x), order_key(v4.y)),
max(order_key(v4.z), order_key(v4.w))));
// A ring is several pixels wide, so consecutive pixels along a row usually fall in the same
// one. Carry a running total for the ring in registers and push it to shared memory only
// when the ring changes - one set of atomics for the run instead of one per pixel, which is
@@ -200,8 +233,10 @@ __global__ void reduce_rings_shared(
const float fv = static_cast<float>(v);
const bool valid = v != INT32_MIN && v != INT32_MAX && b < nbins
&& (clip_k <= 0.0f || clip_keep(fv, mean[b], sigma[b], clip_k));
if (accumulate_corrected)
if (accumulate_corrected) {
count_value(hist, overflow, overflow_n, overflow_cap, valid, b, v);
atomicMax(&word_max[idx / 32], order_key(v));
}
if (!valid) continue;
atomicAdd(&s_sum[b], static_cast<unsigned long long>(static_cast<long long>(v)));
atomicAdd(&s_sum2[b], static_cast<unsigned long long>(static_cast<long long>(v) * v));
@@ -329,57 +364,66 @@ __global__ void finalize_rings(const unsigned long long *__restrict__ sum,
}
}
// Flag strong pixels (value >= ring threshold, or saturated) into the packed bit buffer. Strong
// pixels are sparse, so a plain atomicOr per strong pixel is simpler than warp aggregation and the
// contention is negligible. Mirrors AdaptiveSpotFinderCPU Stage C exactly.
// The range of rings the pixels of each 32-pixel word fall in (x = lowest, y = highest; x > y when
// none of them is in a ring). A function of the geometry alone, filled once.
__global__ void word_ring_range(const uint16_t *__restrict__ pixel_to_bin, ushort2 *__restrict__ range,
size_t npix, int nbins) {
const size_t nwords = (npix + 31) / 32;
for (size_t w = blockIdx.x * blockDim.x + threadIdx.x; w < nwords; w += static_cast<size_t>(blockDim.x) * gridDim.x) {
int lo = nbins, hi = -1;
for (size_t idx = 32 * w; idx < min(32 * w + 32, npix); idx++) {
const int b = pixel_to_bin[idx];
if (b < nbins) {
lo = min(lo, b);
hi = max(hi, b);
}
}
range[w] = (hi < 0) ? make_ushort2(1, 0) : make_ushort2(lo, hi);
}
}
// Flag strong pixels (value >= ring threshold, or saturated) into the packed bit buffer. Mirrors
// AdaptiveSpotFinderCPU Stage C exactly.
// One thread per 32-pixel word, writing the whole word, so the buffer needs neither clearing nor
// atomics. The image is read only for the words that can hold a strong pixel: the plain pass noted
// each word's largest value, and a word whose largest value is below the lowest threshold of the
// rings it touches has none. Strong pixels are sparse, so that is a few thousand words of the image
// instead of all of it. The test reads the same float comparison as the per-pixel one below, and
// the conversion to float keeps the order, so it never drops a word the per-pixel test would flag.
__global__ void flag_strong(const int32_t *__restrict__ image,
const uint16_t *__restrict__ pixel_to_bin,
const uint32_t *__restrict__ word_max,
const ushort2 *__restrict__ word_range,
const float *__restrict__ thr,
uint32_t *__restrict__ strong,
size_t npix, int nbins) {
// Four pixels per thread, read as one 16-byte and one 8-byte transaction instead of four of
// each, exactly as the ring reduction above reads them - and flagged with a single atomicOr,
// because four consecutive pixels always fall in the same word of the bit buffer. The last
// npix % 4 pixels are done one at a time below, so nothing is read past the end.
const size_t stride = static_cast<size_t>(blockDim.x) * gridDim.x;
const size_t nquad = npix / 4;
for (size_t q = blockIdx.x * blockDim.x + threadIdx.x; q < nquad; q += stride) {
const int4 v4 = reinterpret_cast<const int4 *>(image)[q];
const ushort4 b4 = reinterpret_cast<const ushort4 *>(pixel_to_bin)[q];
const int32_t vq[4] = {v4.x, v4.y, v4.z, v4.w};
const uint16_t bq[4] = {b4.x, b4.y, b4.z, b4.w};
const size_t nwords = (npix + 31) / 32;
for (size_t w = blockIdx.x * blockDim.x + threadIdx.x; w < nwords; w += static_cast<size_t>(blockDim.x) * gridDim.x) {
const int32_t vmax = static_cast<int32_t>(word_max[w] ^ 0x80000000u);
bool check = (vmax == INT32_MAX);
if (!check && vmax != INT32_MIN) {
float lowest = INFINITY;
for (int b = word_range[w].x; b <= word_range[w].y; b++)
lowest = fminf(lowest, thr[b]);
check = static_cast<float>(vmax) >= lowest;
}
uint32_t bits = 0;
#pragma unroll
for (int k = 0; k < 4; k++) {
const int32_t v = vq[k];
if (v == INT32_MAX) {
bits |= 1u << k;
} else if (v != INT32_MIN) {
const int b = bq[k];
if (b < nbins && static_cast<float>(v) >= thr[b])
bits |= 1u << k;
if (check) {
for (size_t idx = 32 * w; idx < min(32 * w + 32, npix); idx++) {
const int32_t v = image[idx];
bool s = false;
if (v == INT32_MAX) {
s = true;
} else if (v != INT32_MIN) {
const uint16_t b = pixel_to_bin[idx];
if (b < nbins && static_cast<float>(v) >= thr[b])
s = true;
}
if (s)
bits |= 1u << (idx % 32);
}
}
if (bits) {
const size_t idx = 4 * q;
atomicOr(&strong[idx / 32], bits << (idx % 32));
}
}
for (size_t idx = 4 * nquad + blockIdx.x * blockDim.x + threadIdx.x; idx < npix; idx += stride) {
const int32_t v = image[idx];
bool s = false;
if (v == INT32_MAX) {
s = true;
} else if (v != INT32_MIN) {
const uint16_t b = pixel_to_bin[idx];
if (b < nbins && static_cast<float>(v) >= thr[b])
s = true;
}
if (s)
atomicOr(&strong[idx / 32], 1u << (idx % 32));
strong[w] = bits;
}
}
@@ -401,6 +445,8 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
gpu_sum2_corr(nbins),
gpu_thr(nbins),
gpu_ring(OutputSize()),
gpu_word_max(OutputSize()),
gpu_word_range(OutputSize()),
host_sum(nbins),
host_sum2(nbins),
host_count(nbins),
@@ -423,9 +469,10 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
cudaDeviceProp prop{};
cuda_err(cudaGetDeviceProperties(&prop, device));
reduce_blocks = 8 * prop.multiProcessorCount;
// flag_strong stays at four: it is bandwidth-shaped rather than atomic-bound, and eight
// measured no better (181 vs 175 us/launch).
flag_blocks = 4 * prop.multiProcessorCount;
// flag_strong: a thread per word. Its few loads per word depend on one another (the word's maximum,
// its ring range, their thresholds), so it wants every word in flight at once rather than a
// grid-stride walk.
flag_blocks = static_cast<int>((OutputSize() + flag_threads - 1) / flag_threads);
shared_plain = static_cast<size_t>(nbins) * (4 * sizeof(unsigned long long) + sizeof(uint32_t));
shared_clip = static_cast<size_t>(nbins) * (2 * sizeof(unsigned long long) + sizeof(uint32_t));
@@ -462,6 +509,9 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
gpu_corrections = SharedDeviceTable(mapping.Corrections().data(), npix,
mapping.Corrections().data(),
mapping.GetCorrectionsChecksum(), *stream);
word_ring_range<<<flag_blocks, flag_threads, 0, *stream>>>(gpu_pixel_to_bin->get(), gpu_word_range,
npix, nbins);
cuda_err(cudaGetLastError());
}
void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, float clip_k,
@@ -479,7 +529,7 @@ void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, flo
reduce_rings_shared<<<blocks, reduce_threads, shared, *stream>>>(
gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma,
clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr,
gpu_hist, gpu_overflow, gpu_overflow_n, overflow_cap, npix, nbins);
gpu_hist, gpu_overflow, gpu_overflow_n, overflow_cap, gpu_word_max, npix, nbins);
cuda_err(cudaGetLastError());
} else {
reduce_rings_global<<<reduce_blocks, reduce_threads, 0, *stream>>>(
@@ -555,6 +605,10 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,
if (use_shared) {
cuda_err(cudaMemsetAsync(gpu_hist, 0, sizeof(uint32_t) * nbins * HIST_VALUES, *stream));
cuda_err(cudaMemsetAsync(gpu_overflow_n, 0, sizeof(uint32_t), *stream));
cuda_err(cudaMemsetAsync(gpu_word_max, 0, OutputByteSize(), *stream));
} else {
// The global-atomics pass notes no word maxima: all ones makes flag_strong read every word.
cuda_err(cudaMemsetAsync(gpu_word_max, 0xFF, OutputByteSize(), *stream));
}
ReducePass(image, 0.0f, true); // plain pass also fills the corrected profile accumulators
@@ -602,9 +656,9 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,
// --- Stage C: the ring threshold, intersected with the classic local-box SNR test ---
cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream));
cuda_err(cudaMemsetAsync(gpu_ring, 0, OutputByteSize(), *stream));
flag_strong<<<flag_blocks, flag_threads, 0, *stream>>>(
image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_ring, npix, nbins);
image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_word_max, gpu_word_range, gpu_thr, gpu_ring,
npix, nbins);
cuda_err(cudaGetLastError());
if (settings.signal_to_noise_threshold <= 0.0f) {
image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
+6
View File
@@ -85,6 +85,12 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinderGPU {
CudaDevicePtr<float> gpu_thr;
CudaDevicePtr<uint32_t> gpu_ring;
// Per 32-pixel word of the image: the largest value of the frame (as an order-keeping unsigned
// key, noted by the plain pass) and the range of rings its pixels fall in (geometry, filled once).
// Together they let flag_strong skip the words that cannot hold a strong pixel.
CudaDevicePtr<uint32_t> gpu_word_max;
CudaDevicePtr<ushort2> gpu_word_range;
// Host mirrors of the small per-ring transfers.
std::vector<unsigned long long> host_sum; // clipped raw sum } input to the host threshold computation
std::vector<unsigned long long> host_sum2; // clipped raw sum^2 } (exact integers - see the kernel)
image_analysis/spot_finding/ImageSpotFinderGPU.cu
+77 -44
View File
@@ -86,10 +86,11 @@ __device__ __forceinline__ uint8_t pixel_result(const spot_parameters& params, c
}
// Inclusive prefix sum over the lanes of a warp (modular arithmetic).
__device__ __forceinline__ uint64_t warp_inclusive_scan(uint64_t v) {
template <typename T>
__device__ __forceinline__ T warp_inclusive_scan(T v) {
const int32_t lane = threadIdx.x & (WARP_SIZE - 1);
for (int offset = 1; offset < WARP_SIZE; offset *= 2) {
const uint64_t t = __shfl_up_sync(UINT32_MAX, v, offset);
const T t = __shfl_up_sync(UINT32_MAX, v, offset);
if (lane >= offset)
v += t;
}
@@ -99,13 +100,14 @@ __device__ __forceinline__ uint64_t warp_inclusive_scan(uint64_t v) {
// Sum over the 2 * NBX + 1 = 31 columns of each lane's window, the columns held two per lane: a from
// lane i is column i of the warp's 64, b from lane i is column 32 + i (lane 31 holds no b). Lane l's
// window is columns l + 1 .. l + 31: a of lanes l+1..31 plus b of lanes 0..l-1.
__device__ __forceinline__ uint64_t warp_window_sum(uint64_t a, uint64_t b) {
template <typename T>
__device__ __forceinline__ T warp_window_sum(T a, T b) {
static_assert(2 * ImageSpotFinder::NBX + 1 == WARP_SIZE - 1, "the two-columns-per-lane layout needs a 31-pixel window");
const int32_t lane = threadIdx.x & (WARP_SIZE - 1);
const uint64_t pa = warp_inclusive_scan(a);
const uint64_t pb = warp_inclusive_scan(b);
const uint64_t total_a = __shfl_sync(UINT32_MAX, pa, WARP_SIZE - 1);
uint64_t pb_before = __shfl_up_sync(UINT32_MAX, pb, 1);
const T pa = warp_inclusive_scan(a);
const T pb = warp_inclusive_scan(b);
const T total_a = __shfl_sync(UINT32_MAX, pa, WARP_SIZE - 1);
T pb_before = __shfl_up_sync(UINT32_MAX, pb, 1);
if (lane < 1)
pb_before = 0;
return total_a - pa + pb_before;
@@ -120,6 +122,7 @@ __device__ __forceinline__ uint64_t warp_window_sum(uint64_t a, uint64_t b) {
// looks within NBX pixels of them (analyze_candidates). A warp whose tile holds no pixel that
// any of them can reach leaves its tile unwritten.
// params: spot finding parameters
// rowsPerWave: rows per band
//
// One warp per 32 output columns and a band of rows (a "wave", gridDim-independent: the warps are
// numbered group-fastest). Each lane keeps the vertical sums (sum, sum2, count) of two input columns
@@ -129,19 +132,19 @@ __device__ __forceinline__ uint64_t warp_window_sum(uint64_t a, uint64_t b) {
// sums come from warp prefix scans, so there is no shared memory and no block-wide synchronisation.
// All sums are integer (modular 64-bit), so the result is exact and independent of the order.
__global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint32_t *out,
const uint32_t *needed, const spot_parameters params, int32_t nWaves)
const uint32_t *needed, const spot_parameters params, int32_t rowsPerWave)
{
constexpr int32_t NBX = ImageSpotFinder::NBX;
const int32_t lane = threadIdx.x & (WARP_SIZE - 1);
const int32_t warp = (blockIdx.x * blockDim.x + threadIdx.x) / WARP_SIZE;
const int32_t ngroups = (params.width + WARP_SIZE - 1) / WARP_SIZE;
const int32_t nWaves = (params.height + rowsPerWave - 1) / rowsPerWave;
const int32_t group = warp % ngroups;
const int32_t wave = warp / ngroups;
const int32_t rowsPerWave = (params.height + nWaves - 1) / nWaves;
const int32_t rmin = wave * rowsPerWave;
const int32_t rmax = min(rmin + rowsPerWave, params.height);
// Warp-uniform, so the shuffles and ballots below stay collective.
if (wave >= nWaves || rmin >= params.height)
if (wave >= nWaves)
return;
const int32_t out_col = group * WARP_SIZE + lane; // the column this lane writes
@@ -179,44 +182,51 @@ __global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint3
return sat ? INT32_MAX : in[npixel];
};
uint64_t sum_a = 0, sum2_a = 0, count_a = 0;
uint64_t sum_b = 0, sum2_b = 0, count_b = 0;
// Add (sign +1) or remove (-1) row r of this lane's two columns; sentinels count as nothing.
const auto slide = [&](int32_t row, uint64_t sign) {
if (use_a) {
const int32_t v = value_at(row, col_a);
if (v != INT32_MAX && v != INT32_MIN) {
sum_a += sign * static_cast<uint64_t>(static_cast<int64_t>(v));
sum2_a += sign * static_cast<uint64_t>(static_cast<int64_t>(v) * v);
count_a += sign;
}
}
if (use_b) {
const int32_t v = value_at(row, col_b);
if (v != INT32_MAX && v != INT32_MIN) {
sum_b += sign * static_cast<uint64_t>(static_cast<int64_t>(v));
sum2_b += sign * static_cast<uint64_t>(static_cast<int64_t>(v) * v);
count_b += sign;
}
uint64_t sum_a = 0, sum2_a = 0;
uint64_t sum_b = 0, sum2_b = 0;
// The count of a window is at most (2 * NBX + 1)^2, so 32 bits hold it - and halve its scans.
uint32_t count_a = 0, count_b = 0;
// Row r of column a or b, or the masked sentinel where the lane has no such column - which, like
// the sentinels in the image, counts as nothing.
const auto read_a = [&](int32_t row) { return use_a ? value_at(row, col_a) : INT32_MIN; };
const auto read_b = [&](int32_t row) { return use_b ? value_at(row, col_b) : INT32_MIN; };
// Add (sign +1) or remove (-1) one value of a column; sentinels count as nothing.
const auto slide = [](int32_t v, uint64_t sign, uint64_t &sum, uint64_t &sum2, uint32_t &count) {
if (v != INT32_MAX && v != INT32_MIN) {
sum += sign * static_cast<uint64_t>(static_cast<int64_t>(v));
sum2 += sign * static_cast<uint64_t>(static_cast<int64_t>(v) * v);
count += static_cast<uint32_t>(sign);
}
};
for (int32_t r = max(rmin - NBX, 0); r < min(rmin + NBX, params.height - 1) + 1; r++)
slide(r, 1);
for (int32_t r = max(rmin - NBX, 0); r < min(rmin + NBX, params.height - 1) + 1; r++) {
slide(read_a(r), 1, sum_a, sum2_a, count_a);
slide(read_b(r), 1, sum_b, sum2_b, count_b);
}
for (int32_t row = rmin; row < rmax; row++) {
// The rows leaving and entering the window are read before the window sums are formed, so the
// loads are in flight while the scans run rather than after them.
const bool leave = row - NBX >= 0;
const bool enter = row + NBX + 1 < params.height;
const int32_t leave_a = leave ? read_a(row - NBX) : INT32_MIN;
const int32_t leave_b = leave ? read_b(row - NBX) : INT32_MIN;
const int32_t enter_a = enter ? read_a(row + NBX + 1) : INT32_MIN;
const int32_t enter_b = enter ? read_b(row + NBX + 1) : INT32_MIN;
const int32_t value = (out_col < params.width) ? value_at(row, out_col) : INT32_MIN;
const auto sum = static_cast<int64_t>(warp_window_sum(sum_a, sum_b));
const auto sum2 = static_cast<int64_t>(warp_window_sum(sum2_a, sum2_b));
const auto count = static_cast<int64_t>(warp_window_sum(count_a, count_b));
uint8_t val = 0;
if (out_col < params.width)
val = pixel_result(params, value_at(row, out_col), sum, sum2, count);
val = pixel_result(params, value, sum, sum2, count);
write_result(params, out, row * params.width + out_col, val);
if (row - NBX >= 0)
slide(row - NBX, UINT64_MAX); // -1: the row leaving the window
if (row + NBX + 1 < params.height)
slide(row + NBX + 1, 1);
slide(leave_a, UINT64_MAX, sum_a, sum2_a, count_a); // -1: the row leaving the window
slide(leave_b, UINT64_MAX, sum_b, sum2_b, count_b);
slide(enter_a, 1, sum_a, sum2_a, count_a);
slide(enter_b, 1, sum_b, sum2_b, count_b);
}
}
@@ -267,11 +277,27 @@ __global__ void analyze_candidates(const int32_t *in, const uint32_t *prev_out,
int64_t sum = 0, sum2 = 0;
int32_t count = 0;
if (lane < 2 * NBX + 1 && col >= 0 && col < params.width) {
const int32_t rlo = max(row - NBX, 0);
const int32_t rhi = min(row + NBX, params.height - 1);
for (int32_t r = rlo; r <= rhi; r++) {
const int32_t val = value_at(r * params.width + col);
if (val != INT32_MAX && val != INT32_MIN) {
// Every read of the column first, all of them unconditional (rows past an edge are
// read at the edge and dropped below), so the 2 * (2 * NBX + 1) loads are issued
// together. Read through value_at, each pixel's load waited for its previous-pass
// bit and each row for the one before: a candidate cost a column of memory
// latencies, and a warp walks the candidates of its words one after another.
int32_t vals[2 * NBX + 1];
uint32_t prev[2 * NBX + 1];
#pragma unroll
for (int k = 0; k < 2 * NBX + 1; k++) {
const int32_t r = min(max(row + k - NBX, 0), params.height - 1);
const int32_t pix = r * params.width + col;
vals[k] = in[pix];
prev[k] = prev_out[pix / 32] >> (pix % 32);
}
#pragma unroll
for (int k = 0; k < 2 * NBX + 1; k++) {
const int32_t r = row + k - NBX;
// A previous-pass pixel reads as INT32_MAX, so like the sentinels it is not counted.
const int32_t val = vals[k];
if (r >= 0 && r < params.height && (prev[k] & 1u) == 0
&& val != INT32_MAX && val != INT32_MIN) {
sum += val;
sum2 += static_cast<int64_t>(val) * val;
count += 1;
@@ -347,19 +373,26 @@ void ImageSpotFinderGPU::RunDetect(const ImagePreprocessorBuffer &image, const S
if (windowSizeLimit > spot_params.height)
throw JFJochException(JFJochExceptionCategory::SpotFinderError, "window size limit exceeds number of height");
// With candidates the first pass is wanted only near them, and a warp whose tile is not near any
// returns at once - so its tiles are made short: the work then follows the candidates closely, and
// each warp's serial run of rows (the pass is bound by that chain, not by memory) is short. Dense,
// every tile is worked, and long tiles amortise the 2 * NBX rows each one spends filling its window.
const int32_t rowsPerWave = (gpu_candidates == nullptr) ? (spot_params.height + numberOfWaves - 1) / numberOfWaves
: sparseRowsPerWave;
const int32_t nWaves = (spot_params.height + rowsPerWave - 1) / rowsPerWave;
const int32_t ngroups = (spot_params.width + 31) / 32;
const int32_t warpsPerBlock = numberOfCudaThreads / 32;
const int32_t nBlocks = (ngroups * numberOfWaves + warpsPerBlock - 1) / warpsPerBlock;
const int32_t nBlocks = (ngroups * nWaves + warpsPerBlock - 1) / warpsPerBlock;
cuda_err(cudaMemsetAsync(gpu_out_0, 0, OutputSize() * sizeof(uint32_t), *stream));
cuda_err(cudaMemsetAsync(gpu_out_1, 0, OutputSize() * sizeof(uint32_t), *stream));
// The first pass has no previous one (the classic engine used to hand it the zeroed gpu_out_1).
analyze_pixel<<<nBlocks, numberOfCudaThreads, 0, *stream>>>
(image.getGPUBuffer(), nullptr, gpu_out_0, gpu_candidates, spot_params, numberOfWaves);
(image.getGPUBuffer(), nullptr, gpu_out_0, gpu_candidates, spot_params, rowsPerWave);
cuda_err(cudaGetLastError());
if (gpu_candidates == nullptr)
analyze_pixel<<<nBlocks, numberOfCudaThreads, 0, *stream>>>
(image.getGPUBuffer(), gpu_out_0, gpu_out_1, nullptr, spot_params, numberOfWaves);
(image.getGPUBuffer(), gpu_out_0, gpu_out_1, nullptr, spot_params, rowsPerWave);
else
analyze_candidates<<<candidateBlocks, 256, 0, *stream>>>
(image.getGPUBuffer(), gpu_out_0, gpu_candidates, gpu_out_1, spot_params, OutputSize());
image_analysis/spot_finding/ImageSpotFinderGPU.h
+1
View File
@@ -25,6 +25,7 @@ protected:
private:
const int numberOfCudaThreads = 128; // #threads per block of analyze_pixel (one warp per 32 columns)
const int numberOfWaves = 32; // #row bands of analyze_pixel
const int sparseRowsPerWave = 16; // rows per band of analyze_pixel when only near candidates
const int windowSizeLimit = 32; // limit on the window size (2nby+1, 2nbx+1): a warp holds 2 x 32 columns
int candidateBlocks = 0; // grid of analyze_candidates (256 threads per block)
image_analysis/spot_finding/SpotExtractorGPU.cu
+108 -73
View File
@@ -69,29 +69,41 @@ __global__ void scan_block_counts(const uint32_t *__restrict__ in, uint32_t *__r
if (t == nthreads - 1) *total = shared[t];
}
// One thread per block emits its range. Strong pixels are ~1e-4 of the image, so a block's range
// holds a handful of them and serial emission is both trivially ordered and fast; the parallelism
// comes from the block count.
// One warp per block emits its range, 32 words at a time: each lane takes a word, a warp prefix sum
// over the words' bit counts gives every lane where its pixels go, and the lane writes them in
// increasing order - so the list comes out in exactly the order a single thread walking the range
// would write it. Strong pixels are ~1e-4 of the image, so most passes find nothing at all; doing the
// walk on one thread left the other 31 idle and the reads of the range one after another.
__global__ void scatter_bits(const uint32_t *__restrict__ strong, const uint32_t *__restrict__ res_mask,
const uint32_t *__restrict__ block_offset, const int32_t *__restrict__ image,
uint32_t *__restrict__ out_index, int32_t *__restrict__ out_value,
size_t nwords, uint32_t capacity) {
if (threadIdx.x != 0) return;
const int lane = threadIdx.x % 32;
const size_t per_block = (nwords + gridDim.x - 1) / gridDim.x;
const size_t w0 = static_cast<size_t>(blockIdx.x) * per_block;
const size_t w1 = min(w0 + per_block, nwords);
uint32_t pos = block_offset[blockIdx.x];
for (size_t w = w0; w < w1; ++w) {
uint32_t word = strong[w] & ~res_mask[w];
for (size_t base = w0; base < w1; base += 32) {
const size_t w = base + lane;
uint32_t word = (w < w1) ? (strong[w] & ~res_mask[w]) : 0;
const uint32_t n = __popc(word);
uint32_t before = n; // inclusive prefix sum of n over the lanes
for (int offset = 1; offset < 32; offset *= 2) {
const uint32_t t = __shfl_up_sync(0xffffffff, before, offset);
if (lane >= offset) before += t;
}
const uint32_t total = __shfl_sync(0xffffffff, before, 31);
uint32_t p = pos + before - n;
while (word != 0) {
const uint32_t flat = static_cast<uint32_t>(w * 32 + (__ffs(word) - 1));
word &= word - 1;
if (pos < capacity) {
out_index[pos] = flat;
out_value[pos] = image[flat];
if (p < capacity) {
out_index[p] = flat;
out_value[p] = image[flat];
}
++pos;
++p;
}
pos += total;
}
}
@@ -173,17 +185,19 @@ __global__ void resolve_roots(uint32_t *__restrict__ parent, uint32_t *__restric
root[i] = find_root(parent, static_cast<uint32_t>(i));
}
// Everything after the labelling in ONE block, so a frame needs a single host synchronisation:
// hand out labels, count the members of each component, sum the surviving ones, filter by max-pix
// and compact - all of it order-preserving.
__global__ void finish_components(const uint32_t *__restrict__ index, const int32_t *__restrict__ value,
const uint32_t *__restrict__ root, uint32_t *__restrict__ label,
int32_t *__restrict__ count, SpotExtractorGPUSpot *__restrict__ scratch,
SpotExtractorGPUSpot *__restrict__ out, uint32_t *__restrict__ nout,
const uint32_t *__restrict__ nstrong, uint32_t capacity,
int width, int max_pix, int shape_free_pix, int min_fill_percent) {
// Everything after the labelling, in three launches and still one host synchronisation: hand out
// labels and clear each component's sums (one block), add every pixel to its component (all blocks),
// then filter by max-pix and shape and compact (one block) - each step order-preserving or made of
// integer sums, so the result does not depend on the order anything runs in.
__global__ void label_components(const uint32_t *__restrict__ root, uint32_t *__restrict__ label,
SpotExtractorGPUSpot *__restrict__ scratch, uint32_t *__restrict__ nlabel_out,
uint32_t *__restrict__ nout, const uint32_t *__restrict__ nstrong,
uint32_t capacity) {
const int n = static_cast<int>(min(*nstrong, capacity));
if (threadIdx.x == 0) *nout = 0;
if (threadIdx.x == 0) {
*nout = 0;
*nlabel_out = 0;
}
__syncthreads();
// The same give-up the host makes at StrongPixelLimit - except that here the count is known
// before a single pixel has been written anywhere, so the frame costs nothing to reject.
@@ -194,8 +208,8 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3
const int chunk = (n + nthreads - 1) / nthreads;
const int lo = min(t * chunk, n), hi = min(lo + chunk, n);
// 1) labels, by a prefix sum over the roots in ascending order - the order the host's second
// scan hands them out in, which is what makes the spot ORDER identical.
// Labels, by a prefix sum over the roots in ascending order - the order the host's second scan
// hands them out in, which is what makes the spot ORDER identical.
uint32_t nroot = 0;
for (int i = lo; i < hi; i++) nroot += (root[i] == static_cast<uint32_t>(i)) ? 1u : 0u;
shared[t] = nroot;
@@ -210,54 +224,65 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3
for (int i = lo; i < hi; i++)
if (root[i] == static_cast<uint32_t>(i)) label[i] = next_label++;
const int nlabel = static_cast<int>(shared[nthreads - 1]);
__syncthreads();
if (t == 0) *nlabel_out = nlabel;
// 2) member counts
for (int i = t; i < nlabel; i += nthreads) count[i] = 0;
__syncthreads();
for (int i = t; i < n; i += nthreads) atomicAdd(&count[label[root[i]]], 1);
__syncthreads();
// 3) sums, one thread per component, walking its members in ascending list order so the
// accumulation matches DiffractionSpot::AddPixel term for term. A component bigger than
// max-pix is thrown away below, so it is not summed - which is also what keeps a whole lit
// module or diffraction ring from turning into one thread walking tens of thousands of
// entries.
for (int i = t; i < n; i += nthreads) {
if (root[i] != static_cast<uint32_t>(i)) continue;
const uint32_t l = label[i];
const int want = count[l];
scratch[l].pixel_count = want;
if (want > max_pix) continue;
long long x = 0, y = 0;
long long photons = 0, max_photons = LLONG_MIN;
int min_col = INT_MAX, max_col = INT_MIN, min_line = INT_MAX, max_line = INT_MIN;
int found = 0;
for (int j = i; j < n && found < want; j++) {
if (root[j] != static_cast<uint32_t>(i)) continue;
const long long counts = value[j];
const int col = static_cast<int>(index[j] % width), line = static_cast<int>(index[j] / width);
min_col = min(min_col, col); max_col = max(max_col, col);
min_line = min(min_line, line); max_line = max(max_line, line);
// Integers, exactly as DiffractionSpot::AddPixel does them, so host and device agree by
// construction - no rounding mode to match and nothing for either compiler to contract.
x += static_cast<long long>(index[j] % width) * counts;
y += static_cast<long long>(index[j] / width) * counts;
photons += counts;
max_photons = max(max_photons, counts);
found++;
}
scratch[l].x = x;
scratch[l].y = y;
scratch[l].photons = photons;
scratch[l].max_photons = max_photons;
scratch[l].bbox_side = max(max_col - min_col, max_line - min_line) + 1;
for (int l = t; l < nlabel; l += nthreads) {
SpotExtractorGPUSpot &s = scratch[l];
s.x = 0;
s.y = 0;
s.photons = 0;
s.max_photons = LLONG_MIN;
s.pixel_count = 0;
s.bbox_side = 0;
}
__syncthreads();
}
// 4) size and shape filter, compacted by another prefix sum so the surviving spots keep their
// order. The test is SpotShapeAccepted written out - the constants come in as arguments rather
// than being included here, so there is exactly one definition of them.
// The bounding box of each component (x, y = lowest and highest column, z, w = lowest and highest
// line), kept aside from the sums while they are being accumulated.
__global__ void init_boxes(int4 *__restrict__ box, const uint32_t *__restrict__ nlabel) {
const int n = static_cast<int>(*nlabel);
for (int l = blockIdx.x * blockDim.x + threadIdx.x; l < n; l += blockDim.x * gridDim.x)
box[l] = make_int4(INT_MAX, INT_MIN, INT_MAX, INT_MIN);
}
// Every strong pixel into its component: the member count, the sums exactly as DiffractionSpot::AddPixel
// forms them (integers, so the order the atomics arrive in does not matter and host and device agree by
// construction) and the bounding box.
__global__ void sum_components(const uint32_t *__restrict__ index, const int32_t *__restrict__ value,
const uint32_t *__restrict__ root, const uint32_t *__restrict__ label,
SpotExtractorGPUSpot *__restrict__ scratch, int4 *__restrict__ box,
const uint32_t *__restrict__ nstrong, uint32_t capacity, int width) {
const int n = static_cast<int>(min(*nstrong, capacity));
if (static_cast<uint32_t>(n) >= capacity) return;
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) {
const uint32_t l = label[root[i]];
const long long counts = value[i];
const int col = static_cast<int>(index[i] % width), line = static_cast<int>(index[i] / width);
SpotExtractorGPUSpot &s = scratch[l];
atomicAdd(reinterpret_cast<unsigned long long *>(&s.x), static_cast<unsigned long long>(col * counts));
atomicAdd(reinterpret_cast<unsigned long long *>(&s.y), static_cast<unsigned long long>(line * counts));
atomicAdd(reinterpret_cast<unsigned long long *>(&s.photons), static_cast<unsigned long long>(counts));
atomicMax(reinterpret_cast<long long *>(&s.max_photons), counts);
atomicAdd(&s.pixel_count, 1);
atomicMin(&box[l].x, col);
atomicMax(&box[l].y, col);
atomicMin(&box[l].z, line);
atomicMax(&box[l].w, line);
}
}
// Size and shape filter, compacted by a prefix sum so the surviving spots keep their order. The test
// is SpotShapeAccepted written out - the constants come in as arguments rather than being included
// here, so there is exactly one definition of them.
__global__ void filter_components(SpotExtractorGPUSpot *__restrict__ scratch, const int4 *__restrict__ box,
SpotExtractorGPUSpot *__restrict__ out, uint32_t *__restrict__ nout,
const uint32_t *__restrict__ nlabel_in,
int max_pix, int shape_free_pix, int min_fill_percent) {
const int nlabel = static_cast<int>(*nlabel_in);
if (nlabel == 0) return;
__shared__ uint32_t shared[FINISH_THREADS];
const int t = threadIdx.x, nthreads = blockDim.x;
auto keep = [&](const SpotExtractorGPUSpot &s) {
if (s.pixel_count > max_pix) return false;
if (s.pixel_count <= shape_free_pix) return true;
@@ -267,7 +292,10 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3
const int label_chunk = (nlabel + nthreads - 1) / nthreads;
const int label_lo = min(t * label_chunk, nlabel), label_hi = min(label_lo + label_chunk, nlabel);
uint32_t nkeep = 0;
for (int i = label_lo; i < label_hi; i++) nkeep += keep(scratch[i]) ? 1u : 0u;
for (int i = label_lo; i < label_hi; i++) {
scratch[i].bbox_side = max(box[i].y - box[i].x, box[i].w - box[i].z) + 1;
nkeep += keep(scratch[i]) ? 1u : 0u;
}
shared[t] = nkeep;
__syncthreads();
for (int d = 1; d < nthreads; d <<= 1) {
@@ -296,7 +324,8 @@ SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::sha
gpu_parent(max_strong),
gpu_root(max_strong),
gpu_label(max_strong),
gpu_count(max_strong),
gpu_box(max_strong),
gpu_nlabel(1),
gpu_spot(max_strong),
gpu_spot_out(max_strong),
gpu_nspot(1),
@@ -353,10 +382,16 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im
cuda_err(cudaGetLastError());
resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong);
cuda_err(cudaGetLastError());
finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label,
gpu_count, gpu_spot, gpu_spot_out, gpu_nspot,
gpu_nstrong, max_strong, width, max_pix,
static_cast<int>(SPOT_SHAPE_FREE_PIXELS),
label_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_root, gpu_label, gpu_spot, gpu_nlabel, gpu_nspot,
gpu_nstrong, max_strong);
cuda_err(cudaGetLastError());
init_boxes<<<512, THREADS, 0, *stream>>>(gpu_box, gpu_nlabel);
cuda_err(cudaGetLastError());
sum_components<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_spot, gpu_box,
gpu_nstrong, max_strong, width);
cuda_err(cudaGetLastError());
filter_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_spot, gpu_box, gpu_spot_out, gpu_nspot, gpu_nlabel,
max_pix, static_cast<int>(SPOT_SHAPE_FREE_PIXELS),
static_cast<int>(SPOT_MIN_FILL_PERCENT));
cuda_err(cudaGetLastError());
// Rides along with the spot count on the frame's one synchronisation, so knowing how many strong
image_analysis/spot_finding/SpotExtractorGPU.h
+7 -6
View File
@@ -21,9 +21,9 @@
// * both make a component's root its lowest list index, so both find the same roots;
// * labels are handed out by a prefix sum over the roots in ascending order, which is the order the
// host's second scan hands them out in, so the SPOT ORDER is identical;
// * the centroid sums are accumulated per component in ascending list order, in integers, term for
// term as DiffractionSpot::AddPixel does them, so there is no rounding for the two compilers to
// disagree about.
// * the centroid sums are accumulated per component in integers, term for term as
// DiffractionSpot::AddPixel does them, so there is no rounding for the two compilers to disagree
// about, and the order the device adds them in changes nothing.
// tests/SpotExtractorGPUParityTest.cpp holds the two to each other on realistic, occupancy-swept and
// pathological frames, and checks that repeating a frame gives byte-identical output.
@@ -54,8 +54,8 @@ class SpotExtractorGPU {
const size_t nwords;
// Strong pixels this engine's buffers hold, and above which the extraction gives up on the frame -
// StrongPixelLimit, so it follows the detector rather than standing at a constant. 104 bytes of
// device memory apiece, 29 MB on an 18-megapixel detector.
// StrongPixelLimit, so it follows the detector rather than standing at a constant. 116 bytes of
// device memory apiece, 32 MB on an 18-megapixel detector.
const uint32_t max_strong;
// Spots copied back together with their count in one transfer. A frame with more than this many
// surviving spots - far past anything indexable - simply takes a second copy.
@@ -72,7 +72,8 @@ class SpotExtractorGPU {
CudaDevicePtr<uint32_t> gpu_parent; // union-find parent
CudaDevicePtr<uint32_t> gpu_root;
CudaDevicePtr<uint32_t> gpu_label; // compact label, indexed by root
CudaDevicePtr<int32_t> gpu_count; // pixels per component
CudaDevicePtr<int4> gpu_box; // bounding box per component
CudaDevicePtr<uint32_t> gpu_nlabel; // components found
CudaDevicePtr<SpotExtractorGPUSpot> gpu_spot;
CudaDevicePtr<SpotExtractorGPUSpot> gpu_spot_out;
CudaDevicePtr<uint32_t> gpu_nspot;