From 9f49e5abb7e62ab0a6099465e31421cb814f9df4 Mon Sep 17 00:00:00 2001 From: jungfrau Date: Tue, 18 Aug 2026 01:52:27 -0400 Subject: [PATCH] Keep 16-bit images 16-bit through the GPU pipeline A detector reading out 16 bits had its frame widened to int32 the moment it was decoded, and every per-pixel pass over that frame then moved four bytes a pixel to carry two. Those passes - the ring statistics three times over, the strong-pixel search, spot extraction, the azimuthal and ROI integrators, Bragg integration - are the bulk of the image loop's device traffic, and 16 bits is the mode a fast acquisition runs in, which is exactly where throughput matters. The preprocessed image now keeps the width of its source. Two codes at the top of the 16-bit range carry the two special states, and they cannot collide with a real value: 0xFFFF masked, or the source's own bad-pixel marker. saturation a pixel at or above the saturation limit. 0xFFFE where the limit code leaves room - a 16-bit EIGER declares a count-rate limit of a few thousand, so there is room to spare - and 0xFFFF where the limit is the whole range, in which case the "is error" test has already claimed 0xFFFF, nothing can be saturated, and 0xFFFE stays a real value. Either way a real value is strictly below the saturation limit and so below both codes. Nothing is clipped and nothing is lost, and which code is in force is carried with the image rather than assumed. No engine learns a second convention. PixelView widens on load, so a masked pixel still reads as INT32_MIN and a saturated one as INT32_MAX, and every existing `v != INT32_MIN && v != INT32_MAX` test keeps its meaning. One code path, not two instantiations that can drift apart; the branch is on a pointer that is the same for every thread of every block, on kernels whose time is the loads it selects between. The vector loads are kept - four pixels still arrive in one transaction, 16 bytes wide or 8, whichever the image is. The wide path is unchanged, and is still taken for anything that is not a 16-bit source, and for any caller that wants the preprocessed image copied back to the host - that mirror is int32 and the CPU engines know only that convention. Measured on the one 16-bit dataset in the rotation test set, which is also the smallest detector in it (2.5M pixels, where per-pixel work is a small part of the loop): image loop 1.025 s -> 1.005 s at one GPU, whole run 5.64 s -> 5.52 s. The gain scales with the frame, so a 16M-pixel detector - where six full-frame passes are 86 % of the loop's GPU time - has much more to gain, and nothing here can measure that: every other dataset in the test set is stored 32-bit. Correctness on that dataset is exact where it can be: indexing rate, first-pass validation score and the integrated partial count are identical to the wide path, and its whole battery row - reflections, observations, space group, R_meas, CC1/2, ISa, mosaicity - is unchanged. Battery 6m17s -> 6m16s, 21/24 space groups, no failures. Also logs, once per run, the width the images are stored in, since it decides how much of the frame moves through every pass. Co-Authored-By: Claude Opus 5 --- docs/CHANGELOG.md | 1 + image_analysis/azint/AzIntEngineGPU.cu | 9 +- .../BraggIntegrationEngineGPU.cu | 11 +- .../ImagePreprocessorBuffer.h | 9 ++ .../ImagePreprocessorBufferGPU.cu | 13 +- .../ImagePreprocessorBufferGPU.h | 12 ++ .../ImagePreprocessorGPU.cu | 114 ++++++++++++++---- .../ImagePreprocessorGPU.h | 2 + .../image_preprocessing/PreprocessedPixel.h | 88 ++++++++++++++ image_analysis/roi/ROIIntegrationGPU.cu | 5 +- .../spot_finding/AdaptiveSpotFinderGPU.cu | 23 ++-- .../spot_finding/ImageSpotFinderGPU.cu | 9 +- .../spot_finding/SpotExtractorGPU.cu | 4 +- .../spot_finding/SpotExtractorGPU.h | 3 +- rugnux/Rugnux.cpp | 12 ++ 15 files changed, 260 insertions(+), 55 deletions(-) create mode 100644 image_analysis/image_preprocessing/PreprocessedPixel.h diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 0134e2bb..7d83164f 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### Unreleased +* **rugnux: 16-bit images keep their width all the way through the GPU pipeline** instead of being widened to 32-bit as soon as they are decoded. Every per-pixel pass over a frame - the ring statistics, the strong-pixel search, spot extraction, azimuthal and ROI integration, Bragg integration - then moves half as many bytes. This is the mode a fast acquisition reads out in, and the gain grows with the detector: on the 2.5M-pixel dataset in the test set the image loop drops 2%, and it is worth proportionally more on a large detector, where per-pixel work is most of the loop. Merged results are unchanged. * **rugnux: the first pass of the rotation two-pass no longer merges, reports or writes anything.** It exists to measure the detector geometry, the goniometer rotation scale, the mosaicity and the space group, and the second pass makes the merged result again at the refined geometry - so the first pass now stops once it has those. **The `_01_*` files are no longer written.** A battery of 24 rotation crystals drops from 6m52s to 6m17s and no merged result changes. * rugnux: every run ends with a `Time by phase` table - wall time and mean cores busy for each phase of the whole run, both passes together - so a slow dataset can be diagnosed from its own log without a profiler. * rugnux: **fixed** the `Per-image cost` line, which divided the per-stage means by the thread count while the image loop runs a few workers per GPU. On a 2-GPU machine at `-N 48` every stage was reported 6x too small. It now divides by, and names, the loop's own worker count. diff --git a/image_analysis/azint/AzIntEngineGPU.cu b/image_analysis/azint/AzIntEngineGPU.cu index dce1d43c..645465aa 100644 --- a/image_analysis/azint/AzIntEngineGPU.cu +++ b/image_analysis/azint/AzIntEngineGPU.cu @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "AzIntEngineGPU.h" +#include "../image_preprocessing/PreprocessedPixel.h" inline void cuda_err(cudaError_t val) { if (val != cudaSuccess) @@ -12,7 +13,7 @@ __global__ void gpu_azim_shared( const uint16_t *__restrict__ pixel_to_bin, const float *__restrict__ corrections, - const int32_t *__restrict__ input_buffer, + PixelView input_buffer, float *__restrict__ azint_sum, float *__restrict__ azint_sum2, uint32_t *__restrict__ azint_count, @@ -64,7 +65,7 @@ __global__ void gpu_azim( const uint16_t *__restrict__ pixel_to_bin, const float *__restrict__ corrections, - const int32_t *__restrict__ input_buffer, + PixelView input_buffer, float *__restrict__ azint_sum, float *__restrict__ azint_sum2, uint32_t *__restrict__ azint_count, @@ -125,12 +126,12 @@ void AzIntEngineGPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrat if (shared_needed < shared_size) { gpu_azim_shared<<>>( - gpu_pixel_to_bin->get(),gpu_azint_correction->get(),image.getGPUBuffer(), gpu_sum, gpu_sum2, + gpu_pixel_to_bin->get(),gpu_azint_correction->get(),ViewOf(image), gpu_sum, gpu_sum2, gpu_count, npixel, azint_bins ); } else { gpu_azim<<>>( - gpu_pixel_to_bin->get(),gpu_azint_correction->get(),image.getGPUBuffer(), gpu_sum, gpu_sum2, + gpu_pixel_to_bin->get(),gpu_azint_correction->get(),ViewOf(image), gpu_sum, gpu_sum2, gpu_count, npixel, azint_bins ); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 912e9aa2..6170d99c 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "BraggIntegrationEngineGPU.h" +#include "../image_preprocessing/PreprocessedPixel.h" using namespace bragg_engine; @@ -95,7 +96,7 @@ __device__ __forceinline__ T warp_sum(T v) { // --- Pass A box-sum: rough I / background / centroid / strong flag, one block per reflection. --- __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, - const int32_t *img, const uint8_t *mask, const uint32_t *owner, + PixelView img, const uint8_t *mask, const uint32_t *owner, BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, @@ -372,7 +373,7 @@ __global__ void reset(float *shell_grid, float *global_grid, float *mom, int *sh // the spot's own radial/tangential frame, because a detector-frame stack is azimuthally averaged // and cannot tell a radially smeared spot from a tangentially wide one (see the CPU engine). // One block per reflection. --- -__global__ void learn_profile(const int32_t *img, const uint32_t *owner, +__global__ void learn_profile(PixelView img, const uint32_t *owner, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, @@ -496,7 +497,7 @@ __global__ void radial_correct(const float *rad_sum, const int *rad_cnt, int n_r // reweighting is the Kabsch/Otwinowski iteration: Kabsch, Acta Cryst D66, 133-144 (2010); // Otwinowski & Minor, Methods Enzymol 276, 307-326 (1997). // One block per reflection; the (possibly elongated) profile is built in shared memory. --- -__global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x, const float *px_y, +__global__ void fit(PixelView img, const uint32_t *owner, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_seed, const float *sigma_seed, const float *bkg_a, const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, @@ -785,8 +786,8 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (image.size() != npixel || npredicted == 0) return Finalize(predicted, npredicted, results, image_number); - const int32_t *img = image.getGPUBuffer(); - if (img == nullptr) + const PixelView img = ViewOf(image); + if (img.wide == nullptr && img.narrow == nullptr) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "BraggIntegrationEngineGPU: image buffer is not on the GPU"); diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h index c5730515..8b294c50 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h @@ -40,6 +40,15 @@ public: virtual int32_t *getGPUBuffer(); virtual const int32_t *getGPUBuffer() const; + // Whether the device image kept the 16-bit width of its source instead of being widened to int32 + // (see PreprocessedPixel.h). Always false on the host path, which is int32 whatever came in. + virtual bool IsNarrow() const { return false; } + virtual const uint16_t *getGPUBufferNarrow() const { return nullptr; } + // Which 16-bit code means saturated in this image (see PreprocessedPixel.h). + virtual uint16_t NarrowSatCode() const { return 0; } + // Set by the preprocessor for each image it writes. A no-op on a host buffer, which is int32. + virtual void SetNarrow(bool, uint16_t) {} + // Values of the pixels with the given flat indices. Not just operator[] in a loop: on the GPU the // preprocessed image can live in device memory only, and then this gathers them there. virtual void Gather(const std::vector &npixel, std::vector &values) const; diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu index 54265f5c..eba71dfc 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu @@ -2,11 +2,14 @@ // SPDX-License-Identifier: GPL-3.0-only #include "ImagePreprocessorBufferGPU.h" +#include "PreprocessedPixel.h" -__global__ void gather_kernel(const int32_t *__restrict__ image, +__global__ void gather_kernel(PixelView image, const uint32_t *__restrict__ npixel, int32_t *__restrict__ values, int count) { + // Hands back int32 in the pipeline's convention whatever width the image is stored in, so the + // spot finder that asks for these values sees no difference. for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < count; i += blockDim.x * gridDim.x) values[i] = image[npixel[i]]; } @@ -28,6 +31,12 @@ const int32_t *ImagePreprocessorBufferGPU::getGPUBuffer() const { return gpu_image; } +const uint16_t *ImagePreprocessorBufferGPU::getGPUBufferNarrow() const { + // The same allocation, read as 16-bit. It is sized for the wide form, so the narrow image uses + // its first half and there is nothing to allocate when a run turns out to be 16-bit. + return reinterpret_cast(static_cast(gpu_image)); +} + void ImagePreprocessorBufferGPU::Gather(const std::vector &npixel, std::vector &values) const { values.resize(npixel.size()); if (npixel.empty()) @@ -37,7 +46,7 @@ void ImagePreprocessorBufferGPU::Gather(const std::vector &npixel, std cudaMemcpyAsync(gpu_gather_index.get(), npixel.data(), count * sizeof(uint32_t), cudaMemcpyHostToDevice, gather_stream); gather_kernel<<<(count + 255) / 256, 256, 0, gather_stream>>>( - gpu_image.get(), gpu_gather_index.get(), gpu_gather_value.get(), count); + ViewOf(*this), gpu_gather_index.get(), gpu_gather_value.get(), count); cudaMemcpyAsync(values.data(), gpu_gather_value.get(), count * sizeof(int32_t), cudaMemcpyDeviceToHost, gather_stream); cudaStreamSynchronize(gather_stream); diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h index dd223be1..149c538e 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h @@ -7,10 +7,15 @@ #include "ImagePreprocessorBuffer.h" #include "../indexing/CUDAMemHelpers.h" +#include "PreprocessedPixel.h" class ImagePreprocessorBufferGPU : public ImagePreprocessorBuffer { CudaDevicePtr gpu_image; CudaRegisteredVector buffer_reg; + // Set per image by the preprocessor: true once it has written this frame in its source's 16-bit + // width rather than widening it, with the code that frame uses for a saturated pixel. + bool narrow = false; + uint16_t narrow_sat_code = 0; // Staging for Gather(). Its only caller is ImageSpotFinder::ExtractSpots, which gives up on a frame // with UINT16_MAX or more strong pixels (the connected-component search rejects it anyway), so that @@ -29,5 +34,12 @@ public: int32_t *getGPUBuffer() override; const int32_t *getGPUBuffer() const override; + + // The device image is one allocation either way - the narrow form is the same memory read two + // bytes at a time, so a run that switches width allocates nothing and frees nothing. + bool IsNarrow() const override { return narrow; } + const uint16_t *getGPUBufferNarrow() const override; + uint16_t NarrowSatCode() const override { return narrow_sat_code; } + void SetNarrow(bool v, uint16_t sat_code) override { narrow = v; narrow_sat_code = sat_code; } void Gather(const std::vector &npixel, std::vector &values) const override; }; diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 37c772af..74f6ef10 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -4,6 +4,7 @@ #include #include "ImagePreprocessorGPU.h" +#include "PreprocessedPixel.h" #include "../../common/JFJochException.h" namespace { @@ -13,15 +14,46 @@ namespace { } } -template + +// Store a finished pixel in the width this image uses. The narrow form is only ever chosen when the +// source cannot produce a saturated pixel (see ChooseNarrow), so INT32_MIN - masked or bad - is the +// only special value that can arrive here; it becomes the one reserved 16-bit code. INT32_MAX would +// truncate to that same code, so even a broken invariant degrades to "bad" rather than to a wrong +// intensity. +template __device__ __forceinline__ OutT StorePixel(int32_t v, uint16_t sat_code); +template<> __device__ __forceinline__ int32_t StorePixel(int32_t v, uint16_t) { return v; } +template<> __device__ __forceinline__ uint16_t StorePixel(int32_t v, uint16_t sat_code) { + if (v == INT32_MIN) return preprocessed_pixel::NARROW_BAD; + if (v == INT32_MAX) return sat_code; + return static_cast(v); +} + + +// Eight finished pixels in one aligned store, in whichever width the image is stored in. +template __device__ __forceinline__ void StorePixels8(OutT *dst, const int32_t o[8], uint16_t sat_code); +template<> __device__ __forceinline__ void StorePixels8(int32_t *dst, const int32_t o[8], uint16_t) { + int4 *d = reinterpret_cast(dst); + d[0] = make_int4(o[0], o[1], o[2], o[3]); + d[1] = make_int4(o[4], o[5], o[6], o[7]); +} +template<> __device__ __forceinline__ void StorePixels8(uint16_t *dst, const int32_t o[8], uint16_t sc) { + ushort4 *d = reinterpret_cast(dst); + d[0] = make_ushort4(StorePixel(o[0], sc), StorePixel(o[1], sc), + StorePixel(o[2], sc), StorePixel(o[3], sc)); + d[1] = make_ushort4(StorePixel(o[4], sc), StorePixel(o[5], sc), + StorePixel(o[6], sc), StorePixel(o[7], sc)); +} + +template __global__ void preprocess_kernel( const T *__restrict__ input, const uint8_t *__restrict__ mask, - int32_t *__restrict__ output, + OutT *__restrict__ output, ImageStatistics *__restrict__ stats, T saturation_limit, T err_value, - int npixels) { + int npixels, + uint16_t narrow_sat_code) { // Shared block accumulators __shared__ unsigned long long s_masked; __shared__ unsigned long long s_saturated; @@ -60,8 +92,9 @@ __global__ void preprocess_kernel( bool valid = !(is_masked || is_sat || is_err); // Output - output[i] = - is_masked ? INT32_MIN : is_err ? INT32_MIN : is_sat ? INT32_MAX : (int32_t) v; + output[i] = StorePixel( + is_masked ? INT32_MIN : is_err ? INT32_MIN : is_sat ? INT32_MAX : (int32_t) v, + narrow_sat_code); // Counters local_masked += is_masked; @@ -166,15 +199,16 @@ __device__ __forceinline__ uint64_t transpose8_fused(uint64_t x) { // // The last CUDA block (blockIdx.x == nblocks) finishes the handful of elements bitshuffle stores // verbatim; they are already on the device inside the uploaded chunk. -template +template __global__ __launch_bounds__(256) void untranspose_preprocess_kernel( const uint8_t *__restrict__ shuffled, const BSLZ4BlockDesc *__restrict__ desc, const uint8_t *__restrict__ mask, - int32_t *__restrict__ out, + OutT *__restrict__ out, ImageStatistics *__restrict__ stats, T sat_value, T err_value, int nblocks, - const uint8_t *__restrict__ tail_src, uint32_t tail_elems, uint32_t tail_elem0) { + const uint8_t *__restrict__ tail_src, uint32_t tail_elems, uint32_t tail_elem0, + uint16_t narrow_sat_code) { PreprocessAccum l; // The bytes are assembled in the unsigned counterpart of T - shifting a byte into the top of a @@ -189,7 +223,9 @@ __global__ __launch_bounds__(256) void untranspose_preprocess_kernel( for (int p = 0; p < ES; p++) uv |= (U)((U)tail_src[threadIdx.x * ES + p] << (8 * p)); const T v = (T) uv; - out[tail_elem0 + threadIdx.x] = l.Apply(v, mask[tail_elem0 + threadIdx.x] != 0, sat_value, err_value); + out[tail_elem0 + threadIdx.x] = StorePixel( + l.Apply(v, mask[tail_elem0 + threadIdx.x] != 0, sat_value, err_value), + narrow_sat_code); } FlushStats(l, stats); return; @@ -219,10 +255,9 @@ __global__ __launch_bounds__(256) void untranspose_preprocess_kernel( for (int p = 0; p < ES; p++) uv |= (U)((U)((x[p] >> (8 * k)) & 0xff) << (8 * p)); o[k] = l.Apply((T) uv, mask[elem0 + i * 8 + k] != 0, sat_value, err_value); } - // elem0 is a multiple of 8 (bitshuffle blocks are), so this is 32-byte aligned. - int4 *dst = reinterpret_cast(out + elem0 + i * 8); - dst[0] = make_int4(o[0], o[1], o[2], o[3]); - dst[1] = make_int4(o[4], o[5], o[6], o[7]); + // elem0 is a multiple of 8 (bitshuffle blocks are), so eight pixels are one aligned store of + // 32 bytes wide or 16, whichever this image is. + StorePixels8(out + elem0 + i * 8, o, narrow_sat_code); } FlushStats(l, stats); } @@ -266,6 +301,17 @@ void ImagePreprocessorGPU::PinInputBuffer(std::vector &buffer, size_t s input_reg.rebind(buffer); } + +// Whether this image can keep its 16-bit width instead of being widened to int32. +// +// Whenever the source is 16-bit unsigned: the two reserved codes carry masked/bad and saturated, +// and a real value is always below both (see PreprocessedPixel.h). The one exception is a caller +// that wants the preprocessed image copied back to the host - that mirror is int32, and the CPU +// engines reading it know only that convention. +bool ImagePreprocessorGPU::ChooseNarrow(bool source_is_uint16) const { + return source_is_uint16 && !copy_image_to_host; +} + ImageStatistics ImagePreprocessorGPU::Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *image_ptr, CompressedImageMode image_mode) { switch (image_mode) { @@ -334,20 +380,31 @@ ImageStatistics ImagePreprocessorGPU::UntransposeAndAnalyze(ImagePreprocessorBuf cpu_stats[0] = ImageStatistics{.max_value = INT64_MIN, .min_value = INT64_MAX}; cuda_err(cudaMemcpyAsync(gpu_stats, cpu_stats.data(), sizeof(ImageStatistics), cudaMemcpyHostToDevice, *stream)); + const bool narrow = ChooseNarrow(std::is_same::value); + const uint16_t sat_code = preprocessed_pixel::NarrowSaturatedCode( + std::min(saturation_limit, static_cast(sat_value))); + processed_image.SetNarrow(narrow, sat_code); + // One CUDA block per bitshuffle block, plus one for the verbatim tail when there is one. const int nb = shuffled.nblocks + (shuffled.tail_elems > 0 ? 1 : 0); - untranspose_preprocess_kernel <<< nb, 256, 0, *stream >>>( + if (narrow) + untranspose_preprocess_kernel <<< nb, 256, 0, *stream >>>( shuffled.shuffled, shuffled.desc, gpu_mask->get(), - processed_image.getGPUBuffer(), + reinterpret_cast(processed_image.getGPUBuffer()), gpu_stats, sat_value, err_value, shuffled.nblocks, shuffled.tail_src, shuffled.tail_elems, - shuffled.tail_elem0); + shuffled.tail_elem0, sat_code); + else + untranspose_preprocess_kernel <<< nb, 256, 0, *stream >>>( + shuffled.shuffled, shuffled.desc, gpu_mask->get(), + processed_image.getGPUBuffer(), gpu_stats, sat_value, err_value, + shuffled.nblocks, shuffled.tail_src, shuffled.tail_elems, shuffled.tail_elem0, sat_code); cuda_err(cudaGetLastError()); if (copy_image_to_host) @@ -393,14 +450,23 @@ ImageStatistics ImagePreprocessorGPU::AnalyzeOnDevice(ImagePreprocessorBuffer &p cpu_stats[0] = ImageStatistics{.max_value = INT64_MIN, .min_value = INT64_MAX}; cuda_err(cudaMemcpyAsync(gpu_stats, cpu_stats.data(), sizeof(ImageStatistics), cudaMemcpyHostToDevice, *stream)); - preprocess_kernel <<< blocks, threads, 0, *stream >>>( - reinterpret_cast(gpu_decompressed_image.get()), - gpu_mask->get(), - processed_image.getGPUBuffer(), - gpu_stats, - sat_value, - err_value, - npixels); + + const bool narrow = ChooseNarrow(std::is_same::value); + const uint16_t sat_code = preprocessed_pixel::NarrowSaturatedCode( + std::min(saturation_limit, static_cast(sat_value))); + processed_image.SetNarrow(narrow, sat_code); + // The narrow image is the same allocation written two bytes at a time, so there is nothing to + // allocate and the wide pointer is simply reinterpreted. + if (narrow) + preprocess_kernel <<< blocks, threads, 0, *stream >>>( + reinterpret_cast(gpu_decompressed_image.get()), gpu_mask->get(), + reinterpret_cast(processed_image.getGPUBuffer()), + gpu_stats, sat_value, err_value, npixels, sat_code); + else + preprocess_kernel <<< blocks, threads, 0, *stream >>>( + reinterpret_cast(gpu_decompressed_image.get()), gpu_mask->get(), + processed_image.getGPUBuffer(), + gpu_stats, sat_value, err_value, npixels, sat_code); cuda_err(cudaGetLastError()); // The preprocessed image is 4 bytes per pixel - by far the largest transfer here - and every GPU // engine reads it straight from the device buffer, so it only comes back when a CPU engine needs it. diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h index 60679c7e..ad8f403b 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h @@ -40,6 +40,8 @@ class ImagePreprocessorGPU : public ImagePreprocessor { template ImageStatistics UntransposeAndAnalyze(ImagePreprocessorBuffer &processed_image, const BSLZ4ShuffledImage &shuffled, T err_value, T sat_value); + // Whether an image of this mode can keep its 16-bit width rather than being widened to int32. + bool ChooseNarrow(bool source_is_uint16) const; public: // copy_image_to_host copies the preprocessed image back after every frame. It is only needed when // something on the CPU reads it - the GPU engines all work off the device buffer - and at 4 bytes diff --git a/image_analysis/image_preprocessing/PreprocessedPixel.h b/image_analysis/image_preprocessing/PreprocessedPixel.h new file mode 100644 index 00000000..fb1cff82 --- /dev/null +++ b/image_analysis/image_preprocessing/PreprocessedPixel.h @@ -0,0 +1,88 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// CUDA only: the engines that read a preprocessed image on the device all compile as .cu. The CPU +// engines read ImagePreprocessorBuffer directly, which is int32 whatever the source was. +#include +#include + +#include "ImagePreprocessorBuffer.h" + +// How a preprocessed image is stored, and how every engine reads it. +// +// The preprocessed image is normally int32 per pixel, in the convention the whole pipeline shares: +// INT32_MIN is masked or bad, INT32_MAX is saturated, anything else is a real value. On a detector +// reading out 16 bits - which is what a fast acquisition uses, and what jfjoch writes for DECTRIS +// by default - that spends four bytes a pixel to carry two, and the per-pixel passes over a full +// frame are the largest single cost of the image loop. +// +// So a 16-bit source keeps its width. Two codes at the top of the range carry the two special +// states, and they are unambiguous because a detector's saturation limit bounds what a real value +// can be: +// +// 0xFFFF masked, or the source's own bad-pixel marker +// the saturation a pixel at or above the saturation limit. 0xFFFE where the limit leaves room +// code for it - a 16-bit EIGER declares a count-rate limit of a few thousand, so +// there is room to spare - and 0xFFFF where the limit is the full range, in +// which case the "is error" test has already claimed 0xFFFF and no pixel can +// be saturated at all, so 0xFFFE stays a real value. +// +// Either way a real value is strictly below the saturation limit and so below both codes. Nothing +// is clipped and nothing is lost. +// +// Engines do not learn a second convention. PixelView widens on load, so a masked pixel still reads +// as INT32_MIN and every existing `v != INT32_MIN && v != INT32_MAX` test keeps its meaning. +namespace preprocessed_pixel { + // Masked or bad. Always this value: a 16-bit source uses it as its own bad-pixel marker. + constexpr uint16_t NARROW_BAD = UINT16_MAX; + + // Which code means "saturated", given the saturation limit this image was preprocessed with. + // Falls back to NARROW_BAD when the limit is the whole range, which is exactly the case where + // no pixel can be saturated - so the code is then unreachable rather than ambiguous. + constexpr uint16_t NarrowSaturatedCode(int64_t saturation_limit) { + return saturation_limit <= UINT16_MAX - 2 ? static_cast(UINT16_MAX - 1) : NARROW_BAD; + } +} + +struct PixelView { + const int32_t *wide = nullptr; + const uint16_t *narrow = nullptr; // set instead of `wide` when the image kept its 16-bit width + uint16_t sat_code = preprocessed_pixel::NARROW_BAD; + + __host__ __device__ __forceinline__ int32_t Widen(uint16_t v) const { + // Order matters: where sat_code is NARROW_BAD the first test claims it, which is right - + // that is the case where nothing can be saturated. + if (v == preprocessed_pixel::NARROW_BAD) return INT32_MIN; + if (v == sat_code) return INT32_MAX; + return static_cast(v); + } + + // One pixel. The branch is on a pointer that is the same for every thread of every block, so it + // costs a predicted branch on kernels whose time is entirely the loads it selects between. + __host__ __device__ __forceinline__ int32_t operator[](size_t i) const { + return narrow ? Widen(narrow[i]) : wide[i]; + } + + // Pixels 4q..4q+3 in one transaction - 16 bytes wide or 8, whichever the image is. Callers that + // use this must keep to the same quad indexing for both forms, which they do: the layout is the + // same image either way, only narrower. + __host__ __device__ __forceinline__ void Load4(size_t q, int32_t out[4]) const { + if (narrow) { + const ushort4 v = reinterpret_cast(narrow)[q]; + out[0] = Widen(v.x); out[1] = Widen(v.y); out[2] = Widen(v.z); out[3] = Widen(v.w); + } else { + const int4 v = reinterpret_cast(wide)[q]; + out[0] = v.x; out[1] = v.y; out[2] = v.z; out[3] = v.w; + } + } +}; + +// The only way an engine should get a view of a preprocessed image: it reads the width off the +// buffer, so the pointer and the width cannot be taken from different places and disagree. +inline PixelView ViewOf(const ImagePreprocessorBuffer &image) { + return image.IsNarrow() + ? PixelView{nullptr, image.getGPUBufferNarrow(), image.NarrowSatCode()} + : PixelView{image.getGPUBuffer(), nullptr, preprocessed_pixel::NARROW_BAD}; +} diff --git a/image_analysis/roi/ROIIntegrationGPU.cu b/image_analysis/roi/ROIIntegrationGPU.cu index db215670..9aff2e94 100644 --- a/image_analysis/roi/ROIIntegrationGPU.cu +++ b/image_analysis/roi/ROIIntegrationGPU.cu @@ -4,6 +4,7 @@ #include #include "ROIIntegrationGPU.h" +#include "../image_preprocessing/PreprocessedPixel.h" #include "../../common/DiffractionExperiment.h" inline void cuda_err(cudaError_t val) { @@ -16,7 +17,7 @@ inline void cuda_err(cudaError_t val) { __global__ void gpu_roi( const uint16_t *__restrict__ roi_map, - const int32_t *__restrict__ input_buffer, + PixelView input_buffer, size_t num_pixels, size_t width, int roi_count, @@ -138,7 +139,7 @@ void ROIIntegrationGPU::Run(const ImagePreprocessorBuffer &image, std::map>>( - gpu_roi_map, image.getGPUBuffer(), npixel, width, roi_count, + gpu_roi_map, ViewOf(image), npixel, width, roi_count, gpu_sum, gpu_sum2, gpu_pixels, gpu_x_weighted, gpu_y_weighted, gpu_max); cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 6ba4a8f2..23839d13 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "AdaptiveSpotFinderGPU.h" +#include "../image_preprocessing/PreprocessedPixel.h" #include "AdaptiveThreshold.h" #include "../../common/JFJochException.h" @@ -35,7 +36,7 @@ __device__ __forceinline__ void flush_ring(unsigned long long *s_sum, unsigned l __global__ void reduce_rings_shared( const uint16_t *__restrict__ pixel_to_bin, const float *__restrict__ corrections, - const int32_t *__restrict__ image, + PixelView image, const float *__restrict__ mean, const float *__restrict__ sigma, float clip_k, @@ -80,13 +81,13 @@ __global__ void reduce_rings_shared( 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]; + int32_t vq[4]; + image.Load4(q, vq); const ushort4 b4 = reinterpret_cast(pixel_to_bin)[q]; float4 c4 = make_float4(0.0f, 0.0f, 0.0f, 0.0f); if (accumulate_corrected) c4 = reinterpret_cast(corrections)[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 float cq[4] = {c4.x, c4.y, c4.z, c4.w}; @@ -173,7 +174,7 @@ __global__ void reduce_rings_shared( __global__ void reduce_rings_global( const uint16_t *__restrict__ pixel_to_bin, const float *__restrict__ corrections, - const int32_t *__restrict__ image, + PixelView image, const float *__restrict__ mean, const float *__restrict__ sigma, float clip_k, @@ -225,7 +226,7 @@ __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. -__global__ void flag_strong(const int32_t *__restrict__ image, +__global__ void flag_strong(PixelView image, const uint16_t *__restrict__ pixel_to_bin, const float *__restrict__ thr, uint32_t *__restrict__ strong, @@ -238,9 +239,9 @@ __global__ void flag_strong(const int32_t *__restrict__ image, 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]; + int32_t vq[4]; + image.Load4(q, vq); 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}; uint32_t bits = 0; @@ -358,12 +359,12 @@ void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, flo const size_t shared = accumulate_corrected ? shared_plain : shared_clip; const int blocks = accumulate_corrected ? reduce_blocks_plain : reduce_blocks_clip; reduce_rings_shared<<>>( - gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma, + gpu_pixel_to_bin->get(), gpu_corrections->get(), ViewOf(image), gpu_mean, gpu_sigma, clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr, npix, nbins); } else { reduce_rings_global<<>>( - gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma, + gpu_pixel_to_bin->get(), gpu_corrections->get(), ViewOf(image), gpu_mean, gpu_sigma, clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr, npix, nbins); } @@ -471,7 +472,7 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); flag_strong<<>>( - image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins); + ViewOf(image), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins); // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream, // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host. } @@ -483,6 +484,6 @@ void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector &mask) { const std::vector &AdaptiveSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { - extractor.Extract(gpu_strong, image.getGPUBuffer(), settings, components); + extractor.Extract(gpu_strong, ViewOf(image), settings, components); return components; } diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index 27f3f82a..f5a899bb 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -5,6 +5,7 @@ // Copyright (2019-2023) Paul Scherrer Institute #include "ImageSpotFinderGPU.h" +#include "../image_preprocessing/PreprocessedPixel.h" #include "../../common/JFJochException.h" struct spot_parameters { @@ -98,7 +99,7 @@ __device__ __forceinline__ uint8_t pixel_result(const spot_parameters& params, c // Every cuda thread is associated with a particular column. The thread maintains // the sum/sum2/count values in shared memory for it's column. To do this, the input // pixel values for the hight of the aggregation window are saved in shared memory. -__global__ void analyze_pixel(const int32_t *in, uint32_t *prev_out, uint32_t *out, const spot_parameters params) +__global__ void analyze_pixel(PixelView in, uint32_t *prev_out, uint32_t *out, const spot_parameters params) { // assumption: 2 * params.nby + 1 <= params.rows and 2 * params.nbx + 1 <= params.width const int32_t window = 2 * (int)ImageSpotFinder::NBX + 1; // vertical window @@ -265,7 +266,7 @@ void ImageSpotFinderGPU::SetResolutionMask(const std::vector &mask) { const std::vector &ImageSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { - extractor.Extract(gpu_out_1, image.getGPUBuffer(), settings, components); + extractor.Extract(gpu_out_1, ViewOf(image), settings, components); return components; } @@ -299,9 +300,9 @@ void ImageSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, const Spot cuda_err(cudaMemsetAsync(gpu_out_0, 0, OutputSize() * sizeof(uint32_t), *stream)); cuda_err(cudaMemsetAsync(gpu_out_1, 0, OutputSize() * sizeof(uint32_t), *stream)); analyze_pixel<<>> - (image.getGPUBuffer(), gpu_out_1, gpu_out_0, spot_params); + (ViewOf(image), gpu_out_1, gpu_out_0, spot_params); analyze_pixel<<>> - (image.getGPUBuffer(), gpu_out_0, gpu_out_1, spot_params); + (ViewOf(image), gpu_out_0, gpu_out_1, spot_params); // The bit buffer stays on the device - ExtractComponents reads it there. cuda_err(cudaStreamSynchronize(*stream)); } diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 2a5e9c38..85478cce 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -73,7 +73,7 @@ __global__ void scan_block_counts(const uint32_t *__restrict__ in, uint32_t *__r // holds a handful of them and serial emission is both trivially ordered and fast; the parallelism // comes from the block count. __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, + const uint32_t *__restrict__ block_offset, PixelView image, uint32_t *__restrict__ out_index, int32_t *__restrict__ out_value, size_t nwords, uint32_t capacity) { if (threadIdx.x != 0) return; @@ -318,7 +318,7 @@ void SpotExtractorGPU::SetResolutionMask(const std::vector &packed_mas cuda_err(cudaStreamSynchronize(*stream)); } -void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_image, +void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, PixelView gpu_image, const SpotFindingSettings &settings, std::vector &spots) { const int max_pix = static_cast(settings.max_pix_per_spot); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index e00a86cc..b09f81aa 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -34,6 +34,7 @@ #include "../../common/DiffractionSpot.h" #include "../indexing/CUDAMemHelpers.h" #include "SpotFindingSettings.h" +#include "../image_preprocessing/PreprocessedPixel.h" // Per-component sums, in exactly the form DiffractionSpot holds them: x and y are sum(col*photons) // and sum(line*photons), not a centroid. @@ -86,6 +87,6 @@ public: // gpu_strong is the finder's device bit buffer, gpu_image the preprocessed image it was built // from. Fills spots with every component of at most max-pix pixels, in the same order the host // extractor would. - void Extract(const uint32_t *gpu_strong, const int32_t *gpu_image, + void Extract(const uint32_t *gpu_strong, PixelView gpu_image, const SpotFindingSettings &settings, std::vector &spots); }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 5f725eec..f21ed4db 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1162,6 +1162,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b : calibration ? "powder calibration" : "azimuthal integration"; logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]", images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name); + MarkPhase(observer, full ? "Setting up" : calibration ? "Powder calibration" : "Azimuthal integration"); // Full-analysis shared engines. @@ -1670,6 +1671,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } }; std::atomic total_uncompressed_bytes = 0; + std::once_flag depth_logged; // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of @@ -1765,6 +1767,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (dataset->efficiency.size() > image_idx) msg.image_collection_efficiency = dataset->efficiency[image_idx]; total_uncompressed_bytes += msg.image.GetUncompressedSize(); + // Said once, from the first frame: it is the difference between moving two bytes a pixel + // through every per-pixel pass of the loop and moving four, and a fast acquisition reads + // out 16 bits. + std::call_once(depth_logged, [&] { + const int bits = 8 * static_cast(msg.image.GetByteDepth()); + if (bits == 16) + logger.Info("Images are 16-bit: the preprocessed frame keeps that width"); + else + logger.Info("Images are {}-bit: the preprocessed frame is 32-bit", bits); + }); const auto t0 = std::chrono::steady_clock::now(); try {