From 58910274bf257e8f53bd5ecaaab66a5f7d363073 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Wed, 17 Jun 2026 21:15:54 +0200 Subject: [PATCH] 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 --- image_analysis/CMakeLists.txt | 3 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 9 ++ image_analysis/MXAnalysisWithoutFPGA.h | 2 + image_analysis/roi/CMakeLists.txt | 5 + image_analysis/roi/ROIIntegration.cpp | 34 +++++ image_analysis/roi/ROIIntegration.h | 46 +++++++ image_analysis/roi/ROIIntegrationCPU.cpp | 11 ++ image_analysis/roi/ROIIntegrationCPU.h | 71 +++++++++++ image_analysis/roi/ROIIntegrationGPU.cu | 152 +++++++++++++++++++++++ image_analysis/roi/ROIIntegrationGPU.h | 37 ++++++ tests/CMakeLists.txt | 1 + tests/ROIIntegrationCPUTest.cpp | 95 ++++++++++++++ 12 files changed, 465 insertions(+), 1 deletion(-) create mode 100644 image_analysis/roi/CMakeLists.txt create mode 100644 image_analysis/roi/ROIIntegration.cpp create mode 100644 image_analysis/roi/ROIIntegration.h create mode 100644 image_analysis/roi/ROIIntegrationCPU.cpp create mode 100644 image_analysis/roi/ROIIntegrationCPU.h create mode 100644 image_analysis/roi/ROIIntegrationGPU.cu create mode 100644 image_analysis/roi/ROIIntegrationGPU.h create mode 100644 tests/ROIIntegrationCPUTest.cpp diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 1bf6e448..00996187 100644 --- a/image_analysis/CMakeLists.txt +++ b/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) diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 7bade492..cbe26e6a 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/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(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); azint = std::make_unique(integration); preprocessor = std::make_unique(in_experiment, in_mask); + if (experiment.ROI().size() >= 1) + roi = std::make_unique(experiment); #ifdef JFJOCH_USE_CUDA } else { auto stream = std::make_shared(); @@ -49,6 +53,8 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp preprocessor = std::make_unique(in_experiment, in_mask, stream); spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream); azint = std::make_unique(integration, stream); + if (experiment.ROI().size() >= 1) + roi = std::make_unique(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(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 diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index e5864d41..9d6bd0a9 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/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 azint; + std::unique_ptr roi; std::unique_ptr spotFinder; IndexAndRefine &indexer; std::unique_ptr prediction; diff --git a/image_analysis/roi/CMakeLists.txt b/image_analysis/roi/CMakeLists.txt new file mode 100644 index 00000000..afa1daa6 --- /dev/null +++ b/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() diff --git a/image_analysis/roi/ROIIntegration.cpp b/image_analysis/roi/ROIIntegration.cpp new file mode 100644 index 00000000..0cef8cca --- /dev/null +++ b/image_analysis/roi/ROIIntegration.cpp @@ -0,0 +1,34 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// 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 &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], + }; +} diff --git a/image_analysis/roi/ROIIntegration.h b/image_analysis/roi/ROIIntegration.h new file mode 100644 index 00000000..2d0ddd26 --- /dev/null +++ b/image_analysis/roi/ROIIntegration.h @@ -0,0 +1,46 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#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 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 roi_name; // ROI name indexed by bit position + + // Result accumulators, filled by Run() in the derived class + std::vector roi_sum; + std::vector roi_sum2; + std::vector roi_pixels; + std::vector roi_x_weighted; + std::vector roi_y_weighted; + std::vector roi_max; + + void Export(std::map &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 &out) = 0; +}; diff --git a/image_analysis/roi/ROIIntegrationCPU.cpp b/image_analysis/roi/ROIIntegrationCPU.cpp new file mode 100644 index 00000000..432b2ca4 --- /dev/null +++ b/image_analysis/roi/ROIIntegrationCPU.cpp @@ -0,0 +1,11 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// 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 &out) { + RunROI(image, out); +} diff --git a/image_analysis/roi/ROIIntegrationCPU.h b/image_analysis/roi/ROIIntegrationCPU.h new file mode 100644 index 00000000..0e5fb4eb --- /dev/null +++ b/image_analysis/roi/ROIIntegrationCPU.h @@ -0,0 +1,71 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#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 + void RunROI(const T &image, std::map &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>; + + 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) { + if (v == std::numeric_limits::min()) + continue; + } + // saturated pixels still count towards the max, but not the sums + const bool saturated = (v == std::numeric_limits::max()); + + const int64_t val = static_cast(v); + const int64_t x = static_cast(i % width); + const int64_t y = static_cast(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(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 &out) override; +}; diff --git a/image_analysis/roi/ROIIntegrationGPU.cu b/image_analysis/roi/ROIIntegrationGPU.cu new file mode 100644 index 00000000..66701cfa --- /dev/null +++ b/image_analysis/roi/ROIIntegrationGPU.cu @@ -0,0 +1,152 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#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 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 &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<<>>( + 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(host_sum[r]); + roi_sum2[r] = host_sum2[r]; + roi_pixels[r] = host_pixels[r]; + roi_x_weighted[r] = static_cast(host_x_weighted[r]); + roi_y_weighted[r] = static_cast(host_y_weighted[r]); + roi_max[r] = (host_max[r] == INT_MIN) ? INT64_MIN : static_cast(host_max[r]); + } + Export(out); +} diff --git a/image_analysis/roi/ROIIntegrationGPU.h b/image_analysis/roi/ROIIntegrationGPU.h new file mode 100644 index 00000000..5528b37b --- /dev/null +++ b/image_analysis/roi/ROIIntegrationGPU.h @@ -0,0 +1,37 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "ROIIntegration.h" +#include "../indexing/CUDAMemHelpers.h" + +class ROIIntegrationGPU : public ROIIntegration { + std::shared_ptr stream; + int threads; + int blocks; + size_t shared_needed; + + CudaDevicePtr 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 gpu_sum; + CudaDevicePtr gpu_sum2; + CudaDevicePtr gpu_pixels; + CudaDevicePtr gpu_x_weighted; + CudaDevicePtr gpu_y_weighted; + CudaDevicePtr gpu_max; + + std::vector host_sum; + std::vector host_sum2; + std::vector host_pixels; + std::vector host_x_weighted; + std::vector host_y_weighted; + std::vector host_max; + std::vector max_init; // INT_MIN seed copied into gpu_max each frame +public: + ROIIntegrationGPU(const DiffractionExperiment &experiment, std::shared_ptr stream); + void Run(const ImagePreprocessorBuffer &image, std::map &out) override; +}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 751e8f64..59af6e96 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/ROIIntegrationCPUTest.cpp b/tests/ROIIntegrationCPUTest.cpp new file mode 100644 index 00000000..747d9b16 --- /dev/null +++ b/tests/ROIIntegrationCPUTest.cpp @@ -0,0 +1,95 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#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((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(i); + else if (masked_index < 0 && (roi_map[i] & (1 << 1))) + masked_index = static_cast(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(i % width); + const int64_t y = static_cast(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(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 out; + roi.Run(image, out); + + REQUIRE(out.size() == 2); + REQUIRE(out.contains("roiA")); + REQUIRE(out.contains("roiB")); + + const std::pair 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); +}