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