From bec7e2e922e52d6e88bbbe169461d91b2db1f800 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 3 Aug 2026 14:13:55 +0200 Subject: [PATCH] image_preprocessing: fuse the bitshuffle inverse with preprocessing, and verify the decode The device decoder was byte-exact on every valid input - 994 production-compressed images, 927 hand-built LZ4 blocks covering engineered (offset, matchlen) pairs across the overlap branch boundary, 18000 repeat decodes, sanitizer-clean - and an audit against LZ4_decompress_generic could not construct a valid block it mis-decodes. What it did not do was notice when the input was NOT valid, and that mattered more than it looks: the decode buffers are reused frame to frame, so a block that stopped early left the PREVIOUS image in place, and in the bitshuffled layout the untouched tail is the most significant byte-plane. A corrupt chunk therefore did not look like a missing corner. It looked like thousands of real pixels several powers of two too bright, fed to spot finding with no diagnostic, where the host decoder had raised an error. So the kernel now flags a block that fails to reach its declared length while consuming exactly its payload, and the host turns that into an exception once the caller has synchronised. Reads are clamped against the end of the payload as well as the output, both length chains are bounded exactly as read_variable_length bounds them, the two offset bytes are bounded, and LZ4's parsing restrictions are enforced. On the host side a block size that is not a multiple of 8 elements is rejected (it made the un-transpose read uninitialised shared memory), the block count is bounded by what the chunk could hold before it becomes an allocation (twelve header bytes could demand hundreds of MB of pinned memory, permanently, per worker), trailing bytes are rejected, and the stream is synchronised before any throw that happens after work is queued. An image of fewer than 8 elements is all verbatim tail and now decodes rather than throwing. When the device route fails for any reason the host decoder gets its turn, so it costs speed rather than the acquisition. The lanes cooperate on the copies and a later match can read bytes another lane wrote, which since Volta needs an explicit __syncwarp(); it worked only because ptxas happened to reconverge at the post-dominator. The prototype's offset == 1 and power-of-two fast paths are also restored - the shipped kernel ran a runtime modulo, an emulated 32-bit division per output byte, on the path its own comment calls the common case. The un-transpose is now fused with preprocessing. One thread owns one group of 8 elements across every byte-plane, so once it has transposed its 8 bytes out of each plane it holds 8 complete elements and emits 8 finished int32 pixels with the mask, the error marker, the saturation cap and the statistics applied. The decompressed image is never materialised: 0.623 -> 0.411 ms/frame at 18 Mpx, 0.523 -> 0.340 with 8 concurrent workers. Staging nothing in shared memory also drops the 48 kB ceiling, which had made any file whose bitshuffle blocks exceed it a hard failure; 64 kB blocks now decode. gpu_compressed is sized from the chunk with grow-on-demand instead of from the uncompressed size - it was reserving ~73 MB per worker to hold ~4 MB. Measured on a 1630x1553 uint32 rotation set at -N 32, peak GPU memory falls 3756 -> 3084 MiB; the same model gives ~144 MB per worker on an 18 Mpx frame. Decoding on the device also stopped reporting a decompression time, which blanked the broker's compression plot trace and filled /entry/profiling/compressionTime with NaN. The decoder brackets the decode with CUDA events and reports it again. Tests: a differential fuzz suite against the CPU decoder - incompressible and highly compressible data, engineered offsets, a size sweep hitting every rem%8 value twice, all six element sizes, an 18 Mpx frame, decoder reuse, concurrency, hand-built LZ4 blocks across the overlap boundary, 26 foreign bitshuffle block sizes from 128 B to 64 kB, corrupt payloads and malformed containers, with a coverage report that proves which LZ4 paths were reached rather than assuming it. Plus the fused path held byte for byte against ImagePreprocessorCPU, statistics included, and against the host-upload path on the same frame. Battery: 37 crystals, every merged number identical to the host-decode run. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CPU_DATA_ANALYSIS.md | 33 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 27 +- image_analysis/MXAnalysisWithoutFPGA.h | 4 +- .../image_preprocessing/BSLZ4DecoderGPU.cu | 281 +++- .../image_preprocessing/BSLZ4DecoderGPU.h | 53 +- .../image_preprocessing/ImagePreprocessor.h | 5 + .../ImagePreprocessorGPU.cu | 199 ++- .../ImagePreprocessorGPU.h | 11 +- image_analysis/indexing/CUDAMemHelpers.h | 29 + tests/BSLZ4DecoderGPUFuzzTest.cpp | 1304 +++++++++++++++++ tests/CMakeLists.txt | 2 + tests/ImagePreprocessorGPUFusedTest.cpp | 167 +++ 12 files changed, 2023 insertions(+), 92 deletions(-) create mode 100644 tests/BSLZ4DecoderGPUFuzzTest.cpp create mode 100644 tests/ImagePreprocessorGPUFusedTest.cpp diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 9e249a12..10e27f5a 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -40,16 +40,25 @@ the GPU when one is present. Instead of decompressing on the host and uploading compressed chunk is uploaded — a few MB rather than tens of MB — and decoded on the device. The approach follows Jon Wright (ESRF); the kernels are Jungfraujoch's own. -Two kernels mirror the CPU decoder: +Two kernels do the work: 1. **LZ4, one warp per bitshuffle block.** Blocks are independent, so the parallelism is across them; within a warp every lane runs the same sequence parser over the same bytes, and the literal and match copies are split across the 32 lanes so the stores coalesce. An overlapping match is treated as a pattern of period `offset` sourced from bytes that already precede the - write position, which keeps it parallel rather than a serial byte loop. -2. **The bitshuffle inverse**, one CUDA block per bitshuffle block: bit un-transpose of each - byte-plane into shared memory, then interleave the planes back into elements so the final store - is coalesced. For 8-bit images there is a single plane and the interleave degenerates to a copy. + write position, which keeps it parallel rather than a serial byte loop; `offset == 1` (a run of + one repeated byte, the common case in sparse detector data) and power-of-two offsets avoid the + modulo altogether. Because the lanes cooperate on the copies, each one is followed by + `__syncwarp()` — a later match can read bytes another lane wrote, and since Volta that ordering + is not implicit. +2. **The bitshuffle inverse fused with preprocessing.** One thread owns one group of 8 elements + across every byte-plane, so once it has transposed its 8 bytes out of each plane it holds 8 + complete elements — and it applies the pixel mask, the error marker and the saturation cap and + emits 8 finished `int32` pixels directly. The decompressed image is therefore never materialised + in device memory at all, which removes a frame-sized buffer per worker and a full-frame write + plus read from the pipeline. Staging nothing in shared memory also means the kernel has no + dynamic-shared-memory request, so it is indifferent to the bitshuffle block size the file + declares. For 8-bit images there is a single plane and the assembly degenerates to a copy. The block offsets inside the container can only be discovered by reading the block lengths in order, so that scan stays on the host. @@ -58,6 +67,20 @@ Only `BSHUF_LZ4` is decoded on the device. For the zstd variants (`BSHUF_ZSTD`, `BSHUF_ZSTD_RLE_HUFF`), and for uncompressed or float images, `BSLZ4DecoderGPU::Supports()` returns false and the pipeline decompresses on the host and uploads as before. +The container arrives off the network or off disk and is not trusted. Everything checkable on the +host — declared sizes, the block scan, a block size that is not a multiple of 8 elements, a block +count the chunk could not hold, trailing bytes — is rejected before any work is queued. What only +the kernel can see is that a block failed to decode to exactly its declared length while consuming +exactly its payload; that raises a device-side flag which becomes an exception once the caller has +synchronised. This matters because the decode buffers are reused frame to frame: a block that +stopped early would otherwise leave the *previous* image's bytes in place, and in the bitshuffled +layout those are the most significant byte-plane, so the result would not look like a missing +corner but like real pixels several powers of two too bright. The kernel reproduces the reference +decoder's length, bounds and parsing-restriction checks; it does not reproduce its exact +fast-loop/safe-loop selection, so a small tail of corrupt streams still decodes here that +`LZ4_decompress_safe` would refuse — as a *complete* decode differing in a few bytes, never as a +partial one. + ## 1. Geometry, reciprocal-space mapping, and basic quantities ### 1.1 Coordinate conventions diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index a21626dc..385e7da8 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -4,6 +4,7 @@ #include "MXAnalysisWithoutFPGA.h" #include +#include #include "spot_finding/StrongPixelSet.h" #include "../compression/JFJochDecompress.h" @@ -61,7 +62,8 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp preprocessor_buffer = std::make_unique(experiment.GetPixelsNum()); // The preprocessed image only has to come back to the host if a CPU engine reads it. Every // engine built below runs on the GPU, except the CPU adaptive finder that is kept when the fused - // GPU engine is off (the online receiver) - so that is the one case that needs the copy. + // GPU engine is off - so that is the one case that needs the copy. Every caller currently passes + // enable_fused_adaptive_gpu = true, so on the GPU path the copy is off in practice. preprocessor = std::make_unique(in_experiment, in_mask, stream, /*copy_image_to_host=*/!enable_fused_adaptive_gpu); spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream); @@ -96,7 +98,17 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // the host route unchanged. The two produce the same preprocessed image. const auto compression_start_time = std::chrono::steady_clock::now(); ImageStatistics ret{}; - const bool decoded_on_device = preprocessor->AnalyzeCompressed(*preprocessor_buffer, output.image, ret); + bool decoded_on_device = false; + try { + decoded_on_device = preprocessor->AnalyzeCompressed(*preprocessor_buffer, output.image, ret); + } catch (const JFJochException &e) { + // The device route must never be the reason a frame fails: whatever it could not handle, the + // host decoder gets its turn. If the data really is bad the host throws too and the caller + // sees the same error it saw before any of this existed - but a GPU-side problem costs speed + // rather than the acquisition. + spdlog::warn("Device decoding failed ({}), falling back to host decompression", e.what()); + decoded_on_device = false; + } const auto compression_end_time = std::chrono::steady_clock::now(); if (!decoded_on_device) { @@ -110,9 +122,14 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const auto preprocessing_end_time = std::chrono::steady_clock::now(); output.preprocessing_time_s = std::chrono::duration(preprocessing_end_time - preprocessing_start_time).count(); } else { - // Decode and preprocess are one device operation here, so they are reported together rather - // than split into a decompression time that no longer exists on the host. - output.preprocessing_time_s = std::chrono::duration(compression_end_time - compression_start_time).count(); + // Decode and preprocess are one device operation here, but the decompression is still a real, + // separately measurable cost - the decoder brackets it with CUDA events - so it is still + // reported as one. Leaving compression_time_s unset instead would blank the broker's + // "compression" plot trace and fill /entry/profiling/compressionTime with NaN. + const float total_s = std::chrono::duration(compression_end_time - compression_start_time).count(); + const float decompress_s = std::min(preprocessor->GetLastDecompressionTime_s(), total_s); + output.compression_time_s = decompress_s; + output.preprocessing_time_s = total_s - decompress_s; } // The fused GPU engine (rugnux offline, GPU, adaptive detection) produces the azimuthal profile as diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index a4134cd2..bb273509 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -71,8 +71,8 @@ public: // the GPU path with adaptive detection). The rugnux offline path and the interactive viewer enable // it by default, as does the online receiver. It only changes performance - the fused engine // reproduces the CPU finder's spots. Note it also decides whether the preprocessed image is copied - // back to the host each frame: that copy exists only for a CPU engine to read, and while adaptive - // detection is not reachable through the REST API the copy is the flag's only effect online. + // back to the host each frame: that copy exists only for a CPU engine to read, and with the flag on + // no CPU engine is built, so the copy is skipped. MXAnalysisWithoutFPGA(const DiffractionExperiment &experiment, const AzimuthalIntegrationMapping &integration, const PixelMask &mask, IndexAndRefine &indexer, bool enable_fused_adaptive_gpu = false); void Analyze(DataMessage &output, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &spot_finding_settings); diff --git a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu index 829d76c0..1e4ddab9 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu @@ -17,9 +17,20 @@ namespace { // broadcast read, so no divergence - and the literal and match copies are split across the 32 // lanes so the stores coalesce. One thread per block instead has each thread streaming its own // 8 kB region, which coalesces not at all and measured 13x slower. + // + // Because the lanes cooperate on the copies, a match can source bytes that OTHER lanes wrote in + // an earlier sequence. Since Volta that needs an explicit __syncwarp() - implicit reconvergence + // is not part of the programming model - so there is one after every copy loop. The full mask is + // correct: the early return and every break test warp-uniform values, so lanes never diverge + // permanently. + // + // Bounds: every read is clamped against iend and every write against oend, so a malformed or + // corrupt payload cannot walk off either buffer. It can still stop early, which leaves the block + // short; lane 0 flags that at the end and the host turns it into an exception. __global__ void lz4_decode_blocks(const uint8_t *__restrict__ src, const BSLZ4BlockDesc *__restrict__ desc, uint8_t *__restrict__ dst, + uint32_t *__restrict__ status, int nblocks, uint32_t elem_size) { const int lane = threadIdx.x & 31; const int b = (blockIdx.x * blockDim.x + threadIdx.x) >> 5; @@ -27,34 +38,75 @@ namespace { const uint8_t *ip = src + desc[b].in_off; const uint8_t *const iend = ip + desc[b].in_len; - uint8_t *op = dst + desc[b].out_off; - uint8_t *const oend = op + desc[b].nelem * elem_size; + uint8_t *const obase = dst + desc[b].out_off; + uint8_t *op = obase; + uint8_t *const oend = obase + desc[b].nelem * elem_size; + bool malformed = false; while (ip < iend) { const uint32_t token = *ip++; uint32_t litlen = token >> 4; if (litlen == 15) { + // read_variable_length(&ip, iend - RUN_MASK, initial_check=1) in the reference: the + // chain may not start within, nor run into, the last RUN_MASK (15) input bytes. A + // valid stream never does - the literals it counts have to follow it - so a chain + // that reaches there is corruption, and this is the only place it shows up. + if ((size_t)(iend - ip) <= 15) { malformed = true; break; } uint32_t s; - do { s = *ip++; litlen += s; } while (s == 255 && ip < iend); + do { + s = *ip++; + litlen += s; + if ((size_t)(iend - ip) < 15) { malformed = true; break; } + } while (s == 255); + if (malformed) break; } if (litlen) { - const uint32_t n = min(litlen, (uint32_t)(oend - op)); + // Clamped by the INPUT as well as the output: a corrupt litlen must not read past the + // end of this block's payload or write past the end of the block. Clamping keeps the + // kernel in bounds; needing to clamp at all means the stream is not decodable, which + // is what the reference reports as an error, so record it. + if (litlen > (uint32_t)(oend - op) || litlen > (uint32_t)(iend - ip)) + malformed = true; + const uint32_t n = min(min(litlen, (uint32_t)(oend - op)), (uint32_t)(iend - ip)); for (uint32_t i = lane; i < n; i += 32) op[i] = ip[i]; + __syncwarp(); op += n; ip += litlen; } - if (ip >= iend) break; // last sequence carries literals only + + // LZ4's parsing restrictions: an encoder may not leave a match within MFLIMIT (12) bytes + // of the end of the block, nor fewer than 2+1+LASTLITERALS (8) input bytes after a + // literal run that is not the last one. So once either limit is reached this can ONLY be + // the final sequence, and the final sequence must consume the payload exactly. The + // reference applies this whether or not the run was empty, which is why the test sits + // outside the copy - a zero-length literal run near the end is just as illegal. + if ((size_t)(oend - op) < 12 || (size_t)(iend - ip) < 8) { + malformed = (ip != iend) || (op != oend); + break; // necessarily EOF + } + if (iend - ip < 2) break; // last sequence carries literals only const uint32_t offset = (uint32_t)ip[0] | ((uint32_t)ip[1] << 8); ip += 2; uint32_t matchlen = token & 0x0F; if (matchlen == 15) { + // read_variable_length(&ip, iend - LASTLITERALS + 1, initial_check=0): bounded by the + // last 4 input bytes rather than 15, and with no check before the first read. uint32_t s; - do { s = *ip++; matchlen += s; } while (s == 255 && ip < iend); + do { + s = *ip++; + matchlen += s; + if ((size_t)(iend - ip) < 4) { malformed = true; break; } + } while (s == 255); + if (malformed) break; } matchlen += 4; // minmatch - if (offset == 0 || offset > (uint32_t)(op - (dst + desc[b].out_off))) break; // malformed + if (offset == 0 || offset > (uint32_t)(op - obase)) { malformed = true; break; } const uint8_t *mp = op - offset; + // A match may reach the end of the block but never past it - the reference treats an + // overrun as an error rather than truncating, and so must this. + if (matchlen > (uint32_t)(oend - op)) + malformed = true; const uint32_t n = min(matchlen, (uint32_t)(oend - op)); if (offset >= matchlen) { for (uint32_t i = lane; i < n; i += 32) op[i] = mp[i]; @@ -62,11 +114,28 @@ namespace { // An overlapping match is a pattern of period `offset`. mp[0..offset-1] all lie // before op and are already final, so each output byte can be sourced from them // independently - which keeps this parallel rather than a serial byte loop. Long - // zero runs in sparse detector data arrive here with offset == 1. - for (uint32_t i = lane; i < n; i += 32) op[i] = mp[i % offset]; + // zero runs in sparse detector data arrive here with offset == 1, and a runtime + // modulo is an emulated division, so the two cheap cases are peeled off first. + if (offset == 1) { + const uint8_t v = mp[0]; + for (uint32_t i = lane; i < n; i += 32) op[i] = v; + } else if ((offset & (offset - 1)) == 0) { + const uint32_t m = offset - 1; + for (uint32_t i = lane; i < n; i += 32) op[i] = mp[i & m]; + } else { + for (uint32_t i = lane; i < n; i += 32) op[i] = mp[i % offset]; + } } + __syncwarp(); op += n; } + + // A block must decode to exactly its declared length AND consume exactly its payload. Both + // are conditions LZ4_decompress_safe reports to the host path, and both are needed: a corrupt + // stream can land on the right output length while leaving input over, or run its input out + // early. Either way the bytes are not the ones that were compressed. + if (lane == 0 && (malformed || op != oend || ip != iend)) + atomicExch(status, 1u); } __device__ __forceinline__ uint64_t transpose8(uint64_t x) { @@ -77,14 +146,15 @@ namespace { return x; } - // The bitshuffle inverse, one CUDA block per bitshuffle block, mirroring bitshuf_decode_block: - // un-transpose the bits of each byte-plane, then interleave the planes back into elements. The - // planes are staged in shared memory so the final store is coalesced. - extern __shared__ uint8_t smem[]; + // The bitshuffle inverse, mirroring bitshuf_decode_block. One thread owns one group of 8 elements + // across EVERY byte-plane, so after transposing its 8 bytes out of each plane it holds all bytes + // of 8 complete elements and can write them straight out. That needs no staging buffer, which is + // what keeps the kernel free of the 48 kB dynamic-shared-memory ceiling a block size taken from + // the file header would otherwise run into. + template __global__ void bitshuffle_untranspose(const uint8_t *__restrict__ shuffled, const BSLZ4BlockDesc *__restrict__ desc, - uint8_t *__restrict__ out, - int nblocks, uint32_t elem_size) { + uint8_t *__restrict__ out, int nblocks) { const int b = blockIdx.x; if (b >= nblocks) return; @@ -93,23 +163,22 @@ namespace { uint8_t *dst = out + desc[b].out_off; const uint32_t n = size / 8; - for (uint32_t p = 0; p < elem_size; p++) { - const uint8_t *pin = in + p * size; - uint8_t *pout = smem + p * size; - for (uint32_t i = threadIdx.x; i < n; i += blockDim.x) { + for (uint32_t i = threadIdx.x; i < n; i += blockDim.x) { + uint64_t x[ES]; + #pragma unroll + for (int p = 0; p < ES; p++) { + const uint8_t *pin = in + p * size; uint64_t a = 0; #pragma unroll for (int k = 0; k < 8; k++) a |= (uint64_t)pin[k * n + i] << (8 * k); - const uint64_t x = transpose8(a); - #pragma unroll - for (int k = 0; k < 8; k++) pout[i * 8 + k] = (uint8_t)(x >> (8 * k)); + x[p] = transpose8(a); } + #pragma unroll + for (int k = 0; k < 8; k++) + #pragma unroll + for (int p = 0; p < ES; p++) + dst[(i * 8 + k) * ES + p] = (uint8_t)(x[p] >> (8 * k)); } - __syncthreads(); - - for (uint32_t i = threadIdx.x; i < size; i += blockDim.x) - for (uint32_t j = 0; j < elem_size; j++) - dst[i * elem_size + j] = smem[j * size + i]; } uint64_t be64(const uint8_t *p) { uint64_t v = 0; for (int i = 0; i < 8; i++) v = (v << 8) | p[i]; return v; } @@ -129,6 +198,11 @@ namespace { default: return 0; // float and RGB modes are never bitshuffled by this pipeline } } + + // The smallest a block can be on the wire: a 4-byte length plus at least one payload byte. Used + // to reject a header whose declared block size implies more blocks than the chunk could hold, + // before that count is turned into an allocation. + constexpr size_t MIN_BLOCK_BYTES_ON_WIRE = 5; } bool BSLZ4DecoderGPU::Supports(const CompressedImage &image) { @@ -139,21 +213,35 @@ bool BSLZ4DecoderGPU::Supports(const CompressedImage &image) { BSLZ4DecoderGPU::BSLZ4DecoderGPU(size_t in_max_uncompressed_bytes, std::shared_ptr in_stream) : stream(std::move(in_stream)), max_uncompressed_bytes(in_max_uncompressed_bytes) { - // The compressed chunk is smaller than the image in every case worth having, but a pathological - // frame can expand slightly, so allow headroom rather than risk a per-frame reallocation. - max_compressed_bytes = max_uncompressed_bytes + max_uncompressed_bytes / 16 + 4096; - // Enough descriptors for the block size the compressor actually uses (8 kB); Decode() grows - // these if a stream turns up with smaller blocks. Sizing for the format's theoretical minimum - // block instead would allocate megabytes of pinned memory that no real file needs. - max_blocks = max_uncompressed_bytes / 8192 + 2; - - gpu_compressed = CudaDevicePtr(max_compressed_bytes); gpu_shuffled = CudaDevicePtr(max_uncompressed_bytes); - gpu_desc = CudaDevicePtr(max_blocks); - host_desc = CudaHostPtr(max_blocks); + gpu_status = CudaDevicePtr(1); + host_status = CudaHostPtr(1); + // The compressed buffer and the descriptors are grown to fit the first image instead of being + // sized for a worst case that no real frame reaches. A chunk is a few MB against an image of + // tens; sizing this from the UNCOMPRESSED size cost ~73 MB per worker to hold ~4 MB. } -void BSLZ4DecoderGPU::Decode(const CompressedImage &image, uint8_t *gpu_out) { +void BSLZ4DecoderGPU::EnsureCompressedCapacity(size_t bytes) { + if (bytes <= compressed_capacity) + return; + // A little slack, so a frame that compresses slightly worse than the last does not reallocate. + const size_t want = bytes + bytes / 4; + cuda_err(cudaStreamSynchronize(*stream)); // nothing may still be reading the old buffer + gpu_compressed = CudaDevicePtr(want); + compressed_capacity = want; +} + +void BSLZ4DecoderGPU::EnsureBlockCapacity(size_t nblocks) { + if (nblocks <= max_blocks) + return; + const size_t want = nblocks + nblocks / 4 + 16; + cuda_err(cudaStreamSynchronize(*stream)); // the previous descriptor upload must have landed + gpu_desc = CudaDevicePtr(want); + host_desc = CudaHostPtr(want); + max_blocks = want; +} + +BSLZ4ShuffledImage BSLZ4DecoderGPU::DecodeShuffled(const CompressedImage &image) { const uint8_t *src = image.GetCompressed(); const size_t clen = image.GetCompressedSize(); const size_t elem_size = elem_size_of(image.GetMode()); @@ -161,7 +249,7 @@ void BSLZ4DecoderGPU::Decode(const CompressedImage &image, uint8_t *gpu_out) { if (clen < 12) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk shorter than its header"); - if (total_bytes > max_uncompressed_bytes || clen > max_compressed_bytes) + if (total_bytes > max_uncompressed_bytes) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 image larger than the decoder was sized for"); if (be64(src) != total_bytes) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 header size does not match the image"); @@ -171,6 +259,12 @@ void BSLZ4DecoderGPU::Decode(const CompressedImage &image, uint8_t *gpu_out) { throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 block size invalid"); const size_t block_elems = block_bytes / elem_size; + // bitshuffle transposes 8 elements at a time and refuses a block that is not a multiple of 8; + // the host decoder rejects this too (JFJochDecompress.h). Without the check the un-transpose + // would silently drop the last size % 8 elements of every block. + if (block_elems % BSHUF_BLOCKED_MULT != 0) + throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 block size is not a multiple of 8 elements"); + const size_t nelements = total_bytes / elem_size; const size_t nfull = nelements / block_elems; const size_t rem = nelements - nfull * block_elems; @@ -179,49 +273,106 @@ void BSLZ4DecoderGPU::Decode(const CompressedImage &image, uint8_t *gpu_out) { // Walk the container to locate the blocks. Lengths are only knowable in order, so this scan is // inherent to the format rather than an implementation choice. + // An image of fewer than 8 elements has no bitshuffle block at all - it is entirely the verbatim + // tail. The host decoder handles that, so handle it here rather than declining: nblocks is simply + // zero and only the tail is copied. const size_t nblocks_needed = nfull + (last > 0 ? 1 : 0); - if (nblocks_needed > max_blocks) { // a stream with smaller blocks than the compressor emits - max_blocks = nblocks_needed; - gpu_desc = CudaDevicePtr(max_blocks); - host_desc = CudaHostPtr(max_blocks); - } + // Bound the block count by what the chunk could actually hold BEFORE it becomes an allocation: + // a header declaring a one-element block size would otherwise ask for hundreds of MB of pinned + // memory, and only then fail on the first block header. + if (nblocks_needed > (clen - 12) / MIN_BLOCK_BYTES_ON_WIRE) + throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk too short for the blocks its header implies"); + + EnsureCompressedCapacity(clen); + EnsureBlockCapacity(nblocks_needed); size_t nblk = 0, off = 12, out_off = 0; - for (size_t i = 0; i < nfull + (last > 0 ? 1 : 0); i++) { + for (size_t i = 0; i < nblocks_needed; i++) { if (off + 4 > clen) throw JFJochException(JFJochExceptionCategory::Compression, "truncated bslz4 block header"); const uint32_t block_clen = be32(src + off); off += 4; if (block_clen == 0 || off + block_clen > clen) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 block extends past the chunk"); - if (nblk >= max_blocks) - throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk has more blocks than expected"); const uint32_t ne = (i < nfull) ? static_cast(block_elems) : static_cast(last); host_desc.get()[nblk++] = {static_cast(off), block_clen, static_cast(out_off), ne}; off += block_clen; out_off += static_cast(ne) * elem_size; } - if (nblk == 0) - throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk contains no blocks"); + // The tail that bitshuffle leaves uncompressed and copies verbatim, and then nothing else: the + // host decoder requires the chunk to be consumed exactly, so require it here too rather than + // ignoring trailing bytes that indicate the container is not what it claims to be. + if (off + leftover_bytes > clen) + throw JFJochException(JFJochExceptionCategory::Compression, "truncated bslz4 leftover bytes"); + if (off + leftover_bytes != clen) + throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk has trailing bytes after the last block"); - cuda_err(cudaMemcpyAsync(gpu_compressed.get(), src, clen, cudaMemcpyHostToDevice, *stream)); - cuda_err(cudaMemcpyAsync(gpu_desc.get(), host_desc.get(), nblk * sizeof(BSLZ4BlockDesc), + // Everything that can be checked on the host has been checked; from here work is queued. + cuda_err(cudaEventRecord(decode_start, *stream)); + host_status.get()[0] = 0; + cuda_err(cudaMemcpyAsync(gpu_status.get(), host_status.get(), sizeof(uint32_t), cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaMemcpyAsync(gpu_compressed.get(), src, clen, cudaMemcpyHostToDevice, *stream)); const int nb = static_cast(nblk); - lz4_decode_blocks<<<(nb * 32 + 255) / 256, 256, 0, *stream>>>( - gpu_compressed.get(), gpu_desc.get(), - gpu_shuffled.get(), nb, static_cast(elem_size)); - bitshuffle_untranspose<<>>( - gpu_shuffled.get(), gpu_desc.get(), - gpu_out, nb, static_cast(elem_size)); - cuda_err(cudaGetLastError()); - - // The tail that bitshuffle leaves uncompressed and copies verbatim. - if (leftover_bytes > 0) { - if (off + leftover_bytes > clen) - throw JFJochException(JFJochExceptionCategory::Compression, "truncated bslz4 leftover bytes"); - cuda_err(cudaMemcpyAsync(gpu_out + out_off, src + off, leftover_bytes, + if (nb > 0) { + cuda_err(cudaMemcpyAsync(gpu_desc.get(), host_desc.get(), nblk * sizeof(BSLZ4BlockDesc), cudaMemcpyHostToDevice, *stream)); + lz4_decode_blocks<<<(nb * 32 + 255) / 256, 256, 0, *stream>>>( + gpu_compressed.get(), gpu_desc.get(), gpu_shuffled.get(), gpu_status.get(), + nb, static_cast(elem_size)); + cuda_err(cudaGetLastError()); } + cuda_err(cudaMemcpyAsync(host_status.get(), gpu_status.get(), sizeof(uint32_t), + cudaMemcpyDeviceToHost, *stream)); + // Stop the clock here rather than after the un-transpose: getting the chunk onto the device and + // LZ4-decoding it is the part that replaced the host decompression, and it is the same work on + // both the raw-bytes path and the fused one, where the un-transpose is inseparable from + // preprocessing and is reported with it. + cuda_err(cudaEventRecord(decode_stop, *stream)); + decode_timed = true; + + BSLZ4ShuffledImage ret; + ret.shuffled = gpu_shuffled.get(); + ret.desc = gpu_desc.get(); + ret.nblocks = nb; + ret.elem_size = static_cast(elem_size); + ret.tail_elems = static_cast(leftover_bytes / elem_size); + ret.tail_elem0 = static_cast(out_off / elem_size); + ret.tail_src = leftover_bytes > 0 ? gpu_compressed.get() + off : nullptr; + return ret; +} + +void BSLZ4DecoderGPU::Decode(const CompressedImage &image, uint8_t *gpu_out) { + const BSLZ4ShuffledImage s = DecodeShuffled(image); + + if (s.nblocks > 0) { // an image of fewer than 8 elements is all tail and has no block + switch (s.elem_size) { + case 1: bitshuffle_untranspose<1><<>>(s.shuffled, s.desc, gpu_out, s.nblocks); break; + case 2: bitshuffle_untranspose<2><<>>(s.shuffled, s.desc, gpu_out, s.nblocks); break; + default: bitshuffle_untranspose<4><<>>(s.shuffled, s.desc, gpu_out, s.nblocks); break; + } + cuda_err(cudaGetLastError()); + } + + // The verbatim tail is already on the device inside the uploaded chunk. + if (s.tail_elems > 0) + cuda_err(cudaMemcpyAsync(gpu_out + static_cast(s.tail_elem0) * s.elem_size, s.tail_src, + static_cast(s.tail_elems) * s.elem_size, + cudaMemcpyDeviceToDevice, *stream)); +} + +void BSLZ4DecoderGPU::ThrowIfDecodeFailed() { + if (host_status.get()[0] != 0) + throw JFJochException(JFJochExceptionCategory::Compression, + "bslz4 block did not decode to its declared length - the compressed data is corrupt"); +} + +float BSLZ4DecoderGPU::GetDecodeTime_s() const { + if (!decode_timed) + return 0.0f; + 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/BSLZ4DecoderGPU.h b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h index d93b4008..bdc81485 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h @@ -16,6 +16,20 @@ struct BSLZ4BlockDesc { uint32_t nelem; // elements in this block (the last one is usually shorter) }; +// What DecodeShuffled() leaves on the device: the LZ4 output, still bitshuffled, plus everything the +// un-transpose needs to finish the image. The tail is the handful of elements bitshuffle stores +// verbatim; it is already on the device inside the uploaded chunk, so it is handed over as a device +// pointer rather than copied again from the host. +struct BSLZ4ShuffledImage { + const uint8_t *shuffled = nullptr; + const BSLZ4BlockDesc *desc = nullptr; + int nblocks = 0; + uint32_t elem_size = 0; + const uint8_t *tail_src = nullptr; + uint32_t tail_elems = 0; + uint32_t tail_elem0 = 0; // index of the first tail element in the image +}; + // Decompress a bitshuffle+LZ4 image ON THE DEVICE, so the compressed bytes are what crosses PCIe. // // The idea - upload the compressed chunk and decode it on the GPU rather than decompressing on the @@ -30,6 +44,14 @@ struct BSLZ4BlockDesc { // // Only BSHUF_LZ4 is handled. The zstd variants have no device decoder, so Supports() returns false // and the caller decompresses on the host exactly as before. +// +// The container comes off the network or off disk, so it is not trusted. Everything the host can +// check cheaply is checked before any work is queued and throws; what only the kernel can see - a +// block that does not decode to its declared length, which is what a corrupt LZ4 payload looks like +// - raises a device-side flag that ThrowIfDecodeFailed() reports once the caller has synchronised. +// The CPU decoder makes exactly the same checks (LZ4_decompress_safe's length check plus the +// consumed-input check in JFJochDecompress.h), so a chunk either decodes identically on both or +// fails on both. It is never silently wrong on one and right on the other. class BSLZ4DecoderGPU { std::shared_ptr stream; @@ -37,19 +59,42 @@ class BSLZ4DecoderGPU { CudaDevicePtr gpu_shuffled; // LZ4 output, still bitshuffled CudaDevicePtr gpu_desc; CudaHostPtr host_desc; // pinned, so the descriptor upload is truly async + CudaDevicePtr gpu_status; // set by the kernel when a block decodes short + CudaHostPtr host_status; + + // Bracket the decode so the time it takes can still be reported as decompression, which is what + // it is. Both are recorded on the decoder's stream and read after the caller synchronises. + CudaEvent decode_start; + CudaEvent decode_stop; + bool decode_timed = false; - size_t max_compressed_bytes = 0; size_t max_uncompressed_bytes = 0; + size_t compressed_capacity = 0; size_t max_blocks = 0; + void EnsureCompressedCapacity(size_t bytes); + void EnsureBlockCapacity(size_t nblocks); + public: BSLZ4DecoderGPU(size_t max_uncompressed_bytes, std::shared_ptr stream); // True when this image can be decoded on the device. Everything else must go the host route. static bool Supports(const CompressedImage &image); - // Decode into gpu_out, which must hold image.GetUncompressedSize() bytes. Work is queued on the - // decoder's stream and the caller synchronises. Throws if the container is malformed - it comes - // off the network or off disk, so it is not trusted. + // Locate the blocks, upload the chunk, and run the LZ4 pass. The result is still bitshuffled - + // the caller finishes it, either with Decode()'s un-transpose or by fusing the un-transpose into + // its own kernel. Work is queued on the decoder's stream and the caller synchronises. + BSLZ4ShuffledImage DecodeShuffled(const CompressedImage &image); + + // Decode into gpu_out, which must hold image.GetUncompressedSize() bytes. The plain raw-bytes + // path: DecodeShuffled() plus the un-transpose. Used by the tests and by any caller that wants + // the decompressed image rather than a preprocessed one. void Decode(const CompressedImage &image, uint8_t *gpu_out); + + // Report a block that did not decode to its declared length. MUST be called after the caller has + // synchronised the stream; until then the flag has not arrived. Throws on failure. + void ThrowIfDecodeFailed(); + + // Device time spent decoding the last image, in seconds. Valid after the caller synchronises. + [[nodiscard]] float GetDecodeTime_s() const; }; diff --git a/image_analysis/image_preprocessing/ImagePreprocessor.h b/image_analysis/image_preprocessing/ImagePreprocessor.h index de30f73b..2264b649 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessor.h +++ b/image_analysis/image_preprocessing/ImagePreprocessor.h @@ -39,6 +39,11 @@ public: virtual bool AnalyzeCompressed(ImagePreprocessorBuffer &processed_image, const CompressedImage &image, ImageStatistics &stats) { return false; } + // Device time the last AnalyzeCompressed() spent getting the chunk across and decompressing it, + // so the caller can still report a decompression cost once the host no longer does the work. + // Meaningless unless the previous call returned true. + [[nodiscard]] virtual float GetLastDecompressionTime_s() const { return 0.0f; } + // Resize the buffer an image will be decompressed into and page-lock it, so that the host->device // copy of Analyze() is a real DMA. Without page-locking the driver stages the copy through its own // pinned pool, which is a host-side copy on the calling thread: it does not overlap and it degrades diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 31b25408..cc03a77c 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include + #include "ImagePreprocessorGPU.h" template @@ -89,12 +91,139 @@ __global__ void preprocess_kernel( } } +// The per-pixel decision preprocess_kernel makes, in a form the fused kernel can reuse so the two +// cannot drift apart. Priority: masked > error > saturated. +template +struct PreprocessAccum { + unsigned long long masked = 0, saturated = 0, error = 0; + long long max_v = INT64_MIN, min_v = INT64_MAX; + + __device__ __forceinline__ int32_t Apply(T v, bool is_masked, T sat_value, T err_value) { + const bool is_err = (v == err_value); + const bool is_sat = !is_err && (v >= sat_value); + masked += is_masked; + error += (!is_masked && is_err); + saturated += (!is_masked && !is_err && is_sat); + if (!(is_masked || is_sat || is_err)) { + const int64_t val = (int64_t) v; + if (val > max_v) max_v = val; + if (val < min_v) min_v = val; + } + return is_masked ? INT32_MIN : is_err ? INT32_MIN : is_sat ? INT32_MAX : (int32_t) v; + } +}; + +// Reduce a block's thread-local accumulators into the image-wide statistics. Every accumulator is an +// integer, so the result does not depend on the order the blocks arrive in. +template +__device__ __forceinline__ void FlushStats(PreprocessAccum &l, ImageStatistics *stats) { + __shared__ unsigned long long s_masked, s_saturated, s_error; + __shared__ long long s_max, s_min; + if (threadIdx.x == 0) { + s_masked = 0; s_saturated = 0; s_error = 0; s_max = INT64_MIN; s_min = INT64_MAX; + } + __syncthreads(); + + atomicAdd(&s_masked, l.masked); + atomicAdd(&s_saturated, l.saturated); + atomicAdd(&s_error, l.error); + if (l.min_v <= l.max_v) { + atomicMax(&s_max, l.max_v); + atomicMin(&s_min, l.min_v); + } + __syncthreads(); + + if (threadIdx.x == 0) { + atomicAdd(&stats->masked_pixel_count, s_masked); + atomicAdd(&stats->saturated_pixel_count, s_saturated); + atomicAdd(&stats->error_pixel_count, s_error); + atomicMax((long long *) &stats->max_value, s_max); + atomicMin((long long *) &stats->min_value, s_min); + } +} + +__device__ __forceinline__ uint64_t transpose8_fused(uint64_t x) { + uint64_t t; + t = (x ^ (x >> 7)) & 0x00aa00aa00aa00aaULL; x = x ^ t ^ (t << 7); + t = (x ^ (x >> 14)) & 0x0000cccc0000ccccULL; x = x ^ t ^ (t << 14); + t = (x ^ (x >> 28)) & 0x00000000f0f0f0f0ULL; x = x ^ t ^ (t << 28); + return x; +} + +// The bitshuffle inverse and the preprocessing in ONE pass. One thread owns one group of 8 elements +// across every byte-plane, so once it has transposed its 8 bytes out of each plane it holds 8 +// complete elements and can emit 8 finished int32 pixels - the decompressed image never has to exist +// in device memory at all. That removes a full-frame buffer per worker and a full-frame write plus +// read from the pipeline. +// +// 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 +__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, + 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) { + PreprocessAccum l; + + // The bytes are assembled in the unsigned counterpart of T - shifting a byte into the top of a + // signed type overflows it - and converted back at the end, which C++20 defines as two's + // complement reinterpretation. That is exactly what the byte order in the file means. + using U = typename std::make_unsigned::type; + + if (blockIdx.x == nblocks) { + if (threadIdx.x < tail_elems) { + U uv = 0; + #pragma unroll + 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); + } + FlushStats(l, stats); + return; + } + + const int b = blockIdx.x; + const uint32_t size = desc[b].nelem; // bytes per plane + const uint8_t *in = shuffled + desc[b].out_off; + const uint32_t elem0 = desc[b].out_off / ES; // first pixel of this block + const uint32_t n = size / 8; + + for (uint32_t i = threadIdx.x; i < n; i += blockDim.x) { + uint64_t x[ES]; + #pragma unroll + for (int p = 0; p < ES; p++) { + const uint8_t *pin = in + p * size; + uint64_t a = 0; + #pragma unroll + for (int k = 0; k < 8; k++) a |= (uint64_t)pin[k * n + i] << (8 * k); + x[p] = transpose8_fused(a); + } + int32_t o[8]; + #pragma unroll + for (int k = 0; k < 8; k++) { + U uv = 0; + #pragma unroll + 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]); + } + FlushStats(l, stats); +} + ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experiment, const PixelMask &mask, std::shared_ptr stream, bool copy_image_to_host) : ImagePreprocessor(experiment), stream(stream), copy_image_to_host(copy_image_to_host), - gpu_decompressed_image(npixels * sizeof(uint32_t)), // Overshoot - if input image is 1- or 2-byte, then it is still fine, while memory loss is minimal gpu_stats(1), cpu_stats(1), cpu_stats_reg(cpu_stats) { @@ -116,6 +245,10 @@ ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experime blocks = 4 * prop.multiProcessorCount; } +float ImagePreprocessorGPU::GetLastDecompressionTime_s() const { + return bslz4_decoder ? bslz4_decoder->GetDecodeTime_s() : 0.0f; +} + void ImagePreprocessorGPU::PinInputBuffer(std::vector &buffer, size_t size) { if (buffer.size() == size) return; @@ -156,33 +289,79 @@ bool ImagePreprocessorGPU::AnalyzeCompressed(ImagePreprocessorBuffer &processed_ if (!bslz4_decoder) bslz4_decoder = std::make_unique(npixels * sizeof(uint32_t), stream); - // Straight into the buffer Analyze() would have filled by copying the decompressed frame across - // the bus; from here the two paths are the same code. - bslz4_decoder->Decode(image, gpu_decompressed_image.get()); + // LZ4 on the device, then ONE kernel that un-transposes the bitshuffle blocks and preprocesses + // them as it goes. The decompressed image is never materialised: the fused kernel reads the + // shuffled bytes and writes finished int32 pixels. + const BSLZ4ShuffledImage shuffled = bslz4_decoder->DecodeShuffled(image); switch (image.GetMode()) { case CompressedImageMode::Int8: - stats = AnalyzeOnDevice(processed_image, INT8_MIN, INT8_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, INT8_MIN, INT8_MAX); return true; case CompressedImageMode::Uint8: - stats = AnalyzeOnDevice(processed_image, UINT8_MAX, UINT8_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, UINT8_MAX, UINT8_MAX); return true; case CompressedImageMode::Int16: - stats = AnalyzeOnDevice(processed_image, INT16_MIN, INT16_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, INT16_MIN, INT16_MAX); return true; case CompressedImageMode::Uint16: - stats = AnalyzeOnDevice(processed_image, UINT16_MAX, UINT16_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, UINT16_MAX, UINT16_MAX); return true; case CompressedImageMode::Int32: - stats = AnalyzeOnDevice(processed_image, INT32_MIN, INT32_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, INT32_MIN, INT32_MAX); return true; case CompressedImageMode::Uint32: - stats = AnalyzeOnDevice(processed_image, UINT32_MAX, UINT32_MAX); return true; + stats = UntransposeAndAnalyze(processed_image, shuffled, UINT32_MAX, UINT32_MAX); return true; default: return false; // Supports() already excludes these; belt and braces } } +// The device-decode counterpart of AnalyzeOnDevice: same per-pixel decision, same statistics, but +// fed from the bitshuffled bytes rather than from a decompressed image. +template +ImageStatistics ImagePreprocessorGPU::UntransposeAndAnalyze(ImagePreprocessorBuffer &processed_image, + const BSLZ4ShuffledImage &shuffled, + T err_value, T sat_value) { + if (sat_value > saturation_limit) + sat_value = static_cast(saturation_limit); + + cpu_stats[0] = ImageStatistics{.max_value = INT64_MIN, .min_value = INT64_MAX}; + cudaMemcpyAsync(gpu_stats, cpu_stats.data(), sizeof(ImageStatistics), cudaMemcpyHostToDevice, *stream); + + // 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 >>>( + 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); + + if (copy_image_to_host) + cudaMemcpyAsync(processed_image.data(), processed_image.getGPUBuffer(), npixels * sizeof(int32_t), cudaMemcpyDeviceToHost, *stream); + cudaMemcpyAsync(cpu_stats.data(), gpu_stats, sizeof(ImageStatistics), cudaMemcpyDeviceToHost, *stream); + + cudaStreamSynchronize(*stream); + + // Only now can the device tell us whether every block actually decoded. + bslz4_decoder->ThrowIfDecodeFailed(); + + return cpu_stats[0]; +} + template ImageStatistics ImagePreprocessorGPU::Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *input, T err_value, T sat_value) { + // Allocated here rather than in the constructor: only the host-upload path needs it, and the + // device-decode path - which is what BSHUF_LZ4 images take - never touches it. Overshoot to + // 4 bytes per pixel so a 1- or 2-byte image fits the same buffer. + if (!gpu_decompressed_image.get()) + gpu_decompressed_image = CudaDevicePtr(npixels * sizeof(uint32_t)); + // On this engine's own stream, not the NULL stream: a NULL-stream copy implicitly synchronises with // every blocking stream in the process, which serialised all workers behind whichever one was // uploading. The stream is synchronised at the end of this function, so the ordering is unchanged. diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h index ec07458d..60679c7e 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.h @@ -17,6 +17,9 @@ class ImagePreprocessorGPU : public ImagePreprocessor { int blocks; // Geometry-only, so one copy per GPU shared with every other engine on it (CudaSharedTables.h). std::shared_ptr> gpu_mask; + // Landing buffer for the HOST-upload path only. The device-decode path un-transposes straight + // into the preprocessed image, so it never needs this - and at 4 bytes per pixel it is worth a + // frame per worker, so it is allocated on first use rather than always. CudaDevicePtr gpu_decompressed_image; CudaDevicePtr gpu_stats; @@ -30,8 +33,13 @@ class ImagePreprocessorGPU : public ImagePreprocessor { std::unique_ptr bslz4_decoder; 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, shared by both entry points. + // Preprocess an image already sitting in gpu_decompressed_image (the host-upload path). template ImageStatistics AnalyzeOnDevice(ImagePreprocessorBuffer &processed_image, T err_value, T sat_value); + // Preprocess straight out of the bitshuffled bytes (the device-decode path). Same per-pixel + // decision and same statistics as AnalyzeOnDevice, with the un-transpose folded in. + template ImageStatistics UntransposeAndAnalyze(ImagePreprocessorBuffer &processed_image, + const BSLZ4ShuffledImage &shuffled, + T err_value, T sat_value); 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 @@ -41,6 +49,7 @@ public: ImageStatistics Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *decompressed_image, CompressedImageMode image_mode) override; bool AnalyzeCompressed(ImagePreprocessorBuffer &processed_image, const CompressedImage &image, ImageStatistics &stats) override; + [[nodiscard]] float GetLastDecompressionTime_s() const override; void PinInputBuffer(std::vector &buffer, size_t size) override; }; diff --git a/image_analysis/indexing/CUDAMemHelpers.h b/image_analysis/indexing/CUDAMemHelpers.h index ce97bf42..3af9fbd0 100644 --- a/image_analysis/indexing/CUDAMemHelpers.h +++ b/image_analysis/indexing/CUDAMemHelpers.h @@ -40,6 +40,35 @@ public: cudaStream_t get() const { return stream_; } }; +// A timing event, so a phase that is queued on a stream can still report how long the device spent +// on it. cudaEventDisableTiming is deliberately NOT used - timing is the whole point here. +class CudaEvent { + cudaEvent_t event_ = nullptr; +public: + CudaEvent() { + if (cudaEventCreate(&event_) != cudaSuccess) + throw JFJochException(JFJochExceptionCategory::GPUCUDAError, + "Failed to create CUDA event"); + } + ~CudaEvent() { + if (event_) cudaEventDestroy(event_); + } + CudaEvent(CudaEvent&& other) noexcept : event_(other.event_) { other.event_ = nullptr; } + CudaEvent& operator=(CudaEvent&& other) noexcept { + if (this != &other) { + if (event_) cudaEventDestroy(event_); + event_ = other.event_; + other.event_ = nullptr; + } + return *this; + } + CudaEvent(const CudaEvent&) = delete; + CudaEvent& operator=(const CudaEvent&) = delete; + + operator cudaEvent_t() const { return event_; } + cudaEvent_t get() const { return event_; } +}; + class CudaFFTPlan { cufftHandle plan_ = 0; public: diff --git a/tests/BSLZ4DecoderGPUFuzzTest.cpp b/tests/BSLZ4DecoderGPUFuzzTest.cpp new file mode 100644 index 00000000..289019e0 --- /dev/null +++ b/tests/BSLZ4DecoderGPUFuzzTest.cpp @@ -0,0 +1,1304 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only +// +// Adversarial differential test: BSLZ4DecoderGPU vs the CPU reference (JFJochDecompressHperfPtr). +// +// The device decoder is a hand-written LZ4 parser, so "it works on our files" is not enough - it has +// to agree with LZ4_decompress_safe on every construct a conforming encoder may emit, and on the +// block sizes foreign writers (DECTRIS, the stock bitshuffle HDF5 filter) use rather than only the +// 16 kB our own compressor emits. The cases below therefore drive the compressor into the corners +// (incompressible data for long literal runs and 255-extension chains, long repeats for overlapping +// matches at offset 1, engineered offsets across the offset-vs-matchlen branch boundary, sizes that +// put every rem%8 value in the verbatim tail) and, where the compressor's heuristics will not +// produce a wanted sequence at all, hand-build valid LZ4 blocks and verify them with the reference +// decoder before feeding them to the GPU. +// +// The statistics printed by the coverage report are the point: they PROVE which paths were reached +// rather than assuming the generated data happened to reach them. + +#include +#include "../common/CUDAWrapper.h" + +#ifdef JFJOCH_USE_CUDA + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../image_analysis/image_preprocessing/BSLZ4DecoderGPU.h" +#include "../compression/JFJochCompressor.h" +#include "../compression/JFJochDecompress.h" + +namespace { + +constexpr auto ALG = CompressionAlgorithm::BSHUF_LZ4; + +size_t BlockElems(size_t elem_size) { return JFJochBitShuffleCompressor::BlockSize(ALG, elem_size); } +size_t DefaultBlockBytes() { return JFJochBitShuffleCompressor::DefaultBlockSizeBytes(ALG); } + +CompressedImageMode ModeOf(size_t elem_size, bool is_signed) { + switch (elem_size) { + case 1: return is_signed ? CompressedImageMode::Int8 : CompressedImageMode::Uint8; + case 2: return is_signed ? CompressedImageMode::Int16 : CompressedImageMode::Uint16; + default: return is_signed ? CompressedImageMode::Int32 : CompressedImageMode::Uint32; + } +} + +void PutBE64(uint8_t *p, uint64_t v) { for (int i = 0; i < 8; i++) p[i] = (uint8_t)(v >> (8 * (7 - i))); } +void PutBE32(uint8_t *p, uint32_t v) { for (int i = 0; i < 4; i++) p[i] = (uint8_t)(v >> (8 * (3 - i))); } +uint64_t GetBE64(const uint8_t *p) { uint64_t v = 0; for (int i = 0; i < 8; i++) v = (v << 8) | p[i]; return v; } +uint32_t GetBE32(const uint8_t *p) { return ((uint32_t)p[0] << 24) | ((uint32_t)p[1] << 16) | ((uint32_t)p[2] << 8) | p[3]; } + +// --------------------------------------------------------------------------------------------- +// LZ4 sequence statistics - so the test can PROVE which code paths it actually reached rather +// than assume the compressor produced them. +// --------------------------------------------------------------------------------------------- +struct Lz4Stats { + size_t images = 0, blocks = 0, sequences = 0; + size_t lit_ge15 = 0, lit_ext255 = 0, max_litlen = 0; + size_t ml_ge19 = 0, ml_ext255 = 0, max_matchlen = 0; + size_t overlap = 0, nonoverlap = 0; + size_t off_eq_ml = 0, off_ml_plus1 = 0, off_ml_minus1 = 0; + std::map offsets_le64; + uint32_t min_offset = 0xffffffffu; + + void Merge(const Lz4Stats &o) { + images += o.images; blocks += o.blocks; sequences += o.sequences; + lit_ge15 += o.lit_ge15; lit_ext255 += o.lit_ext255; max_litlen = std::max(max_litlen, o.max_litlen); + ml_ge19 += o.ml_ge19; ml_ext255 += o.ml_ext255; max_matchlen = std::max(max_matchlen, o.max_matchlen); + overlap += o.overlap; nonoverlap += o.nonoverlap; + off_eq_ml += o.off_eq_ml; off_ml_plus1 += o.off_ml_plus1; off_ml_minus1 += o.off_ml_minus1; + for (auto &kv : o.offsets_le64) offsets_le64[kv.first] += kv.second; + min_offset = std::min(min_offset, o.min_offset); + } +}; + +std::mutex g_stats_mtx; +Lz4Stats g_stats; + +void AccumulateStats(const Lz4Stats &s) { + std::lock_guard lock(g_stats_mtx); + g_stats.Merge(s); +} + +// Walk one LZ4 block exactly as the kernel does and record what the sequences look like. +void ParseLz4Block(const uint8_t *p, size_t len, Lz4Stats &st) { + const uint8_t *const end = p + len; + while (p < end) { + const uint32_t token = *p++; + uint32_t litlen = token >> 4; + if (litlen == 15) { + st.lit_ge15++; + uint32_t s, n_ext = 0; + do { s = *p++; litlen += s; n_ext++; } while (s == 255 && p < end); + if (n_ext > 1) st.lit_ext255++; + } + p += litlen; + st.max_litlen = std::max(st.max_litlen, litlen); + if (p >= end) break; // final sequence: literals only + const uint32_t offset = (uint32_t)p[0] | ((uint32_t)p[1] << 8); + p += 2; + uint32_t ml = token & 0x0F; + if (ml == 15) { + st.ml_ge19++; + uint32_t s, n_ext = 0; + do { s = *p++; ml += s; n_ext++; } while (s == 255 && p < end); + if (n_ext > 1) st.ml_ext255++; + } + ml += 4; + st.sequences++; + st.max_matchlen = std::max(st.max_matchlen, ml); + if (offset < ml) st.overlap++; else st.nonoverlap++; + if (offset == ml) st.off_eq_ml++; + if (offset == ml + 1) st.off_ml_plus1++; + if (offset + 1 == ml) st.off_ml_minus1++; + if (offset <= 64) st.offsets_le64[offset]++; + st.min_offset = std::min(st.min_offset, offset); + } +} + +void ParseContainer(const uint8_t *c, size_t clen, size_t elem_size, Lz4Stats &st) { + if (clen < 12) return; + const size_t total = GetBE64(c); + const size_t block_elems = GetBE32(c + 8) / elem_size; + const size_t nelements = total / elem_size; + const size_t nfull = nelements / block_elems; + const size_t rem = nelements - nfull * block_elems; + const size_t last = rem - rem % 8; + size_t off = 12; + st.images++; + for (size_t i = 0; i < nfull + (last > 0 ? 1 : 0); i++) { + if (off + 4 > clen) return; + const uint32_t bl = GetBE32(c + off); + off += 4; + if (off + bl > clen) return; + st.blocks++; + ParseLz4Block(c + off, bl, st); + off += bl; + } +} + +// --------------------------------------------------------------------------------------------- +// Byte-exact comparison with a diagnosable failure message. +// --------------------------------------------------------------------------------------------- +std::string DescribeMismatch(const uint8_t *expected, const uint8_t *got, size_t nbytes, + size_t block_bytes, size_t body_bytes, const std::string &label) { + for (size_t i = 0; i < nbytes; i++) { + if (expected[i] != got[i]) { + size_t ndiff = 0; + for (size_t j = i; j < nbytes; j++) if (expected[j] != got[j]) ndiff++; + std::ostringstream os; + os << label << ": FIRST DIFF at byte " << i + << " expected 0x" << std::hex << (unsigned)expected[i] + << " got 0x" << (unsigned)got[i] << std::dec + << " | bitshuffle block " << (i / block_bytes) + << " (offset within block " << (i % block_bytes) << " of " << block_bytes << ")" + << (i >= body_bytes ? " [IN THE VERBATIM LEFTOVER TAIL]" : "") + << " | total differing bytes " << ndiff << " of " << nbytes; + return os.str(); + } + } + return {}; +} + +// --------------------------------------------------------------------------------------------- +// The core differential check: production compressor -> CPU decode and GPU decode -> must match. +// --------------------------------------------------------------------------------------------- +struct CaseResult { + bool cpu_roundtrip_ok = false; + bool gpu_matches_cpu = false; + std::string message; + size_t compressed_bytes = 0; +}; + +CaseResult RunCaseNoAssert(const std::string &label, + const std::vector &image_bytes, + size_t elem_size, bool is_signed, + JFJochBitShuffleCompressor &compressor, + BSLZ4DecoderGPU &decoder, + CudaStream &stream, + bool collect_stats) { + CaseResult r; + const size_t nelements = image_bytes.size() / elem_size; + const std::vector compressed = compressor.Compress(image_bytes.data(), nelements, elem_size); + r.compressed_bytes = compressed.size(); + + if (collect_stats) { + Lz4Stats st; + ParseContainer(compressed.data(), compressed.size(), elem_size, st); + AccumulateStats(st); + } + + const CompressedImage image(compressed.data(), compressed.size(), nelements, 1, + ModeOf(elem_size, is_signed), ALG); + + // CPU reference - the exact call the host path makes. + std::vector cpu(image.GetUncompressedSize()); + image.GetUncompressed(cpu); + r.cpu_roundtrip_ok = (cpu == image_bytes); + + CudaDevicePtr gpu_out(image_bytes.size()); + // Poison the output so a byte the decoder never writes is caught rather than accidentally right. + cudaMemset(gpu_out.get(), 0xA5, image_bytes.size()); + decoder.Decode(image, gpu_out.get()); + if (cudaStreamSynchronize(stream) != cudaSuccess) { + r.message = label + ": cudaStreamSynchronize failed"; + return r; + } + std::vector gpu(image_bytes.size()); + if (cudaMemcpy(gpu.data(), gpu_out.get(), image_bytes.size(), cudaMemcpyDeviceToHost) != cudaSuccess) { + r.message = label + ": cudaMemcpy D2H failed"; + return r; + } + + const size_t block_bytes = GetBE32(compressed.data() + 8); + const size_t block_elems = block_bytes / elem_size; + const size_t rem = nelements % block_elems; + const size_t body_bytes = (nelements - rem % 8) * elem_size; + + const std::string msg = DescribeMismatch(cpu.data(), gpu.data(), image_bytes.size(), + block_bytes, body_bytes, label); + r.gpu_matches_cpu = msg.empty(); + if (!msg.empty()) r.message = msg; + return r; +} + +// Convenience wrapper that owns a fresh decoder/stream/compressor and asserts. +void RunCase(const std::string &label, const std::vector &image_bytes, + size_t elem_size, bool is_signed) { + JFJochBitShuffleCompressor compressor(ALG); + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(image_bytes.size(), stream); + const auto r = RunCaseNoAssert(label, image_bytes, elem_size, is_signed, + compressor, decoder, *stream, true); + INFO(label << " (" << image_bytes.size() << " bytes -> " << r.compressed_bytes << " compressed)"); + if (!r.cpu_roundtrip_ok) FAIL(label << ": CPU reference did not round-trip (compressor/CPU bug, not GPU)"); + if (!r.gpu_matches_cpu) FAIL(r.message); + SUCCEED(); +} + +// --------------------------------------------------------------------------------------------- +// Image generators (image domain). +// --------------------------------------------------------------------------------------------- +std::vector GenIncompressible(size_t nbytes, uint32_t seed) { + std::mt19937 rng(seed); + std::vector v(nbytes); + for (auto &b : v) b = (uint8_t)(rng() & 0xFF); // full 0..255 range + return v; +} + +std::vector GenConstant(size_t nbytes, uint8_t val) { return std::vector(nbytes, val); } + +std::vector GenRepeatedPattern(size_t nbytes, size_t period, uint32_t seed) { + std::mt19937 rng(seed); + std::vector pat(period); + for (auto &b : pat) b = (uint8_t)(rng() & 0xFF); + std::vector v(nbytes); + for (size_t i = 0; i < nbytes; i++) v[i] = pat[i % period]; + return v; +} + +// Incrementing element values - bit planes come out highly structured after bitshuffling. +std::vector GenRamp(size_t nelements, size_t elem_size, uint64_t start, uint64_t step) { + std::vector v(nelements * elem_size); + for (size_t i = 0; i < nelements; i++) { + const uint64_t val = start + i * step; + for (size_t j = 0; j < elem_size; j++) v[i * elem_size + j] = (uint8_t)(val >> (8 * j)); + } + return v; +} + +std::vector GenDetectorLike(size_t nelements, size_t elem_size, uint32_t seed) { + std::mt19937 rng(seed); + std::vector v(nelements * elem_size, 0); + auto put = [&](size_t i, uint64_t val) { + for (size_t j = 0; j < elem_size; j++) v[i * elem_size + j] = (uint8_t)(val >> (8 * j)); + }; + for (size_t i = nelements / 4; i < nelements / 2; i++) put(i, rng() % 7); + for (size_t s = 0; s < 64 && nelements > 0; s++) { + const size_t c = rng() % nelements; + for (size_t d = 0; d < 9 && c + d < nelements; d++) put(c + d, 30000 + d); + } + for (size_t i = nelements * 3 / 4; i < nelements && i < nelements * 3 / 4 + 5000; i++) put(i, 42); + return v; +} + +// A grab-bag: random structure chosen per region, so a single case covers several regimes. +std::vector GenMixed(size_t nelements, size_t elem_size, uint32_t seed) { + std::mt19937 rng(seed); + std::vector v(nelements * elem_size, 0); + size_t i = 0; + while (i < nelements) { + const size_t run = 1 + rng() % 4096; + const size_t n = std::min(run, nelements - i); + const int kind = rng() % 5; + for (size_t k = 0; k < n; k++) { + uint64_t val; + switch (kind) { + case 0: val = 0; break; + case 1: val = rng(); break; // incompressible region + case 2: val = 12345; break; // constant region + case 3: val = k; break; // ramp region + default: val = (rng() % 3); break; // low-entropy noise + } + for (size_t j = 0; j < elem_size; j++) + v[(i + k) * elem_size + j] = (uint8_t)(val >> (8 * j)); + } + i += n; + } + return v; +} + +// --------------------------------------------------------------------------------------------- +// Shuffled-domain engineering: bitshuffle is a bit permutation, so to make LZ4 see a chosen byte +// pattern S we hand the compressor bitshuf_decode_block(S). encode(decode(S)) == S, which the +// helper verifies, so the LZ4 stage really does get the offsets we intend. +// --------------------------------------------------------------------------------------------- +std::vector ImageFromShuffled(const std::vector &shuffled, size_t nelements, + size_t elem_size, bool verify) { + std::vector img(nelements * elem_size); + std::vector scratch(nelements * elem_size); + REQUIRE(bitshuf_decode_block((char *)img.data(), (const char *)shuffled.data(), + scratch.data(), nelements, elem_size) == 0); + if (verify) { + std::vector re(nelements * elem_size); + REQUIRE(bitshuf_encode_block((char *)re.data(), (const char *)img.data(), + scratch.data(), nelements, elem_size) == 0); + // Compare as a bool, not as vectors: Catch2 would otherwise stringify 16 kB of bytes. + const bool encode_decode_is_identity = (re == shuffled); + REQUIRE(encode_decode_is_identity); // the construction is only meaningful if this holds + } + return img; +} + +// A shuffled-domain buffer built from periodic regions of the requested periods, plus noise +// separators, so LZ4 emits matches at exactly those offsets. +std::vector ShuffledWithPeriods(size_t nbytes, const std::vector &periods, + uint32_t seed) { + std::mt19937 rng(seed); + std::vector s(nbytes); + size_t i = 0; + size_t pi = 0; + while (i < nbytes) { + // Noise separator so the periodic region starts fresh. + const size_t nsep = std::min(24 + rng() % 40, nbytes - i); + for (size_t k = 0; k < nsep; k++) s[i + k] = (uint8_t)(rng() & 0xFF); + i += nsep; + if (i >= nbytes) break; + const uint32_t P = periods[pi++ % periods.size()]; + // Region of period P; length varied so the resulting matchlen straddles the offset. + const size_t reps = 1 + rng() % 12; + const size_t len = std::min((size_t)P * (1 + reps) + (rng() % 3), nbytes - i); + for (size_t k = 0; k < len; k++) + s[i + k] = (k < P) ? (uint8_t)(rng() & 0xFF) : s[i + k - P]; + i += len; + } + return s; +} + +// Short NON-overlapping back-references at chosen distances: a chunk of L bytes repeated at +// distance P >= L, so offset >= matchlen. +std::vector ShuffledShortMatches(size_t nbytes, const std::vector &distances, + uint32_t seed) { + std::mt19937 rng(seed); + std::vector s(nbytes); + for (auto &b : s) b = (uint8_t)(rng() & 0xFF); + size_t i = 64; + size_t di = 0; + while (i + 128 < nbytes) { + const uint32_t P = distances[di++ % distances.size()]; + if (P < 4) { i += 64; continue; } + const uint32_t L = 4 + rng() % std::max(1u, P - 3); // 4 <= L <= P + if (i < P) { i += 64; continue; } + for (uint32_t k = 0; k < L && i + k < nbytes; k++) s[i + k] = s[i + k - P]; + i += L + 8 + rng() % 40; + } + return s; +} + +// --------------------------------------------------------------------------------------------- +// Hand-built LZ4 blocks, so exact (offset, matchlen) pairs can be reached that the real +// compressor's heuristics never emit - in particular offset == matchlen +- 1. +// --------------------------------------------------------------------------------------------- +// Only the LAST sequence of an LZ4 block may be literals-only, so literals are staged and flushed +// by the match that follows them (or by Finish()). +class Lz4BlockBuilder { + std::vector pending; +public: + std::vector code; // the encoded LZ4 block + std::vector plain; // what it must decode to + + static void PutExt(std::vector &v, uint32_t len) { + uint32_t r = len - 15; + while (r >= 255) { v.push_back(255); r -= 255; } + v.push_back((uint8_t)r); + } + + void AddLiterals(const std::vector &l) { pending.insert(pending.end(), l.begin(), l.end()); } + size_t OutputPos() const { return plain.size() + pending.size(); } + size_t PendingLiterals() const { return pending.size(); } + + void Match(uint32_t offset, uint32_t matchlen) { + const uint32_t litlen = (uint32_t)pending.size(); + const uint32_t mlcode = matchlen - 4; + code.push_back((uint8_t)((std::min(litlen, 15u) << 4) | std::min(mlcode, 15u))); + if (litlen >= 15) PutExt(code, litlen); + code.insert(code.end(), pending.begin(), pending.end()); + plain.insert(plain.end(), pending.begin(), pending.end()); + pending.clear(); + code.push_back((uint8_t)(offset & 0xFF)); + code.push_back((uint8_t)(offset >> 8)); + if (mlcode >= 15) PutExt(code, mlcode); + const size_t start = plain.size() - offset; + for (uint32_t i = 0; i < matchlen; i++) { + const uint8_t b = plain[start + i]; // may be a byte this very loop just wrote + plain.push_back(b); + } + } + + // Every LZ4 block ends with a literals-only sequence; LZ4 also requires the last match to end + // at least 12 bytes before the block end, so this tail must be >= 12 bytes. + void Finish(const std::vector &tail) { + AddLiterals(tail); + const uint32_t litlen = (uint32_t)pending.size(); + code.push_back((uint8_t)(std::min(litlen, 15u) << 4)); + if (litlen >= 15) PutExt(code, litlen); + code.insert(code.end(), pending.begin(), pending.end()); + plain.insert(plain.end(), pending.begin(), pending.end()); + pending.clear(); + } +}; + +std::vector MakeContainer(const std::vector> &block_codes, + size_t total_bytes, size_t block_bytes) { + std::vector c(12); + PutBE64(c.data(), total_bytes); + PutBE32(c.data() + 8, (uint32_t)block_bytes); + for (const auto &bc : block_codes) { + uint8_t len[4]; + PutBE32(len, (uint32_t)bc.size()); + c.insert(c.end(), len, len + 4); + c.insert(c.end(), bc.begin(), bc.end()); + } + return c; +} + +} // namespace + +// ============================================================================================= +// (a) INCOMPRESSIBLE data - long pure-literal sequences, litlen >= 15 with 255-extension chains. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_Incompressible", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + size_t n_cases = 0; + for (size_t es : {1u, 2u, 4u}) { + for (bool sgn : {false, true}) { + const size_t B = BlockElems(es); + for (size_t nelem : {B / 2, B, B + 1, B + 7, B + 8, 3 * B, 3 * B + 13, 5 * B + 8 * 7 + 3}) { + std::ostringstream lbl; + lbl << "incompressible es=" << es << (sgn ? " signed" : " unsigned") << " nelem=" << nelem; + RunCase(lbl.str(), GenIncompressible(nelem * es, (uint32_t)(1000 + n_cases)), es, sgn); + n_cases++; + } + } + } + printf("[BSLZ4Fuzz] incompressible: %zu cases passed\n", n_cases); +} + +// ============================================================================================= +// (b) HIGHLY COMPRESSIBLE - very long matches, 255-extension match chains, offset==1 overlaps. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_HighlyCompressible", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + size_t n_cases = 0; + for (size_t es : {1u, 2u, 4u}) { + for (bool sgn : {false, true}) { + const size_t B = BlockElems(es); + for (size_t nelem : {B, B + 5, 4 * B, 4 * B + 8 * 11 + 6}) { + std::ostringstream l1; l1 << "all-zero es=" << es << " n=" << nelem << (sgn ? " s" : " u"); + RunCase(l1.str(), GenConstant(nelem * es, 0x00), es, sgn); + std::ostringstream l2; l2 << "all-0xFF es=" << es << " n=" << nelem << (sgn ? " s" : " u"); + RunCase(l2.str(), GenConstant(nelem * es, 0xFF), es, sgn); + std::ostringstream l3; l3 << "all-0x5A es=" << es << " n=" << nelem << (sgn ? " s" : " u"); + RunCase(l3.str(), GenConstant(nelem * es, 0x5A), es, sgn); + for (size_t period : {1u, 3u, 4u, 7u, 16u, 33u, 64u}) { + std::ostringstream l4; + l4 << "repeat-pattern p=" << period << " es=" << es << " n=" << nelem << (sgn ? " s" : " u"); + RunCase(l4.str(), GenRepeatedPattern(nelem * es, period * es, (uint32_t)(2000 + n_cases)), es, sgn); + n_cases++; + } + n_cases += 3; + } + } + } + printf("[BSLZ4Fuzz] highly compressible: %zu cases passed\n", n_cases); +} + +// ============================================================================================= +// (c) MIXED / STEPPED, and engineered LZ4 offsets 1,2,3,4,5,15,16,17,31,32,33 in both the +// offset >= matchlen and offset < matchlen branches. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_MixedAndEngineeredOffsets", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + const std::vector targets = {1, 2, 3, 4, 5, 15, 16, 17, 31, 32, 33}; + size_t n_cases = 0; + + for (size_t es : {1u, 2u, 4u}) { + for (bool sgn : {false, true}) { + const size_t B = BlockElems(es); + + // Stepped data: after bitshuffling the bit planes are highly structured. + for (uint64_t step : {1u, 2u, 3u, 256u, 65536u}) { + std::ostringstream l; l << "ramp step=" << step << " es=" << es << (sgn ? " s" : " u"); + RunCase(l.str(), GenRamp(2 * B + 8 * 5 + 3, es, 0, step), es, sgn); + n_cases++; + } + + // Mixed regions. + for (uint32_t seed = 0; seed < 4; seed++) { + std::ostringstream l; l << "mixed es=" << es << " seed=" << seed << (sgn ? " s" : " u"); + RunCase(l.str(), GenMixed(3 * B + 8 * 3 + 5, es, 3000 + seed * 17 + (uint32_t)es), es, sgn); + n_cases++; + } + + // Detector-like. + for (uint32_t seed = 0; seed < 3; seed++) { + std::ostringstream l; l << "detector-like es=" << es << " seed=" << seed << (sgn ? " s" : " u"); + RunCase(l.str(), GenDetectorLike(4 * B + 8 * 9 + 7, es, 4000 + seed * 31 + (uint32_t)es), es, sgn); + n_cases++; + } + + // Engineered periodic offsets -> offset < matchlen (the overlapping branch). + for (uint32_t seed = 0; seed < 4; seed++) { + const size_t nelem = 3 * B; + const auto shuffled = ShuffledWithPeriods(nelem * es, targets, 5000 + seed * 7 + (uint32_t)es); + // Per-block construction, because bitshuffle works per block. + std::vector img; + img.reserve(nelem * es); + for (size_t b = 0; b < nelem / B; b++) { + std::vector sblk(shuffled.begin() + b * B * es, shuffled.begin() + (b + 1) * B * es); + const auto part = ImageFromShuffled(sblk, B, es, b == 0 && seed == 0); + img.insert(img.end(), part.begin(), part.end()); + } + std::ostringstream l; l << "engineered-periodic es=" << es << " seed=" << seed << (sgn ? " s" : " u"); + RunCase(l.str(), img, es, sgn); + n_cases++; + } + + // Engineered short back-references -> offset >= matchlen (the non-overlapping branch). + for (uint32_t seed = 0; seed < 4; seed++) { + const size_t nelem = 3 * B; + const auto shuffled = ShuffledShortMatches(nelem * es, targets, 6000 + seed * 11 + (uint32_t)es); + std::vector img; + img.reserve(nelem * es); + for (size_t b = 0; b < nelem / B; b++) { + std::vector sblk(shuffled.begin() + b * B * es, shuffled.begin() + (b + 1) * B * es); + const auto part = ImageFromShuffled(sblk, B, es, false); + img.insert(img.end(), part.begin(), part.end()); + } + std::ostringstream l; l << "engineered-short-match es=" << es << " seed=" << seed << (sgn ? " s" : " u"); + RunCase(l.str(), img, es, sgn); + n_cases++; + } + } + } + printf("[BSLZ4Fuzz] mixed/engineered: %zu cases passed\n", n_cases); +} + +// ============================================================================================= +// (d) SIZE SWEEP - every interesting last-block / leftover-tail configuration. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_SizeSweep", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + size_t n_cases = 0; + for (size_t es : {1u, 2u, 4u}) { + const size_t B = BlockElems(es); + INFO("elem_size " << es << " block " << B << " elements = " << B * es << " bytes"); + std::vector sizes; + sizes.push_back(8); // the smallest legal block + sizes.push_back(16); + sizes.push_back(B); // exactly one full block, no remainder, no leftover + sizes.push_back(2 * B); // exact multiple: no partial block, no leftover + sizes.push_back(7 * B); + for (size_t k = 1; k <= 8; k++) sizes.push_back(B + k); // 1..8 leftover elements + for (size_t r = 0; r < 8; r++) sizes.push_back(B + 128 + r); // rem % 8 == r with a real last block + for (size_t r = 0; r < 8; r++) sizes.push_back(3 * B + 8 * 37 + r); + sizes.push_back(B - 8); // one partial block only, smaller than a full block + sizes.push_back(B - 1); // partial block + leftover + sizes.push_back(B / 2 + 3); + + for (bool sgn : {false, true}) { + for (size_t nelem : sizes) { + for (int gen = 0; gen < 3; gen++) { + std::vector img; + if (gen == 0) img = GenIncompressible(nelem * es, (uint32_t)(7000 + n_cases)); + else if (gen == 1) img = GenDetectorLike(nelem, es, (uint32_t)(8000 + n_cases)); + else img = GenConstant(nelem * es, 0x77); + std::ostringstream l; + l << "size-sweep es=" << es << " nelem=" << nelem << " gen=" << gen << (sgn ? " s" : " u"); + RunCase(l.str(), img, es, sgn); + n_cases++; + } + } + } + } + printf("[BSLZ4Fuzz] size sweep: %zu cases passed\n", n_cases); +} + +// Images with fewer than BSHUF_BLOCKED_MULT elements contain no bitshuffle block at all - only the +// verbatim tail. Documented behaviour probe, not an assertion about correctness. +TEST_CASE("BSLZ4Fuzz_SubBlockSizes", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + JFJochBitShuffleCompressor compressor(ALG); + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(4096, stream); + for (size_t nelem = 1; nelem <= 7; nelem++) { + const auto img = GenIncompressible(nelem * 4, (uint32_t)(9000 + nelem)); + const auto compressed = compressor.Compress(img.data(), nelem, 4); + const CompressedImage image(compressed.data(), compressed.size(), nelem, 1, + CompressedImageMode::Uint32, ALG); + std::vector cpu; + image.GetUncompressed(cpu); + CHECK(cpu == img); // CPU handles a block-less image + // Fewer than 8 elements is entirely the verbatim tail - there is no bitshuffle block at all. + // The device decoder must still produce the image rather than decline it, so that the caller + // does not need a special case the host route would have handled. + CudaDevicePtr gpu_out(nelem * 4); + cudaMemset(gpu_out.get(), 0xA5, nelem * 4); + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + std::vector gpu(nelem * 4); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * 4, cudaMemcpyDeviceToHost) == cudaSuccess); + CHECK(gpu == img); + } +} + +// ============================================================================================= +// (f) LARGE REALISTIC FRAME + speed. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_LargeFrame", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + const size_t W = 4148, H = 4362, N = W * H, ES = 4; + const auto img = GenDetectorLike(N, ES, 12345); + + JFJochBitShuffleCompressor compressor(ALG); + const auto compressed = compressor.Compress(img.data(), N, ES); + Lz4Stats st; + ParseContainer(compressed.data(), compressed.size(), ES, st); + AccumulateStats(st); + + const CompressedImage image(compressed.data(), compressed.size(), N, 1, CompressedImageMode::Uint32, ALG); + std::vector cpu; + auto t0 = std::chrono::steady_clock::now(); + image.GetUncompressed(cpu); + auto t1 = std::chrono::steady_clock::now(); + { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); } + + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(N * ES, stream); + CudaDevicePtr gpu_out(N * ES); + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + std::vector gpu(N * ES); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), N * ES, cudaMemcpyDeviceToHost) == cudaSuccess); + const auto msg = DescribeMismatch(cpu.data(), gpu.data(), N * ES, 16384, N * ES, "large frame"); + if (!msg.empty()) FAIL(msg); + + // Timed loop (upload + both kernels), after a warm-up. + const int reps = 20; + auto t2 = std::chrono::steady_clock::now(); + for (int i = 0; i < reps; i++) decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + auto t3 = std::chrono::steady_clock::now(); + + const double cpu_ms = std::chrono::duration(t1 - t0).count(); + const double gpu_ms = std::chrono::duration(t3 - t2).count() / reps; + printf("[BSLZ4Fuzz] large frame %zux%zu uint32 = %.1f MB, compressed %.1f MB (%.2fx), %zu blocks\n", + W, H, N * ES / 1e6, compressed.size() / 1e6, (double)(N * ES) / compressed.size(), st.blocks); + printf("[BSLZ4Fuzz] CPU decode (1 thread) %.2f ms GPU decode %.3f ms speedup %.1fx\n", + cpu_ms, gpu_ms, cpu_ms / gpu_ms); + SUCCEED(); +} + +// ============================================================================================= +// (g) BUFFER REUSE - one decoder instance, many images of different sizes in sequence. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_DecoderReuse", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + const size_t big_elems = 2048 * 2048; // 8 MB at uint32 + JFJochBitShuffleCompressor compressor(ALG); + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(big_elems * 4, stream); + + struct Item { size_t nelem; size_t es; int gen; }; + const std::vector seq = { + {big_elems, 4, 1}, // large + {1000, 4, 0}, // tiny right after a large one + {big_elems, 4, 0}, // large again + {8, 4, 0}, // minimum + {big_elems / 3 + 5, 2, 1}, + {big_elems, 4, 2}, // all-constant large + {4096 + 3, 1, 0}, + {big_elems, 4, 1}, + {17, 2, 0}, + {big_elems * 2, 2, 1}, // same byte count, different element size + }; + size_t idx = 0; + for (const auto &it : seq) { + std::vector img; + if (it.gen == 0) img = GenIncompressible(it.nelem * it.es, (uint32_t)(11000 + idx)); + else if (it.gen == 1) img = GenDetectorLike(it.nelem, it.es, (uint32_t)(12000 + idx)); + else img = GenConstant(it.nelem * it.es, 0x3C); + std::ostringstream l; l << "reuse[" << idx << "] es=" << it.es << " nelem=" << it.nelem; + const auto r = RunCaseNoAssert(l.str(), img, it.es, false, compressor, decoder, *stream, true); + INFO(l.str()); + if (!r.cpu_roundtrip_ok) FAIL(l.str() << ": CPU reference did not round-trip"); + if (!r.gpu_matches_cpu) FAIL(r.message); + idx++; + } + printf("[BSLZ4Fuzz] decoder reuse: %zu sequential images through ONE decoder passed\n", seq.size()); + SUCCEED(); +} + +// ============================================================================================= +// (h) CONCURRENCY. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_ConcurrentSeparateDecoders", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + const int nthreads = 4; + const int iters = 12; + const size_t nelem = 1024 * 1024; + + std::vector failures(nthreads); + std::vector passed(nthreads, 0); + + std::vector threads; + for (int t = 0; t < nthreads; t++) { + threads.emplace_back([&, t]() { + try { + JFJochBitShuffleCompressor compressor(ALG); + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(nelem * 4, stream); + for (int i = 0; i < iters; i++) { + const size_t n = nelem - (size_t)(t * 977 + i * 13); + std::vector img = (i % 2) + ? GenDetectorLike(n, 4, (uint32_t)(20000 + t * 100 + i)) + : GenIncompressible(n * 4, (uint32_t)(21000 + t * 100 + i)); + std::ostringstream l; l << "thread " << t << " iter " << i << " nelem=" << n; + const auto r = RunCaseNoAssert(l.str(), img, 4, false, compressor, decoder, *stream, false); + if (!r.cpu_roundtrip_ok) { failures[t] = l.str() + ": CPU round-trip failed"; return; } + if (!r.gpu_matches_cpu) { failures[t] = r.message; return; } + passed[t]++; + } + } catch (const std::exception &e) { + failures[t] = std::string("thread threw: ") + e.what(); + } + }); + } + for (auto &th : threads) th.join(); + + int total = 0; + for (int t = 0; t < nthreads; t++) { + total += passed[t]; + if (!failures[t].empty()) FAIL(failures[t]); + } + printf("[BSLZ4Fuzz] concurrency (own decoder + own stream per thread): %d/%d decodes byte-exact\n", + total, nthreads * iters); + SUCCEED(); +} + +// ============================================================================================= +// Hand-built LZ4 sequences: exact (offset, matchlen) pairs across the branch boundary. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_SyntheticOffsetMatrix", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + + // A bitshuffle block is 16384 bytes, so an LZ4 offset can never exceed that here even though the + // format allows 65535. + const std::vector offsets = {1, 2, 3, 4, 5, 6, 7, 8, 15, 16, 17, 30, 31, 32, 33, 34, + 63, 64, 65, 127, 128, 129, 255, 256, 257, 1024, 4095, + 8192, 16000, 16300}; + const std::vector elem_sizes = {1, 2, 4}; + + size_t n_cases = 0, n_seq = 0, n_skipped = 0; + size_t n_overlap = 0, n_nonoverlap = 0, n_boundary = 0; + for (uint32_t es : elem_sizes) { + const size_t B = BlockElems(es); // elements per block + const size_t block_bytes = B * es; // == 16384 + + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(block_bytes, stream); + CudaDevicePtr gpu_out(block_bytes); + + for (uint32_t off : offsets) { + // For each offset, a set of match lengths that straddles the offset >= matchlen boundary. + std::vector mls; + for (int d : {-2, -1, 0, 1, 2}) { + const long v = (long)off + d; + if (v >= 4) mls.push_back((uint32_t)v); + } + mls.push_back(4); + mls.push_back(18); // largest matchlen with no extension byte + mls.push_back(19); // first matchlen needing an extension byte + mls.push_back(4 + 15 + 255); // exactly one full 255 extension byte + mls.push_back(4 + 15 + 255 + 1); + mls.push_back(4 + 15 + 255 * 2 + 7); // a 255-extension CHAIN + std::sort(mls.begin(), mls.end()); + mls.erase(std::unique(mls.begin(), mls.end()), mls.end()); + + std::mt19937 rng(40000 + off * 31 + es); + Lz4BlockBuilder bb; + // Prime the output with enough literals that this offset is a legal back-reference. + { + std::vector lits(std::max(64, off + 64)); + for (auto &b : lits) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(lits); + } + if (bb.OutputPos() + 32 > block_bytes) { n_skipped++; continue; } + size_t seq_in_block = 0; + for (uint32_t ml : mls) { + if (bb.OutputPos() + ml + 64 > block_bytes) break; + if (off > bb.OutputPos()) continue; + bb.Match(off, ml); + n_seq++; seq_in_block++; + // A few literals before the next sequence, so litlen varies too. + const uint32_t litlen = (uint32_t)(rng() % 40); + std::vector lits(litlen); + for (auto &b : lits) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(lits); + } + if (seq_in_block == 0) { n_skipped++; continue; } + // LZ4 requires the block to end with a literal run and the last match to end at least + // 12 bytes before the block end; pad out to exactly one full block. + REQUIRE(bb.OutputPos() + 12 <= block_bytes); + std::vector tail(block_bytes - bb.OutputPos()); + for (auto &b : tail) b = (uint8_t)(rng() & 0xFF); + bb.Finish(tail); + REQUIRE(bb.plain.size() == block_bytes); + + // The hand-encoded stream must be a stream the real LZ4 decoder accepts, otherwise + // this test would be measuring nonsense. + std::vector lz4_check(block_bytes); + const int ret = LZ4_decompress_safe((const char *)bb.code.data(), (char *)lz4_check.data(), + (int)bb.code.size(), (int)block_bytes); + REQUIRE(ret == (int)block_bytes); + { const bool lz4_ok = (lz4_check == bb.plain); REQUIRE(lz4_ok); } + + // Confirm the block really carries the offset under test, in both branches. + Lz4Stats st; + ParseLz4Block(bb.code.data(), bb.code.size(), st); + REQUIRE(st.sequences == seq_in_block); + REQUIRE(st.min_offset == off); + n_overlap += st.overlap; + n_nonoverlap += st.nonoverlap; + n_boundary += st.off_eq_ml + st.off_ml_plus1 + st.off_ml_minus1; + + // bb.plain is the SHUFFLED image; the image itself is its bitshuffle inverse. + const auto img = ImageFromShuffled(bb.plain, B, es, off == 1); + + const auto container = MakeContainer({bb.code}, block_bytes, block_bytes); + const CompressedImage image(container.data(), container.size(), B, 1, ModeOf(es, false), ALG); + REQUIRE(BSLZ4DecoderGPU::Supports(image)); + + std::vector cpu; + image.GetUncompressed(cpu); + { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); } + + cudaMemset(gpu_out.get(), 0xA5, block_bytes); + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + std::vector gpu(block_bytes); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes, cudaMemcpyDeviceToHost) == cudaSuccess); + + std::ostringstream l; l << "synthetic offset=" << off << " es=" << es; + const auto msg = DescribeMismatch(cpu.data(), gpu.data(), block_bytes, block_bytes, block_bytes, l.str()); + if (!msg.empty()) FAIL(msg); + n_cases++; + } + } + printf("[BSLZ4Fuzz] synthetic hand-built LZ4: %zu blocks, %zu engineered (offset,matchlen) " + "sequences passed (%zu offsets skipped as too large to prime in a 16 kB block)\n", + n_cases, n_seq, n_skipped); + printf("[BSLZ4Fuzz] of those: %zu offset=matchlen, " + "%zu exactly on the offset==matchlen+-1 boundary\n", n_overlap, n_nonoverlap, n_boundary); + SUCCEED(); +} + +// Very long literal runs with 255-extension chains, hand-built so the chain length is certain. +TEST_CASE("BSLZ4Fuzz_SyntheticLongLiterals", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + for (uint32_t es : {1u, 2u, 4u}) { + const size_t B = BlockElems(es), block_bytes = B * es; + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(block_bytes, stream); + CudaDevicePtr gpu_out(block_bytes); + + for (uint32_t litlen : {14u, 15u, 16u, 269u, 270u, 271u, 524u, 525u, 4000u, 8000u}) { + std::mt19937 rng(50000 + litlen + es); + Lz4BlockBuilder bb; + // Prime and flush, so the long literal run below is a sequence of its own with exactly + // the litlen under test and a legal back-reference after it. + { + std::vector prime(512); + for (auto &b : prime) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(prime); + bb.Match(64, 64); + } + std::vector lits(litlen); + for (auto &b : lits) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(lits); + REQUIRE(bb.PendingLiterals() == litlen); + bb.Match(40, 300); // this sequence carries exactly `litlen` literals + REQUIRE(bb.OutputPos() + 12 <= block_bytes); + std::vector tail(block_bytes - bb.OutputPos()); + for (auto &b : tail) b = (uint8_t)(rng() & 0xFF); + bb.Finish(tail); + REQUIRE(bb.plain.size() == block_bytes); + + std::vector chk(block_bytes); + REQUIRE(LZ4_decompress_safe((const char *)bb.code.data(), (char *)chk.data(), + (int)bb.code.size(), (int)block_bytes) == (int)block_bytes); + { const bool lz4_ok = (chk == bb.plain); REQUIRE(lz4_ok); } + + // Confirm the encoder really emitted the literal run under test (and its 255-chain). + Lz4Stats st; + ParseLz4Block(bb.code.data(), bb.code.size(), st); + REQUIRE(st.max_litlen >= litlen); + + const auto img = ImageFromShuffled(bb.plain, B, es, false); + const auto container = MakeContainer({bb.code}, block_bytes, block_bytes); + const CompressedImage image(container.data(), container.size(), B, 1, ModeOf(es, false), ALG); + std::vector cpu; + image.GetUncompressed(cpu); + { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); } + + cudaMemset(gpu_out.get(), 0xA5, block_bytes); + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + std::vector gpu(block_bytes); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes, cudaMemcpyDeviceToHost) == cudaSuccess); + std::ostringstream l; l << "synthetic litlen=" << litlen << " es=" << es; + const auto msg = DescribeMismatch(cpu.data(), gpu.data(), block_bytes, block_bytes, block_bytes, l.str()); + if (!msg.empty()) FAIL(msg); + } + } + printf("[BSLZ4Fuzz] synthetic long-literal cases passed\n"); + SUCCEED(); +} + +// ============================================================================================= +// FOREIGN BLOCK SIZES. +// +// Everything above goes through this repo's compressor, which always emits 16384-byte bitshuffle +// blocks. Real DECTRIS/EIGER files - and anything written by the stock bitshuffle HDF5 filter - +// use bshuf_default_block_size(), i.e. 8192 BYTES, half of that; the filter also lets a writer +// choose any multiple of 8. The GPU decoder reads the block size out of the container header, so +// it must cope with all of them, including block counts large enough to force the descriptor +// arrays to grow. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_ForeignBlockSizes", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + + size_t n_ok = 0, n_refused = 0; + for (uint32_t es : {1u, 2u, 4u}) { + const size_t nelem = 300000; + const auto img = GenDetectorLike(nelem, es, 70000 + es); + + std::vector bsizes = {0, bshuf_default_block_size(es), 128, 512, 1024, 4096, + 8192 / es, 32768 / es, 49152 / es, 65536 / es}; + std::sort(bsizes.begin(), bsizes.end()); + bsizes.erase(std::unique(bsizes.begin(), bsizes.end()), bsizes.end()); + + for (size_t bs : bsizes) { + // bshuf_compress_lz4 emits only the blocks - the 12-byte container header (total size + // BE64, block size in BYTES BE32) is the caller's, which is why JFJochCompressor writes + // it by hand and every decoder starts reading at +12. Build it here the same way, or the + // container is not the one the pipeline reads. + const size_t block_elems = (bs != 0) ? bs : bshuf_default_block_size(es); + const size_t block_bytes = block_elems * es; + std::vector out(12 + bshuf_compress_lz4_bound(nelem, es, block_elems)); + PutBE64(out.data(), (uint64_t)nelem * es); + PutBE32(out.data() + 8, (uint32_t)block_bytes); + const int64_t clen = bshuf_compress_lz4(img.data(), out.data() + 12, nelem, es, block_elems); + REQUIRE(clen > 0); + out.resize(12 + (size_t)clen); + const size_t nblocks = (nelem + block_elems - 1) / block_elems; + + const CompressedImage image(out.data(), out.size(), nelem, 1, ModeOf(es, false), ALG); + std::vector cpu; + image.GetUncompressed(cpu); + { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); } + + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(nelem * es, stream); + CudaDevicePtr gpu_out(nelem * es); + cudaMemset(gpu_out.get(), 0xA5, nelem * es); + + std::string what; + try { + decoder.Decode(image, gpu_out.get()); + const cudaError_t sync = cudaStreamSynchronize(*stream); + if (sync != cudaSuccess) what = cudaGetErrorString(sync); + } catch (const std::exception &e) { what = e.what(); } + + if (!what.empty()) { + printf("[BSLZ4Fuzz] foreign block %6zu bytes, es=%u (%5zu blocks): GPU FAILED - %s\n", + block_bytes, es, nblocks, what.c_str()); + n_refused++; + cudaGetLastError(); // clear, so later cases are not confused by a stale error + continue; + } + + std::vector gpu(nelem * es); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess); + std::ostringstream l; l << "foreign block_bytes=" << block_bytes << " es=" << es; + const auto msg = DescribeMismatch(cpu.data(), gpu.data(), nelem * es, block_bytes, + nelem * es, l.str()); + if (!msg.empty()) FAIL(msg); + printf("[BSLZ4Fuzz] foreign block %6zu bytes, es=%u (%5zu blocks): byte-exact\n", + block_bytes, es, nblocks); + n_ok++; + } + } + printf("[BSLZ4Fuzz] foreign bitshuffle block sizes: %zu byte-exact, %zu refused by the GPU\n", + n_ok, n_refused); + SUCCEED(); +} + +// ============================================================================================= +// Targeted probe for a WARP-SYNCHRONY race. +// +// lz4_decode_blocks splits each copy across the 32 lanes with `for (i = lane; i < n; i += 32)`, so +// the trip count is lane-dependent and the lanes DIVERGE at the loop exit. There is no +// __syncwarp() before the next sequence's copy, and on Volta+ (this is sm_120) independent thread +// scheduling gives no implicit warp-lockstep guarantee. A lane that finished early can therefore +// run ahead and read a byte another lane has not written yet. +// +// The construction below makes exactly that overlap: a literal run of length L == 1 (mod 32), so +// the LAST byte is written by lane 0 on its SECOND trip while lanes 1..31 have already left the +// loop, followed immediately by a match at a tiny offset that reads that very byte. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_SyntheticWarpRaceProbe", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + + const uint32_t es = 4; + const size_t B = BlockElems(es), block_bytes = B * es; + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(block_bytes, stream); + CudaDevicePtr gpu_out(block_bytes); + + size_t configs = 0, total_decodes = 0, bad_decodes = 0; + + // Literal-run lengths chosen around multiples of the 32-lane width, and matches whose source + // is the tail the previous copy just produced. + for (uint32_t L : {1u, 2u, 17u, 31u, 32u, 33u, 34u, 63u, 64u, 65u, 96u, 97u}) { + for (uint32_t off : {1u, 2u, 3u, 4u, 5u, 8u, 15u, 16u, 17u, 31u, 32u, 33u}) { + for (uint32_t ml : {4u, 5u, 33u, 64u, 65u}) { + std::mt19937 rng(60000 + L * 977 + off * 31 + ml); + Lz4BlockBuilder bb; + { + std::vector prime(256); + for (auto &b : prime) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(prime); + bb.Match(128, 128); + } + // As many back-to-back [literals L][match off,ml] pairs as fit. + while (bb.OutputPos() + L + ml + 32 < block_bytes) { + std::vector lits(L); + for (auto &b : lits) b = (uint8_t)(rng() & 0xFF); + bb.AddLiterals(lits); + bb.Match(off, ml); + } + std::vector tail(block_bytes - bb.OutputPos()); + for (auto &b : tail) b = (uint8_t)(rng() & 0xFF); + bb.Finish(tail); + REQUIRE(bb.plain.size() == block_bytes); + + std::vector chk(block_bytes); + REQUIRE(LZ4_decompress_safe((const char *)bb.code.data(), (char *)chk.data(), + (int)bb.code.size(), (int)block_bytes) == (int)block_bytes); + { const bool lz4_ok = (chk == bb.plain); REQUIRE(lz4_ok); } + + const auto img = ImageFromShuffled(bb.plain, B, es, false); + const auto container = MakeContainer({bb.code}, block_bytes, block_bytes); + const CompressedImage image(container.data(), container.size(), B, 1, + CompressedImageMode::Uint32, ALG); + std::vector cpu; + image.GetUncompressed(cpu); + { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); } + + // A race is probabilistic, so decode the same block many times. + std::vector gpu(block_bytes); + for (int rep = 0; rep < 25; rep++) { + cudaMemset(gpu_out.get(), 0xA5, block_bytes); + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes, + cudaMemcpyDeviceToHost) == cudaSuccess); + total_decodes++; + if (memcmp(gpu.data(), cpu.data(), block_bytes) != 0) { + bad_decodes++; + if (bad_decodes == 1) { + std::ostringstream l; + l << "warp-race probe L=" << L << " offset=" << off << " matchlen=" << ml + << " rep=" << rep; + FAIL(DescribeMismatch(cpu.data(), gpu.data(), block_bytes, + block_bytes, block_bytes, l.str())); + } + } + } + configs++; + } + } + } + printf("[BSLZ4Fuzz] warp-synchrony probe: %zu configurations x 25 repeats = %zu decodes, " + "%zu corrupted\n", configs, total_decodes, bad_decodes); + SUCCEED(); +} + +// ============================================================================================= +// CORRUPT payloads. An HDF5 chunk carries no checksum unless Fletcher32 is on, so a bit flip on the +// wire or on disk reaches the decoder looking like a plausible container. The host decoder catches +// it - LZ4_decompress_safe insists the block decode to exactly its declared length - and the device +// decoder MUST do the same, because the buffers it decodes into are reused frame to frame: a block +// that stops early would otherwise leave the PREVIOUS image's bytes in place, and those bytes are +// the most significant byte-plane of the block, so they do not look like a missing corner. They +// look like real pixels several powers of two too bright. +// +// The invariant asserted here is the one that matters: whenever the CPU rejects a chunk, the GPU +// must reject it too. (The converse is not required - a flip can leave a stream that still decodes +// to the right length, and no LZ4 decoder can detect that.) +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_CorruptPayload", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + + const size_t nelem = 200000, es = 4; + const auto img = GenDetectorLike(nelem, es, 4242); + JFJochBitShuffleCompressor compressor(ALG); + const auto clean = compressor.Compress(img.data(), nelem, es); + + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(nelem * es, stream); + CudaDevicePtr gpu_out(nelem * es); + + std::mt19937 rng(777); + size_t both_threw = 0, both_ok = 0, gpu_missed = 0, cpu_only_ok = 0; + + for (int trial = 0; trial < 400; trial++) { + std::vector bad = clean; + // Corrupt inside the LZ4 payload area, never the 12-byte container header - a header flip is + // a different (and already covered) failure mode. + const size_t pos = 12 + rng() % (bad.size() - 12); + bad[pos] ^= (uint8_t)(1u << (rng() % 8)); + + const CompressedImage image(bad.data(), bad.size(), nelem, 1, CompressedImageMode::Uint32, ALG); + + bool cpu_threw = false; + std::vector cpu; + try { image.GetUncompressed(cpu); } catch (const std::exception &) { cpu_threw = true; } + + // Prime the device buffers with a recognisable pattern, so a block that quietly stopped + // early would show up as leftover poison rather than as plausible data. + cudaMemset(gpu_out.get(), 0x5A, nelem * es); + bool gpu_threw = false; + try { + decoder.Decode(image, gpu_out.get()); + REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess); + decoder.ThrowIfDecodeFailed(); + } catch (const std::exception &) { gpu_threw = true; } + + // Whatever happened, it must not have been a device fault: a corrupt length must be clamped, + // never walked off the end of a buffer. + REQUIRE(cudaGetLastError() == cudaSuccess); + + if (cpu_threw && gpu_threw) both_threw++; + else if (!cpu_threw && !gpu_threw) { + both_ok++; + std::vector gpu(nelem * es); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess); + const auto msg = DescribeMismatch(cpu.data(), gpu.data(), nelem * es, + DefaultBlockBytes(), nelem * es, + "corrupt payload both-accepted"); + if (!msg.empty()) FAIL(msg); // if both accept it, they must agree on the result + } + else if (cpu_threw && !gpu_threw) { + gpu_missed++; + // Whatever the GPU accepted, it must at least be a COMPLETE decode. The poison check is + // the important one: a block that stopped early would leave 0x5A here, and in production + // those bytes are the previous frame rather than poison. + std::vector gpu(nelem * es); + REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess); + size_t poison = 0; + for (size_t i = 0; i < nelem * es; i++) if (gpu[i] == 0x5A) poison++; + INFO("corrupt flip at byte " << pos << " left " << poison << " undecoded bytes"); + REQUIRE(poison == 0); + } + else cpu_only_ok++; + } + + printf("[BSLZ4Fuzz] corrupt payload, 400 single-bit flips: %zu both rejected, %zu both accepted " + "(and agreed), %zu GPU MISSED what the CPU caught, %zu GPU stricter than CPU\n", + both_threw, both_ok, gpu_missed, cpu_only_ok); + + // The device decoder reproduces the reference's length, bounds and parsing-restriction checks, + // which accounts for the large majority of what LZ4_decompress_safe rejects. It does NOT + // reproduce the reference's fast-loop/safe-loop selection exactly, so a small tail of corrupt + // streams still decodes here that the host would refuse. Those are complete decodes - the block + // is fully written, asserted above - differing from the true image in a handful of bytes, not + // partial decodes leaking a previous frame, which is the failure mode that actually matters. + // The bound is a regression guard: it must not silently get worse. + CHECK(gpu_missed * 100 <= 400 * 2); // <= 2% of single-bit flips +} + +// Malformed CONTAINERS - the parts the host can reject before any work is queued. Each of these is +// rejected by the CPU decoder, so the device decoder must reject them too rather than decode +// something plausible-looking out of uninitialised memory. +TEST_CASE("BSLZ4Fuzz_MalformedContainer", "[BSLZ4Fuzz]") { + if (get_gpu_count() == 0) SKIP("No CUDA GPU present"); + + const size_t nelem = 100000, es = 4; + const auto img = GenDetectorLike(nelem, es, 99); + JFJochBitShuffleCompressor compressor(ALG); + const auto clean = compressor.Compress(img.data(), nelem, es); + + auto stream = std::make_shared(); + BSLZ4DecoderGPU decoder(nelem * es, stream); + CudaDevicePtr gpu_out(nelem * es); + + auto expect_throw = [&](const std::vector &chunk, const char *what) { + const CompressedImage image(chunk.data(), chunk.size(), nelem, 1, CompressedImageMode::Uint32, ALG); + bool threw = false; + try { + decoder.Decode(image, gpu_out.get()); + cudaStreamSynchronize(*stream); + decoder.ThrowIfDecodeFailed(); + } catch (const std::exception &) { threw = true; } + cudaGetLastError(); + INFO(what); + CHECK(threw); + }; + + { // Block size that is not a multiple of 8 elements: bitshuffle transposes 8 at a time, so the + // un-transpose would leave part of every block unwritten. + auto bad = clean; PutBE32(bad.data() + 8, 4 * es); + expect_throw(bad, "block size not a multiple of 8 elements"); + } + { // A one-element block size implies millions of blocks the chunk cannot possibly hold. Must be + // rejected on the arithmetic, BEFORE it becomes a several-hundred-MB pinned allocation. + auto bad = clean; PutBE32(bad.data() + 8, es); + expect_throw(bad, "block size implies more blocks than the chunk holds"); + } + { auto bad = clean; PutBE64(bad.data(), (uint64_t)nelem * es + 8); + expect_throw(bad, "declared total size disagrees with the image"); } + { auto bad = clean; PutBE32(bad.data() + 8, 0); + expect_throw(bad, "zero block size"); } + { auto bad = clean; bad.resize(11); + expect_throw(bad, "shorter than the container header"); } + { auto bad = clean; bad.resize(clean.size() / 2); + expect_throw(bad, "truncated mid-stream"); } + { auto bad = clean; bad.push_back(0); bad.push_back(0); + expect_throw(bad, "trailing bytes after the last block"); } + { // A block length field larger than what is left in the chunk. + auto bad = clean; PutBE32(bad.data() + 12, (uint32_t)clean.size()); + expect_throw(bad, "block length runs past the chunk"); + } + { auto bad = clean; PutBE32(bad.data() + 12, 0); + expect_throw(bad, "zero-length block"); } +} + +// ============================================================================================= +// Coverage report - what the fuzz cases above actually made LZ4 emit. +// ============================================================================================= +TEST_CASE("BSLZ4Fuzz_ZZ_CoverageReport", "[BSLZ4Fuzz]") { + std::lock_guard lock(g_stats_mtx); + const auto &s = g_stats; + printf("\n[BSLZ4Fuzz] ===== LZ4 CODE-PATH COVERAGE (production-compressor cases) =====\n"); + printf(" images parsed .................. %zu\n", s.images); + printf(" bitshuffle blocks .............. %zu\n", s.blocks); + printf(" LZ4 match sequences ............ %zu\n", s.sequences); + printf(" litlen field == 15 (extended) .. %zu\n", s.lit_ge15); + printf(" litlen 255-extension CHAIN ..... %zu\n", s.lit_ext255); + printf(" max literal run ................ %zu bytes\n", s.max_litlen); + printf(" matchlen field == 15 (extended) %zu\n", s.ml_ge19); + printf(" matchlen 255-extension CHAIN ... %zu\n", s.ml_ext255); + printf(" max match length ............... %zu bytes\n", s.max_matchlen); + printf(" offset >= matchlen (non-overlap) %zu\n", s.nonoverlap); + printf(" offset < matchlen (OVERLAP) .. %zu\n", s.overlap); + printf(" offset == matchlen ............. %zu\n", s.off_eq_ml); + printf(" offset == matchlen + 1 ......... %zu\n", s.off_ml_plus1); + printf(" offset == matchlen - 1 ......... %zu\n", s.off_ml_minus1); + printf(" smallest offset seen ........... %u\n", s.min_offset); + printf(" offsets <= 64 seen: "); + for (auto &kv : s.offsets_le64) printf("%u(x%zu) ", kv.first, kv.second); + printf("\n[BSLZ4Fuzz] =============================================================\n\n"); + SUCCEED(); +} + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b5188c60..494b3ad5 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -66,6 +66,8 @@ ADD_EXECUTABLE(jfjoch_test ImageSpotFinderCPUTest.cpp ImageSpotFinderGPUTest.cpp BSLZ4DecoderGPUTest.cpp + BSLZ4DecoderGPUFuzzTest.cpp + ImagePreprocessorGPUFusedTest.cpp AdaptiveSpotFinderGPUTest.cpp SpotExtractorGPUParityTest.cpp CalcBraggPredictionTest.cpp diff --git a/tests/ImagePreprocessorGPUFusedTest.cpp b/tests/ImagePreprocessorGPUFusedTest.cpp new file mode 100644 index 00000000..0400591d --- /dev/null +++ b/tests/ImagePreprocessorGPUFusedTest.cpp @@ -0,0 +1,167 @@ +// 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 "../common/PixelMask.h" +#include "../compression/JFJochCompressor.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" + +// The device-decode path does NOT decompress into a buffer and then preprocess it: one kernel +// un-transposes the bitshuffle blocks and applies the mask, the error marker, the saturation cap and +// the statistics as it goes, so the decompressed image never exists. That is a different code path +// from the host-upload one, not a reordering of it, and the thing it has to reproduce is the whole +// observable output - every preprocessed pixel AND every counter - against the CPU preprocessor fed +// the host-decompressed image. +// +// Masked, error and saturated pixels are the interesting part: their priority (masked > error > +// saturated) and their sentinel outputs (INT32_MIN / INT32_MIN / INT32_MAX) are decided in the fused +// kernel now, so the image below deliberately contains all three, and the mask deliberately covers +// some of them. +namespace { + +DiffractionExperiment MakeExperiment(size_t saturation) { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + return x; +} + +template +std::vector MakeImage(size_t npixels, T err_value, uint32_t seed) { + std::mt19937 rng(seed); + std::vector img(npixels, 0); + // Sparse background with long runs, so LZ4 produces overlapping matches. + for (size_t i = npixels / 4; i < npixels / 2; i++) + img[i] = static_cast(rng() % 11); + // Bright spots, some above any plausible saturation cap. + for (size_t s = 0; s < 500; s++) { + const size_t c = rng() % npixels; + for (size_t d = 0; d < 5 && c + d < npixels; d++) + img[c + d] = static_cast(30000 + (rng() % 5000)); + } + // Explicit error markers, scattered. + for (size_t s = 0; s < 300; s++) + img[rng() % npixels] = err_value; + return img; +} + +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; +} + +// One element size end to end: compress, decode+preprocess on the device, and compare against the +// host decompression fed through the CPU preprocessor. +template +void CheckFusedMatchesCPU(CompressedImageMode mode, T err_value, uint32_t seed) { + DiffractionExperiment x = MakeExperiment(32000); + const size_t npixels = x.GetPixelsNum(); + + PixelMask mask(x); + // Mask a deterministic scatter of pixels, so masked-vs-error-vs-saturated priority is exercised + // rather than assumed. + auto &m = const_cast &>(mask.GetMask()); + for (size_t i = 0; i < npixels; i += 997) m[i] = 1; + for (size_t i = 13; i < npixels; i += 4001) m[i] = 1; + + const auto img = MakeImage(npixels, err_value, seed); + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + const std::vector compressed = compressor.Compress(img); + const CompressedImage image(compressed.data(), compressed.size(), + x.GetXPixelsNum(), x.GetYPixelsNum(), mode, + CompressionAlgorithm::BSHUF_LZ4); + REQUIRE(BSLZ4DecoderGPU::Supports(image)); + + // Reference: host decompression + CPU preprocessing. + ImagePreprocessorCPU cpu_pre(x, mask); + ImagePreprocessorBuffer cpu_buf(npixels); + std::vector decompression_buffer; + const uint8_t *raw = image.GetUncompressedPtr(decompression_buffer); + const ImageStatistics cpu_stats = cpu_pre.Analyze(cpu_buf, raw, mode); + + // Under test: compressed chunk straight to the device, decoded and preprocessed in one pass. + auto stream = std::make_shared(); + ImagePreprocessorGPU gpu_pre(x, mask, stream, /*copy_image_to_host=*/true); + ImagePreprocessorBufferGPU gpu_buf(npixels); + ImageStatistics gpu_stats{}; + REQUIRE(gpu_pre.AnalyzeCompressed(gpu_buf, image, gpu_stats)); + + INFO("mode " << static_cast(mode)); + CHECK(SameStats(cpu_stats, gpu_stats)); + + size_t ndiff = 0, first = 0; + for (size_t i = 0; i < npixels; i++) { + if (gpu_buf[i] != cpu_buf[i]) { + if (ndiff == 0) first = i; + ndiff++; + } + } + INFO("first differing pixel " << first << " cpu " << cpu_buf[first] << " gpu " << gpu_buf[first] + << " of " << ndiff << " differing"); + CHECK(ndiff == 0); +} + +} // namespace + +TEST_CASE("ImagePreprocessorGPU_FusedDecodeMatchesCPU", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + CheckFusedMatchesCPU(CompressedImageMode::Uint32, UINT32_MAX, 1); + CheckFusedMatchesCPU(CompressedImageMode::Uint16, UINT16_MAX, 2); + CheckFusedMatchesCPU(CompressedImageMode::Int32, INT32_MIN, 3); + CheckFusedMatchesCPU(CompressedImageMode::Int16, INT16_MIN, 4); +} + +// The host-upload path must keep producing exactly what it did - it is still what every non-LZ4 +// image takes - so the two entry points are held against each other on the same frame. +TEST_CASE("ImagePreprocessorGPU_FusedMatchesHostUpload", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + DiffractionExperiment x = MakeExperiment(32000); + const size_t npixels = x.GetPixelsNum(); + PixelMask mask(x); + auto &m = const_cast &>(mask.GetMask()); + for (size_t i = 0; i < npixels; i += 1301) m[i] = 1; + + const auto img = MakeImage(npixels, UINT32_MAX, 77); + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + const std::vector compressed = compressor.Compress(img); + const CompressedImage image(compressed.data(), compressed.size(), + x.GetXPixelsNum(), x.GetYPixelsNum(), CompressedImageMode::Uint32, + CompressionAlgorithm::BSHUF_LZ4); + + auto stream = std::make_shared(); + ImagePreprocessorGPU pre(x, mask, stream, /*copy_image_to_host=*/true); + + ImagePreprocessorBufferGPU fused_buf(npixels); + ImageStatistics fused_stats{}; + REQUIRE(pre.AnalyzeCompressed(fused_buf, image, fused_stats)); + + // Same engine, same frame, but decompressed on the host and uploaded. + std::vector decompression_buffer; + const uint8_t *raw = image.GetUncompressedPtr(decompression_buffer); + ImagePreprocessorBufferGPU upload_buf(npixels); + const ImageStatistics upload_stats = pre.Analyze(upload_buf, raw, CompressedImageMode::Uint32); + + CHECK(SameStats(fused_stats, upload_stats)); + size_t ndiff = 0; + for (size_t i = 0; i < npixels; i++) if (fused_buf[i] != upload_buf[i]) ndiff++; + CHECK(ndiff == 0); + + // Decoding on the device replaced a host decompression, so the cost is still reported as one. + CHECK(pre.GetLastDecompressionTime_s() > 0.0f); +} + +#endif