From f36cd88795c56db983540d7c2c4d549022897dad Mon Sep 17 00:00:00 2001 From: jungfrau Date: Tue, 18 Aug 2026 18:41:53 -0400 Subject: [PATCH] Give a 16-bit image a saturation code when its limit is 65534 The narrow encoding picked NARROW_BAD as the saturation code whenever the saturation limit was above UINT16_MAX - 2. That is one too strict. A real value is STRICTLY below the limit, so a limit of UINT16_MAX - 1 still leaves UINT16_MAX - 1 free to be the code; only a limit of the whole range has nothing to spare, and that is exactly the case where no pixel can be saturated, because the "is error" test claims UINT16_MAX first. At a limit of exactly 65534 the old condition therefore stored a saturated pixel as the masked code, and it widened back to INT32_MIN instead of INT32_MAX. Masked and saturated are not interchangeable: the strong-pixel search flags a saturated pixel unconditionally and a masked one never, so the overloaded core of the strongest spots dropped out of the strong-pixel mask. Bragg integration treats the two alike and the image statistics are taken from the raw value, so nothing downstream of those moved - which is why a battery over 24 crystals showed nothing. That value is not a corner case. GetByteDepthImage()-driven writing stores saturation_value = GetSaturationLimit() - 1 and the readers take it back as-is, so a 16-bit acquisition whose detector cutoff is at or above the full range comes back with a limit of exactly 65534. The one 16-bit dataset in the rotation test set declares 11963, which is safe, so it could not have caught this. Found by review, not by testing, because nothing tested the narrow path at all: every GPU test writes into the wide buffer directly and never asks for the narrow one, and every preprocessor test copies the image back to the host, which forces the wide path. So the test comes with the fix. It runs the wide and the narrow preprocessor over the same synthetic frame - values around each boundary, masked on a stride coprime with the value cycle so every value appears both masked and unmasked - across saturation limits of none, 5000, 0xFFFD, 0xFFFE and 0xFFFF, on both the host-upload and the device-decode entry point, and compares the statistics and every pixel. Against the old condition it fails with 625860 differing pixels; against this one it passes. Also makes the header self-contained: it uses __host__/__device__ and the CUDA vector types and only compiled because every includer happened to pull in cuda_runtime.h first. Co-Authored-By: Claude Opus 5 --- docs/CHANGELOG.md | 1 + .../image_preprocessing/PreprocessedPixel.h | 10 +- tests/CMakeLists.txt | 1 + tests/ImagePreprocessorGPUNarrowTest.cpp | 141 ++++++++++++++++++ 4 files changed, 150 insertions(+), 3 deletions(-) create mode 100644 tests/ImagePreprocessorGPUNarrowTest.cpp diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 81fc3068..a9d92885 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### Unreleased +* rugnux: **fixed** a case in the new 16-bit image path where a saturated pixel was recorded as masked. It needed a detector saturation limit of exactly 65534 - which is what a 16-bit file written by jfjoch and read back reports, because the saturation value is written one count low and read back as-is - and it removed the core of the strongest spots from the strong-pixel mask. Integrated intensities and image statistics were unaffected. * jfjoch_broker: a DECTRIS detector no longer claims a 16-bit image and readout depth before anything has asked the detector. Both are now unset until the SIMPLON configuration or the incoming stream supplies them. **This makes `image_format_settings.bit_depth_image` work on DECTRIS detectors** - it was documented as controlling the output depth, but the detector's assumed default took priority and it was silently discarded. Where nothing supplies a depth the result is unchanged (16-bit). * **rugnux: 16-bit images keep their width all the way through the GPU pipeline** instead of being widened to 32-bit as soon as they are decoded. Every per-pixel pass over a frame - the ring statistics, the strong-pixel search, spot extraction, azimuthal and ROI integration, Bragg integration - then moves half as many bytes. This is the mode a fast acquisition reads out in, and the gain grows with the detector: on the 2.5M-pixel dataset in the test set the image loop drops 2%, and it is worth proportionally more on a large detector, where per-pixel work is most of the loop. Merged results are unchanged. * **rugnux: the first pass of the rotation two-pass no longer merges, reports or writes anything.** It exists to measure the detector geometry, the goniometer rotation scale, the mosaicity and the space group, and the second pass makes the merged result again at the refined geometry - so the first pass now stops once it has those. **The `_01_*` files are no longer written.** A battery of 24 rotation crystals drops from 6m52s to 6m17s and no merged result changes. diff --git a/image_analysis/image_preprocessing/PreprocessedPixel.h b/image_analysis/image_preprocessing/PreprocessedPixel.h index fb1cff82..2e53b826 100644 --- a/image_analysis/image_preprocessing/PreprocessedPixel.h +++ b/image_analysis/image_preprocessing/PreprocessedPixel.h @@ -5,6 +5,8 @@ // CUDA only: the engines that read a preprocessed image on the device all compile as .cu. The CPU // engines read ImagePreprocessorBuffer directly, which is int32 whatever the source was. +#include + #include #include @@ -39,10 +41,12 @@ namespace preprocessed_pixel { constexpr uint16_t NARROW_BAD = UINT16_MAX; // Which code means "saturated", given the saturation limit this image was preprocessed with. - // Falls back to NARROW_BAD when the limit is the whole range, which is exactly the case where - // no pixel can be saturated - so the code is then unreachable rather than ambiguous. + // A real value is strictly below the limit, so the limit may be as high as UINT16_MAX - 1 and + // still leave UINT16_MAX - 1 free to be the code. Only a limit of the whole range falls back to + // NARROW_BAD, and that is exactly the case where no pixel can be saturated - the "is error" test + // claims UINT16_MAX first - so the code is then unreachable rather than ambiguous. constexpr uint16_t NarrowSaturatedCode(int64_t saturation_limit) { - return saturation_limit <= UINT16_MAX - 2 ? static_cast(UINT16_MAX - 1) : NARROW_BAD; + return saturation_limit <= UINT16_MAX - 1 ? static_cast(UINT16_MAX - 1) : NARROW_BAD; } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4439900f..1ebe0a50 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -72,6 +72,7 @@ ADD_EXECUTABLE(jfjoch_test BSLZ4DecoderGPUTest.cpp BSLZ4DecoderGPUFuzzTest.cpp ImagePreprocessorGPUFusedTest.cpp + ImagePreprocessorGPUNarrowTest.cpp AdaptiveSpotFinderGPUTest.cpp SpotExtractorGPUParityTest.cpp CalcBraggPredictionTest.cpp diff --git a/tests/ImagePreprocessorGPUNarrowTest.cpp b/tests/ImagePreprocessorGPUNarrowTest.cpp new file mode 100644 index 00000000..72d280bb --- /dev/null +++ b/tests/ImagePreprocessorGPUNarrowTest.cpp @@ -0,0 +1,141 @@ +// 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 "../common/PixelMask.h" +#include "../compression/JFJochCompressor.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" + +// A 16-bit unsigned source keeps its width through the GPU pipeline (PreprocessedPixel.h): masked +// and saturated pixels ride on two reserved codes at the top of the range and are widened again on +// load. That has to give back exactly what the 32-bit path gives, for every value a 16-bit detector +// can produce and for every saturation limit it can declare - so both are run over the same frame +// and compared pixel for pixel, at both entry points (host upload and device decode). +namespace { + +// The values the encoding turns on: the two reserved codes, their neighbours, and the pixels either +// side of the saturation limit. +std::vector InterestingValues(int64_t saturation) { + std::vector v{0, 1, 2, 100, 0xFFFC, 0xFFFD, 0xFFFE, 0xFFFF}; + for (int64_t d = -2; d <= 2; d++) { + const int64_t s = saturation + d; + if (s >= 0 && s <= UINT16_MAX) + v.push_back(static_cast(s)); + } + return v; +} + +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; +} + +// Gather() reads through PixelView, so it reports the image in the pipeline's int32 convention +// whatever width it is stored in - which is the only way to read a narrow image back. +std::vector ReadBack(const ImagePreprocessorBuffer &buffer) { + std::vector out(buffer.size()); + std::vector index; + std::vector values; + for (size_t i = 0; i < buffer.size(); i += UINT16_MAX) { + const size_t n = std::min(UINT16_MAX, buffer.size() - i); + index.resize(n); + for (size_t k = 0; k < n; k++) + index[k] = static_cast(i + k); + buffer.Gather(index, values); + std::copy(values.begin(), values.end(), out.begin() + i); + } + return out; +} + +void CheckNarrowMatchesWide(std::optional saturation, bool device_decode) { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + // A 16-bit unsigned container, which is the only source the narrow path is taken for. + x.BitDepthImage(16).PixelSigned(false); + x.Detector().SaturationLimit(saturation); + const size_t npixels = x.GetPixelsNum(); + + const std::vector values = InterestingValues(x.GetSaturationLimit()); + std::vector img(npixels); + for (size_t i = 0; i < npixels; i++) + img[i] = values[i % values.size()]; + + // 11 is coprime with the number of values, so every value appears both masked and unmasked. + PixelMask mask(x); + auto &m = const_cast &>(mask.GetMask()); + for (size_t i = 3; i < npixels; i += 11) + m[i] = 1; + + auto stream = std::make_shared(); + // The 32-bit path, which is what every consumer saw before, is the reference. + ImagePreprocessorGPU wide_pre(x, mask, stream, /*copy_image_to_host=*/true); + ImagePreprocessorBufferGPU wide_buffer(npixels); + // Telling the engine nothing has to come back to the host is what lets the image keep its width. + ImagePreprocessorGPU narrow_pre(x, mask, stream, /*copy_image_to_host=*/false); + ImagePreprocessorBufferGPU narrow_buffer(npixels, /*host_mirror=*/false); + + ImageStatistics wide_stats{}; + ImageStatistics narrow_stats{}; + if (device_decode) { + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + const std::vector compressed = compressor.Compress(img); + const CompressedImage image(compressed.data(), compressed.size(), + x.GetXPixelsNum(), x.GetYPixelsNum(), + CompressedImageMode::Uint16, CompressionAlgorithm::BSHUF_LZ4); + REQUIRE(wide_pre.AnalyzeCompressed(wide_buffer, image, wide_stats)); + REQUIRE(narrow_pre.AnalyzeCompressed(narrow_buffer, image, narrow_stats)); + } else { + const auto *raw = reinterpret_cast(img.data()); + wide_stats = wide_pre.Analyze(wide_buffer, raw, CompressedImageMode::Uint16); + narrow_stats = narrow_pre.Analyze(narrow_buffer, raw, CompressedImageMode::Uint16); + } + + INFO("saturation limit " << x.GetSaturationLimit() << (device_decode ? " device decode" : " host upload")); + REQUIRE(narrow_buffer.IsNarrow()); + REQUIRE_FALSE(wide_buffer.IsNarrow()); + CHECK(SameStats(wide_stats, narrow_stats)); + + const std::vector narrow_read = ReadBack(narrow_buffer); + size_t ndiff = 0; + size_t first = 0; + for (size_t i = 0; i < npixels; i++) { + if (narrow_read[i] != wide_buffer[i]) { + if (ndiff == 0) first = i; + ndiff++; + } + } + INFO("first differing pixel " << first << " source value " << img[first] + << (m[first] ? " (masked)" : "") + << " wide " << wide_buffer[first] << " narrow " << narrow_read[first] + << " of " << ndiff << " differing"); + CHECK(ndiff == 0); +} + +} // namespace + +TEST_CASE("ImagePreprocessorGPU_NarrowMatchesWide", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + for (bool device_decode : {false, true}) { + // No declared limit: the limit is the whole 16-bit range, and nothing can saturate. + CheckNarrowMatchesWide(std::nullopt, device_decode); + // A limit with plenty of room below the reserved codes - the ordinary case. + CheckNarrowMatchesWide(5000, device_decode); + // The largest limit that still leaves a code of its own, and the two above it. + CheckNarrowMatchesWide(UINT16_MAX - 2, device_decode); + CheckNarrowMatchesWide(UINT16_MAX - 1, device_decode); + CheckNarrowMatchesWide(UINT16_MAX, device_decode); + } +} + +#endif