From f4252ad39d108c5363f962a936e80ee47c8459af Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 22:37:37 +0200 Subject: [PATCH] miniCBF frames are decoded on the device The CBF reader decoded every frame's x-CBF_BYTE_OFFSET stream on the worker thread - 19-30 % of the image-loop worker CPU on CBF sets (7qij, 4mtk) - and the 4-byte image was then uploaded. ReadRawImage now hands over the stream as the file stores it, in the image's own buffer, as a new read-only CompressionAlgorithm::CBF_BYTE_OFFSET (mode Int32, what it decodes to), and ImagePreprocessorGPU decodes it on the card into the same int32 image and runs the same preprocessing kernel on it. The device decode (CBFDecoderGPU) cuts the stream into 256-byte tiles. A token is 1, 3, 7 or 15 bytes, so a tile can be entered at 15 offsets; for each one thread works out where the decode leaves the tile, how many pixels it emitted and their sum modulo 2^32 (the walks from different entries meet within a few bytes, after which they are one walk). An exclusive scan composing these maps gives each tile its true entry, first pixel and running value, and a second walk writes the pixels. A stream that ends before the image is refused as on the host. The host decoder moves to compression/ (DecodeCBFByteOffset) and is reached through JFJochDecompressPtr, so every host consumer - the CPU build, the azimuthal-integration worker, the pre-scan's shadow and gain passes - decodes exactly as before, now from the stream. It is the CPU path, not a fallback: a frame the device refuses is not retried on the host. A process file written from CBF input stores the decoded int32 image, as it always did. .cbf.gz is still inflated on the host. Tests: CBFDecoderGPU_MatchesHost (every token size, 0x80 inside payloads, escape runs across tiles, trailing CIF text), CBFDecoderGPU_ShortStreamIsRefused, ImagePreprocessorGPU_StoredEncodingsMatchCPU. Checked locally against the host decoder on 36 real frames of 9 open-arm CBF sets (signed 32-bit, signed and unsigned 16-bit, plain and gzip): identical. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- compression/CBFByteOffset.cpp | 89 +++++++ compression/CBFByteOffset.h | 12 + compression/CMakeLists.txt | 2 + compression/CompressionAlgorithmEnum.h | 6 +- compression/JFJochDecompress.h | 7 + image_analysis/MXAnalysisWithoutFPGA.cpp | 4 + .../image_preprocessing/CBFDecoderGPU.cu | 221 ++++++++++++++++++ .../image_preprocessing/CBFDecoderGPU.h | 67 ++++++ .../image_preprocessing/CMakeLists.txt | 1 + .../ImagePreprocessorGPU.cu | 25 +- .../ImagePreprocessorGPU.h | 9 + reader/JFJochCBFReader.cpp | 12 +- reader/MiniCBF.cpp | 97 ++------ reader/MiniCBF.h | 11 +- rugnux/RugnuxPipeline.h | 9 + tests/CMakeLists.txt | 1 + tests/StoredImageDecodeGPUTest.cpp | 181 ++++++++++++++ 17 files changed, 662 insertions(+), 92 deletions(-) create mode 100644 compression/CBFByteOffset.cpp create mode 100644 compression/CBFByteOffset.h create mode 100644 image_analysis/image_preprocessing/CBFDecoderGPU.cu create mode 100644 image_analysis/image_preprocessing/CBFDecoderGPU.h create mode 100644 tests/StoredImageDecodeGPUTest.cpp diff --git a/compression/CBFByteOffset.cpp b/compression/CBFByteOffset.cpp new file mode 100644 index 000000000..f618e4e02 --- /dev/null +++ b/compression/CBFByteOffset.cpp @@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "CBFByteOffset.h" + +#include +#include + +#if defined(__SSE2__) || defined(_M_X64) || defined(_M_AMD64) +# include +# define MINICBF_STREAM_STORE 1 +#endif + +#include "../common/JFJochException.h" + +namespace { + +int16_t ReadI16(const uint8_t *p) { int16_t v; std::memcpy(&v, p, 2); return v; } +int32_t ReadI32(const uint8_t *p) { int32_t v; std::memcpy(&v, p, 4); return v; } +int64_t ReadI64(const uint8_t *p) { int64_t v; std::memcpy(&v, p, 8); return v; } + +// Nothing on this side reads the decoded image back: it goes to the GPU by DMA, or a CPU engine +// walks it once and drops it. An ordinary store therefore pays twice over - it reads every cache +// line before overwriting it, and it evicts 72 MB of live cache to make room - so the pixels go +// straight past the cache. Measured on an 18 Mpx frame: 27 ms of ordinary stores against 11 ms of +// these. The bytes written are the same either way; the fence is what makes them visible to the +// copy that follows. +void StorePixel(int32_t *out, int32_t value) { +#ifdef MINICBF_STREAM_STORE + _mm_stream_si32(out, value); +#else + *out = value; +#endif +} + +void StoreFence() { +#ifdef MINICBF_STREAM_STORE + _mm_sfence(); +#endif +} + +} // namespace + +// The x-CBF_BYTE_OFFSET scheme: a running value, each pixel stored as a delta in the smallest +// container that holds it, escaping to the next size with that container's most negative value. +// Following Bernstein & Hammersley (2006) Int. Tables Cryst. G, 37-43 +void DecodeCBFByteOffset(const uint8_t *data, size_t size, int32_t *out, size_t n_pixels) { + int64_t value = 0; + size_t pos = 0; + size_t written = 0; + + while (written < n_pixels) { + if (pos >= size) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "miniCBF byte-offset stream ended after " + std::to_string(written) + + " of " + std::to_string(n_pixels) + " pixels"); + + int64_t delta = static_cast(data[pos]); + pos += 1; + + if (delta == -128) { + if (pos + 2 > size) break; + delta = ReadI16(data + pos); + pos += 2; + if (delta == -32768) { + if (pos + 4 > size) break; + delta = ReadI32(data + pos); + pos += 4; + if (delta == INT32_MIN) { + if (pos + 8 > size) break; + delta = ReadI64(data + pos); + pos += 8; + } + } + } + + value += delta; + StorePixel(out + written, static_cast(value)); + written++; + } + + StoreFence(); + + if (written != n_pixels) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "miniCBF byte-offset stream ended after " + std::to_string(written) + + " of " + std::to_string(n_pixels) + " pixels"); +} + diff --git a/compression/CBFByteOffset.h b/compression/CBFByteOffset.h new file mode 100644 index 000000000..b2c8665c4 --- /dev/null +++ b/compression/CBFByteOffset.h @@ -0,0 +1,12 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// x-CBF_BYTE_OFFSET -> int32. Deltas against a running value, smallest container first: int8, +// escaping to int16 via -128, to int32 via -32768, to int64 via INT32_MIN. Little-endian, packed. +// Throws if the stream ends before n_pixels are produced. out must hold n_pixels. +void DecodeCBFByteOffset(const uint8_t *data, size_t size, int32_t *out, size_t n_pixels); diff --git a/compression/CMakeLists.txt b/compression/CMakeLists.txt index 385615cc3..e69f3c462 100644 --- a/compression/CMakeLists.txt +++ b/compression/CMakeLists.txt @@ -16,6 +16,8 @@ ADD_LIBRARY(Compression STATIC JFJochCompressor.h JFJochDecompress.h BitShuffleBlock.h + CBFByteOffset.cpp + CBFByteOffset.h MaxCompressedSize.cpp MaxCompressedSize.h) set_target_properties(Compression PROPERTIES POSITION_INDEPENDENT_CODE ON) diff --git a/compression/CompressionAlgorithmEnum.h b/compression/CompressionAlgorithmEnum.h index b8dbded6d..56406b6dc 100644 --- a/compression/CompressionAlgorithmEnum.h +++ b/compression/CompressionAlgorithmEnum.h @@ -6,5 +6,9 @@ // LZ4_NO_SHUFFLE is READ-ONLY: HDF5 filter 32004, plain LZ4 with no bit shuffle, written by DECTRIS // Eiger firmware 1.x. Nothing here produces it. It is last so the values of the others do not move - // the enum crosses the CBOR stream. +// +// CBF_BYTE_OFFSET is READ-ONLY too, and never leaves the process: it is how the miniCBF reader hands +// over a frame as the file stores it, so that the GPU path can upload those bytes and decode them on +// the device. The image's mode is what it decodes to (Int32). enum class CompressionAlgorithm {BSHUF_LZ4, BSHUF_ZSTD, BSHUF_ZSTD_RLE, BSHUF_ZSTD_RLE_HUFF, NO_COMPRESSION, - LZ4_NO_SHUFFLE}; + LZ4_NO_SHUFFLE, CBF_BYTE_OFFSET}; diff --git a/compression/JFJochDecompress.h b/compression/JFJochDecompress.h index c5a90df63..cc027d330 100644 --- a/compression/JFJochDecompress.h +++ b/compression/JFJochDecompress.h @@ -12,6 +12,7 @@ #include #include "BitShuffleBlock.h" +#include "CBFByteOffset.h" #include "../compression/CompressionAlgorithmEnum.h" #include "../common/JFJochException.h" @@ -233,6 +234,12 @@ inline void JFJochDecompressPtr(uint8_t *output, JFJochDecompressLZ4Ptr(output, source, source_size, nelements, elem_size); return; } + if (algorithm == CompressionAlgorithm::CBF_BYTE_OFFSET) { + if (elem_size != sizeof(int32_t)) + throw JFJochException(JFJochExceptionCategory::Compression, "CBF byte offset decodes to int32"); + DecodeCBFByteOffset(source, source_size, reinterpret_cast(output), nelements); + return; + } size_t block_size = 0; if (algorithm != CompressionAlgorithm::NO_COMPRESSION) diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 9f6354a7c..3409013e9 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -121,6 +121,10 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // sees the same error it saw before any of this existed - but a GPU-side problem costs speed // rather than the acquisition. Unless the context is lost, which the host route cannot survive. cuda_throw_if_context_lost(); + // The miniCBF encoding is decoded on the device whenever there is one; the host decoder is the + // CPU path, not a slower second try. A stream the device refuses, the host refuses as well. + if (output.image.GetCompressionAlgorithm() == CompressionAlgorithm::CBF_BYTE_OFFSET) + throw; spdlog::warn("Device decoding failed ({}), falling back to host decompression", e.what()); decoded_on_device = false; // Handled here, so it must not stay behind as the thread's last error: the host route below diff --git a/image_analysis/image_preprocessing/CBFDecoderGPU.cu b/image_analysis/image_preprocessing/CBFDecoderGPU.cu new file mode 100644 index 000000000..13cf93845 --- /dev/null +++ b/image_analysis/image_preprocessing/CBFDecoderGPU.cu @@ -0,0 +1,221 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include + +#include "CBFDecoderGPU.h" +#include "../../common/JFJochException.h" + +namespace { + void cuda_err(cudaError_t val) { + if (val != cudaSuccess) + throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val)); + } + + constexpr size_t TILE_BYTES = 256; + + // The token starting at byte p: its length, and its delta (low 32 bits) in `delta`. 0 when the + // token runs past the end of the stream - where the host decoder stops too. + // Following Bernstein & Hammersley (2006) Int. Tables Cryst. G, 37-43 + __device__ inline int TokenAt(const uint8_t *__restrict__ d, size_t size, size_t p, uint32_t &delta) { + const int8_t b = static_cast(d[p]); + if (b != -128) { + delta = static_cast(static_cast(b)); + return 1; + } + if (p + 3 > size) return 0; + const int16_t v16 = static_cast(d[p + 1] | (d[p + 2] << 8)); + if (v16 != -32768) { + delta = static_cast(static_cast(v16)); + return 3; + } + if (p + 7 > size) return 0; + const uint32_t v32 = d[p + 3] | (d[p + 4] << 8) | (d[p + 5] << 16) | (static_cast(d[p + 6]) << 24); + if (v32 != 0x80000000u) { + delta = v32; + return 7; + } + if (p + 15 > size) return 0; + // An int64 delta: only its low 32 bits reach an int32 pixel. + delta = d[p + 7] | (d[p + 8] << 8) | (d[p + 9] << 16) | (static_cast(d[p + 10]) << 24); + return 15; + } + + // One step of a walk at position pos: false when the stream has ended there. + __device__ inline bool Step(const uint8_t *__restrict__ d, size_t size, size_t &pos, + uint32_t &count, uint32_t &sum) { + uint32_t delta = 0; + const int len = TokenAt(d, size, pos, delta); + if (len == 0) return false; + pos += len; + count++; + sum += delta; + return true; + } + + __device__ inline uint8_t ExitOf(size_t pos, size_t tile_end, size_t size, bool ended) { + return (ended || pos >= size) ? CBFTileMap::DEAD : static_cast(pos - tile_end); + } + + __global__ void tile_map_kernel(const uint8_t *__restrict__ d, size_t size, size_t ntiles, + CBFTileMap *__restrict__ maps) { + const size_t t = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; + if (t >= ntiles) return; + const size_t start = t * TILE_BYTES; + const size_t end = start + TILE_BYTES; + const size_t limit = min(end, size); + CBFTileMap m; + + // The walk from the tile's first byte, in full. + size_t pos0 = start; + uint32_t count0 = 0, sum0 = 0; + bool ended0 = false; + while (pos0 < limit && !ended0) + ended0 = !Step(d, size, pos0, count0, sum0); + m.count[0] = count0; + m.sum[0] = sum0; + m.exit[0] = ExitOf(pos0, end, size, ended0); + + // The walks from the other entries. Each is stepped beside the first walk, whichever is behind + // moving, until both stand on the same byte - from there they are one walk, and the rest is + // the first walk's remainder. + for (int e = 1; e < CBFTileMap::DEAD; e++) { + size_t a = start, b = start + e; + uint32_t count_a = 0, sum_a = 0, count_b = 0, sum_b = 0; + bool a_ended = false, b_ended = false, merged = false; + while (b < limit && !b_ended) { + if (a == b && !a_ended) { + merged = true; + break; + } + if (a < b && !a_ended) + a_ended = !Step(d, size, a, count_a, sum_a); + else + b_ended = !Step(d, size, b, count_b, sum_b); + } + if (merged) { + m.count[e] = count_b + (count0 - count_a); + m.sum[e] = sum_b + (sum0 - sum_a); + m.exit[e] = m.exit[0]; + } else { + m.count[e] = count_b; + m.sum[e] = sum_b; + m.exit[e] = ExitOf(b, end, size, b_ended); + } + } + m.count[CBFTileMap::DEAD] = 0; + m.sum[CBFTileMap::DEAD] = 0; + m.exit[CBFTileMap::DEAD] = CBFTileMap::DEAD; + maps[t] = m; + } + + // a, then b. The scan's look-back also composes the uninitialised state of tiles before the first + // one and discards the result, so an exit is masked to the table rather than trusted to be in it. + struct ComposeTileMaps { + __host__ __device__ CBFTileMap operator()(const CBFTileMap &a, const CBFTileMap &b) const { + CBFTileMap r; + for (int e = 0; e < 16; e++) { + const int mid = a.exit[e] & 15; + r.exit[e] = b.exit[mid]; + r.count[e] = a.count[e] + b.count[mid]; + r.sum[e] = a.sum[e] + b.sum[mid]; + } + return r; + } + }; + + __global__ void tile_write_kernel(const uint8_t *__restrict__ d, size_t size, size_t ntiles, + const CBFTileMap *__restrict__ entries, const CBFTileMap *__restrict__ maps, + int32_t *__restrict__ out, uint32_t npixels, uint32_t *__restrict__ total) { + const size_t t = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; + if (t >= ntiles) return; + const CBFTileMap &in = entries[t]; + if (t == ntiles - 1) + *total = in.count[0] + maps[t].count[in.exit[0]]; + if (in.exit[0] == CBFTileMap::DEAD) + return; + const size_t end = min((t + 1) * TILE_BYTES, size); + size_t pos = t * TILE_BYTES + in.exit[0]; + uint32_t pixel = in.count[0]; + uint32_t value = in.sum[0]; + while (pos < end && pixel < npixels) { + uint32_t delta = 0; + const int len = TokenAt(d, size, pos, delta); + if (len == 0) break; + pos += len; + value += delta; + out[pixel++] = static_cast(value); + } + } +} + +CBFDecoderGPU::CBFDecoderGPU(std::shared_ptr in_stream) + : stream(std::move(in_stream)), gpu_total(1), host_total(1) {} + +void CBFDecoderGPU::Decode(const CompressedImage &image, int32_t *gpu_out) { + const size_t size = image.GetCompressedSize(); + last_npixels = image.GetWidth() * image.GetHeight(); + if (size == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "miniCBF byte-offset stream ended after 0 of " + std::to_string(last_npixels) + " pixels"); + const size_t ntiles = (size + TILE_BYTES - 1) / TILE_BYTES; + + if (size > compressed_capacity || ntiles > tile_capacity) { + // A frame's stream is a little longer or shorter than the last one's: grow with slack. + cuda_err(cudaStreamSynchronize(*stream)); + compressed_capacity = size + size / 4; + tile_capacity = (compressed_capacity + TILE_BYTES - 1) / TILE_BYTES; + host_compressed = CudaHostPtr(compressed_capacity); + gpu_compressed = CudaDevicePtr(compressed_capacity); + gpu_maps = CudaDevicePtr(tile_capacity); + gpu_entries = CudaDevicePtr(tile_capacity); + CBFTileMap identity{}; + cuda_err(cub::DeviceScan::ExclusiveScan(nullptr, scan_temp_bytes, gpu_maps.get(), gpu_entries.get(), + ComposeTileMaps(), identity, tile_capacity, *stream)); + gpu_scan_temp = CudaDevicePtr(scan_temp_bytes); + } + + // The upload waits for the previous frame's to have left the staging buffer. + cuda_err(cudaStreamSynchronize(*stream)); + std::memcpy(host_compressed.get(), image.GetCompressed(), size); + cuda_err(cudaEventRecord(decode_start, *stream)); + cuda_err(cudaMemcpyAsync(gpu_compressed, host_compressed.get(), size, cudaMemcpyHostToDevice, *stream)); + + const int threads = 128; + const int blocks = static_cast((ntiles + threads - 1) / threads); + tile_map_kernel<<>>(gpu_compressed, size, ntiles, gpu_maps); + cuda_err(cudaGetLastError()); + + // Entering tile 0 at its first byte: every entry maps to itself, nothing emitted yet. + CBFTileMap identity{}; + for (int e = 0; e < 16; e++) + identity.exit[e] = static_cast(e); + size_t temp_bytes = scan_temp_bytes; + cuda_err(cub::DeviceScan::ExclusiveScan(gpu_scan_temp.get(), temp_bytes, gpu_maps.get(), gpu_entries.get(), + ComposeTileMaps(), identity, ntiles, *stream)); + + tile_write_kernel<<>>(gpu_compressed, size, ntiles, gpu_entries, gpu_maps, + gpu_out, static_cast(last_npixels), gpu_total); + cuda_err(cudaGetLastError()); + cuda_err(cudaEventRecord(decode_stop, *stream)); + cuda_err(cudaMemcpyAsync(host_total.get(), gpu_total.get(), sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); +} + +void CBFDecoderGPU::ThrowIfDecodeFailed() const { + const size_t total = host_total.get()[0]; + if (total < last_npixels) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "miniCBF byte-offset stream ended after " + std::to_string(total) + + " of " + std::to_string(last_npixels) + " pixels"); +} + +float CBFDecoderGPU::GetDecodeTime_s() const { + float ms = 0.0f; + if (cudaEventElapsedTime(&ms, decode_start, decode_stop) != cudaSuccess) + return 0.0f; + return ms * 1e-3f; +} diff --git a/image_analysis/image_preprocessing/CBFDecoderGPU.h b/image_analysis/image_preprocessing/CBFDecoderGPU.h new file mode 100644 index 000000000..b264f5090 --- /dev/null +++ b/image_analysis/image_preprocessing/CBFDecoderGPU.h @@ -0,0 +1,67 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../../common/CompressedImage.h" +#include "../indexing/CUDAMemHelpers.h" + +// What one tile of the byte-offset stream does to a decode that enters it at each of the 15 byte +// positions a token can straddle the tile boundary by: where the decode leaves the tile (the offset +// into the next one), how many pixels it emitted, and what they added to the running value. Entry and +// exit 15 is "the stream has ended", which every later tile passes on. +struct CBFTileMap { + static constexpr int DEAD = 15; + uint32_t count[16]; + uint32_t sum[16]; // modulo 2^32, which is all an int32 pixel keeps of the int64 running value + uint8_t exit[16]; +}; + +// x-CBF_BYTE_OFFSET decoded ON THE DEVICE, exactly as DecodeCBFByteOffset does on the host: the same +// int32 pixels, and the same refusal of a stream that ends before the image does. +// +// The format is a running sum of deltas whose tokens are 1, 3, 7 or 15 bytes long, so where a token +// starts depends on every token before it. The stream is cut into tiles; for each tile and each of the +// 15 offsets a decode can enter it at, one thread works out where it leaves, how many pixels it emits +// and what they sum to (the walks from different offsets meet within a few bytes, after which they +// are the same walk). Composing those maps across tiles is a scan, which gives every tile its true +// entry offset, its first pixel and the running value there; a second walk of each tile then writes +// the pixels. +class CBFDecoderGPU { + std::shared_ptr stream; + + // The reader's buffer is pageable; copied here first, the upload is one DMA (see BSLZ4DecoderGPU). + CudaHostPtr host_compressed; + CudaDevicePtr gpu_compressed; + size_t compressed_capacity = 0; + + CudaDevicePtr gpu_maps; + CudaDevicePtr gpu_entries; // exclusive scan of gpu_maps: the state each tile is entered in + CudaDevicePtr gpu_scan_temp; + size_t scan_temp_bytes = 0; + size_t tile_capacity = 0; + + CudaDevicePtr gpu_total; // pixels the whole stream decodes to + CudaHostPtr host_total; + size_t last_npixels = 0; + + CudaEvent decode_start; + CudaEvent decode_stop; + +public: + explicit CBFDecoderGPU(std::shared_ptr stream); + ~CBFDecoderGPU() { if (stream) cudaStreamSynchronize(*stream); } + + // Upload the stream and queue its decode into gpu_out (image width x height int32 pixels) on the + // decoder's stream. The caller synchronises, then calls ThrowIfDecodeFailed(). + void Decode(const CompressedImage &image, int32_t *gpu_out); + + // Report a stream that ended before the image did. MUST be called after the caller has + // synchronised the stream. + void ThrowIfDecodeFailed() const; + + // Device time of the last upload and decode, in seconds. Valid after the caller synchronises. + [[nodiscard]] float GetDecodeTime_s() const; +}; diff --git a/image_analysis/image_preprocessing/CMakeLists.txt b/image_analysis/image_preprocessing/CMakeLists.txt index 0b958acc5..fe6dd3ea9 100644 --- a/image_analysis/image_preprocessing/CMakeLists.txt +++ b/image_analysis/image_preprocessing/CMakeLists.txt @@ -11,6 +11,7 @@ IF (JFJOCH_CUDA_AVAILABLE) ../indexing/CUDAMemHelpers.h ImagePreprocessorGPU.cu ImagePreprocessorGPU.h BSLZ4DecoderGPU.cu BSLZ4DecoderGPU.h BSLZ4DecodeWarp.h + CBFDecoderGPU.cu CBFDecoderGPU.h ImagePreprocessorBufferGPU.cu ImagePreprocessorBufferGPU.h) ENDIF() \ No newline at end of file diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 079525caf..5ec3cf53d 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -340,7 +340,11 @@ ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experime } float ImagePreprocessorGPU::GetLastDecompressionTime_s() const { - return bslz4_decoder ? bslz4_decoder->GetDecodeTime_s() : 0.0f; + switch (last_decode) { + case LastDecode::BSLZ4: return bslz4_decoder->GetDecodeTime_s(); + case LastDecode::CBF: return cbf_decoder->GetDecodeTime_s(); + default: return 0.0f; + } } ImagePreprocessorGPU::~ImagePreprocessorGPU() { @@ -406,8 +410,13 @@ ImageStatistics ImagePreprocessorGPU::Analyze(ImagePreprocessorBuffer &processed bool ImagePreprocessorGPU::AnalyzeCompressed(ImagePreprocessorBuffer &processed_image, const CompressedImage &image, ImageStatistics &stats) { + if (image.GetCompressionAlgorithm() == CompressionAlgorithm::CBF_BYTE_OFFSET) { + stats = AnalyzeStoredInt32(processed_image, image); + return true; + } if (!BSLZ4DecoderGPU::Supports(image)) return false; // caller decompresses on the host and uses Analyze() + last_decode = LastDecode::BSLZ4; if (image.GetUncompressedSize() != npixels * image.GetByteDepth()) return false; @@ -529,6 +538,20 @@ ImageStatistics ImagePreprocessorGPU::Analyze(ImagePreprocessorBuffer &processed return AnalyzeOnDevice(processed_image, err_value, sat_value); } +ImageStatistics ImagePreprocessorGPU::AnalyzeStoredInt32(ImagePreprocessorBuffer &processed_image, + const CompressedImage &image) { + if (!gpu_decompressed_image.get()) + gpu_decompressed_image = CudaDevicePtr(npixels * sizeof(uint32_t)); + + last_decode = LastDecode::CBF; + if (!cbf_decoder) + cbf_decoder = std::make_unique(stream); + cbf_decoder->Decode(image, reinterpret_cast(gpu_decompressed_image.get())); + const ImageStatistics ret = AnalyzeOnDevice(processed_image, INT32_MIN, INT32_MAX); + cbf_decoder->ThrowIfDecodeFailed(); + return ret; +} + // Everything after the image is on the device, shared by the host-upload and the device-decode // entry points so the two cannot drift apart. template diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h index 9a056de92..4e7c8e640 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h @@ -7,6 +7,7 @@ #include "ImagePreprocessor.h" #include "BSLZ4DecoderGPU.h" +#include "CBFDecoderGPU.h" #include "../indexing/CUDAMemHelpers.h" #include "../indexing/CudaSharedTables.h" @@ -33,6 +34,14 @@ class ImagePreprocessorGPU : public ImagePreprocessor { // Built on first use: a decoder that can serve this engine's images, sized to the frame. std::unique_ptr bslz4_decoder; + std::unique_ptr cbf_decoder; + // Which decoder the last AnalyzeCompressed used, for GetLastDecompressionTime_s. + enum class LastDecode { None, BSLZ4, CBF } last_decode = LastDecode::None; + + // The file encoding of the miniCBF reader (CompressionAlgorithmEnum.h): the stored bytes are + // uploaded, decoded into gpu_decompressed_image as the int32 image the host would have produced, + // and preprocessed by the same kernel as an uploaded int32 image. + ImageStatistics AnalyzeStoredInt32(ImagePreprocessorBuffer &processed_image, const CompressedImage &image); template ImageStatistics Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *input, T err_value, T sat_value); // Preprocess an image already sitting in gpu_decompressed_image (the host-upload path). diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index 8c06459d3..8cec01a4c 100644 --- a/reader/JFJochCBFReader.cpp +++ b/reader/JFJochCBFReader.cpp @@ -394,7 +394,17 @@ bool JFJochCBFReader::HasImage(int64_t image_number) const { bool JFJochCBFReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { if (!HasImage(image_number)) return false; - image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer); + // The file as it is stored: the consumer decodes the byte-offset stream, on the device where there + // is one, so the host neither decodes nor uploads a four-byte image. The stream sits in the image's + // own buffer, which the image keeps (read_buffer is not part of it). + size_t start = 0; + const auto h = minicbf::ReadCompressed(files_[image_number], image.image_buffer, start); + if (h.nx != header0_.nx || h.ny != header0_.ny) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "CBF image size differs from the first image of the sweep"); + image.image = CompressedImage(image.image_buffer.data() + start, image.image_buffer.size() - start, + static_cast(header0_.nx), static_cast(header0_.ny), + CompressedImageMode::Int32, CompressionAlgorithm::CBF_BYTE_OFFSET); return true; } diff --git a/reader/MiniCBF.cpp b/reader/MiniCBF.cpp index aa48f0aba..3793b0cc1 100644 --- a/reader/MiniCBF.cpp +++ b/reader/MiniCBF.cpp @@ -15,12 +15,8 @@ #include #include -#if defined(__SSE2__) || defined(_M_X64) || defined(_M_AMD64) -# include -# define MINICBF_STREAM_STORE 1 -#endif - #include "../common/JFJochException.h" +#include "../compression/CBFByteOffset.h" namespace minicbf { @@ -225,10 +221,6 @@ void ParseAxisTable(const std::string &text, Header &h) { } } -int16_t ReadI16(const uint8_t *p) { int16_t v; std::memcpy(&v, p, 2); return v; } -int32_t ReadI32(const uint8_t *p) { int32_t v; std::memcpy(&v, p, 4); return v; } -int64_t ReadI64(const uint8_t *p) { int64_t v; std::memcpy(&v, p, 8); return v; } - } // namespace bool ScansPhi(const Header &h) { @@ -306,76 +298,6 @@ Header ParseHeader(const char *data, size_t size) { namespace { -// Nothing on this side reads the decoded image back: it goes to the GPU by DMA, or a CPU engine -// walks it once and drops it. An ordinary store therefore pays twice over - it reads every cache -// line before overwriting it, and it evicts 72 MB of live cache to make room - so the pixels go -// straight past the cache. Measured on an 18 Mpx frame: 27 ms of ordinary stores against 11 ms of -// these. The bytes written are the same either way; the fence is what makes them visible to the -// copy that follows. -void StorePixel(int32_t *out, int32_t value) { -#ifdef MINICBF_STREAM_STORE - _mm_stream_si32(out, value); -#else - *out = value; -#endif -} - -void StoreFence() { -#ifdef MINICBF_STREAM_STORE - _mm_sfence(); -#endif -} - -} // namespace - -// The x-CBF_BYTE_OFFSET scheme: a running value, each pixel stored as a delta in the smallest -// container that holds it, escaping to the next size with that container's most negative value. -// Following Bernstein & Hammersley (2006) Int. Tables Cryst. G, 37-43 -void DecodeByteOffset(const uint8_t *data, size_t size, int32_t *out, size_t n_pixels) { - int64_t value = 0; - size_t pos = 0; - size_t written = 0; - - while (written < n_pixels) { - if (pos >= size) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "miniCBF byte-offset stream ended after " + std::to_string(written) - + " of " + std::to_string(n_pixels) + " pixels"); - - int64_t delta = static_cast(data[pos]); - pos += 1; - - if (delta == -128) { - if (pos + 2 > size) break; - delta = ReadI16(data + pos); - pos += 2; - if (delta == -32768) { - if (pos + 4 > size) break; - delta = ReadI32(data + pos); - pos += 4; - if (delta == INT32_MIN) { - if (pos + 8 > size) break; - delta = ReadI64(data + pos); - pos += 8; - } - } - } - - value += delta; - StorePixel(out + written, static_cast(value)); - written++; - } - - StoreFence(); - - if (written != n_pixels) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "miniCBF byte-offset stream ended after " + std::to_string(written) - + " of " + std::to_string(n_pixels) + " pixels"); -} - -namespace { - // Whether the file is gzip-wrapped, from its first two bytes rather than its name - a sweep is // named consistently but the magic is what decides how to read it. bool IsGzip(const std::string &path) { @@ -392,7 +314,8 @@ bool IsGzip(const std::string &path) { // up front the way the plain path can. It grows only past what it already holds: a buffer reused // frame to frame is the previous frame's size, and every resize above that zero-fills the bytes // it adds - tens of megabytes a frame of a 16 Mpx sweep, all of them overwritten by the read. -void SlurpGz(const std::string &path, size_t max_bytes, std::vector &buf) { +template +void SlurpGz(const std::string &path, size_t max_bytes, Buffer &buf) { gzFile g = gzopen(path.c_str(), "rb"); if (g == nullptr) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, @@ -419,7 +342,8 @@ void SlurpGz(const std::string &path, size_t max_bytes, std::vector &bu gzclose(g); } -void Slurp(const std::string &path, size_t max_bytes, std::vector &buf) { +template +void Slurp(const std::string &path, size_t max_bytes, Buffer &buf) { // The plain path is kept exactly as it was rather than reading everything through zlib: a // gzFile over an uncompressed file works, but it copies every byte through zlib's own buffer, // and the overwhelming majority of sweeps are not compressed. @@ -466,7 +390,8 @@ Header ReadHeader(const std::string &path) { namespace { // One read of the file, header parsed, dimensions checked. The caller supplies where the pixels go. -Header ReadCommon(const std::string &path, std::vector &buf, size_t &binary_start) { +template +Header ReadCommon(const std::string &path, Buffer &buf, size_t &binary_start) { Slurp(path, std::numeric_limits::max(), buf); const auto start = FindBinarySection(buf.data(), buf.size()); if (!start.has_value()) @@ -502,18 +427,22 @@ Header Read(const std::string &path, std::vector &out) { size_t start = 0; const Header h = ReadCommon(path, buf, start); out.resize(static_cast(h.nelem)); - DecodeByteOffset(buf.data() + start, buf.size() - start, out.data(), + DecodeCBFByteOffset(buf.data() + start, buf.size() - start, out.data(), static_cast(h.nelem)); return h; } +Header ReadCompressed(const std::string &path, RawByteBuffer &buf, size_t &binary_start) { + return ReadCommon(path, buf, binary_start); +} + Header ReadInto(const std::string &path, int32_t *out, size_t capacity, std::vector &scratch) { size_t start = 0; const Header h = ReadCommon(path, scratch, start); if (static_cast(h.nelem) > capacity) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "CBF image does not fit the supplied buffer: " + path); - DecodeByteOffset(scratch.data() + start, scratch.size() - start, out, static_cast(h.nelem)); + DecodeCBFByteOffset(scratch.data() + start, scratch.size() - start, out, static_cast(h.nelem)); return h; } diff --git a/reader/MiniCBF.h b/reader/MiniCBF.h index 9cb3fd96d..6b33be351 100644 --- a/reader/MiniCBF.h +++ b/reader/MiniCBF.h @@ -9,6 +9,8 @@ #include #include +#include "../common/CompressedImage.h" + // PILATUS miniCBF: an ASCII header, four separator bytes, then one byte-offset compressed image. // No CIF parser and no libcbf are needed - every value is on a "# " comment line or a MIME line. namespace minicbf { @@ -80,11 +82,6 @@ std::optional FindBinarySection(const uint8_t *data, size_t size); // Diamond I24 set carries 6335 bytes, well past a 4 kB guess), so never parse a fixed prefix. Header ParseHeader(const char *data, size_t size); -// x-CBF_BYTE_OFFSET -> int32. Deltas against a running value, smallest container first: int8, -// escaping to int16 via -128, to int32 via -32768, to int64 via INT32_MIN. Little-endian, packed. -// Throws if the stream ends before n_pixels are produced. out must hold n_pixels. -void DecodeByteOffset(const uint8_t *data, size_t size, int32_t *out, size_t n_pixels); - // Header + pixels of one file, read from disk. Header Read(const std::string &path, std::vector &out); @@ -94,6 +91,10 @@ Header Read(const std::string &path, std::vector &out); // Throws if the image does not fit in capacity pixels. Header ReadInto(const std::string &path, int32_t *out, size_t capacity, std::vector &scratch); +// The whole file read into buf, header parsed and checked, without decoding the pixels: the byte-offset +// stream starts at binary_start. For a consumer that decodes it itself - on the device, where there is one. +Header ReadCompressed(const std::string &path, RawByteBuffer &buf, size_t &binary_start); + // Header only - reads just enough of the file to reach the separator. Cheap enough to call per // frame for the rotation angles. Header ReadHeader(const std::string &path); diff --git a/rugnux/RugnuxPipeline.h b/rugnux/RugnuxPipeline.h index 353c07b4c..73fe2158d 100644 --- a/rugnux/RugnuxPipeline.h +++ b/rugnux/RugnuxPipeline.h @@ -106,6 +106,14 @@ namespace rugnux_internal { } space_.notify_one(); try { + // A frame its reader handed over as the file stores it - miniCBF - is written as the + // int32 image it decodes to: that encoding never leaves the process. + if (job.msg.image.GetCompressionAlgorithm() == CompressionAlgorithm::CBF_BYTE_OFFSET) { + const CompressedImage stored = job.msg.image; + stored.GetUncompressed(decoded_); + job.msg.image = CompressedImage(decoded_.data(), decoded_.size(), stored.GetWidth(), + stored.GetHeight(), stored.GetMode()); + } writer_.Write(job.msg); } catch (...) { std::lock_guard lock(m_); @@ -124,6 +132,7 @@ namespace rugnux_internal { std::condition_variable work_; std::condition_variable space_; std::deque queue_; + std::vector decoded_; bool done_ = false; bool failed_ = false; std::exception_ptr error_; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1ff1356f4..c88e132c7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -127,6 +127,7 @@ ADD_EXECUTABLE(jfjoch_test BSLZ4DecoderGPUTest.cpp BSLZ4DecoderGPUFuzzTest.cpp ImagePreprocessorGPUFusedTest.cpp + StoredImageDecodeGPUTest.cpp AdaptiveSpotFinderGPUTest.cpp SpotExtractorGPUParityTest.cpp CUDAMemHelpersTest.cpp diff --git a/tests/StoredImageDecodeGPUTest.cpp b/tests/StoredImageDecodeGPUTest.cpp new file mode 100644 index 000000000..6e133845c --- /dev/null +++ b/tests/StoredImageDecodeGPUTest.cpp @@ -0,0 +1,181 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "../common/CUDAWrapper.h" + +#ifdef JFJOCH_USE_CUDA + +#include +#include +#include + +#include "../common/PixelMask.h" +#include "../compression/CBFByteOffset.h" +#include "../image_analysis/image_preprocessing/CBFDecoderGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" + +// The miniCBF frames reach the GPU as their files store them, and the device has to turn them into +// exactly the int32 image the host decoder would have: the same pixels, the same refusal of a stream +// that ends early. +namespace { + +// x-CBF_BYTE_OFFSET with the smallest container for each delta; force64 puts a delta in the 15-byte +// form whatever its size, which no real writer does but the format allows. +void AppendDelta(std::vector &out, int64_t delta, bool force64) { + auto put = [&out](int64_t v, int bytes) { + for (int i = 0; i < bytes; i++) out.push_back(static_cast(v >> (8 * i))); + }; + if (!force64 && delta > -128 && delta <= 127) { put(delta, 1); return; } + put(-128, 1); + if (!force64 && delta > -32768 && delta <= 32767) { put(delta, 2); return; } + put(-32768, 2); + if (!force64 && delta > INT32_MIN && delta <= INT32_MAX) { put(delta, 4); return; } + put(INT32_MIN, 4); + put(delta, 8); +} + +std::vector EncodeByteOffset(const std::vector &pixels, std::mt19937 &rng) { + std::vector out; + int64_t prev = 0; + for (const int32_t p : pixels) { + AppendDelta(out, p - prev, rng() % 5000 == 0); + prev = p; + } + return out; +} + +// Mostly small deltas, with runs of every escape: int16 and int32 deltas whose payload bytes include +// 0x80 (the escape byte, which must not be taken for one inside a payload), and negative pixels. +std::vector MakePixels(size_t n, uint32_t seed) { + std::mt19937 rng(seed); + std::vector p(n); + for (size_t i = 0; i < n; i++) { + const uint32_t r = rng() % 100; + if (r < 80) p[i] = static_cast(rng() % 40); + else if (r < 90) p[i] = 128 + static_cast(rng() % 1000) * 256; // payload bytes 0x80 + else if (r < 95) p[i] = static_cast(rng() % 2000000) - 1000000; + else if (r < 98) p[i] = -1; + else p[i] = static_cast(rng()); + } + // A long run of escapes back to back, so whole tiles hold nothing but multi-byte tokens. + for (size_t i = n / 3; i < n / 3 + 5000 && i < n; i++) + p[i] = (i % 2) ? 1000000 : -1000000; + return p; +} + +std::vector DecodeOnDevice(const std::vector &stream_bytes, size_t nx, size_t ny, + CBFDecoderGPU &decoder, CudaStream &stream) { + const CompressedImage image(stream_bytes.data(), stream_bytes.size(), nx, ny, CompressedImageMode::Int32, + CompressionAlgorithm::CBF_BYTE_OFFSET); + CudaDevicePtr out(nx * ny); + decoder.Decode(image, out); + std::vector ret(nx * ny); + REQUIRE(cudaMemcpyAsync(ret.data(), out.get(), ret.size() * sizeof(int32_t), cudaMemcpyDeviceToHost, stream) + == cudaSuccess); + REQUIRE(cudaStreamSynchronize(stream) == cudaSuccess); + decoder.ThrowIfDecodeFailed(); + return ret; +} + +bool SameStats(const ImageStatistics &a, const ImageStatistics &b) { + return a.max_value == b.max_value && a.min_value == b.min_value + && a.masked_pixel_count == b.masked_pixel_count + && a.error_pixel_count == b.error_pixel_count + && a.saturated_pixel_count == b.saturated_pixel_count; +} + +} // namespace + +TEST_CASE("CBFDecoderGPU_MatchesHost", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + auto stream = std::make_shared(); + CBFDecoderGPU decoder(stream); + std::mt19937 rng(11); + for (const auto &[nx, ny] : std::vector>{{1, 1}, {7, 3}, {1000, 1013}, {2463, 1000}}) { + const auto pixels = MakePixels(nx * ny, static_cast(nx * 31 + ny)); + auto bytes = EncodeByteOffset(pixels, rng); + // A real file carries the CIF trailer after the stream; nothing past the last pixel is read. + for (const char c : std::string("\r\n--CIF-BINARY-FORMAT-SECTION----\r\n;\r\n\x80\x80\x80")) + bytes.push_back(static_cast(c)); + + std::vector host(nx * ny); + DecodeCBFByteOffset(bytes.data(), bytes.size(), host.data(), host.size()); + REQUIRE(host == pixels); + INFO(nx << " x " << ny); + CHECK(DecodeOnDevice(bytes, nx, ny, decoder, *stream) == host); + } +} + +TEST_CASE("CBFDecoderGPU_ShortStreamIsRefused", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + auto stream = std::make_shared(); + CBFDecoderGPU decoder(stream); + std::mt19937 rng(5); + const size_t nx = 300, ny = 200; + auto pixels = MakePixels(nx * ny, 3); + pixels.back() = pixels[pixels.size() - 2] + 1000000; // the last delta takes the 7-byte form + const auto bytes = EncodeByteOffset(pixels, rng); + + // One token short, and cut inside the last (multi-byte) token: both decoders refuse both. + const std::vector short_stream(bytes.begin(), bytes.end() - 2); + for (const auto &s : {std::vector(bytes.begin(), bytes.end() - 7), short_stream}) { + std::vector host(nx * ny); + CHECK_THROWS(DecodeCBFByteOffset(s.data(), s.size(), host.data(), host.size())); + CHECK_THROWS(DecodeOnDevice(s, nx, ny, decoder, *stream)); + } +} + +// Through the preprocessor, as the image loop calls it: the stored encoding against the CPU +// preprocessor fed what the host decodes it to. +TEST_CASE("ImagePreprocessorGPU_StoredEncodingsMatchCPU", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + DetectorSetup detector = DetJF4M(); + detector.SaturationLimit(65535); + DiffractionExperiment x(detector); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + const size_t npixels = x.GetPixelsNum(); + PixelMask mask(x); + std::vector user_mask(mask.GetMask().size(), 0); + for (size_t i = 0; i < npixels; i += 997) user_mask[i] = 1; + mask.LoadUserMask(x, user_mask); + + std::mt19937 rng(17); + const auto pixels = MakePixels(npixels, 9); + const auto cbf = EncodeByteOffset(pixels, rng); + + const std::vector images = { + CompressedImage(cbf.data(), cbf.size(), x.GetXPixelsNum(), x.GetYPixelsNum(), CompressedImageMode::Int32, + CompressionAlgorithm::CBF_BYTE_OFFSET) + }; + auto stream = std::make_shared(); + ImagePreprocessorGPU gpu_pre(x, mask, stream, /*copy_image_to_host=*/true); + ImagePreprocessorCPU cpu_pre(x, mask); + for (const auto &image : images) { + std::vector decoded; + const uint8_t *raw = image.GetUncompressedPtr(decoded); + ImagePreprocessorBuffer cpu_buf(npixels); + const ImageStatistics cpu_stats = cpu_pre.Analyze(cpu_buf, raw, CompressedImageMode::Int32); + + ImagePreprocessorBufferGPU gpu_buf(npixels); + ImageStatistics gpu_stats{}; + REQUIRE(gpu_pre.AnalyzeCompressed(gpu_buf, image, gpu_stats)); + INFO("algorithm " << static_cast(image.GetCompressionAlgorithm())); + CHECK(SameStats(cpu_stats, gpu_stats)); + CHECK(cpu_stats.saturated_pixel_count > 0); + size_t ndiff = 0; + for (size_t i = 0; i < npixels; i++) + ndiff += gpu_buf[i] != cpu_buf[i]; + CHECK(ndiff == 0); + } +} + +#endif