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 {