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 <noreply@anthropic.com>

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<uint32_t>(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<uint32_t>(binary_mask[i] != 0) << (i % 32);
+        packed_run_mask[i / 32] |= static_cast<uint32_t>(binary_mask[i] == BinaryMaskRunOnly) << (i % 32);
+    }
 }
 
 void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) {
@@ -134,6 +137,10 @@ const std::vector<uint32_t> &PixelMask::GetPackedMask() const {
     return packed_mask;
 }
 
+const std::vector<uint32_t> &PixelMask::GetPackedRunMask() const {
+    return packed_run_mask;
+}
+
 uint64_t PixelMask::GetBinaryMaskChecksum() const {
     return binary_mask_checksum;
 }

common/PixelMask.h

@@ -20,8 +20,8 @@ struct PixelMaskStatistics {
 class PixelMask {
     std::vector<uint32_t> mask;
     std::vector<uint32_t> 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<uint8_t> 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<uint32_t> packed_mask;
+    // Packed the same way: the pixels masked only by what a run found on its own frames.
+    std::vector<uint32_t> packed_run_mask;
 
     uint32_t LoadMask(const std::vector<uint32_t>& 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<uint32_t> &GetMask() const;
     [[nodiscard]] const std::vector<uint8_t> &GetBinaryMask() const;
     [[nodiscard]] const std::vector<uint32_t> &GetPackedMask() const;
+    [[nodiscard]] const std::vector<uint32_t> &GetPackedRunMask() const;
     [[nodiscard]] uint64_t GetBinaryMaskChecksum() const;
 
     [[nodiscard]] std::vector<uint32_t> GetUserMask(const DiffractionExperiment& experiment) const;

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 `<file_prefix>.mtz`. An MTZ table (layout below) | |

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).

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

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
 

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();

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 {

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<T>;
     constexpr U flip = std::is_signed_v<T> ? static_cast<U>(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<U>::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<U>(-static_cast<int>(!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<U>(-static_cast<int>(!masked & !is_error & !is_saturated));   // all bits, or none
+        const U summed = static_cast<U>(-static_cast<int>((!masked | run_masked) & !is_error & !is_saturated));
         out[i] = (masked | error) ? INT32_MIN : saturated ? INT32_MAX : static_cast<int32_t>(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<T>(static_cast<U>(v) & summed);
         const U key = static_cast<U>(v) ^ flip;
         key_max = std::max<U>(key_max, key & counted);
         key_min = std::min<U>(key_min, key | static_cast<U>(~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<long long>(ret.max_value, static_cast<T>(key_max ^ flip));
         ret.min_value = std::min<long long>(ret.min_value, static_cast<T>(key_min ^ flip));

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<uint32_t> mask_bits;
+    std::vector<uint32_t> run_mask_bits;
 
     template <class T> void AnalyzeBlock(ImagePreprocessorBuffer &processed_image, size_t first, const uint8_t *input,
                                          size_t n, T err_value, T sat_value, ImageStatistics &stats);

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<class T>
 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<T> &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<int4 *>(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);
     }
 }
 

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<int64_t>("/entry/image/max_value");
+            dataset->pixel_sum = master_file->ReadOptVector<int64_t>("/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)

reader/JFJochMultiReader.cpp

@@ -99,6 +99,7 @@ JFJochMultiReader::JFJochMultiReader(std::vector<std::unique_ptr<JFJochReader>>
     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();

reader/JFJochReaderDataset.h

@@ -79,6 +79,7 @@ struct JFJochReaderDataset {
     std::vector<float> image_scale_factor;
     std::vector<float> image_scale_cc;
     std::vector<int64_t> max_value;
+    std::vector<int64_t> 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

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

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;

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<class T>
+    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<int16_t> frame(RAW_MODULE_SIZE);
+        for (auto &v: frame)
+            v = static_cast<int16_t>(static_cast<int>(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<uint32_t> module_mask(RAW_MODULE_SIZE, 0);
+        module_mask[1024 + 5] = 1;
+        module_mask[1024 + 6] = 1;
+
+        HLSSimulatedDevice test(0, 64);
+        std::vector<uint32_t> 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<T> 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<const uint8_t *>(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<CudaStream>(), /*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<int16_t>(1, CompressedImageMode::Int16);
+    CheckPixelSumParity<int32_t>(4, CompressedImageMode::Int32);
+}

tests/ImagePreprocessorGPUFusedTest.cpp

@@ -54,7 +54,7 @@ std::vector<T> 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;

tests/ImagePreprocessorPixelSumTest.cpp

@@ -0,0 +1,266 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <limits>
+#include <random>
+#include <vector>
+
+#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<class T>
+    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<T>);
+
+        // Pixels 0-63 hold a masked pixel, so they are tested one by one; the words after them are whole.
+        std::vector<uint32_t> 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<T> 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<T>) {
+            image[71] = static_cast<T>(-30);
+            expected -= 30;
+        }
+
+        ImagePreprocessorCPU preprocessor(x, mask);
+        ImagePreprocessorBuffer buffer(x.GetPixelsNum());
+        const auto input = reinterpret_cast<const uint8_t *>(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<int8_t>(CompressedImageMode::Int8, INT8_MIN);
+    CheckEveryPixelType<uint8_t>(CompressedImageMode::Uint8, UINT8_MAX);
+    CheckEveryPixelType<int16_t>(CompressedImageMode::Int16, INT16_MIN);
+    CheckEveryPixelType<uint16_t>(CompressedImageMode::Uint16, UINT16_MAX);
+    CheckEveryPixelType<int32_t>(CompressedImageMode::Int32, INT32_MIN);
+    CheckEveryPixelType<uint32_t>(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<int32_t>(x.GetSaturationLimit()) - 1;
+    const std::vector<int32_t> image(x.GetPixelsNum(), value);
+
+    ImagePreprocessorCPU preprocessor(x, mask);
+    ImagePreprocessorBuffer buffer(x.GetPixelsNum());
+    const ImageStatistics stats = preprocessor.Analyze(buffer, reinterpret_cast<const uint8_t *>(image.data()),
+                                                       CompressedImageMode::Int32);
+    CHECK(stats.pixel_sum == static_cast<int64_t>(value) * static_cast<int64_t>(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<uint32_t> 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<int32_t> 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<const uint8_t *>(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<uint32_t> 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<int32_t> 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<uint8_t> 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<int64_t>{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<uint32_t> user_mask(npixels, 0);
+    std::vector<uint32_t> 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<int32_t>(x.GetSaturationLimit());
+    std::vector<int32_t> image(npixels);
+    for (auto &v: image) {
+        const uint32_t r = rng() % 1000;
+        v = r == 0 ? INT32_MIN : r == 1 ? limit : r < 100 ? -static_cast<int32_t>(rng() % 100000) : limit - 1 - static_cast<int32_t>(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<const uint8_t *>(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<CudaStream>();
+    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<uint8_t> 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<uint8_t> 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<uint8_t>(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<uint8_t>(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

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));

tests/JFJochReceiverLiteTest.cpp

@@ -115,6 +115,13 @@ TEST_CASE("JFJochReceiverLite", "[JFJochReceiver]") {
     CHECK(!master.GetBool("/entry/instrument/detector/pixel_mask_applied"));
     CHECK(master.ReadVector<uint32_t>("/entry/instrument/detector/countrate_correction_lookup_table")
           == std::vector<uint32_t>{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<int64_t>("/entry/image/pixel_sum") == std::vector<int64_t>(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());
-}
+}
+// 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<int16_t> image_conv(file_space.GetDimensions()[1] * file_space.GetDimensions()[2]);
+    std::vector<hsize_t> start = {4, 0, 0};
+    std::vector<hsize_t> 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<TestImagePuller>();
+    StartMessage start_msg;
+    experiment.FillMessage(start_msg);
+    puller->Put(ImagePullerOutput{.cbor = std::make_shared<CBORStream2DeserializerOutput>(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<CBORStream2DeserializerOutput>(data_msg)});
+    }
+    puller->Put(ImagePullerOutput{.cbor = std::make_shared<CBORStream2DeserializerOutput>(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<int64_t>("/entry/image/pixel_sum") == std::vector<int64_t>(20, expected_sum));
+}

tests/StoredImageDecodeGPUTest.cpp

@@ -83,7 +83,7 @@ std::vector<int32_t> DecodeOnDevice(const std::vector<uint8_t> &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;