diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 2c2d169cf..eae6cc795 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -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(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(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(static_cast(v))); atomicAdd(&s_sum2[b], static_cast(static_cast(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(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(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(image)[q]; - const ushort4 b4 = reinterpret_cast(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(blockDim.x) * gridDim.x) { + const int32_t vmax = static_cast(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(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(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(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(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((OutputSize() + flag_threads - 1) / flag_threads); shared_plain = static_cast(nbins) * (4 * sizeof(unsigned long long) + sizeof(uint32_t)); shared_clip = static_cast(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<<>>(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<<>>( 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<<>>( @@ -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<<>>( - 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) { diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index e6ae453cc..2e0f077f5 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -85,6 +85,12 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinderGPU { CudaDevicePtr gpu_thr; CudaDevicePtr 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 gpu_word_max; + CudaDevicePtr gpu_word_range; + // 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 } (exact integers - see the kernel) diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index ca542f180..72abcfe2d 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -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 +__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 +__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(static_cast(v)); - sum2_a += sign * static_cast(static_cast(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(static_cast(v)); - sum2_b += sign * static_cast(static_cast(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(static_cast(v)); + sum2 += sign * static_cast(static_cast(v) * v); + count += static_cast(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(warp_window_sum(sum_a, sum_b)); const auto sum2 = static_cast(warp_window_sum(sum2_a, sum2_b)); const auto count = static_cast(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(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<<>> - (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<<>> - (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<<>> (image.getGPUBuffer(), gpu_out_0, gpu_candidates, gpu_out_1, spot_params, OutputSize()); diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index 97c2812cf..f363fb5ca 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -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) diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index d823ec084..b91e97205 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -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(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(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(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(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(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(i)) label[i] = next_label++; const int nlabel = static_cast(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(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(i)) continue; - const long long counts = value[j]; - const int col = static_cast(index[j] % width), line = static_cast(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(index[j] % width) * counts; - y += static_cast(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(*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(min(*nstrong, capacity)); + if (static_cast(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(index[i] % width), line = static_cast(index[i] / width); + SpotExtractorGPUSpot &s = scratch[l]; + atomicAdd(reinterpret_cast(&s.x), static_cast(col * counts)); + atomicAdd(reinterpret_cast(&s.y), static_cast(line * counts)); + atomicAdd(reinterpret_cast(&s.photons), static_cast(counts)); + atomicMax(reinterpret_cast(&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(*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(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(SPOT_SHAPE_FREE_PIXELS), static_cast(SPOT_MIN_FILL_PERCENT)); cuda_err(cudaGetLastError()); // Rides along with the spot count on the frame's one synchronisation, so knowing how many strong diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index 289582abf..8639a46f9 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -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 gpu_parent; // union-find parent CudaDevicePtr gpu_root; CudaDevicePtr gpu_label; // compact label, indexed by root - CudaDevicePtr gpu_count; // pixels per component + CudaDevicePtr gpu_box; // bounding box per component + CudaDevicePtr gpu_nlabel; // components found CudaDevicePtr gpu_spot; CudaDevicePtr gpu_spot_out; CudaDevicePtr gpu_nspot;