From 15af14357d2cef1d2955788aa0a76e0ce497fedb Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 7 Nov 2026 13:41:17 +0100 Subject: [PATCH] pixel_sum on the CPU and GPU path, under the device mask, and read back from HDF5 pixel_sum existed only on the FPGA path, so a DECTRIS master file carried an all-zero /entry/image/pixel_sum. The CPU and GPU image preprocessors now compute it in the pass that takes the minimum and maximum, as a signed 64-bit sum, and MXAnalysisWithoutFPGA reports it: Lite, live stills, live rotation and rugnux. On the simulated card it equals the card's own sum, for 16-bit and for summed 32-bit frames. The sum excludes error and saturated pixels and the device mask (detector and calibration). The hot pixels and the beam-stop shadow a run finds on its own frames do not change it, though they stay masked in the image and in the other statistics, so the number does not depend on what a run decided. Those two have bits of their own (PixelMask::RunMaskBits); the byte-per-pixel mask keeps the difference, 1 for the device mask and 2 for the run's alone, and the CPU preprocessor gets a packed mask of the second. Live stills reports this sum of the assembled image in place of the card's, and with the run's mask out of it the two agree on the test frames. Frames a live rotation run does not reach stay 0 in the master file, as before, and have no dataset in the data files. The reader now reads pixel_sum into the dataset and the per-image message; a file without it opens as before. Co-Authored-By: Claude Opus 5.5 --- common/PixelMask.cpp | 13 +- common/PixelMask.h | 12 +- docs/CBOR.md | 4 +- docs/FPGA_DATA_ANALYSIS.md | 5 + docs/GPU_DATA_ANALYSIS.md | 6 +- docs/HDF5.md | 4 + image_analysis/MXAnalysisWithoutFPGA.cpp | 1 + .../image_preprocessing/ImagePreprocessor.h | 3 + .../ImagePreprocessorCPU.cpp | 23 +- .../ImagePreprocessorCPU.h | 1 + .../ImagePreprocessorGPU.cu | 25 +- reader/HDF5MetadataSource.cpp | 9 +- reader/JFJochMultiReader.cpp | 1 + reader/JFJochReaderDataset.h | 1 + tests/CMakeLists.txt | 1 + tests/DirectHostReadGPUTest.cpp | 2 +- tests/FPGAIntegrationTest.cpp | 78 +++++ tests/ImagePreprocessorGPUFusedTest.cpp | 2 +- tests/ImagePreprocessorPixelSumTest.cpp | 266 ++++++++++++++++++ tests/JFJochReaderTest.cpp | 4 + tests/JFJochReceiverLiteTest.cpp | 61 +++- tests/StoredImageDecodeGPUTest.cpp | 2 +- 22 files changed, 496 insertions(+), 28 deletions(-) create mode 100644 tests/ImagePreprocessorPixelSumTest.cpp diff --git a/common/PixelMask.cpp b/common/PixelMask.cpp index 5c413b4b3..34db40f6b 100644 --- a/common/PixelMask.cpp +++ b/common/PixelMask.cpp @@ -57,11 +57,14 @@ void PixelMask::UpdateDerived(const DiffractionExperiment &experiment) { void PixelMask::UpdateBinaryMask() { binary_mask.resize(mask.size()); for (size_t i = 0; i < mask.size(); i++) - binary_mask[i] = (mask[i] != 0); + binary_mask[i] = mask[i] == 0 ? 0 : (mask[i] & ~RunMaskBits) != 0 ? BinaryMaskDevice : BinaryMaskRunOnly; binary_mask_checksum = TableChecksum(binary_mask.data(), binary_mask.size()); packed_mask.assign((mask.size() + 31) / 32, 0); - for (size_t i = 0; i < mask.size(); i++) - packed_mask[i / 32] |= static_cast(binary_mask[i]) << (i % 32); + packed_run_mask.assign((mask.size() + 31) / 32, 0); + for (size_t i = 0; i < mask.size(); i++) { + packed_mask[i / 32] |= static_cast(binary_mask[i] != 0) << (i % 32); + packed_run_mask[i / 32] |= static_cast(binary_mask[i] == BinaryMaskRunOnly) << (i % 32); + } } void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) { @@ -134,6 +137,10 @@ const std::vector &PixelMask::GetPackedMask() const { return packed_mask; } +const std::vector &PixelMask::GetPackedRunMask() const { + return packed_run_mask; +} + uint64_t PixelMask::GetBinaryMaskChecksum() const { return binary_mask_checksum; } diff --git a/common/PixelMask.h b/common/PixelMask.h index 7e4dff850..2b2df286a 100644 --- a/common/PixelMask.h +++ b/common/PixelMask.h @@ -20,8 +20,8 @@ struct PixelMaskStatistics { class PixelMask { std::vector mask; std::vector raw_mask; - // One byte per pixel, 1 where the pixel is masked at all - the form the GPU image preprocessor - // uploads - and the checksum the shared device-table cache keys on (CudaSharedTables.h). Both are + // One byte per pixel, non-zero where the pixel is masked at all (see BinaryMaskDevice) - the form + // the GPU image preprocessor uploads - and the checksum the shared device-table cache keys on (CudaSharedTables.h). Both are // pure functions of `mask`, and an analysis engine is built per worker per pass, so they are // derived here once instead of in every one of those engines. std::vector binary_mask; @@ -29,6 +29,8 @@ class PixelMask { // The same flag packed 32 pixels to a word (pixel i is bit i % 32 of word i / 32) - the form the CPU // image preprocessor reads, derived here for the same reason. std::vector packed_mask; + // Packed the same way: the pixels masked only by what a run found on its own frames. + std::vector packed_run_mask; uint32_t LoadMask(const std::vector& mask, uint8_t bit); // Everything that follows from the mask, recomputed wherever the mask changes. @@ -49,6 +51,11 @@ public: constexpr static const uint8_t HotPixelBit = 10; constexpr static const uint8_t ChipGapPixelBit = 31; constexpr static const uint8_t ModuleEdgePixelBit = 30; + // What a run determines on its own frames, as against the detector's and the calibration's mask. + constexpr static const uint32_t RunMaskBits = (1u << BeamStopPixelBit) | (1u << HotPixelBit); + // The values of GetBinaryMask(): masked by the device mask, or only by what the run determined. + constexpr static const uint8_t BinaryMaskDevice = 1; + constexpr static const uint8_t BinaryMaskRunOnly = 2; PixelMask(); explicit PixelMask(size_t width, size_t height); @@ -74,6 +81,7 @@ public: [[nodiscard]] const std::vector &GetMask() const; [[nodiscard]] const std::vector &GetBinaryMask() const; [[nodiscard]] const std::vector &GetPackedMask() const; + [[nodiscard]] const std::vector &GetPackedRunMask() const; [[nodiscard]] uint64_t GetBinaryMaskChecksum() const; [[nodiscard]] std::vector GetUserMask(const DiffractionExperiment& experiment) const; diff --git a/docs/CBOR.md b/docs/CBOR.md index c3950eb56..852bc6f52 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -220,7 +220,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | receiver_buf_in_preparation | uint64 | Receiver internal number of buffer locations currently in processing | | | | storage_cell | uint64 | Storage cell number | | | | saturated_pixel_count | uint64 | Saturated pixel count | | | -| pixel_sum | uint64 | Sum of all pixels, excl. error and saturation | | | +| pixel_sum | int64 | Sum of pixels, excl. device-masked, error and saturated (negative values included) | | | | error_pixel_count | uint64 | Error pixel count | | | | strong_pixel_count | uint64 | Strong pixel count (first stage of spot finding) | | | | min_viable_pixel_value | int64 | Minimal pixel value, excl. error and saturation | | | @@ -345,7 +345,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | integrated_reflections | Array(int32) | Per-image count of integrated reflections | | | indexed_lattice_count | Array(int32) | Per-image count of indexed lattices | | | niggli_class | Array(uint8) | Per-image Niggli class identifier for indexed images; 0 if unavailable | | -| pixel_sum | Array(int64) | Per-image sum of all valid pixels, excluding error/saturated pixels | | +| pixel_sum | Array(int64) | Per-image sum of pixels, excl. device-masked, error and saturated (0 for an image with no analysis) | | | image_scale_mosaicity | Array(float) | Scaling result: Image scale mosaicity \[deg\] | | | image_scale_cc | Array(float) | Scaling result: Image scale CC | | | merged | object (optional) | Merged reflections of an in-broker analysis run, written by the writer as `.mtz`. An MTZ table (layout below) | | diff --git a/docs/FPGA_DATA_ANALYSIS.md b/docs/FPGA_DATA_ANALYSIS.md index 2afd2c48c..fcb0ba202 100644 --- a/docs/FPGA_DATA_ANALYSIS.md +++ b/docs/FPGA_DATA_ANALYSIS.md @@ -77,6 +77,11 @@ The following statistics are collected for each module: Valid pixels are not masked, not saturated, not error pixels. +These are not specific to the FPGA path: the CPU and GPU image preprocessors of the shared +`image_analysis/` library collect the same statistics per image, for detectors driven through +DECTRIS SIMPLON (e.g. EIGER) and for rugnux. There the sum excludes only the device mask (detector and +calibration), not the hot pixels and the beam-stop shadow a run finds on its own frames. + ## Square root compression Jungfraujoch FPGA includes lossy compression preserving counting statistic properties of X-ray image, while reducing bit width of an image. Scheme was described in [Wakonig et al., J. Appl. Cryst., 53, 574-586, 2020](https://doi.org/10.1107/S1600576720001776). diff --git a/docs/GPU_DATA_ANALYSIS.md b/docs/GPU_DATA_ANALYSIS.md index 9c0e93706..bec91a4af 100644 --- a/docs/GPU_DATA_ANALYSIS.md +++ b/docs/GPU_DATA_ANALYSIS.md @@ -30,9 +30,9 @@ run of a GPU build. and widened to int32 on the card, inside preprocessing. [§0](CPU_DATA_ANALYSIS_IMAGE.md#0-getting-the-image-onto-the-gpu-device-side-decoding) - **Preprocessing** (`ImagePreprocessorGPU.cu`, `ImagePreprocessorBufferGPU.cu`) - applies the pixel - mask and flags error and saturated pixels, counting each, and keeps the preprocessed image on the - device for every engine below; on request it gathers the pixels of the strong spots, or small - squares around given positions, back to the host. Integer, so identical to the host. + mask and flags error and saturated pixels, counting each and summing the rest, and keeps the + preprocessed image on the device for every engine below; on request it gathers the pixels of the + strong spots, or small squares around given positions, back to the host. Integer, so identical to the host. [§1.6](CPU_DATA_ANALYSIS_IMAGE.md#16-defective-pixels) ## Per image: profiles and spots diff --git a/docs/HDF5.md b/docs/HDF5.md index 3e2c00ffe..40f7f221a 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -542,6 +542,10 @@ every pixel belonging to that ROI, so an ROI's footprint can be recovered exactl `[n_images]` vectors: `max_value`, `min_value` (viable min/max, excluding error/saturated pixels), `error_pixels`, `saturated_pixels`, `pixel_sum`. Surfaced in the master file under `/entry/image`. +`pixel_sum` is the signed sum of the pixels that are not error, not saturated and not under the +device mask (detector and calibration); a mask rugnux determines on the run's own frames (hot pixels, +beam-stop shadow) does not change it. The FPGA computes it per module, and CPU/GPU analysis (e.g. +EIGER) per image, on the assembled image, to the same definition. ### 4.6 `/entry/profiling` — per-image timing diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index f6bdee7e1..3c228fba2 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -380,6 +380,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, output.min_viable_pixel_value = ret.min_value; output.error_pixel_count = ret.error_pixel_count; output.saturated_pixel_count = ret.saturated_pixel_count; + output.pixel_sum = ret.pixel_sum; output.az_int_profile = profile.GetResult(); output.az_int_profile_count = profile.GetPixelCount(); output.az_int_profile_std = profile.GetStd(); diff --git a/image_analysis/image_preprocessing/ImagePreprocessor.h b/image_analysis/image_preprocessing/ImagePreprocessor.h index c9a60b97b..0c552d309 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessor.h +++ b/image_analysis/image_preprocessing/ImagePreprocessor.h @@ -17,6 +17,9 @@ struct alignas(8) ImageStatistics { unsigned long long masked_pixel_count = 0; long long max_value = INT64_MIN; long long min_value = INT64_MAX; + // Of the pixels that are not error, saturated or under the device mask, as on the FPGA. A pixel masked + // only by what a run determined (PixelMask::RunMaskBits) counts, so the sum does not depend on it. + long long pixel_sum = 0; }; class ImagePreprocessor { diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp index 8a24acb8a..86c1953f1 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp @@ -9,7 +9,8 @@ ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &experiment, const PixelMask &mask) : ImagePreprocessor(experiment), - mask_bits(mask.GetPackedMask()) { + mask_bits(mask.GetPackedMask()), + run_mask_bits(mask.GetPackedRunMask()) { if (mask.GetMask().size() != npixels) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Pixel mask size doesn't match the image"); @@ -60,19 +61,25 @@ void ImagePreprocessorCPU::AnalyzeBlock(ImagePreprocessorBuffer &processed_image using U = std::make_unsigned_t; constexpr U flip = std::is_signed_v ? static_cast(U{1} << (8 * sizeof(T) - 1)) : U{0}; uint64_t n_masked = 0, n_error = 0, n_saturated = 0; + int64_t sum = 0; U key_max = 0, key_min = std::numeric_limits::max(); - const auto pixel = [&](size_t i, bool masked) { + const auto pixel = [&](size_t i, bool masked, bool run_masked) { const T v = image[i]; // Error/invalid marker = the pixel type's extreme value (0xFFFFFFFF for EIGER uint32). Tested // before saturation, since for unsigned types the marker also exceeds sat_pixel_val (which is // clipped above to the HDF5 saturation_value). - const bool error = !masked & (v == err_pixel_val); - const bool saturated = !masked & !error & (v >= sat_pixel_val); - const U counted = static_cast(-static_cast(!masked & !error & !saturated)); // all bits, or none + const bool is_error = v == err_pixel_val; + const bool is_saturated = !is_error & (v >= sat_pixel_val); + const bool error = !masked & is_error; + const bool saturated = !masked & is_saturated; + const U counted = static_cast(-static_cast(!masked & !is_error & !is_saturated)); // all bits, or none + const U summed = static_cast(-static_cast((!masked | run_masked) & !is_error & !is_saturated)); out[i] = (masked | error) ? INT32_MIN : saturated ? INT32_MAX : static_cast(v); n_masked += masked; n_error += error; n_saturated += saturated; + // Masked in the pixel's own width and only then widened, which is what keeps the loop vectorised. + sum += static_cast(static_cast(v) & summed); const U key = static_cast(v) ^ flip; key_max = std::max(key_max, key & counted); key_min = std::min(key_min, key | static_cast(~counted)); @@ -86,15 +93,17 @@ void ImagePreprocessorCPU::AnalyzeBlock(ImagePreprocessorBuffer &processed_image end += 32; if (end > i) { for (; i < end; i++) - pixel(i, false); + pixel(i, false, false); } else { - pixel(i, (mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U); + const bool masked = (mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U; + pixel(i, masked, masked && ((run_mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U)); i++; } } ret.masked_pixel_count += n_masked; ret.error_pixel_count += n_error; ret.saturated_pixel_count += n_saturated; + ret.pixel_sum += sum; if (n_masked + n_error + n_saturated < n) { // some pixel counted ret.max_value = std::max(ret.max_value, static_cast(key_max ^ flip)); ret.min_value = std::min(ret.min_value, static_cast(key_min ^ flip)); diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.h b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h index c0e1b80d6..dd7aa68d6 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h @@ -10,6 +10,7 @@ class ImagePreprocessorCPU : public ImagePreprocessor { // The pixel mask one bit per pixel, 32 to a word - PixelMask derives it once, so each engine only // copies it. std::vector mask_bits; + std::vector run_mask_bits; template void AnalyzeBlock(ImagePreprocessorBuffer &processed_image, size_t first, const uint8_t *input, size_t n, T err_value, T sat_value, ImageStatistics &stats); diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index c7ccb3793..d21de562d 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -38,6 +38,7 @@ __global__ void preprocess_kernel( __shared__ unsigned long long s_error; __shared__ long long s_max; __shared__ long long s_min; + __shared__ long long s_sum; if (threadIdx.x == 0) { s_masked = 0; @@ -45,6 +46,7 @@ __global__ void preprocess_kernel( s_error = 0; s_max = INT64_MIN; s_min = INT64_MAX; + s_sum = 0; } __syncthreads(); @@ -54,13 +56,15 @@ __global__ void preprocess_kernel( unsigned long long local_error = 0; long long local_max = INT64_MIN; long long local_min = INT64_MAX; + long long local_sum = 0; for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < npixels; i += blockDim.x * gridDim.x) { T v = input[i]; - bool is_masked = mask[i]; + const uint8_t mask_value = mask[i]; + bool is_masked = mask_value != 0; // Error/invalid marker = the pixel type's extreme value (0xFFFFFFFF for EIGER uint32); tested // before saturation, since for unsigned types the marker also exceeds saturation_limit (which is // clipped to the HDF5 saturation_value). Priority: masked > error > saturated. @@ -84,12 +88,15 @@ __global__ void preprocess_kernel( if (val > local_max) local_max = val; if (val < local_min) local_min = val; } + if (mask_value != PixelMask::BinaryMaskDevice && !is_err && !is_sat) + local_sum += (int64_t) v; } - // Reduce to shared memory + // Reduce to shared memory. There is no signed 64-bit atomicAdd; two's complement adds the same bits. atomicAdd(&s_masked, local_masked); atomicAdd(&s_saturated, local_saturated); atomicAdd(&s_error, local_error); + atomicAdd((unsigned long long *) &s_sum, (unsigned long long) local_sum); if (local_min <= local_max) { atomicMax((long long *) &s_max, (long long) local_max); @@ -103,6 +110,7 @@ __global__ void preprocess_kernel( atomicAdd(&stats->masked_pixel_count, s_masked); atomicAdd(&stats->saturated_pixel_count, s_saturated); atomicAdd(&stats->error_pixel_count, s_error); + atomicAdd((unsigned long long *) &stats->pixel_sum, (unsigned long long) s_sum); atomicMax((long long *) &stats->max_value, (long long) s_max); atomicMin((long long *) &stats->min_value, (long long) s_min); @@ -114,9 +122,10 @@ __global__ void preprocess_kernel( template struct PreprocessAccum { unsigned long long masked = 0, saturated = 0, error = 0; - long long max_v = INT64_MIN, min_v = INT64_MAX; + long long max_v = INT64_MIN, min_v = INT64_MAX, sum = 0; - __device__ __forceinline__ int32_t Apply(T v, bool is_masked, T sat_value, T err_value) { + __device__ __forceinline__ int32_t Apply(T v, uint8_t mask_value, T sat_value, T err_value) { + const bool is_masked = mask_value != 0; const bool is_err = (v == err_value); const bool is_sat = !is_err && (v >= sat_value); masked += is_masked; @@ -127,6 +136,8 @@ struct PreprocessAccum { if (val > max_v) max_v = val; if (val < min_v) min_v = val; } + if (mask_value != PixelMask::BinaryMaskDevice && !is_err && !is_sat) + sum += (int64_t) v; return is_masked ? INT32_MIN : is_err ? INT32_MIN : is_sat ? INT32_MAX : (int32_t) v; } }; @@ -147,12 +158,14 @@ __device__ __forceinline__ void FlushStats(PreprocessAccum &l, ImageStatistic l.error += __shfl_down_sync(0xffffffff, l.error, d); l.max_v = max(l.max_v, __shfl_down_sync(0xffffffff, l.max_v, d)); l.min_v = min(l.min_v, __shfl_down_sync(0xffffffff, l.min_v, d)); + l.sum += __shfl_down_sync(0xffffffff, l.sum, d); } if ((threadIdx.x & 31) == 0) { atomicAdd(&stats->masked_pixel_count, l.masked); atomicAdd(&stats->saturated_pixel_count, l.saturated); atomicAdd(&stats->error_pixel_count, l.error); + atomicAdd((unsigned long long *) &stats->pixel_sum, (unsigned long long) l.sum); atomicMax((long long *) &stats->max_value, l.max_v); atomicMin((long long *) &stats->min_value, l.min_v); } @@ -203,7 +216,7 @@ __device__ __forceinline__ void untranspose_preprocess_block(const uint8_t *in, 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); + o[k] = l.Apply((T) uv, mask[elem0 + i * 8 + k], 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); @@ -226,7 +239,7 @@ __device__ __forceinline__ void preprocess_tail(const uint8_t *__restrict__ tail 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); + out[tail_elem0 + threadIdx.x] = l.Apply(v, mask[tail_elem0 + threadIdx.x], sat_value, err_value); } } diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 4d3dd97fc..19c278dfe 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -631,8 +631,10 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen .value_or("")); } } - if (master_file->Exists("/entry/image")) + if (master_file->Exists("/entry/image")) { dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); + dataset->pixel_sum = master_file->ReadOptVector("/entry/image/pixel_sum"); + } }; if (master_file->Exists("/entry/data/data")) { @@ -809,6 +811,9 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen ReadVector(dataset->max_value, data_file, "/entry/image/max_value", number_of_images, fimages); + ReadVector(dataset->pixel_sum, + data_file, "/entry/image/pixel_sum", + number_of_images, fimages); } } catch (JFJochException &e) { // Say why. Everything read here is optional per-image metadata except the @@ -1476,6 +1481,8 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im message.indexing_lattice_count = dataset->indexing_lattice_count[image_number]; if (dataset->bkg_estimate.size() > image_number) message.bkg_estimate = dataset->bkg_estimate[image_number]; + if (dataset->pixel_sum.size() > image_number) + message.pixel_sum = dataset->pixel_sum[image_number]; if (dataset->xfel_pulse_id.size() > image_number) message.xfel_pulse_id = dataset->xfel_pulse_id[image_number]; if (dataset->xfel_event_code.size() > image_number) diff --git a/reader/JFJochMultiReader.cpp b/reader/JFJochMultiReader.cpp index ba91a4f6a..3a65e93df 100644 --- a/reader/JFJochMultiReader.cpp +++ b/reader/JFJochMultiReader.cpp @@ -99,6 +99,7 @@ JFJochMultiReader::JFJochMultiReader(std::vector> dataset->indexing_result.clear(); dataset->indexing_lattice_count.clear(); dataset->bkg_estimate.clear(); + dataset->pixel_sum.clear(); dataset->ice_ring_score.clear(); dataset->resolution_estimate.clear(); dataset->image_scale_factor.clear(); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 5c1cc154e..7dd906bea 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -79,6 +79,7 @@ struct JFJochReaderDataset { std::vector image_scale_factor; std::vector image_scale_cc; std::vector max_value; + std::vector pixel_sum; // Per-image sweep-quality code from /entry/MX/sweepQuality: 0 = the image is in no flagged range, // otherwise a 1-based index into sweep_quality_reasons (the vocabulary stored beside it, so the diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 559864c22..bd4271384 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -167,6 +167,7 @@ ADD_EXECUTABLE(jfjoch_test BSLZ4DecoderGPUTest.cpp BSLZ4DecoderGPUFuzzTest.cpp ImagePreprocessorGPUFusedTest.cpp + ImagePreprocessorPixelSumTest.cpp StoredImageDecodeGPUTest.cpp DirectHostReadGPUTest.cpp AdaptiveSpotFinderGPUTest.cpp diff --git a/tests/DirectHostReadGPUTest.cpp b/tests/DirectHostReadGPUTest.cpp index 2a150aac7..a07a9a98c 100644 --- a/tests/DirectHostReadGPUTest.cpp +++ b/tests/DirectHostReadGPUTest.cpp @@ -24,7 +24,7 @@ namespace { bool SameStats(const ImageStatistics &a, const ImageStatistics &b) { - return a.max_value == b.max_value && a.min_value == b.min_value + return a.max_value == b.max_value && a.min_value == b.min_value && a.pixel_sum == b.pixel_sum && a.masked_pixel_count == b.masked_pixel_count && a.error_pixel_count == b.error_pixel_count && a.saturated_pixel_count == b.saturated_pixel_count; diff --git a/tests/FPGAIntegrationTest.cpp b/tests/FPGAIntegrationTest.cpp index 2fd6358dd..09ad5b2d5 100644 --- a/tests/FPGAIntegrationTest.cpp +++ b/tests/FPGAIntegrationTest.cpp @@ -6,6 +6,11 @@ #include "../jungfrau/JFPedestalCalc.h" #include "../acquisition_device/HLSSimulatedDevice.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#endif #include "CheckImageOutput.h" using namespace std::literals::chrono_literals; @@ -2749,3 +2754,76 @@ TEST_CASE("HLS_C_Simulation_internal_packet_generator_4_images", "[FPGA][Full]") } } } + +// The claim behind ImageStatistics::pixel_sum: the preprocessors, reading the frame the card wrote, give the +// card's own sum over the same pixels. +namespace { + template + void CheckPixelSumParity(int64_t nsummation, CompressedImageMode mode) { + const uint16_t nmodules = 2; + + DiffractionExperiment x(DetJF(nmodules)); + x.Raw(); + x.UseInternalPacketGenerator(true).ImagesPerTrigger(1).PedestalG0Frames(0).Summation(nsummation) + .AutoSummation(true).ApplyPixelMask(true).BitDepthImage(8 * sizeof(T)).PixelSigned(true); + + std::mt19937 rng(5); + std::vector frame(RAW_MODULE_SIZE); + for (auto &v: frame) + v = static_cast(static_cast(rng() % 401) - 200); + frame[1024 + 1] = 2000; + frame[1024 + 2] = -1500; + frame[1024 + 3] = INT16_MAX; + frame[1024 + 4] = INT16_MIN; + frame[1024 + 5] = 1234; + + std::vector module_mask(RAW_MODULE_SIZE, 0); + module_mask[1024 + 5] = 1; + module_mask[1024 + 6] = 1; + + HLSSimulatedDevice test(0, 64); + std::vector user_mask; + for (int m = 0; m < x.GetModulesNum(); m++) { + test.InitializePixelMask(module_mask.data(), m); + test.SetInternalGeneratorFrame((uint16_t *) frame.data(), m); + user_mask.insert(user_mask.end(), module_mask.begin(), module_mask.end()); + } + REQUIRE_NOTHROW(test.StartAction(x)); + REQUIRE_NOTHROW(test.WaitForActionComplete()); + + std::vector image(nmodules * RAW_MODULE_SIZE); + int64_t card_sum = 0; + for (int m = 0; m < nmodules; m++) { + const auto *output = test.GetDeviceOutput(0, m); + memcpy(image.data() + m * RAW_MODULE_SIZE, output->pixels, RAW_MODULE_SIZE * sizeof(T)); + card_sum += output->module_statistics.pixel_sum; + } + CHECK(card_sum != 0); + + // The modules one under the other, as a detector without a card delivers its image. + DiffractionExperiment x_image(DetDECTRIS(RAW_MODULE_COLS, nmodules * RAW_MODULE_LINES, "Test", {}) + .BitDepthImage(8 * sizeof(T))); + x_image.PixelSigned(true); + PixelMask mask(x_image); + mask.LoadUserMask(x_image, user_mask); + const auto input = reinterpret_cast(image.data()); + + ImagePreprocessorCPU cpu(x_image, mask); + ImagePreprocessorBuffer cpu_buffer(x_image.GetPixelsNum()); + CHECK(cpu.Analyze(cpu_buffer, input, mode).pixel_sum == card_sum); + +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() > 0) { + ImagePreprocessorGPU gpu(x_image, mask, std::make_shared(), /*copy_image_to_host=*/false); + ImagePreprocessorBufferGPU gpu_buffer(x_image.GetPixelsNum()); + CHECK(gpu.Analyze(gpu_buffer, input, mode).pixel_sum == card_sum); + } +#endif + } +} + +TEST_CASE("HLS_C_Simulation_pixel_sum_matches_preprocessor", "[FPGA][Full]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + CheckPixelSumParity(1, CompressedImageMode::Int16); + CheckPixelSumParity(4, CompressedImageMode::Int32); +} diff --git a/tests/ImagePreprocessorGPUFusedTest.cpp b/tests/ImagePreprocessorGPUFusedTest.cpp index fcd122b83..37f937377 100644 --- a/tests/ImagePreprocessorGPUFusedTest.cpp +++ b/tests/ImagePreprocessorGPUFusedTest.cpp @@ -54,7 +54,7 @@ std::vector MakeImage(size_t npixels, T err_value, uint32_t seed) { } bool SameStats(const ImageStatistics &a, const ImageStatistics &b) { - return a.max_value == b.max_value && a.min_value == b.min_value + return a.max_value == b.max_value && a.min_value == b.min_value && a.pixel_sum == b.pixel_sum && a.masked_pixel_count == b.masked_pixel_count && a.error_pixel_count == b.error_pixel_count && a.saturated_pixel_count == b.saturated_pixel_count; diff --git a/tests/ImagePreprocessorPixelSumTest.cpp b/tests/ImagePreprocessorPixelSumTest.cpp new file mode 100644 index 000000000..c723a159c --- /dev/null +++ b/tests/ImagePreprocessorPixelSumTest.cpp @@ -0,0 +1,266 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include + +#include "../common/CUDAWrapper.h" +#include "../common/PixelMask.h" +#include "../common/ScanResultGenerator.h" +#include "../compression/JFJochCompressor.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" + +#ifdef JFJOCH_USE_CUDA +#include "../compression/bitshuffle/bitshuffle.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#endif + +// pixel_sum on the CPU and GPU preprocessors, to the card's definition (fpga/hls/pixel_calc.cpp): signed, +// and over the pixels that are not masked, not error and not saturated. +namespace { + constexpr int64_t SATURATION = 110; + + template + void CheckEveryPixelType(CompressedImageMode mode, T err_value) { + DiffractionExperiment x(DetDECTRIS(64, 4, "Test", {}).BitDepthImage(8 * sizeof(T)).SaturationLimit(SATURATION)); + x.PixelSigned(std::is_signed_v); + + // Pixels 0-63 hold a masked pixel, so they are tested one by one; the words after them are whole. + std::vector user_mask(x.GetPixelsNum(), 0); + user_mask[3] = 1; + user_mask[4] = 1; + user_mask[5] = 1; + PixelMask mask(x); + mask.LoadUserMask(x, user_mask); + + std::vector image(x.GetPixelsNum(), 0); + image[1] = 100; + image[3] = 50; // masked + image[4] = err_value; // masked before it is an error + image[5] = SATURATION; // masked before it is saturated + image[70] = 20; + image[72] = SATURATION - 1; + image[64] = err_value; + image[65] = SATURATION; + image[66] = SATURATION + 1; + + int64_t expected = 100 + 20 + SATURATION - 1; + if (std::is_signed_v) { + image[71] = static_cast(-30); + expected -= 30; + } + + ImagePreprocessorCPU preprocessor(x, mask); + ImagePreprocessorBuffer buffer(x.GetPixelsNum()); + const auto input = reinterpret_cast(image.data()); + const ImageStatistics stats = preprocessor.Analyze(buffer, input, mode); + CHECK(stats.pixel_sum == expected); + CHECK(stats.masked_pixel_count == 3); + CHECK(stats.error_pixel_count == 1); + CHECK(stats.saturated_pixel_count == 2); + + // In blocks, as a bitshuffled image is decoded. + ImageStatistics blocks{}; + preprocessor.AnalyzeBlock(buffer, 0, input, 40, mode, blocks); + preprocessor.AnalyzeBlock(buffer, 40, input + 40 * sizeof(T), x.GetPixelsNum() - 40, mode, blocks); + CHECK(blocks.pixel_sum == expected); + } +} + +TEST_CASE("ImagePreprocessorCPU_PixelSum", "[ImagePreprocessor]") { + CheckEveryPixelType(CompressedImageMode::Int8, INT8_MIN); + CheckEveryPixelType(CompressedImageMode::Uint8, UINT8_MAX); + CheckEveryPixelType(CompressedImageMode::Int16, INT16_MIN); + CheckEveryPixelType(CompressedImageMode::Uint16, UINT16_MAX); + CheckEveryPixelType(CompressedImageMode::Int32, INT32_MIN); + CheckEveryPixelType(CompressedImageMode::Uint32, UINT32_MAX); +} + +// Every pixel just under the limit of a large detector: far past 32 bits, in a sum of 32-bit pixels. +TEST_CASE("ImagePreprocessorCPU_PixelSum_NoOverflow", "[ImagePreprocessor]") { + DiffractionExperiment x(DetDECTRIS(4096, 4096, "Test", {}).BitDepthImage(32)); + x.PixelSigned(true); + const PixelMask mask(x); + + const int32_t value = static_cast(x.GetSaturationLimit()) - 1; + const std::vector image(x.GetPixelsNum(), value); + + ImagePreprocessorCPU preprocessor(x, mask); + ImagePreprocessorBuffer buffer(x.GetPixelsNum()); + const ImageStatistics stats = preprocessor.Analyze(buffer, reinterpret_cast(image.data()), + CompressedImageMode::Int32); + CHECK(stats.pixel_sum == static_cast(value) * static_cast(image.size())); + CHECK(stats.pixel_sum > INT32_MAX); +} + +// A mask a run determines (a beam-stop shadow, hot pixels) takes pixels out of the image, not out of the sum. +TEST_CASE("ImagePreprocessorCPU_PixelSum_RunMask", "[ImagePreprocessor]") { + DiffractionExperiment x(DetDECTRIS(64, 4, "Test", {}).BitDepthImage(32)); + x.PixelSigned(true); + + std::vector user_mask(x.GetPixelsNum(), 0), shadow(x.GetPixelsNum(), 0), hot(x.GetPixelsNum(), 0); + user_mask[1] = 1; // the device's + shadow[1] = 1; // and the run's: the device mask decides + shadow[2] = 1; + shadow[3] = 1; // error pixel under the shadow + hot[4] = 1; + hot[5] = 1; // saturated pixel, hot + + PixelMask device(x); + device.LoadUserMask(x, user_mask); + PixelMask with_run = device; + with_run.LoadBeamStopMask(x, shadow); + with_run.LoadHotPixelMask(x, hot); + + std::vector image(x.GetPixelsNum(), 0); + image[0] = 11; + image[1] = 1000; + image[2] = 22; + image[3] = INT32_MIN; + image[4] = -33; + image[5] = INT32_MAX; + const auto input = reinterpret_cast(image.data()); + + ImagePreprocessorBuffer buffer(x.GetPixelsNum()); + ImagePreprocessorCPU device_only(x, device); + CHECK(device_only.Analyze(buffer, input, CompressedImageMode::Int32).pixel_sum == 11 + 22 - 33); + + ImagePreprocessorCPU run(x, with_run); + const ImageStatistics stats = run.Analyze(buffer, input, CompressedImageMode::Int32); + CHECK(stats.pixel_sum == 11 + 22 - 33); + CHECK(stats.masked_pixel_count == 5); + CHECK(stats.max_value == 11); + CHECK(buffer[2] == INT32_MIN); +} + +// CPU or GPU, whichever this machine runs, from the image to the vector the master file is written from. +TEST_CASE("MXAnalysisWithoutFPGA_PixelSum", "[ImagePreprocessor]") { + constexpr size_t width = 1030, height = 1064; + + DiffractionExperiment experiment(DetDECTRIS(width, height, "Test", {}).BitDepthImage(32)); + experiment.ImagesPerTrigger(1).NumTriggers(1).PixelSigned(true) + .DetectorDistance_mm(75).BeamX_pxl(515).BeamY_pxl(532).IncidentEnergy_keV(12.4); + + std::vector user_mask(width * height, 0), hot(width * height, 0); + user_mask[1000] = 1; + hot[10] = 1; + PixelMask pixel_mask(experiment); + pixel_mask.LoadUserMask(experiment, user_mask); + pixel_mask.LoadHotPixelMask(experiment, hot); + + AzimuthalIntegrationMapping mapping(experiment, pixel_mask); + IndexAndRefine index_and_refine(experiment, nullptr); + MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, index_and_refine); + + std::vector pixels(width * height, 0); + pixels[10] = 3'000'000; + pixels[11] = 2'000'000'000; + pixels[12] = -40; + pixels[1000] = 777; + pixels[2000] = INT32_MIN; + pixels[3000] = INT32_MAX; + const int64_t expected = int64_t{3'000'000} + 2'000'000'000 - 40; + + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + const std::vector compressed = compressor.Compress(pixels); + + DataMessage message{}; + message.number = 0; + message.image = CompressedImage(compressed.data(), compressed.size(), width, height, + CompressedImageMode::Int32, CompressionAlgorithm::BSHUF_LZ4); + + SpotFindingSettings settings = DiffractionExperiment::DefaultDataProcessingSettings(); + settings.enable = false; + + AzimuthalIntegrationProfile profile(mapping); + analysis.Analyze(message, profile, settings); + REQUIRE(message.pixel_sum); + CHECK(message.pixel_sum.value() == expected); + + ScanResultGenerator scan_result(experiment); + scan_result.Add(message); + EndMessage end{}; + scan_result.FillEndMessage(end); + CHECK(end.pixel_sum == std::vector{expected}); +} + +#ifdef JFJOCH_USE_CUDA +// All three routes through the kernels against the CPU, with a run's shadow over the device mask: the upload, the decode fused with the preprocessing, +// and the two-kernel route of a bitshuffle block too large for the fused one. +TEST_CASE("ImagePreprocessorGPU_PixelSum", "[ImagePreprocessorGPU]") { + if (get_gpu_count() == 0) + SKIP("No CUDA GPU present"); + + DiffractionExperiment x(DetDECTRIS(4096, 4096, "Test", {}).BitDepthImage(32)); + x.PixelSigned(true); + const size_t npixels = x.GetPixelsNum(); + + std::vector user_mask(npixels, 0); + std::vector shadow(npixels, 0); + for (size_t i = 0; i < npixels; i += 997) + user_mask[i] = 1; + for (size_t i = 0; i < npixels; i += 501) + shadow[i] = 1; + PixelMask mask(x); + mask.LoadUserMask(x, user_mask); + mask.LoadBeamStopMask(x, shadow); + + // Mostly near the limit, so the sum needs far more than 32 bits; negative pixels and markers among them. + std::mt19937 rng(11); + const int32_t limit = static_cast(x.GetSaturationLimit()); + std::vector image(npixels); + for (auto &v: image) { + const uint32_t r = rng() % 1000; + v = r == 0 ? INT32_MIN : r == 1 ? limit : r < 100 ? -static_cast(rng() % 100000) : limit - 1 - static_cast(r); + } + + int64_t expected = 0; + for (size_t i = 0; i < npixels; i++) + if (user_mask[i] == 0 && image[i] != INT32_MIN && image[i] < limit) + expected += image[i]; + REQUIRE(expected > INT32_MAX); + + const auto input = reinterpret_cast(image.data()); + ImagePreprocessorCPU cpu(x, mask); + ImagePreprocessorBuffer cpu_buffer(npixels); + CHECK(cpu.Analyze(cpu_buffer, input, CompressedImageMode::Int32).pixel_sum == expected); + + auto stream = std::make_shared(); + ImagePreprocessorGPU gpu(x, mask, stream, /*copy_image_to_host=*/false); + ImagePreprocessorBufferGPU gpu_buffer(npixels); + + CHECK(gpu.Analyze(gpu_buffer, input, CompressedImageMode::Int32).pixel_sum == expected); + + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + const std::vector compressed = compressor.Compress(image); + ImageStatistics fused{}; + REQUIRE(gpu.AnalyzeCompressed(gpu_buffer, CompressedImage(compressed.data(), compressed.size(), x.GetXPixelsNum(), + x.GetYPixelsNum(), CompressedImageMode::Int32, + CompressionAlgorithm::BSHUF_LZ4), fused)); + CHECK(fused.pixel_sum == expected); + + // The container by hand, with 32 kB blocks: 64-bit size and 32-bit block size in bytes, big-endian. + constexpr size_t block_elements = 8192; + std::vector large_blocks(12 + bshuf_compress_lz4_bound(npixels, sizeof(int32_t), block_elements)); + const uint64_t total_bytes = npixels * sizeof(int32_t); + for (int i = 0; i < 8; i++) + large_blocks[i] = static_cast(total_bytes >> (56 - 8 * i)); + const uint32_t block_bytes = block_elements * sizeof(int32_t); + for (int i = 0; i < 4; i++) + large_blocks[8 + i] = static_cast(block_bytes >> (24 - 8 * i)); + const int64_t written = bshuf_compress_lz4(image.data(), large_blocks.data() + 12, npixels, sizeof(int32_t), + block_elements); + REQUIRE(written > 0); + ImageStatistics two_kernels{}; + REQUIRE(gpu.AnalyzeCompressed(gpu_buffer, CompressedImage(large_blocks.data(), 12 + written, x.GetXPixelsNum(), + x.GetYPixelsNum(), CompressedImageMode::Int32, + CompressionAlgorithm::BSHUF_LZ4), two_kernels)); + CHECK(two_kernels.pixel_sum == expected); +} +#endif diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index b26b3d350..79eed81b1 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -737,6 +737,7 @@ TEST_CASE("JFJochReader_DataI16", "[HDF5][Full][portable]") { message.spots = spots; message.indexing_result = (i % 2 == 0); message.bkg_estimate = i * 345.6; + message.pixel_sum = -5'000'000'000LL + i; // Only one frame has a spindle severity: the others must come back ABSENT (stored as // NaN), because no value is the CANNOT-SAY trigger state and a zero is not. if (i == 2) @@ -762,7 +763,9 @@ TEST_CASE("JFJochReader_DataI16", "[HDF5][Full][portable]") { REQUIRE(dataset->bkg_estimate.size() == 4); REQUIRE(dataset->profile_radius.size() == 4); REQUIRE(dataset->spindle_blind_fraction.size() == 4); + REQUIRE(dataset->pixel_sum.size() == 4); for (int i = 0; i < 4; i++) { + CHECK(dataset->pixel_sum[i] == -5'000'000'000LL + i); if (i == 2) CHECK(dataset->spindle_blind_fraction[i] == Catch::Approx(0.75f)); else @@ -780,6 +783,7 @@ TEST_CASE("JFJochReader_DataI16", "[HDF5][Full][portable]") { CHECK(reader_image->Image()[2] == image[2]); CHECK(reader_image->Image()[3] == image[3]); CHECK(reader_image->Image()[5678] == i); + CHECK(reader_image->ImageData().pixel_sum == -5'000'000'000LL + i); CHECK(dataset->indexing_result[i] == (i % 2 == 0)); CHECK(dataset->bkg_estimate[i] == Catch::Approx(i * 345.6)); CHECK(dataset->profile_radius[i] == Catch::Approx(123.09)); diff --git a/tests/JFJochReceiverLiteTest.cpp b/tests/JFJochReceiverLiteTest.cpp index bcfe1167d..f9a009729 100644 --- a/tests/JFJochReceiverLiteTest.cpp +++ b/tests/JFJochReceiverLiteTest.cpp @@ -115,6 +115,13 @@ TEST_CASE("JFJochReceiverLite", "[JFJochReceiver]") { CHECK(!master.GetBool("/entry/instrument/detector/pixel_mask_applied")); CHECK(master.ReadVector("/entry/instrument/detector/countrate_correction_lookup_table") == std::vector{0, 1, 3, 7}); + + int64_t expected_sum = 0; + for (size_t i = 0; i < image_conv.size(); i++) + if (pixel_mask.GetMask()[i] == 0 && image_conv[i] != INT16_MIN && image_conv[i] < experiment.GetSaturationLimit()) + expected_sum += image_conv[i]; + CHECK(expected_sum != 0); + CHECK(master.ReadVector("/entry/image/pixel_sum") == std::vector(5, expected_sum)); } TEST_CASE("JFJochReceiverLite_Cancel", "[JFJochReceiver]") { @@ -237,4 +244,56 @@ TEST_CASE("JFJochReceiverLite_DarkMask", "[JFJochReceiver]") { // No progress value at the end of the measurement REQUIRE(!service.GetProgress().has_value()); -} \ No newline at end of file +} +// Identical frames, so the run masks a large part of the detector as a shadow and as hot pixels. That mask +// is the run's own: the master file holds the sum under the device mask. +TEST_CASE("JFJochReceiverRugnux_PixelSumUnderDeviceMask", "[JFJochReceiver]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + DiffractionExperiment experiment(DetDECTRIS(2068, 2164, "Test", {})); + experiment.ImagesPerTrigger(20).NumTriggers(1).UseInternalPacketGenerator(true).ImagesPerFile(10) + .FilePrefix("rugnux_live_pixel_sum").JungfrauConvPhotonCnt(false).PixelSigned(true) + .SetFileWriterFormat(FileWriterFormat::NXmxVDS).OverwriteExistingFiles(true) + .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4) + .Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {})); + PixelMask pixel_mask(experiment); + + HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); + HDF5DataSet dataset(data, "/entry/data/data"); + HDF5DataSpace file_space(dataset); + std::vector image_conv(file_space.GetDimensions()[1] * file_space.GetDimensions()[2]); + std::vector start = {4, 0, 0}; + std::vector file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]}; + dataset.ReadVector(image_conv, start, file_size); + + HDF5FilePusher pusher; + auto puller = std::make_shared(); + StartMessage start_msg; + experiment.FillMessage(start_msg); + puller->Put(ImagePullerOutput{.cbor = std::make_shared(start_msg)}); + DataMessage data_msg; + data_msg.image = CompressedImage(image_conv, file_space.GetDimensions()[2], file_space.GetDimensions()[1]); + for (int i = 0; i < experiment.GetImageNum(); i++) { + data_msg.number = i; + puller->Put(ImagePullerOutput{.cbor = std::make_shared(data_msg)}); + } + puller->Put(ImagePullerOutput{.cbor = std::make_shared(EndMessage{})}); + + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.RugnuxLive(true, nullptr, 4, 60); + service.SetSpotFindingSettings(DiffractionExperiment::DefaultDataProcessingSettings()); + service.Indexing(experiment.GetIndexingSettings()); + service.Start(experiment, pixel_mask, nullptr, puller); + service.Stop(); + + int64_t expected_sum = 0; + for (size_t i = 0; i < image_conv.size(); i++) + if (pixel_mask.GetMask()[i] == 0 && image_conv[i] != INT16_MIN && image_conv[i] < experiment.GetSaturationLimit()) + expected_sum += image_conv[i]; + CHECK(expected_sum != 0); + + HDF5ReadOnlyFile master("rugnux_live_pixel_sum_master.h5"); + CHECK(master.ReadVector("/entry/image/pixel_sum") == std::vector(20, expected_sum)); +} diff --git a/tests/StoredImageDecodeGPUTest.cpp b/tests/StoredImageDecodeGPUTest.cpp index 8663a8079..5db2d3052 100644 --- a/tests/StoredImageDecodeGPUTest.cpp +++ b/tests/StoredImageDecodeGPUTest.cpp @@ -83,7 +83,7 @@ std::vector DecodeOnDevice(const std::vector &stream_bytes, si } bool SameStats(const ImageStatistics &a, const ImageStatistics &b) { - return a.max_value == b.max_value && a.min_value == b.min_value + return a.max_value == b.max_value && a.min_value == b.min_value && a.pixel_sum == b.pixel_sum && a.masked_pixel_count == b.masked_pixel_count && a.error_pixel_count == b.error_pixel_count && a.saturated_pixel_count == b.saturated_pixel_count;