image_analysis: compute ROI statistics in the non-FPGA path

MXAnalysisWithoutFPGA never filled DataMessage.roi, so ROI integrals were
only available on the FPGA path. Add a software ROI engine that mirrors the
FPGA roi_calc kernel: per-ROI sum, sum of squares, good-pixel count, max and
intensity-weighted centre of mass, with each pixel carrying a 16-bit mask so
it can contribute to any subset of up to 16 ROIs.

New image_analysis/roi/ library (JFJochROIIntegration), structured like azint:
a base that precomputes the per-pixel mask and names, a templated CPU engine
(generic over pixel type for a future 16-bit path), and a GPU kernel using
per-block shared-memory atomics for the STXM case (half-detector ROIs).

Masked pixels are excluded entirely; saturated pixels are excluded from the
sums but still count towards the max, matching roi_calc exactly. The engine is
only constructed when at least one ROI is defined. Downstream CBOR/HDF5 already
consume message.roi, so no further changes are needed.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

image_analysis/CMakeLists.txt

@@ -50,6 +50,7 @@ ADD_SUBDIRECTORY(lattice_search)
 ADD_SUBDIRECTORY(scale_merge)
 ADD_SUBDIRECTORY(image_preprocessing)
 ADD_SUBDIRECTORY(azint)
+ADD_SUBDIRECTORY(roi)
 ADD_SUBDIRECTORY(pixel_refinement)
 
-TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochPixelRefine gemmi)
+TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochROIIntegration JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochPixelRefine gemmi)

image_analysis/MXAnalysisWithoutFPGA.cpp

@@ -11,9 +11,11 @@
 #include "image_preprocessing/ImagePreprocessorCPU.h"
 
 #include "azint/AzIntEngineCPU.h"
+#include "roi/ROIIntegrationCPU.h"
 #include "spot_finding/ImageSpotFinderCPU.h"
 #ifdef JFJOCH_USE_CUDA
 #include "azint/AzIntEngineGPU.h"
+#include "roi/ROIIntegrationGPU.h"
 #include "spot_finding/ImageSpotFinderGPU.h"
 #include "image_preprocessing/ImagePreprocessorGPU.h"
 #include "image_preprocessing/ImagePreprocessorBufferGPU.h"
@@ -42,6 +44,8 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp
         spotFinder = std::make_unique<ImageSpotFinderCPU>(experiment.GetXPixelsNum(), experiment.GetYPixelsNum());
         azint = std::make_unique<AzIntEngineCPU>(integration);
         preprocessor = std::make_unique<ImagePreprocessorCPU>(in_experiment, in_mask);
+        if (experiment.ROI().size() >= 1)
+            roi = std::make_unique<ROIIntegrationCPU>(experiment);
 #ifdef JFJOCH_USE_CUDA
     } else {
         auto stream = std::make_shared<CudaStream>();
@@ -49,6 +53,8 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp
         preprocessor = std::make_unique<ImagePreprocessorGPU>(in_experiment, in_mask, stream);
         spotFinder = std::make_unique<ImageSpotFinderGPU>(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream);
         azint = std::make_unique<AzIntEngineGPU>(integration, stream);
+        if (experiment.ROI().size() >= 1)
+            roi = std::make_unique<ROIIntegrationGPU>(experiment, stream);
     }
 #endif
 }
@@ -77,6 +83,9 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
     const auto azint_end_time = std::chrono::steady_clock::now();
     output.azint_time_s = std::chrono::duration<float>(azint_end_time - azint_start_time).count();
 
