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