+    if (roi)
+        roi->Run(*preprocessor_buffer, output.roi);
+
     if (spot_finding_settings.enable) {
         // Update resolution mask
         if (mask_high_res != spot_finding_settings.high_resolution_limit

image_analysis/MXAnalysisWithoutFPGA.h

@@ -14,6 +14,7 @@
 #include "spot_finding/ImageSpotFinder.h"
 #include "indexing/IndexerThreadPool.h"
 #include "azint/AzIntEngine.h"
+#include "roi/ROIIntegration.h"
 #include "IndexAndRefine.h"
 #include "image_preprocessing/ImagePreprocessor.h"
 #include "image_preprocessing/ImagePreprocessorBuffer.h"
@@ -31,6 +32,7 @@ class MXAnalysisWithoutFPGA {
     size_t xpixels;
 
     std::unique_ptr<AzIntEngine> azint;
+    std::unique_ptr<ROIIntegration> roi;
     std::unique_ptr<ImageSpotFinder> spotFinder;
     IndexAndRefine &indexer;
     std::unique_ptr<BraggPrediction> prediction;

image_analysis/roi/CMakeLists.txt

@@ -0,0 +1,5 @@
+ADD_LIBRARY(JFJochROIIntegration STATIC ROIIntegration.cpp ROIIntegration.h ROIIntegrationCPU.cpp ROIIntegrationCPU.h)
+TARGET_LINK_LIBRARIES(JFJochROIIntegration JFJochCommon)
+IF (JFJOCH_CUDA_AVAILABLE)
+    TARGET_SOURCES(JFJochROIIntegration PRIVATE ../indexing/CUDAMemHelpers.h ROIIntegrationGPU.cu ROIIntegrationGPU.h)
+ENDIF()

image_analysis/roi/ROIIntegration.cpp

@@ -0,0 +1,34 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ROIIntegration.h"
+#include "../../common/DiffractionExperiment.h"
+
+ROIIntegration::ROIIntegration(const DiffractionExperiment &experiment)
+    : roi_map(experiment.ExportROIMap()),
+      width(experiment.GetXPixelsNumConv()),
+      npixel(roi_map.size()),
+      roi_count(experiment.ROI().size()),
+      roi_name(roi_count),
+      roi_sum(roi_count),
+      roi_sum2(roi_count),
+      roi_pixels(roi_count),
+      roi_x_weighted(roi_count),
+      roi_y_weighted(roi_count),
+      roi_max(roi_count) {
+    for (const auto &[name, id] : experiment.ROI().GetROINameMap())
+        roi_name[id] = name;
+}
+
+void ROIIntegration::Export(std::map<std::string, ROIMessage> &out) const {
+    for (uint16_t r = 0; r < roi_count; r++)
+        out[roi_name[r]] = ROIMessage{
+            .sum           = roi_sum[r],
+            .sum_square    = roi_sum2[r],
+            .max_count     = roi_max[r],
+            .pixels        = roi_pixels[r],
+            .pixels_masked = 0,
+            .x_weighted    = roi_x_weighted[r],
+            .y_weighted    = roi_y_weighted[r],
+        };
+}

image_analysis/roi/ROIIntegration.h

@@ -0,0 +1,46 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstdint>
+#include <map>
+#include <string>
+#include <vector>
+
+#include "../../common/JFJochMessages.h"
+#include "../image_preprocessing/ImagePreprocessorBuffer.h"
+
+class DiffractionExperiment;
+
+// Computes per-ROI statistics over an assembled image: sum, sum of squares,
+// good-pixel count, maximum and intensity-weighted centre of mass. Each pixel
+// carries a 16-bit mask selecting which of up to 16 ROIs it belongs to, so a
+// single pixel can contribute to none, one, or several ROIs at once. This is
+// the software counterpart of the FPGA roi_calc kernel.
+class ROIIntegration {
+protected:
+    std::vector<uint16_t> roi_map;   // per-pixel ROI bitmask, assembled (converted) geometry
+    size_t width;                    // assembled image width, to recover x/y from a pixel index
+    size_t npixel;
+    uint16_t roi_count;
+    std::vector<std::string> roi_name;  // ROI name indexed by bit position
+
+    // Result accumulators, filled by Run() in the derived class
+    std::vector<int64_t>  roi_sum;
+    std::vector<uint64_t> roi_sum2;
+    std::vector<uint64_t> roi_pixels;
+    std::vector<int64_t>  roi_x_weighted;
+    std::vector<int64_t>  roi_y_weighted;
+    std::vector<int64_t>  roi_max;
+
+    void Export(std::map<std::string, ROIMessage> &out) const;
+public:
+    explicit ROIIntegration(const DiffractionExperiment &experiment);
+    virtual ~ROIIntegration() = default;
+
+    [[nodiscard]] uint16_t Count() const { return roi_count; }
+    [[nodiscard]] bool empty() const { return roi_count == 0; }
+
+    virtual void Run(const ImagePreprocessorBuffer &image, std::map<std::string, ROIMessage> &out) = 0;
+};

image_analysis/roi/ROIIntegrationCPU.cpp

@@ -0,0 +1,11 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ROIIntegrationCPU.h"
+
+ROIIntegrationCPU::ROIIntegrationCPU(const DiffractionExperiment &experiment)
+    : ROIIntegration(experiment) {}
+
+void ROIIntegrationCPU::Run(const ImagePreprocessorBuffer &image, std::map<std::string, ROIMessage> &out) {
+    RunROI(image, out);
+}

image_analysis/roi/ROIIntegrationCPU.h

@@ -0,0 +1,71 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <limits>
+#include <type_traits>
+
+#include "ROIIntegration.h"
+#include "../../common/JFJochException.h"
+
+class ROIIntegrationCPU : public ROIIntegration {
+public:
+    explicit ROIIntegrationCPU(const DiffractionExperiment &experiment);
+
+    // image is anything indexable with operator[] and size(). Templated on the
+    // pixel type so a future narrow-integer (e.g. 16-bit) path works as well.
+    template <class T>
+    void RunROI(const T &image, std::map<std::string, ROIMessage> &out) {
+        if (image.size() != npixel)
+            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                                  "ROIIntegration: mismatch in image size");
+
+        for (uint16_t r = 0; r < roi_count; r++) {
+            roi_sum[r] = 0;
+            roi_sum2[r] = 0;
+            roi_pixels[r] = 0;
+            roi_x_weighted[r] = 0;
+            roi_y_weighted[r] = 0;
+            roi_max[r] = INT64_MIN;
+        }
+
+        using pixel_t = std::remove_cv_t<std::remove_reference_t<decltype(image[0])>>;
+
+        for (size_t i = 0; i < npixel; i++) {
+            const uint16_t mask = roi_map[i];
+            if (mask == 0)
+                continue;
+
+            const pixel_t v = image[i];
+            // masked/bad pixels (signed types only) are excluded entirely
+            if constexpr (std::is_signed_v<pixel_t>) {
+                if (v == std::numeric_limits<pixel_t>::min())
+                    continue;
+            }
+            // saturated pixels still count towards the max, but not the sums
+            const bool saturated = (v == std::numeric_limits<pixel_t>::max());
+
+            const int64_t val = static_cast<int64_t>(v);
+            const int64_t x = static_cast<int64_t>(i % width);
+            const int64_t y = static_cast<int64_t>(i / width);
+
+            for (uint16_t r = 0; r < roi_count; r++) {
+                if (!(mask & (1u << r)))
+                    continue;
+                if (!saturated) {
+                    roi_sum[r] += val;
+                    roi_sum2[r] += static_cast<uint64_t>(val * val);
+                    roi_pixels[r] += 1;
+                    roi_x_weighted[r] += val * x;
+                    roi_y_weighted[r] += val * y;
+                }
+                if (val > roi_max[r])
+                    roi_max[r] = val;
+            }
+        }
+        Export(out);
+    }
+
+    void Run(const ImagePreprocessorBuffer &image, std::map<std::string, ROIMessage> &out) override;
+};

image_analysis/roi/ROIIntegrationGPU.cu

@@ -0,0 +1,152 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <climits>
+
+#include "ROIIntegrationGPU.h"
+#include "../../common/DiffractionExperiment.h"
+
+inline void cuda_err(cudaError_t val) {
+    if (val != cudaSuccess)
+        throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val));
+}
+
+// One pixel carries a 16-bit mask, so it can feed any subset of the ROIs.
+// Each block reduces into shared memory first to keep global atomics low.
+__global__
+void gpu_roi(
+    const uint16_t *__restrict__ roi_map,
+    const int32_t  *__restrict__ input_buffer,
+    size_t num_pixels,
+    size_t width,
+    int roi_count,
+    unsigned long long *__restrict__ roi_sum,
+    unsigned long long *__restrict__ roi_sum2,
+    unsigned long long *__restrict__ roi_pixels,
+    unsigned long long *__restrict__ roi_x_weighted,
+    unsigned long long *__restrict__ roi_y_weighted,
+    int *__restrict__ roi_max) {
+    extern __shared__ unsigned long long shared[];
+    unsigned long long *s_sum    = shared;
+    unsigned long long *s_sum2   = &s_sum[roi_count];
+    unsigned long long *s_pixels = &s_sum2[roi_count];
+    unsigned long long *s_xw     = &s_pixels[roi_count];
+    unsigned long long *s_yw     = &s_xw[roi_count];
+    int *s_max = (int *) &s_yw[roi_count];
+
+    for (int r = threadIdx.x; r < roi_count; r += blockDim.x) {
+        s_sum[r] = 0;
+        s_sum2[r] = 0;
+        s_pixels[r] = 0;
+        s_xw[r] = 0;
+        s_yw[r] = 0;
+        s_max[r] = INT_MIN;
+    }
+    __syncthreads();
+
+    for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
+         idx < num_pixels;
+         idx += blockDim.x * gridDim.x) {
+        const uint16_t mask = roi_map[idx];
+        if (mask == 0)
+            continue;
+
+        const int32_t v = input_buffer[idx];
+        if (v == INT32_MIN)        // masked/bad pixel
+            continue;
+        const bool saturated = (v == INT32_MAX);
+
+        const long long val = v;
+        const long long x = idx % width;
+        const long long y = idx / width;
+        const unsigned long long val_u  = (unsigned long long) val;
+        const unsigned long long val2_u = (unsigned long long) (val * val);
+        const unsigned long long vx_u   = (unsigned long long) (val * x);
+        const unsigned long long vy_u   = (unsigned long long) (val * y);
+
+        for (int r = 0; r < roi_count; r++) {
+            if (!(mask & (1u << r)))
+                continue;
+            if (!saturated) {
+                atomicAdd(&s_sum[r], val_u);
+                atomicAdd(&s_sum2[r], val2_u);
+                atomicAdd(&s_pixels[r], 1ULL);
+                atomicAdd(&s_xw[r], vx_u);
+                atomicAdd(&s_yw[r], vy_u);
+            }
+            atomicMax(&s_max[r], v);
+        }
+    }
+    __syncthreads();
+
+    for (int r = threadIdx.x; r < roi_count; r += blockDim.x) {
+        atomicAdd(&roi_sum[r], s_sum[r]);
+        atomicAdd(&roi_sum2[r], s_sum2[r]);
+        atomicAdd(&roi_pixels[r], s_pixels[r]);
+        atomicAdd(&roi_x_weighted[r], s_xw[r]);
+        atomicAdd(&roi_y_weighted[r], s_yw[r]);
+        atomicMax(&roi_max[r], s_max[r]);
+    }
+}
+
+ROIIntegrationGPU::ROIIntegrationGPU(const DiffractionExperiment &experiment, std::shared_ptr<CudaStream> stream)
+    : ROIIntegration(experiment),
+      stream(stream),
+      gpu_roi_map(npixel),
+      gpu_sum(roi_count),
+      gpu_sum2(roi_count),
+      gpu_pixels(roi_count),
+      gpu_x_weighted(roi_count),
+      gpu_y_weighted(roi_count),
+      gpu_max(roi_count),
+      host_sum(roi_count),
+      host_sum2(roi_count),
+      host_pixels(roi_count),
+      host_x_weighted(roi_count),
+      host_y_weighted(roi_count),
+      host_max(roi_count),
+      max_init(roi_count, INT_MIN) {
+
+    cudaDeviceProp prop{};
+    cuda_err(cudaGetDeviceProperties(&prop, 0));
+    threads = 128;
+    blocks = 4 * prop.multiProcessorCount;
+    shared_needed = roi_count * (5 * sizeof(unsigned long long) + sizeof(int));
+
+    cuda_err(cudaMemcpy(gpu_roi_map, roi_map.data(), sizeof(uint16_t) * npixel, cudaMemcpyHostToDevice));
+}
+
+void ROIIntegrationGPU::Run(const ImagePreprocessorBuffer &image, std::map<std::string, ROIMessage> &out) {
+    if (image.size() != npixel)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "ROIIntegration: mismatch in image size");
+
+    cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * roi_count, *stream));
+    cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * roi_count, *stream));
+    cuda_err(cudaMemsetAsync(gpu_pixels, 0, sizeof(unsigned long long) * roi_count, *stream));
+    cuda_err(cudaMemsetAsync(gpu_x_weighted, 0, sizeof(unsigned long long) * roi_count, *stream));
+    cuda_err(cudaMemsetAsync(gpu_y_weighted, 0, sizeof(unsigned long long) * roi_count, *stream));
+    cuda_err(cudaMemcpyAsync(gpu_max, max_init.data(), sizeof(int) * roi_count, cudaMemcpyHostToDevice, *stream));
+
+    gpu_roi<<<blocks, threads, shared_needed, *stream>>>(
+        gpu_roi_map, image.getGPUBuffer(), npixel, width, roi_count,
+        gpu_sum, gpu_sum2, gpu_pixels, gpu_x_weighted, gpu_y_weighted, gpu_max);
+
+    cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cudaMemcpyAsync(host_pixels.data(), gpu_pixels, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cudaMemcpyAsync(host_x_weighted.data(), gpu_x_weighted, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cudaMemcpyAsync(host_y_weighted.data(), gpu_y_weighted, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cudaMemcpyAsync(host_max.data(), gpu_max, sizeof(int) * roi_count, cudaMemcpyDeviceToHost, *stream);
+    cuda_err(cudaStreamSynchronize(*stream));
+
+    for (uint16_t r = 0; r < roi_count; r++) {
+        roi_sum[r] = static_cast<int64_t>(host_sum[r]);
+        roi_sum2[r] = host_sum2[r];
+        roi_pixels[r] = host_pixels[r];
+        roi_x_weighted[r] = static_cast<int64_t>(host_x_weighted[r]);
+        roi_y_weighted[r] = static_cast<int64_t>(host_y_weighted[r]);
+        roi_max[r] = (host_max[r] == INT_MIN) ? INT64_MIN : static_cast<int64_t>(host_max[r]);
+    }
+    Export(out);
+}

image_analysis/roi/ROIIntegrationGPU.h

@@ -0,0 +1,37 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <memory>
+
+#include "ROIIntegration.h"
+#include "../indexing/CUDAMemHelpers.h"
+
+class ROIIntegrationGPU : public ROIIntegration {
+    std::shared_ptr<CudaStream> stream;
+    int threads;
+    int blocks;
+    size_t shared_needed;
+
+    CudaDevicePtr<uint16_t> gpu_roi_map;
+    // 64-bit sums are accumulated as unsigned long long (two's-complement bit
+    // pattern) because CUDA atomicAdd has no signed 64-bit overload.
+    CudaDevicePtr<unsigned long long> gpu_sum;
+    CudaDevicePtr<unsigned long long> gpu_sum2;
+    CudaDevicePtr<unsigned long long> gpu_pixels;
+    CudaDevicePtr<unsigned long long> gpu_x_weighted;
+    CudaDevicePtr<unsigned long long> gpu_y_weighted;
+    CudaDevicePtr<int> gpu_max;
+
+    std::vector<unsigned long long> host_sum;
+    std::vector<unsigned long long> host_sum2;
+    std::vector<unsigned long long> host_pixels;
+    std::vector<unsigned long long> host_x_weighted;
+    std::vector<unsigned long long> host_y_weighted;
+    std::vector<int> host_max;
+    std::vector<int> max_init;   // INT_MIN seed copied into gpu_max each frame
+public:
+    ROIIntegrationGPU(const DiffractionExperiment &experiment, std::shared_ptr<CudaStream> stream);
+    void Run(const ImagePreprocessorBuffer &image, std::map<std::string, ROIMessage> &out) override;
+};

tests/CMakeLists.txt

@@ -32,6 +32,7 @@ ADD_EXECUTABLE(jfjoch_test
         JPEGTest.cpp
         HistogramTest.cpp
         ROIMapTest.cpp
+        ROIIntegrationCPUTest.cpp
         LossyFilterTest.cpp
         ImageBufferTest.cpp
         PixelMaskTest.cpp

tests/ROIIntegrationCPUTest.cpp

@@ -0,0 +1,95 @@
+// 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 "../image_analysis/roi/ROIIntegrationCPU.h"
+#include "../common/DiffractionExperiment.h"
+
+TEST_CASE("ROIIntegrationCPU") {
+    DiffractionExperiment experiment(DetJF(1));
+    experiment.ROI().SetROI(ROIDefinition{.boxes = {
+        ROIBox("roiA", 10, 20, 30, 40),
+        ROIBox("roiB", 15, 25, 35, 45)
+    }});
+
+    ROIIntegrationCPU roi(experiment);
+    REQUIRE(roi.Count() == 2);
+    REQUIRE_FALSE(roi.empty());
+
+    const auto roi_map = experiment.ExportROIMap();
+    const size_t width = experiment.GetXPixelsNumConv();
+    const size_t npixel = roi_map.size();
+
+    ImagePreprocessorBuffer image(npixel);
+    for (size_t i = 0; i < npixel; i++)
+        image[i] = static_cast<int32_t>((i * 7 + 3) % 251); // deterministic, well below INT32_MAX
+
+    // Inject one saturated pixel into roiA (bit 0) and one masked pixel into roiB (bit 1)
+    int64_t saturated_index = -1;
+    int64_t masked_index = -1;
+    for (size_t i = 0; i < npixel && (saturated_index < 0 || masked_index < 0); i++) {
+        if (saturated_index < 0 && (roi_map[i] & (1 << 0)))
+            saturated_index = static_cast<int64_t>(i);
+        else if (masked_index < 0 && (roi_map[i] & (1 << 1)))
+            masked_index = static_cast<int64_t>(i);
+    }
+    REQUIRE(saturated_index >= 0);
+    REQUIRE(masked_index >= 0);
+    image[saturated_index] = INT32_MAX;
+    image[masked_index] = INT32_MIN;
+
+    // Independent reference matching the documented semantics (masked excluded
+    // entirely; saturated excluded from sums but counted towards the max)
+    struct Ref { int64_t sum = 0; uint64_t sum2 = 0; uint64_t pixels = 0;
+                 int64_t xw = 0; int64_t yw = 0; int64_t max = INT64_MIN; };
+    Ref ref[2];
+    for (size_t i = 0; i < npixel; i++) {
+        const uint16_t mask = roi_map[i];
+        if (mask == 0)
+            continue;
+        const int32_t v = image[i];
+        if (v == INT32_MIN)
+            continue;
+        const bool saturated = (v == INT32_MAX);
+        const int64_t val = v;
+        const int64_t x = static_cast<int64_t>(i % width);
+        const int64_t y = static_cast<int64_t>(i / width);
+        for (int r = 0; r < 2; r++) {
+            if (!(mask & (1 << r)))
+                continue;
+            if (!saturated) {
+                ref[r].sum += val;
+                ref[r].sum2 += static_cast<uint64_t>(val * val);
+                ref[r].pixels += 1;
+                ref[r].xw += val * x;
+                ref[r].yw += val * y;
+            }
+            if (val > ref[r].max)
+                ref[r].max = val;
+        }
+    }
+
+    std::map<std::string, ROIMessage> out;
+    roi.Run(image, out);
+
+    REQUIRE(out.size() == 2);
+    REQUIRE(out.contains("roiA"));
+    REQUIRE(out.contains("roiB"));
+
+    const std::pair<std::string, int> roi_ids[2] = {{"roiA", 0}, {"roiB", 1}};
+    for (const auto &[name, r] : roi_ids) {
+        const auto &msg = out.at(name);
+        CHECK(msg.sum == ref[r].sum);
+        CHECK(msg.sum_square == ref[r].sum2);
+        CHECK(msg.pixels == ref[r].pixels);
+        CHECK(msg.x_weighted == ref[r].xw);
+        CHECK(msg.y_weighted == ref[r].yw);
+        CHECK(msg.max_count == ref[r].max);
+        CHECK(msg.pixels_masked == 0);
+    }
+
+    // Targeted checks: saturated pixel sets the max for roiA but is not summed
+    CHECK(out.at("roiA").max_count == INT32_MAX);
+    CHECK(ref[0].pixels > 0);
+}