diff --git a/common/DiffractionSpot.cpp b/common/DiffractionSpot.cpp index 5e6a8857..bb0c473f 100644 --- a/common/DiffractionSpot.cpp +++ b/common/DiffractionSpot.cpp @@ -23,6 +23,10 @@ DiffractionSpot& DiffractionSpot::operator+=(const DiffractionSpot &other) { return *this; } +DiffractionSpot::DiffractionSpot(float x_sum, float y_sum, int64_t in_pixel_count, + int64_t in_photons, int64_t in_max_photons) + : x(x_sum), y(y_sum), pixel_count(in_pixel_count), photons(in_photons), max_photons(in_max_photons) {} + DiffractionSpot::DiffractionSpot(const SpotToSave &save) { x = save.x * static_cast(save.intensity); y = save.y * static_cast(save.intensity); diff --git a/common/DiffractionSpot.h b/common/DiffractionSpot.h index 2614b4c9..e98fc34b 100644 --- a/common/DiffractionSpot.h +++ b/common/DiffractionSpot.h @@ -19,6 +19,10 @@ public: DiffractionSpot() = default; DiffractionSpot(uint32_t col, uint32_t line, int64_t photons); DiffractionSpot(const SpotToSave &save); + // From already-summed quantities. x_sum/y_sum are sum(col * photons) / sum(line * photons), i.e. + // the members as AddPixel leaves them, NOT a centroid. Used by the GPU spot extractor, which + // computes those sums on the device. + DiffractionSpot(float x_sum, float y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); DiffractionSpot& operator+=(const DiffractionSpot& spot); int64_t PixelCount() const; int64_t Count() const; diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index f934ffc9..feafef4e 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -223,13 +223,26 @@ Strong pixels are grouped into connected components (adjacent strong pixels) usi Spot-level filters include minimum/maximum pixel count and resolution limits. +The host implementation (`StrongPixelSet::sparseccl`) is the SparseCCL of the ACTS/traccc project: it +runs over the strong pixels sorted row-major, uses a sliding window over the previous line and a +union-find whose root is each component's lowest index. On the GPU the same labelling runs **on the +device** (`SpotExtractorGPU`): the packed strong-pixel bitmask is compacted into that same sorted +list without leaving the card, each pixel finds its at most four earlier 8-neighbours by binary +search, and a lock-free union-find with path halving labels them. Only the finished spot list — a few +hundred entries — comes back to the host, instead of the whole bitmask (2.26 MB per frame at 18 MP). +The two implementations produce the same components, in the same order, with the same pixel counts +and intensities; `tests/SpotExtractorGPUParityTest.cpp` holds them to it. The device version is also +insensitive to frame content: the host sliding window becomes quadratic when many pixels light up in +one detector line — a hot module, or a diffraction ring where it runs tangent to a row — which costs +tens to hundreds of milliseconds on such a frame, while the device version stays under a millisecond. + ### 3.5 Adaptive per-image minimum spot size The minimum-pixels-per-spot filter (§3.4) trades sensitivity against noise: a small value keeps faint one- or two-pixel spots — real signal on strong data, but detector noise on high-background frames — while a larger value keeps only well-formed spots. The best value is dataset-dependent, so for serial-stills indexing it can be chosen **per image** rather than fixed. The frame is indexed three times, at min-pix 3, 2 and 1, and the setting that maximises $$ \frac{n_\mathrm{indexed}^2}{n_\mathrm{total}} \quad\text{(indexed-spot count weighted by indexed fraction)} $$ -is kept; the frame is then integrated once at that min-pix. The fraction factor discounts the extra spots a smaller min-pix admits *unless the lattice actually explains them*, so strong frames keep their real weak spots (extending resolution) while noise-flooded frames stay strict. Because min-pix filters the connected components *after* detection, strong-pixel detection runs **once** per frame and the three attempts only repeat the cheap CCL and spot-level filter, not the pixel reduction; the azimuthal profile is the one that single detection pass computed. The winning attempt's spot list is kept rather than re-extracted, so the frame that is integrated is exactly the frame that was scored. This is a **stills-only, indexing-path** option — rotation indexing builds one global lattice from all frames and keeps a fixed min-pix. In `rugnux` it is the default; giving an explicit `--min-pix-per-spot` pins a fixed value instead. +is kept; the frame is then integrated once at that min-pix. The fraction factor discounts the extra spots a smaller min-pix admits *unless the lattice actually explains them*, so strong frames keep their real weak spots (extending resolution) while noise-flooded frames stay strict. Because min-pix filters the connected components *after* detection, strong-pixel detection AND the connected-component labelling both run **once** per frame, and the three attempts only repeat the spot-level filter; the azimuthal profile is the one that single detection pass computed. The winning attempt's spot list is kept rather than re-extracted, so the frame that is integrated is exactly the frame that was scored. This is a **stills-only, indexing-path** option — rotation indexing builds one global lattice from all frames and keeps a fixed min-pix. In `rugnux` it is the default; giving an explicit `--min-pix-per-spot` pins a fixed value instead. --- diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index e4a54303..e300a187 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -139,14 +139,15 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, && !experiment.IsRotationIndexing(); if (adaptive_min_pix) { // Choose the per-image min-pix adaptively instead of a fixed one. min-pix filters connected - // components AFTER detection, so detection (the expensive per-pixel pass) runs ONCE and only - // the cheap CCL + spot filter is repeated; the azimuthal profile is the one Detect() computed. + // components AFTER detection, so BOTH the detection (the expensive per-pixel pass) and the + // connected-component search run ONCE and only the filter is repeated; the azimuthal + // profile is the one Detect() computed. // Index at 3/2/1 (index-only, no integration/accumulation), keep whichever maximises // n_indexed^2 / n_total (indexed count weighted by indexed fraction) together with its spot - // list, and integrate that one. Keeping the list also means the frame integrated is exactly - // the frame scored, which re-extracting could not guarantee on the GPU (atomic-order sums). + // list, and integrate that one. const auto detect_start_time = std::chrono::steady_clock::now(); finder.Detect(*preprocessor_buffer, spot_finding_settings); + const auto &components = finder.ExtractComponents(*preprocessor_buffer, spot_finding_settings); float spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - detect_start_time).count(); float indexing_time_s = 0.0f; @@ -158,7 +159,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, for (int mp : {3, 2, 1}) { s.min_pix_per_spot = mp; const auto extract_start_time = std::chrono::steady_clock::now(); - std::vector spots = finder.ExtractSpots(*preprocessor_buffer, s, mask_resolution); + std::vector spots = ImageSpotFinder::Filter(components, s); spot_finding_time_s += std::chrono::duration(std::chrono::steady_clock::now() - extract_start_time).count(); SpotAnalyze(experiment, s, spots, output); @@ -187,7 +188,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, output.indexing_time_s = indexing_time_s; } else { const auto spot_finding_start_time = std::chrono::steady_clock::now(); - const std::vector spots = finder.Run(*preprocessor_buffer, spot_finding_settings, mask_resolution); + const std::vector spots = finder.Run(*preprocessor_buffer, spot_finding_settings); SpotAnalyze(experiment, spot_finding_settings, spots, output); output.spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - spot_finding_start_time).count(); if (spot_finding_settings.indexing) @@ -266,4 +267,9 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett auto const &resolution_map = integration.Resolution(); for (int i = 0; i < mask_resolution.size(); i++) mask_resolution[i] = (resolution_map[i] > mask_low_res) || (resolution_map[i] < high_res); + + // The finders keep their own copy (the GPU ones a bit-packed device copy), so the mask is handed + // over here - when the limits change - rather than with every image. + spotFinder->SetResolutionMask(mask_resolution); + adaptiveSpotFinder->SetResolutionMask(mask_resolution); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 78379d2b..2ff14c9f 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -165,9 +165,9 @@ __global__ void flag_strong(const int32_t *__restrict__ image, AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping, std::shared_ptr in_stream) : ImageSpotFinder(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight())), + static_cast(in_mapping.GetHeight()), false), mapping(in_mapping), - stream(std::move(in_stream)), + stream(in_stream), nbins(in_mapping.GetBinNumber()), npix(in_mapping.GetPixelToBin().size()), gpu_sum(nbins), @@ -185,7 +185,8 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & prof_sum(nbins), prof_sum2(nbins), prof_count(nbins), - output_buffer_reg(output_buffer), + extractor(static_cast(in_mapping.GetWidth()), + static_cast(in_mapping.GetHeight()), std::move(in_stream)), last_profile(in_mapping) { // The current device, not device 0: callers round-robin engines across GPUs, so device 0's shared @@ -314,8 +315,9 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, last_profile.Add(prof_sum, prof_sum2, prof_count); if (host_thr.empty()) { - // Nothing valid to threshold against: leave no strong pixels for ExtractSpots to build on. - std::fill(output_buffer.begin(), output_buffer.end(), 0); + // Nothing valid to threshold against: leave no strong pixels for the extractor to build on. + cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); + cuda_err(cudaStreamSynchronize(*stream)); return; } @@ -324,6 +326,17 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); flag_strong<<>>( image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins); - cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_strong, OutputByteSize(), cudaMemcpyDeviceToHost, *stream)); + // The bit buffer stays on the device - ExtractComponents reads it there. cuda_err(cudaStreamSynchronize(*stream)); } + +void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector &mask) { + ImageSpotFinder::SetResolutionMask(mask); + extractor.SetResolutionMask(res_mask_bits); +} + +const std::vector &AdaptiveSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { + extractor.Extract(gpu_strong, image.getGPUBuffer(), settings, components); + return components; +} diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index 999327ba..7b5d2700 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -12,8 +12,8 @@ // - the azimuthal-integration profile (mean intensity per ring, in flat-field-corrected space), and // - the per-ring background (mean, sigma, peak-excluded via two sigma-clip passes) that sets the // self-calibrating spot-detection threshold (in raw photon counts). -// It then flags strong pixels (value >= ring threshold) into a packed bit buffer and hands it to the -// shared host connected-component extractor (ImageSpotFinder::ExtractSpots). +// It then flags strong pixels (value >= ring threshold) into a packed bit buffer and hands that +// buffer - still on the device - to SpotExtractorGPU, which builds the spots there. // // Numerically it reproduces AdaptiveSpotFinderCPU: the same three-pass robust background, the same // per-ring threshold formula (shared via AdaptiveThreshold.h, computed on the host once per frame), @@ -26,6 +26,7 @@ #include #include "ImageSpotFinder.h" +#include "SpotExtractorGPU.h" #include "SpotFindingSettings.h" #include "../../common/AzimuthalIntegrationProfile.h" #include "../../common/AzimuthalIntegrationMapping.h" @@ -78,7 +79,7 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { std::vector prof_sum2; // plain corrected sum^2 } std::vector prof_count; // plain pixel count } - CudaRegisteredVector output_buffer_reg; // pins the base-class bit buffer for fast D2H + SpotExtractorGPU extractor; // builds the spots from gpu_strong without it leaving the device AzimuthalIntegrationProfile last_profile; // filled every Run(), retrievable via GetProfile() @@ -98,6 +99,9 @@ public: AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete; void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; + void SetResolutionMask(const std::vector &mask) override; + const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) override; // The azimuthal profile computed as a byproduct of the last Detect() - lets this engine replace the // separate azint pass in the analysis pipeline. diff --git a/image_analysis/spot_finding/CMakeLists.txt b/image_analysis/spot_finding/CMakeLists.txt index 8164f10d..0d0d6035 100644 --- a/image_analysis/spot_finding/CMakeLists.txt +++ b/image_analysis/spot_finding/CMakeLists.txt @@ -18,6 +18,7 @@ TARGET_LINK_LIBRARIES(JFJochSpotFinding JFJochCommon) IF (JFJOCH_CUDA_AVAILABLE) TARGET_SOURCES(JFJochSpotFinding PRIVATE + SpotExtractorGPU.cu SpotExtractorGPU.h ImageSpotFinderGPU.cu ImageSpotFinderGPU.h AdaptiveSpotFinderGPU.cu AdaptiveSpotFinderGPU.h) ENDIF() \ No newline at end of file diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 8296d325..74f57d5e 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -1,15 +1,23 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#include +#include +#include +#include "../../common/JFJochException.h" #include "ImageSpotFinder.h" #include "StrongPixelSet.h" -ImageSpotFinder::ImageSpotFinder(int32_t width, int32_t height) +ImageSpotFinder::ImageSpotFinder(int32_t width, int32_t height, bool host_bit_buffer) : width(width), height(height), - output_buffer(width * height / 32 + 1) { + output_buffer(host_bit_buffer ? width * height / 32 + 1 : 0), + res_mask_bits(OutputSize(), 0) { + // Exclude the padding bits of the last word up front, so neither the host scan nor the GPU + // compaction needs a separate "is this bit still inside the image?" test. + const size_t npixel = static_cast(width) * height; + if (npixel % 32 != 0) + res_mask_bits.back() = ~((1u << (npixel % 32)) - 1u); } size_t ImageSpotFinder::OutputSize() const { @@ -20,44 +28,72 @@ size_t ImageSpotFinder::OutputByteSize() const { return OutputSize() * sizeof(uint32_t); } -std::vector ImageSpotFinder::Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) { - Detect(image, settings); - return ExtractSpots(image, settings, res_mask); +void ImageSpotFinder::SetResolutionMask(const std::vector &mask) { + const size_t npixel = static_cast(width) * height; + if (mask.size() != npixel) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ImageSpotFinder::SetResolutionMask: mask size mismatch"); + std::fill(res_mask_bits.begin(), res_mask_bits.end(), 0); + for (size_t i = 0; i < npixel; i++) + if (mask[i]) + res_mask_bits[i / 32] |= 1u << (i % 32); + if (npixel % 32 != 0) + res_mask_bits.back() |= ~((1u << (npixel % 32)) - 1u); } -std::vector ImageSpotFinder::ExtractSpots(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) { +void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { // Collect the strong pixels first and read their values afterwards, instead of reading the image // pixel by pixel: on the GPU that read is a device gather, which is what lets the preprocessed // image stay on the device instead of being copied back in full for every frame. strong_pixel.clear(); - for (int i = 0; i < OutputSize(); i++) { - if (output_buffer[i]) { - std::bitset<32> bset = output_buffer[i]; - for (int bit = 0; bit < 32; bit++) { - if (bset.test(bit)) { - uint32_t npixel = i * 32 + bit; - if (npixel / width < height && res_mask[npixel] == 0) - strong_pixel.push_back(npixel); - } - } + for (size_t i = 0; i < OutputSize(); i++) { + // The resolution mask is packed like the bit buffer, so a whole word of it is excluded at + // once instead of testing 32 bits one at a time. + uint32_t word = output_buffer[i] & ~res_mask_bits[i]; + while (word != 0) { + strong_pixel.push_back(static_cast(i * 32 + std::countr_zero(word))); + word &= word - 1; } } - std::vector vec; + components.clear(); // The connected-component search rejects a frame with this many strong pixels, so their values are // of no use. if (strong_pixel.size() >= UINT16_MAX) - return vec; + return; image.Gather(strong_pixel, strong_pixel_value); StrongPixelSet pixel_set; for (size_t i = 0; i < strong_pixel.size(); i++) pixel_set.AddStrongPixel(strong_pixel[i] % width, strong_pixel[i] / width, strong_pixel_value[i]); - pixel_set.FindSpotsImage(settings, vec); - return vec; + pixel_set.FindComponentsImage(settings, components); +} + +const std::vector &ImageSpotFinder::ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { + ExtractComponentsHost(image, settings); + return components; +} + +std::vector ImageSpotFinder::Filter(const std::vector &in, + const SpotFindingSettings &settings) { + std::vector out; + const int64_t min_pix = settings.min_pix_per_spot.value_or(2); + for (const auto &spot: in) + if (spot.PixelCount() >= min_pix) + out.push_back(spot); + return out; +} + +std::vector ImageSpotFinder::ExtractSpots(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { + return Filter(ExtractComponents(image, settings), settings); +} + +std::vector ImageSpotFinder::Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { + Detect(image, settings); + return ExtractSpots(image, settings); } diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index 41344c72..d03dc9db 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -18,10 +18,22 @@ class ImageSpotFinder { protected: const int32_t width, height; std::vector output_buffer; + // Pixels excluded from spot finding, packed the same way as output_buffer (bit set = excluded). + // The bits past the last image pixel are set at construction, so the padding of the last word + // needs no separate guard. Default: nothing excluded. + std::vector res_mask_bits; + // The connected components of the last extraction. A member so ExtractComponents can hand out a + // reference and reuse the allocation from image to image. + std::vector components; - ImageSpotFinder(int32_t width, int32_t height); + // host_bit_buffer = false leaves output_buffer empty: the GPU finders extract on the device and + // never read the bit buffer on the host, so allocating and pinning 2.26 MB per engine (at 18 MP) + // would be pure waste for them. + ImageSpotFinder(int32_t width, int32_t height, bool host_bit_buffer = true); size_t OutputSize() const; size_t OutputByteSize() const; + // Host extraction: scan the bit buffer, gather the values, run the connected-component search. + void ExtractComponentsHost(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings); public: constexpr static int32_t MIN_VALID_PIXELS = 100; constexpr static int NBX = 15; @@ -29,13 +41,24 @@ public: virtual ~ImageSpotFinder() = default; // Detect flags the image's strong pixels into the internal bit buffer - the expensive step (local - // box or per-ring background over every pixel). ExtractSpots then builds the spots from those - // pixels; min/max-pix and the resolution mask enter only there, so several min-pix values can be - // tried on ONE detection pass (MXAnalysisWithoutFPGA does that when min-pix is chosen per image). + // box or per-ring background over every pixel). ExtractComponents then builds the connected + // components from those pixels. virtual void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) = 0; - std::vector ExtractSpots(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); - std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); + + // Pixels to ignore, one bool per pixel (true = ignore). Set when the resolution limits change, + // not per image: the GPU finders keep a bit-packed device copy of it, and re-uploading that for + // every image would cost more than the extraction it feeds. + virtual void SetResolutionMask(const std::vector &mask); + + // Every connected component of the last Detect() with at most max-pix pixels. min-pix is NOT + // applied here on purpose - it is the only spot setting that changes between the passes of the + // per-image min-pix search, so ONE extraction serves all three of them. + virtual const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings); + // The components that also pass min-pix. + static std::vector Filter(const std::vector &in, + const SpotFindingSettings &settings); + + std::vector ExtractSpots(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings); + std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings); }; - - - diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index b535d8e1..4219bc96 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -238,14 +238,25 @@ __global__ void analyze_pixel(const int32_t *in, uint32_t *prev_out, uint32_t *o } ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height, - std::shared_ptr stream) - : ImageSpotFinder(in_width, in_height), - stream(std::move(stream)), - output_buffer_reg(output_buffer) { + std::shared_ptr in_stream) + : ImageSpotFinder(in_width, in_height, false), + stream(in_stream), + extractor(in_width, in_height, std::move(in_stream)) { gpu_out_0 = CudaDevicePtr(OutputSize()); gpu_out_1 = CudaDevicePtr(OutputSize()); } +void ImageSpotFinderGPU::SetResolutionMask(const std::vector &mask) { + ImageSpotFinder::SetResolutionMask(mask); + extractor.SetResolutionMask(res_mask_bits); +} + +const std::vector &ImageSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { + extractor.Extract(gpu_out_1, image.getGPUBuffer(), settings, components); + return components; +} + void ImageSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { spot_parameters spot_params{}; spot_params.height = height; @@ -279,7 +290,6 @@ void ImageSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, const Spot (image.getGPUBuffer(), gpu_out_1, gpu_out_0, spot_params); analyze_pixel<<>> (image.getGPUBuffer(), gpu_out_0, gpu_out_1, spot_params); - cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_out_1, OutputSize() * sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); - + // The bit buffer stays on the device - ExtractComponents reads it there. cuda_err(cudaStreamSynchronize(*stream)); } diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index 0dab3a74..5f13f238 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -7,14 +7,15 @@ #include "SpotFindingSettings.h" #include "ImageSpotFinder.h" +#include "SpotExtractorGPU.h" #include "../indexing/CUDAMemHelpers.h" class ImageSpotFinderGPU : public ImageSpotFinder { std::shared_ptr stream; CudaDevicePtr gpu_out_0; - CudaDevicePtr gpu_out_1; - CudaRegisteredVector output_buffer_reg; + CudaDevicePtr gpu_out_1; // holds the strong-pixel bits after Detect() + SpotExtractorGPU extractor; const int numberOfCudaThreads = 128; // #threads per block that should work well for Nvidia L4 const int numberOfWaves = 32; // #waves that should work well for Nvidia L4 @@ -24,6 +25,9 @@ public: ~ImageSpotFinderGPU() override = default; void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; + void SetResolutionMask(const std::vector &mask) override; + const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) override; }; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu new file mode 100644 index 00000000..09a1ac7f --- /dev/null +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -0,0 +1,351 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// Sparse connected-component labelling adapted from the design the ACTS/traccc project arrived at +// for sparse silicon-detector clusterization (backward-neighbour graph over a sorted hit list, then +// a parallel union-find), which is itself the GPU counterpart of the SparseCCL that +// StrongPixelSet.cpp uses on the host. + +#include + +#include "SpotExtractorGPU.h" +#include "../../common/JFJochException.h" + +namespace { + +inline void cuda_err(cudaError_t val) { + if (val != cudaSuccess) + throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val)); +} + +constexpr int THREADS = 256; +constexpr int FINISH_THREADS = 1024; + +// --- compaction: packed bit buffer -> strong-pixel list, sorted by flat index ------------------- +// Each block owns a CONTIGUOUS range of words. Pass 1 counts its strong bits, a single-block scan +// turns the counts into offsets, and pass 2 walks the same range in increasing order and writes at +// that offset. No atomics anywhere, which is what keeps the output sorted. +__global__ void count_bits(const uint32_t *__restrict__ strong, const uint32_t *__restrict__ res_mask, + uint32_t *__restrict__ block_count, size_t nwords) { + const size_t per_block = (nwords + gridDim.x - 1) / gridDim.x; + const size_t w0 = static_cast(blockIdx.x) * per_block; + const size_t w1 = min(w0 + per_block, nwords); + uint32_t local = 0; + for (size_t w = w0 + threadIdx.x; w < w1; w += blockDim.x) + local += __popc(strong[w] & ~res_mask[w]); + + __shared__ uint32_t block_total; + if (threadIdx.x == 0) block_total = 0; + __syncthreads(); + atomicAdd(&block_total, local); + __syncthreads(); + if (threadIdx.x == 0) block_count[blockIdx.x] = block_total; +} + +__global__ void scan_block_counts(const uint32_t *__restrict__ in, uint32_t *__restrict__ out, + uint32_t *__restrict__ total, int n) { + __shared__ uint32_t shared[FINISH_THREADS]; + const int t = threadIdx.x, nthreads = blockDim.x; + const int chunk = (n + nthreads - 1) / nthreads; + const int lo = min(t * chunk, n), hi = min(lo + chunk, n); + uint32_t sum = 0; + for (int i = lo; i < hi; i++) sum += in[i]; + shared[t] = sum; + __syncthreads(); + for (int d = 1; d < nthreads; d <<= 1) { + const uint32_t v = (t >= d) ? shared[t - d] : 0u; + __syncthreads(); + shared[t] += v; + __syncthreads(); + } + uint32_t acc = shared[t] - sum; + for (int i = lo; i < hi; i++) { out[i] = acc; acc += in[i]; } + if (t == nthreads - 1) *total = shared[t]; +} + +// One thread per block emits its range. Strong pixels are ~1e-4 of the image, so a block's range +// holds a handful of them and serial emission is both trivially ordered and fast; the parallelism +// comes from the block count. +__global__ void scatter_bits(const uint32_t *__restrict__ strong, const uint32_t *__restrict__ res_mask, + const uint32_t *__restrict__ block_offset, const int32_t *__restrict__ image, + uint32_t *__restrict__ out_index, int32_t *__restrict__ out_value, + size_t nwords, uint32_t capacity) { + if (threadIdx.x != 0) return; + const size_t per_block = (nwords + gridDim.x - 1) / gridDim.x; + const size_t w0 = static_cast(blockIdx.x) * per_block; + const size_t w1 = min(w0 + per_block, nwords); + uint32_t pos = block_offset[blockIdx.x]; + for (size_t w = w0; w < w1; ++w) { + uint32_t word = strong[w] & ~res_mask[w]; + while (word != 0) { + const uint32_t flat = static_cast(w * 32 + (__ffs(word) - 1)); + word &= word - 1; + if (pos < capacity) { + out_index[pos] = flat; + out_value[pos] = image[flat]; + } + ++pos; + } + } +} + +// --- connected components ---------------------------------------------------------------------- +__device__ __forceinline__ int lower_bound_device(const uint32_t *a, int n, uint32_t key) { + int lo = 0, hi = n; + while (lo < hi) { + const int mid = (lo + hi) >> 1; + if (a[mid] < key) lo = mid + 1; else hi = mid; + } + return lo; +} + +// Walk to the root, halving the path on the way. The plain store is safe: a parent only ever +// decreases and the grandparent is always still an ancestor, so a concurrent writer can only make +// the chain shorter. Without halving, a long thin feature - a diffraction ring is exactly one - +// builds a parent chain as long as the feature itself and every later merge walks all of it. +__device__ __forceinline__ uint32_t find_root(uint32_t *parent, uint32_t a) { + uint32_t p = parent[a]; + while (p != a) { + const uint32_t gp = parent[p]; + if (gp == p) return p; + parent[a] = gp; + a = gp; + p = parent[a]; + } + return a; +} + +// Lock-free union. Terminates because max(a,b) strictly decreases; converges on the component's +// LOWEST index as its root, which is what the host's make_union does too. +__device__ __forceinline__ void merge_roots(uint32_t *parent, uint32_t a, uint32_t b) { + a = find_root(parent, a); + b = find_root(parent, b); + while (a != b) { + if (a < b) { const uint32_t t = a; a = b; b = t; } + const uint32_t old = atomicMin(&parent[a], b); + if (old == a) return; // a was a root and now points at b: done + a = find_root(parent, old); // someone re-parented a; carry on from its root + b = find_root(parent, b); + } +} + +__global__ void init_parent(uint32_t *__restrict__ parent, const uint32_t *__restrict__ nstrong, + uint32_t capacity) { + const int n = static_cast(min(*nstrong, capacity)); + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) + parent[i] = i; +} + +// The list is sorted by flat index, so a pixel's 8-neighbours that come EARLIER in it are exactly +// four: (line, col-1), (line-1, col-1), (line-1, col) and (line-1, col+1). Each is one binary +// search away, which is what makes the sparse formulation cheap. +__global__ void union_neighbours(const uint32_t *__restrict__ index, uint32_t *__restrict__ parent, + const uint32_t *__restrict__ nstrong, uint32_t capacity, int width) { + const int n = static_cast(min(*nstrong, capacity)); + if (static_cast(n) >= capacity) return; + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) { + const uint32_t flat = index[i]; + const uint32_t col = flat % width; + const uint32_t line = flat / width; + uint32_t candidate[4]; + int ncandidate = 0; + if (col > 0) candidate[ncandidate++] = flat - 1; + if (line > 0 && col > 0) candidate[ncandidate++] = flat - width - 1; + if (line > 0) candidate[ncandidate++] = flat - width; + if (line > 0 && col + 1 < static_cast(width)) candidate[ncandidate++] = flat - width + 1; + for (int k = 0; k < ncandidate; k++) { + const int j = lower_bound_device(index, i, candidate[k]); + if (j < i && index[j] == candidate[k]) merge_roots(parent, static_cast(i), static_cast(j)); + } + } +} + +__global__ void resolve_roots(uint32_t *__restrict__ parent, uint32_t *__restrict__ root, + const uint32_t *__restrict__ nstrong, uint32_t capacity) { + const int n = static_cast(min(*nstrong, capacity)); + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) + root[i] = find_root(parent, static_cast(i)); +} + +// Everything after the labelling in ONE block, so a frame needs a single host synchronisation: +// hand out labels, count the members of each component, sum the surviving ones, filter by max-pix +// and compact - all of it order-preserving. +__global__ void finish_components(const uint32_t *__restrict__ index, const int32_t *__restrict__ value, + const uint32_t *__restrict__ root, uint32_t *__restrict__ label, + int32_t *__restrict__ count, SpotExtractorGPUSpot *__restrict__ scratch, + SpotExtractorGPUSpot *__restrict__ out, uint32_t *__restrict__ nout, + const uint32_t *__restrict__ nstrong, uint32_t capacity, + int width, int max_pix) { + const int n = static_cast(min(*nstrong, capacity)); + if (threadIdx.x == 0) *nout = 0; + __syncthreads(); + // Same give-up as StrongPixelSet::FindComponentsImage - except that here the count is known + // before a single pixel has been written anywhere, so the frame costs nothing to reject. + if (n == 0 || static_cast(n) >= capacity) return; + + __shared__ uint32_t shared[FINISH_THREADS]; + const int t = threadIdx.x, nthreads = blockDim.x; + const int chunk = (n + nthreads - 1) / nthreads; + const int lo = min(t * chunk, n), hi = min(lo + chunk, n); + + // 1) labels, by a prefix sum over the roots in ascending order - the order the host's second + // scan hands them out in, which is what makes the spot ORDER identical. + uint32_t nroot = 0; + for (int i = lo; i < hi; i++) nroot += (root[i] == static_cast(i)) ? 1u : 0u; + shared[t] = nroot; + __syncthreads(); + for (int d = 1; d < nthreads; d <<= 1) { + const uint32_t v = (t >= d) ? shared[t - d] : 0u; + __syncthreads(); + shared[t] += v; + __syncthreads(); + } + uint32_t next_label = shared[t] - nroot; + for (int i = lo; i < hi; i++) + if (root[i] == static_cast(i)) label[i] = next_label++; + const int nlabel = static_cast(shared[nthreads - 1]); + __syncthreads(); + + // 2) member counts + for (int i = t; i < nlabel; i += nthreads) count[i] = 0; + __syncthreads(); + for (int i = t; i < n; i += nthreads) atomicAdd(&count[label[root[i]]], 1); + __syncthreads(); + + // 3) sums, one thread per component, walking its members in ascending list order so the float + // accumulation matches DiffractionSpot::AddPixel term for term. A component bigger than + // max-pix is thrown away below, so it is not summed - which is also what keeps a whole lit + // module or diffraction ring from turning into one thread walking tens of thousands of + // entries. + for (int i = t; i < n; i += nthreads) { + if (root[i] != static_cast(i)) continue; + const uint32_t l = label[i]; + const int want = count[l]; + scratch[l].pixel_count = want; + if (want > max_pix) continue; + float x = 0.0f, y = 0.0f; + long long photons = 0, max_photons = LLONG_MIN; + int found = 0; + for (int j = i; j < n && found < want; j++) { + if (root[j] != static_cast(i)) continue; + const long long counts = value[j]; + // Spelled out rather than left as "x += col * counts", because the rounding has to match + // DiffractionSpot::AddPixel term for term and the two compilers do not contract alike: + // gcc fuses AddPixel into a vfmadd whenever the build enables FMA (the CI flags do), so + // __fmaf_rn is the counterpart there. Where the host cannot fuse (a baseline -march, or + // MSVC, which does not contract by default) this leaves the last bit of the centroid + // differing by an ulp - see SpotExtractorGPUParityTest.cpp. + x = __fmaf_rn(static_cast(index[j] % width), static_cast(counts), x); + y = __fmaf_rn(static_cast(index[j] / width), static_cast(counts), y); + photons += counts; + max_photons = max(max_photons, counts); + found++; + } + scratch[l].x = x; + scratch[l].y = y; + scratch[l].photons = photons; + scratch[l].max_photons = max_photons; + } + __syncthreads(); + + // 4) max-pix filter, compacted by another prefix sum so the surviving spots keep their order + const int label_chunk = (nlabel + nthreads - 1) / nthreads; + const int label_lo = min(t * label_chunk, nlabel), label_hi = min(label_lo + label_chunk, nlabel); + uint32_t nkeep = 0; + for (int i = label_lo; i < label_hi; i++) nkeep += (scratch[i].pixel_count <= max_pix) ? 1u : 0u; + shared[t] = nkeep; + __syncthreads(); + for (int d = 1; d < nthreads; d <<= 1) { + const uint32_t v = (t >= d) ? shared[t - d] : 0u; + __syncthreads(); + shared[t] += v; + __syncthreads(); + } + uint32_t pos = shared[t] - nkeep; + for (int i = label_lo; i < label_hi; i++) + if (scratch[i].pixel_count <= max_pix) out[pos++] = scratch[i]; + if (t == nthreads - 1) *nout = shared[t]; +} + +} // namespace + +SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::shared_ptr in_stream) + : stream(std::move(in_stream)), + width(in_width), + nwords((static_cast(in_width) * in_height + 31) / 32), + gpu_res_mask(nwords), + gpu_nstrong(1), + gpu_index(MAX_STRONG), + gpu_value(MAX_STRONG), + gpu_parent(MAX_STRONG), + gpu_root(MAX_STRONG), + gpu_label(MAX_STRONG), + gpu_count(MAX_STRONG), + gpu_spot(MAX_STRONG), + gpu_spot_out(MAX_STRONG), + gpu_nspot(1), + host_nspot(1), + host_spot(SPOT_PREFIX) { + // One block per few hundred words: enough blocks to fill the device, few enough that the serial + // emission inside a block stays short even when a whole detector row lights up. + compact_blocks = static_cast((nwords + 255) / 256); + if (compact_blocks > 4096) compact_blocks = 4096; + if (compact_blocks < 1) compact_blocks = 1; + gpu_block_count = CudaDevicePtr(compact_blocks); + gpu_block_offset = CudaDevicePtr(compact_blocks); + + // Nothing excluded except the padding bits of the last word - the same starting point as + // ImageSpotFinder's own mask, so the two agree even when no resolution mask is ever set. + std::vector mask(nwords, 0); + const size_t npixel = static_cast(in_width) * in_height; + if (npixel % 32 != 0) + mask.back() = ~((1u << (npixel % 32)) - 1u); + cuda_err(cudaMemcpy(gpu_res_mask, mask.data(), nwords * sizeof(uint32_t), cudaMemcpyHostToDevice)); +} + +void SpotExtractorGPU::SetResolutionMask(const std::vector &packed_mask) { + if (packed_mask.size() != nwords) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "SpotExtractorGPU::SetResolutionMask: mask size mismatch"); + cuda_err(cudaMemcpyAsync(gpu_res_mask, packed_mask.data(), nwords * sizeof(uint32_t), + cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); +} + +void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_image, + const SpotFindingSettings &settings, std::vector &spots) { + const int max_pix = static_cast(settings.max_pix_per_spot); + + count_bits<<>>(gpu_strong, gpu_res_mask, gpu_block_count, nwords); + scan_block_counts<<<1, FINISH_THREADS, 0, *stream>>>(gpu_block_count, gpu_block_offset, gpu_nstrong, + compact_blocks); + scatter_bits<<>>(gpu_strong, gpu_res_mask, gpu_block_offset, gpu_image, + gpu_index, gpu_value, nwords, MAX_STRONG); + // Fixed grids reading the strong-pixel count from device memory: the host never learns it, so it + // never has to synchronise in the middle of the frame. + init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, MAX_STRONG); + union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, MAX_STRONG, width); + resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, MAX_STRONG); + finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, + gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, + gpu_nstrong, MAX_STRONG, width, max_pix); + cuda_err(cudaMemcpyAsync(host_nspot, gpu_nspot, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_spot, gpu_spot_out, SPOT_PREFIX * sizeof(SpotExtractorGPUSpot), + cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); // the only synchronisation in the frame + + const uint32_t nspot = *host_nspot.get(); + const SpotExtractorGPUSpot *s = host_spot.get(); + if (nspot > SPOT_PREFIX) { + overflow_spot.resize(nspot); + cuda_err(cudaMemcpyAsync(overflow_spot.data(), gpu_spot_out, nspot * sizeof(SpotExtractorGPUSpot), + cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); + s = overflow_spot.data(); + } + + spots.clear(); + spots.reserve(nspot); + for (uint32_t i = 0; i < nspot; i++) + spots.emplace_back(s[i].x, s[i].y, s[i].pixel_count, s[i].photons, s[i].max_photons); +} diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h new file mode 100644 index 00000000..b3eafa2e --- /dev/null +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -0,0 +1,91 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// Device-side connected-component extraction for the GPU spot finders. +// +// The GPU finders flag strong pixels into a packed bit buffer ON THE DEVICE. Reading spots out of it +// used to mean copying that whole buffer back (2.26 MB per frame at 18 MP) and scanning it bit by bit +// on the host. This does the whole extraction where the data already is, so nothing about the image +// comes back - only the finished spot list, a few hundred entries. +// +// The algorithm is the sparse formulation the ACTS/traccc project settled on for the same problem +// (sparse silicon-detector hits): the strong pixels are compacted into a list that is sorted by flat +// index, each pixel finds its at most FOUR backward 8-neighbours by binary search in that list, and +// the resulting graph is labelled with a lock-free union-find. A dense image-wide labelling +// (Playne-equivalence, BUF/BKE, nppiLabelMarkers, cv::cuda::connectedComponents) would label 18 +// million pixels to find five hundred. +// +// It reproduces the host StrongPixelSet::sparseccl EXACTLY, not just equivalently: +// * both make a component's root its lowest list index, so both find the same roots; +// * labels are handed out by a prefix sum over the roots in ascending order, which is the order the +// host's second scan hands them out in, so the SPOT ORDER is identical; +// * the centroid sums are accumulated per component in ascending list order, in float, term for +// term as DiffractionSpot::AddPixel does them, with the rounding spelled out (see the comment at +// the sum itself - the two compilers do not contract a multiply-add alike). +// tests/SpotExtractorGPUParityTest.cpp holds the two to each other on realistic, occupancy-swept and +// pathological frames, and checks that repeating a frame gives byte-identical output. + +#include +#include +#include + +#include "../../common/DiffractionSpot.h" +#include "../indexing/CUDAMemHelpers.h" +#include "SpotFindingSettings.h" + +// Per-component sums, in exactly the form DiffractionSpot holds them: x and y are sum(col*photons) +// and sum(line*photons), not a centroid. +struct SpotExtractorGPUSpot { + float x; + float y; + int64_t photons; + int64_t max_photons; + int32_t pixel_count; + int32_t padding; +}; + +class SpotExtractorGPU { + std::shared_ptr stream; + const int32_t width; + const size_t nwords; + + // The connected-component search gives up above this many strong pixels (see + // StrongPixelSet::FindComponentsImage), so nothing larger is ever built. + static constexpr uint32_t MAX_STRONG = UINT16_MAX; + // Spots copied back together with their count in one transfer. A frame with more than this many + // surviving spots - far past anything indexable - simply takes a second copy. + static constexpr uint32_t SPOT_PREFIX = 4096; + + int compact_blocks = 0; + + CudaDevicePtr gpu_res_mask; // packed, bit set = pixel excluded + CudaDevicePtr gpu_block_count; + CudaDevicePtr gpu_block_offset; + CudaDevicePtr gpu_nstrong; + CudaDevicePtr gpu_index; // strong pixels, sorted by flat index + CudaDevicePtr gpu_value; + CudaDevicePtr gpu_parent; // union-find parent + CudaDevicePtr gpu_root; + CudaDevicePtr gpu_label; // compact label, indexed by root + CudaDevicePtr gpu_count; // pixels per component + CudaDevicePtr gpu_spot; + CudaDevicePtr gpu_spot_out; + CudaDevicePtr gpu_nspot; + + CudaHostPtr host_nspot; + CudaHostPtr host_spot; // SPOT_PREFIX entries, pinned + std::vector overflow_spot; // only for a frame with more spots than that + +public: + SpotExtractorGPU(int32_t width, int32_t height, std::shared_ptr stream); + + void SetResolutionMask(const std::vector &packed_mask); + + // gpu_strong is the finder's device bit buffer, gpu_image the preprocessed image it was built + // from. Fills spots with every component of at most max-pix pixels, in the same order the host + // extractor would. + void Extract(const uint32_t *gpu_strong, const int32_t *gpu_image, + const SpotFindingSettings &settings, std::vector &spots); +}; diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 3a6ac524..7af8acfb 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -84,14 +84,12 @@ std::vector StrongPixelSet::sparseccl() { } -void StrongPixelSet::FindSpotsImage(const SpotFindingSettings &settings, std::vector &spots) { +void StrongPixelSet::FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots) { // Avoid spot finding, when more than 65536 strong pixel count (will be super slow) if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { for (const auto &spot: sparseccl()) { - if ((spot.PixelCount() <= settings.max_pix_per_spot) - && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { + if (spot.PixelCount() <= settings.max_pix_per_spot) spots.push_back(spot); - } } } } diff --git a/image_analysis/spot_finding/StrongPixelSet.h b/image_analysis/spot_finding/StrongPixelSet.h index be242cec..2e53c28c 100644 --- a/image_analysis/spot_finding/StrongPixelSet.h +++ b/image_analysis/spot_finding/StrongPixelSet.h @@ -35,6 +35,9 @@ public: void AddStrongPixel(uint16_t col, uint16_t line, int32_t photons = 1); void FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, std::vector &spots, uint16_t module_number); - void FindSpotsImage(const SpotFindingSettings &settings, std::vector &spots); + // Every connected component of at most max-pix pixels. min-pix is deliberately NOT applied: it is + // the only spot setting that changes between the passes of the per-image min-pix search, so one + // connected-component search serves all of them (see ImageSpotFinder::ExtractComponents). + void FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots); uint32_t GetStrongPixelCount() const; }; diff --git a/tests/AdaptiveSpotFinderCPUTest.cpp b/tests/AdaptiveSpotFinderCPUTest.cpp index f7819812..321a3776 100644 --- a/tests/AdaptiveSpotFinderCPUTest.cpp +++ b/tests/AdaptiveSpotFinderCPUTest.cpp @@ -67,7 +67,8 @@ TEST_CASE("AdaptiveSpotFinderCPU_RawGeometry", "[AdaptiveSpotFinder]") { std::vector res_mask(x.GetPixelsNum(), false); AdaptiveSpotFinderCPU finder(mapping); - const auto spots = finder.Run(buffer, AdaptiveSettings(), res_mask); + finder.SetResolutionMask(res_mask); + const auto spots = finder.Run(buffer, AdaptiveSettings()); REQUIRE(spots.size() == 1); CHECK(std::lround(spots[0].RawCoord().x) == static_cast(spot_col)); @@ -113,7 +114,8 @@ TEST_CASE("AdaptiveSpotFinderCPU_ThresholdTracksBackground", "[AdaptiveSpotFinde buffer[planted[1]] = scale * (10 + 30); // mean + 28 S: well above it std::vector res_mask(x.GetPixelsNum(), false); AdaptiveSpotFinderCPU finder(mapping); - return finder.Run(buffer, AdaptiveSettings(), res_mask); + finder.SetResolutionMask(res_mask); + return finder.Run(buffer, AdaptiveSettings()); }; const auto plain = run_at_scale(1); diff --git a/tests/AdaptiveSpotFinderGPUTest.cpp b/tests/AdaptiveSpotFinderGPUTest.cpp index 74ace77e..2e6f6f47 100644 --- a/tests/AdaptiveSpotFinderGPUTest.cpp +++ b/tests/AdaptiveSpotFinderGPUTest.cpp @@ -98,8 +98,10 @@ TEST_CASE("AdaptiveSpotFinderGPU_SpotFindingParity", "[AdaptiveSpotFinderGPU]") auto stream = std::make_shared(); AdaptiveSpotFinderGPU gpu(mapping, stream); - const auto cpu_spots = cpu.Run(buffer, settings, res_mask); - const auto gpu_spots = gpu.Run(buffer, settings, res_mask); + cpu.SetResolutionMask(res_mask); + gpu.SetResolutionMask(res_mask); + const auto cpu_spots = cpu.Run(buffer, settings); + const auto gpu_spots = gpu.Run(buffer, settings); INFO("cpu spots=" << cpu_spots.size() << " gpu spots=" << gpu_spots.size()); REQUIRE(cpu_spots.size() > 0); @@ -132,7 +134,7 @@ TEST_CASE("AdaptiveSpotFinderGPU_AzimuthalIntegration", "[AdaptiveSpotFinderGPU] auto stream = std::make_shared(); AdaptiveSpotFinderGPU gpu(mapping, stream); - gpu.Run(buffer, settings, res_mask); + gpu.Run(buffer, settings); AzIntEngineGPU azint(mapping, stream); AzimuthalIntegrationProfile ref_profile(mapping); @@ -181,10 +183,10 @@ TEST_CASE("AdaptiveSpotFinderGPU_RunToRunReproducible", "[AdaptiveSpotFinderGPU] auto stream = std::make_shared(); AdaptiveSpotFinderGPU gpu(mapping, stream); - const auto first = gpu.Run(buffer, settings, res_mask); + const auto first = gpu.Run(buffer, settings); REQUIRE(first.size() > 0); for (int repeat = 0; repeat < 4; repeat++) { - const auto again = gpu.Run(buffer, settings, res_mask); + const auto again = gpu.Run(buffer, settings); REQUIRE(again.size() == first.size()); REQUIRE(SortedCoords(again) == SortedCoords(first)); } @@ -230,9 +232,9 @@ TEST_CASE("AdaptiveSpotFinderGPU_Speed", "[AdaptiveSpotFinderGPU][.benchmark]") }; const double t_azint = bench("GPU azint (standalone)", [&] { azint.Run(buffer, profile); }); - const double t_cpu = bench("CPU adaptive spot finding", [&] { cpu.Run(buffer, settings, res_mask); }); - const double t_classic = bench("GPU classic spot finding (local-box)", [&] { gpu_classic.Run(buffer, settings, res_mask); }); - const double t_fused = bench("GPU adaptive FUSED (azint + spot finding)", [&] { gpu_fused.Run(buffer, settings, res_mask); }); + const double t_cpu = bench("CPU adaptive spot finding", [&] { cpu.Run(buffer, settings); }); + const double t_classic = bench("GPU classic spot finding (local-box)", [&] { gpu_classic.Run(buffer, settings); }); + const double t_fused = bench("GPU adaptive FUSED (azint + spot finding)", [&] { gpu_fused.Run(buffer, settings); }); WARN("standard adaptive path (GPU azint + CPU adaptive) = " << (t_azint + t_cpu) << " ms/frame vs fused GPU = " << t_fused << " ms/frame (speedup " diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1e3959e9..e2ff8066 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -66,6 +66,7 @@ ADD_EXECUTABLE(jfjoch_test ImageSpotFinderCPUTest.cpp ImageSpotFinderGPUTest.cpp AdaptiveSpotFinderGPUTest.cpp + SpotExtractorGPUParityTest.cpp CalcBraggPredictionTest.cpp SpotUtilsTest.cpp LatticeSearchTest.cpp diff --git a/tests/ImageSpotFinderCPUTest.cpp b/tests/ImageSpotFinderCPUTest.cpp index 5d1ff417..fe680a6c 100644 --- a/tests/ImageSpotFinderCPUTest.cpp +++ b/tests/ImageSpotFinderCPUTest.cpp @@ -30,7 +30,8 @@ TEST_CASE("ImageSpotFinderCPU_SignalToNoise") { }; std::vector mask_resolution(width * height, false); - auto spots = s.Run(buffer, settings, mask_resolution); + s.SetResolutionMask(mask_resolution); + auto spots = s.Run(buffer, settings); REQUIRE(spots.size() == 2); REQUIRE(spots[0].RawCoord().y == 25); @@ -64,7 +65,8 @@ TEST_CASE("ImageSpotFinderCPU_SignalToNoise_Resolution") { std::vector mask_resolution(width * height, false); mask_resolution[width * 50 + 50] = true; - auto spots = s.Run(buffer, settings, mask_resolution); + s.SetResolutionMask(mask_resolution); + auto spots = s.Run(buffer, settings); REQUIRE(spots.size() == 1); REQUIRE(spots[0].RawCoord().x == 26); @@ -98,7 +100,8 @@ TEST_CASE("ImageSpotFinderCPU_CountThreshold_Resolution") { std::vector mask_resolution(width * height, false); mask_resolution[width * 50 + 50] = true; - auto spots = s.Run(buffer, settings, mask_resolution); + s.SetResolutionMask(mask_resolution); + auto spots = s.Run(buffer, settings); REQUIRE(spots.size() == 2); REQUIRE(spots[0].RawCoord().y == 25); @@ -132,7 +135,8 @@ TEST_CASE("ImageSpotFinderCPU_CountThreshold_Mask") { }; std::vector mask_resolution(width * height, false); - auto spots = s.Run(buffer, settings, mask_resolution); + s.SetResolutionMask(mask_resolution); + auto spots = s.Run(buffer, settings); REQUIRE(spots.size() == 3); REQUIRE(spots[0].RawCoord().x == 26); @@ -171,7 +175,8 @@ TEST_CASE("ImageSpotFinderCPU_SignalToNoise_Mask") { }; std::vector mask_resolution(width * height, false); - auto spots = s.Run(buffer, settings, mask_resolution); + s.SetResolutionMask(mask_resolution); + auto spots = s.Run(buffer, settings); REQUIRE(spots.size() == 3); REQUIRE(spots[0].RawCoord().x == 26); diff --git a/tests/ImageSpotFinderGPUTest.cpp b/tests/ImageSpotFinderGPUTest.cpp index 41de3efc..84226cca 100644 --- a/tests/ImageSpotFinderGPUTest.cpp +++ b/tests/ImageSpotFinderGPUTest.cpp @@ -18,7 +18,7 @@ static void fill_test_image(ImagePreprocessorBuffer& buffer, size_t width, size_ buffer[width * 75 + 25] = 12; } -// Helper to run GPU and get DiffractionSpot list via StrongPixelSet -> FindSpotsImage +// Helper to run the GPU finder and collect its spot list static std::vector run_gpu_and_collect_spots(ImagePreprocessorBufferGPU &buffer, size_t width, size_t height, const SpotFindingSettings &settings, @@ -26,6 +26,7 @@ static std::vector run_gpu_and_collect_spots(ImagePreprocessorB auto stream = std::make_shared(); ImageSpotFinderGPU gpu(static_cast(width), static_cast(height), stream); REQUIRE(get_gpu_count() > 0); + gpu.SetResolutionMask(res_mask); REQUIRE(cudaMemcpyAsync(buffer.getGPUBuffer(), buffer.getBuffer().data(), @@ -33,7 +34,7 @@ static std::vector run_gpu_and_collect_spots(ImagePreprocessorB cudaMemcpyHostToDevice, *stream) == cudaSuccess); - return gpu.Run(buffer, settings, res_mask); + return gpu.Run(buffer, settings); } // Mirror of ImageSpotFinder_SignalToNoise @@ -59,8 +60,8 @@ TEST_CASE("ImageSpotFinderGPU_SignalToNoise") { .low_resolution_limit = 3.0, }; - // GPU produces strong pixels; FindSpotsImage uses mask/resolution implicit in StrongPixelSet. - // StrongPixelSet doesn't carry resolution/mask by itself, but FindSpotsImage(settings, vec) + // GPU produces strong pixels; the resolution mask is handed to the finder separately and the + // connected-component search then runs on the device (SpotExtractorGPU). The spot-level filter // matches CPU ImageSpotFinder test behavior for these synthetic inputs. auto spots = run_gpu_and_collect_spots(buffer, width, height, settings, res_mask); @@ -91,8 +92,8 @@ TEST_CASE("ImageSpotFinderGPU_CountThreshold") { .low_resolution_limit = 3.0, }; - // GPU produces strong pixels; FindSpotsImage uses mask/resolution implicit in StrongPixelSet. - // StrongPixelSet doesn't carry resolution/mask by itself, but FindSpotsImage(settings, vec) + // GPU produces strong pixels; the resolution mask is handed to the finder separately and the + // connected-component search then runs on the device (SpotExtractorGPU). The spot-level filter // matches CPU ImageSpotFinder test behavior for these synthetic inputs. auto spots = run_gpu_and_collect_spots(buffer, width, height, settings, res_mask); @@ -124,8 +125,8 @@ TEST_CASE("ImageSpotFinderGPU_20M") { .low_resolution_limit = 3.0, }; - // GPU produces strong pixels; FindSpotsImage uses mask/resolution implicit in StrongPixelSet. - // StrongPixelSet doesn't carry resolution/mask by itself, but FindSpotsImage(settings, vec) + // GPU produces strong pixels; the resolution mask is handed to the finder separately and the + // connected-component search then runs on the device (SpotExtractorGPU). The spot-level filter // matches CPU ImageSpotFinder test behavior for these synthetic inputs. auto spots = run_gpu_and_collect_spots(buffer, width, height, settings, res_mask); @@ -187,7 +188,8 @@ TEST_CASE("ImageSpotFinder_CPU_GPU_Parity", "[ImageSpotFinder]") { const std::vector res_mask(width * height, false); ImageSpotFinderCPU cpu(static_cast(width), static_cast(height)); - const auto cpu_spots = cpu.Run(cpu_buffer, settings, res_mask); + cpu.SetResolutionMask(res_mask); + const auto cpu_spots = cpu.Run(cpu_buffer, settings); const auto gpu_spots = run_gpu_and_collect_spots(gpu_buffer, width, height, settings, res_mask); diff --git a/tests/SpotExtractorGPUParityTest.cpp b/tests/SpotExtractorGPUParityTest.cpp new file mode 100644 index 00000000..622d30e4 --- /dev/null +++ b/tests/SpotExtractorGPUParityTest.cpp @@ -0,0 +1,370 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// The GPU spot extractor (SpotExtractorGPU) has to reproduce the host connected-component search +// (StrongPixelSet::sparseccl, driven through ImageSpotFinder::ExtractComponentsHost) EXACTLY - the +// same partition of the strong pixels into components, the same number of spots, the same spot +// ORDER, and the same pixel count, photon sum, maximum and centroid sums for each of them. If the +// two ever drift apart, every downstream result (indexing, integration, merging) moves with them +// and nothing else in the suite would notice. +// +// Both extractors are driven from the SAME strong-pixel bit buffer here, so this compares only the +// extraction; the detection stages have their own parity cases +// (ImageSpotFinder_CPU_GPU_Parity, AdaptiveSpotFinderGPU_SpotFindingParity). +// +// The frames are generated, not recorded: dataset-derived fixtures are not allowed in the tree. The +// blob-size and density distributions were chosen to match what real 18 MP frames produce +// (~1e3 - 1.5e4 strong pixels, spots of a few to a few tens of pixels). + +#include +#include "../common/CUDAWrapper.h" + +#ifdef JFJOCH_USE_CUDA + +#include +#include +#include +#include +#include +#include +#include + +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/ImageSpotFinder.h" +#include "../image_analysis/spot_finding/SpotExtractorGPU.h" + +namespace { + +// A detector large enough for the geometries below to be representative, small enough that a case +// costs milliseconds. The pathological shapes are what set the size, not the spot count. +constexpr int32_t WIDTH = 2048; +constexpr int32_t HEIGHT = 2048; + +SpotFindingSettings ExtractionSettings(int64_t max_pix = 50) { + SpotFindingSettings s{}; + s.min_pix_per_spot = 1; + s.max_pix_per_spot = max_pix; + return s; +} + +// A finder whose detection stage is "use the bit buffer I was handed". It exists so both extractors +// can be given byte-identical input. +class PlantedSpotFinder : public ImageSpotFinder { +public: + PlantedSpotFinder() : ImageSpotFinder(WIDTH, HEIGHT) {} + void Detect(const ImagePreprocessorBuffer &, const SpotFindingSettings &) override {} + void Plant(const std::vector &strong_index) { + std::fill(output_buffer.begin(), output_buffer.end(), 0); + for (uint32_t flat : strong_index) + output_buffer[flat / 32] |= 1u << (flat % 32); + } + const std::vector &BitBuffer() const { return output_buffer; } + size_t Words() const { return OutputSize(); } +}; + +// --------------------------------------------------------------------------------------------- +// Frame geometries +// --------------------------------------------------------------------------------------------- +struct Frame { + std::string name; + std::vector index; // flat pixel index, ascending + std::vector value; +}; + +Frame Finish(std::string name, std::vector index, uint64_t seed) { + std::sort(index.begin(), index.end()); + index.erase(std::unique(index.begin(), index.end()), index.end()); + Frame f{std::move(name), std::move(index), {}}; + // Values spanning the range a preprocessed image really carries, including negatives (a + // pedestal-subtracted pixel can go below zero) - the centroid sums must match on those too. + std::mt19937_64 rng(seed); + f.value.reserve(f.index.size()); + for (size_t i = 0; i < f.index.size(); i++) + f.value.push_back(static_cast(rng() % 40000) - 200); + return f; +} + +// Compact blobs scattered over the detector: what a diffraction pattern looks like to the extractor. +Frame BlobFrame(const std::string &name, int nblob, uint64_t seed) { + std::mt19937_64 rng(seed); + std::vector index; + for (int b = 0; b < nblob; b++) { + const int col = 20 + static_cast(rng() % (WIDTH - 40)); + const int line = 20 + static_cast(rng() % (HEIGHT - 40)); + const int radius = 1 + static_cast(rng() % 3); // spots of ~4 to ~30 pixels + for (int dl = -radius; dl <= radius; dl++) + for (int dc = -radius; dc <= radius; dc++) + if (dl * dl + dc * dc <= radius * radius) + index.push_back(static_cast((line + dl) * WIDTH + (col + dc))); + } + return Finish(name, std::move(index), seed ^ 0x9e37u); +} + +// A lit band of whole detector rows - a hot module, a panel edge, and the limit a diffraction ring +// approaches where it runs tangent to a row. This is the shape the host's two-row sliding window +// turns quadratic on. +Frame BandFrame(const std::string &name, int nrows) { + std::vector index; + for (int line = 500; line < 500 + nrows; line++) + for (int col = 0; col < WIDTH; col++) + index.push_back(static_cast(line * WIDTH + col)); + return Finish(name, std::move(index), 0x51ed3u + nrows); +} + +// Concentric diffraction rings. Unlike a hot column, a ring is SIGNAL and cannot be masked away, and +// at its top and bottom it produces exactly the long horizontal run that the band case models. +Frame RingFrame(const std::string &name, int radius, int half_width, int nring) { + const double cx = WIDTH * 0.5, cy = HEIGHT * 0.5; + std::vector index; + for (int line = 0; line < HEIGHT; line++) { + for (int col = 0; col < WIDTH; col++) { + const double d = std::sqrt((col - cx) * (col - cx) + (line - cy) * (line - cy)); + for (int k = 0; k < nring; k++) { + if (std::fabs(d - (radius + k * 130)) <= half_width) { + index.push_back(static_cast(line * WIDTH + col)); + break; + } + } + } + } + return Finish(name, std::move(index), 0x21f9u + radius * 31 + nring); +} + +Frame SolidBlockFrame(const std::string &name, int side) { + std::vector index; + for (int line = 300; line < 300 + side; line++) + for (int col = 300; col < 300 + side; col++) + index.push_back(static_cast(line * WIDTH + col)); + return Finish(name, std::move(index), 0x7c1u + side); +} + +// Isolated single pixels at the density a very noisy frame reaches: many components, all tiny. +Frame SparseNoiseFrame(const std::string &name, int npixel, uint64_t seed) { + std::mt19937_64 rng(seed); + std::vector index; + for (int i = 0; i < npixel; i++) + index.push_back(static_cast(rng() % (static_cast(WIDTH) * HEIGHT))); + return Finish(name, std::move(index), seed ^ 0x1234u); +} + +// --------------------------------------------------------------------------------------------- +// Harness: run both extractors on one frame +// --------------------------------------------------------------------------------------------- +class ParityRunner { + PlantedSpotFinder finder; + ImagePreprocessorBufferGPU image; + std::shared_ptr stream; + SpotExtractorGPU extractor; + CudaDevicePtr gpu_strong; + +public: + ParityRunner() + : image(static_cast(WIDTH) * HEIGHT), + stream(std::make_shared()), + extractor(WIDTH, HEIGHT, stream), + gpu_strong(finder.Words()) {} + + void SetResolutionMask(const std::vector &mask) { + finder.SetResolutionMask(mask); + std::vector packed((static_cast(WIDTH) * HEIGHT + 31) / 32, 0); + for (size_t i = 0; i < mask.size(); i++) + if (mask[i]) packed[i / 32] |= 1u << (i % 32); + extractor.SetResolutionMask(packed); + } + + void Load(const Frame &f) { + auto &buffer = image.getBuffer(); + std::fill(buffer.begin(), buffer.end(), 0); + for (size_t i = 0; i < f.index.size(); i++) + buffer[f.index[i]] = f.value[i]; + REQUIRE(cudaMemcpy(image.getGPUBuffer(), buffer.data(), buffer.size() * sizeof(int32_t), + cudaMemcpyHostToDevice) == cudaSuccess); + finder.Plant(f.index); + REQUIRE(cudaMemcpy(gpu_strong.get(), finder.BitBuffer().data(), + finder.Words() * sizeof(uint32_t), cudaMemcpyHostToDevice) == cudaSuccess); + REQUIRE(cudaDeviceSynchronize() == cudaSuccess); + } + + std::vector Cpu(const SpotFindingSettings &settings) { + return finder.ExtractComponents(image, settings); + } + std::vector Gpu(const SpotFindingSettings &settings) { + std::vector spots; + extractor.Extract(gpu_strong.get(), image.getGPUBuffer(), settings, spots); + return spots; + } +}; + +// How many representable floats apart two values are. +int64_t UlpDistance(float a, float b) { + int32_t ia, ib; + memcpy(&ia, &a, sizeof(ia)); + memcpy(&ib, &b, sizeof(ib)); + if (ia < 0) ia = INT32_MIN - ia; // map to a monotone ordering across the sign + if (ib < 0) ib = INT32_MIN - ib; + return std::abs(static_cast(ia) - static_cast(ib)); +} + +// Everything that decides which spots exist and what they weigh is compared EXACTLY: the number of +// spots, their order, and each one's pixel count, photon sum and maximum. Those are integers, and a +// difference in any of them is a difference in the partition. +// +// The centroid is a float sum, and its last bit is a property of the BUILD rather than of either +// implementation: gcc contracts DiffractionSpot::AddPixel into an FMA whenever the flags allow it +// (the CI -march=x86-64-v3 does), while a baseline -march, or MSVC with its default /fp:precise, +// cannot. The extractor uses __fmaf_rn, so it is bit-exact against a host that fuses and one ulp +// off one that does not. An ulp bound catches a real divergence - which moves a centroid by pixels, +// not by 1e-4 of one - while staying true whichever way the host was built. +void RequireIdentical(const std::string &what, + const std::vector &cpu, + const std::vector &gpu) { + INFO(what); + REQUIRE(gpu.size() == cpu.size()); + for (size_t i = 0; i < cpu.size(); i++) { + INFO("spot " << i << " of " << cpu.size()); + REQUIRE(gpu[i].PixelCount() == cpu[i].PixelCount()); + REQUIRE(gpu[i].Count() == cpu[i].Count()); + REQUIRE(gpu[i].MaxCount() == cpu[i].MaxCount()); + // RawCoord divides the sums by the photon count, so comparing it compares the sums; a spot + // whose photons sum to zero reports (0,0) on both sides by the same branch. + const Coord c = cpu[i].RawCoord(), g = gpu[i].RawCoord(); +#ifdef __FMA__ + REQUIRE(memcmp(&g.x, &c.x, sizeof(float)) == 0); + REQUIRE(memcmp(&g.y, &c.y, sizeof(float)) == 0); +#else + REQUIRE(UlpDistance(g.x, c.x) <= 2); + REQUIRE(UlpDistance(g.y, c.y) <= 2); +#endif + } +} + +std::vector RealisticFrames() { + std::vector frames; + // Spot counts a real 18 MP frame produces, scaled to this detector: a weak still, a good + // rotation frame, a heavily-diffracting one. + frames.push_back(BlobFrame("blobs_sparse", 60, 11)); + frames.push_back(BlobFrame("blobs_typical", 600, 22)); + frames.push_back(BlobFrame("blobs_dense", 2500, 33)); + return frames; +} + +std::vector PathologicalFrames() { + std::vector frames; + frames.push_back(BandFrame("band_2_rows", 2)); + frames.push_back(BandFrame("band_8_rows", 8)); + frames.push_back(BandFrame("band_31_rows", 31)); // just under the 65534-strong-pixel cap + frames.push_back(SolidBlockFrame("solid_block_120", 120)); + frames.push_back(SolidBlockFrame("solid_block_240", 240)); + frames.push_back(RingFrame("ring_r700_w3", 700, 1, 1)); + frames.push_back(RingFrame("ring_r700_w9", 700, 4, 1)); + frames.push_back(RingFrame("ring_r300_w3_x4", 300, 1, 4)); + frames.push_back(RingFrame("ring_r950_w5", 950, 2, 1)); + frames.push_back(SparseNoiseFrame("noise_20k", 20000, 44)); + return frames; +} + +} // namespace + +TEST_CASE("SpotExtractorGPU_Parity_RealisticFrames", "[SpotExtractorGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping SpotExtractorGPU_Parity_RealisticFrames"); + return; + } + ParityRunner runner; + const auto settings = ExtractionSettings(); + for (const auto &f : RealisticFrames()) { + runner.Load(f); + RequireIdentical(f.name + " (" + std::to_string(f.index.size()) + " strong pixels)", + runner.Cpu(settings), runner.Gpu(settings)); + } +} + +// The occupancy range the extractor has to hold up over, from a nearly empty frame to one pixel +// short of the give-up threshold, plus the first frame past it (where both sides must return +// nothing rather than disagree). +TEST_CASE("SpotExtractorGPU_Parity_OccupancySweep", "[SpotExtractorGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping SpotExtractorGPU_Parity_OccupancySweep"); + return; + } + ParityRunner runner; + const auto settings = ExtractionSettings(); + for (int target : {100, 300, 1000, 3000, 10000, 30000, 65000, 70000}) { + // Isolated pixels give the most components per strong pixel, so this is also the worst case + // for the labelling and the compaction. + Frame f = SparseNoiseFrame("occupancy_" + std::to_string(target), target, 900u + target); + runner.Load(f); + RequireIdentical(f.name + " (" + std::to_string(f.index.size()) + " strong pixels)", + runner.Cpu(settings), runner.Gpu(settings)); + } +} + +TEST_CASE("SpotExtractorGPU_Parity_Pathological", "[SpotExtractorGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping SpotExtractorGPU_Parity_Pathological"); + return; + } + ParityRunner runner; + for (const auto &f : PathologicalFrames()) { + runner.Load(f); + // Both the production max-pix and an effectively unlimited one: the second forces the big + // components to be summed as well, which the first deliberately skips. + for (int64_t max_pix : {int64_t(50), int64_t(1000000)}) { + const auto settings = ExtractionSettings(max_pix); + RequireIdentical(f.name + " (" + std::to_string(f.index.size()) + " strong pixels, max_pix " + + std::to_string(max_pix) + ")", + runner.Cpu(settings), runner.Gpu(settings)); + } + } +} + +TEST_CASE("SpotExtractorGPU_Parity_ResolutionMask", "[SpotExtractorGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping SpotExtractorGPU_Parity_ResolutionMask"); + return; + } + ParityRunner runner; + const auto settings = ExtractionSettings(); + // A resolution band, the shape UpdateMaskResolution really produces: excluding it has to split + // and drop components identically on both sides. + std::vector mask(static_cast(WIDTH) * HEIGHT, false); + const double cx = WIDTH * 0.5, cy = HEIGHT * 0.5; + for (int line = 0; line < HEIGHT; line++) + for (int col = 0; col < WIDTH; col++) { + const double d = std::sqrt((col - cx) * (col - cx) + (line - cy) * (line - cy)); + if (d < 250.0 || d > 900.0) + mask[static_cast(line) * WIDTH + col] = true; + } + runner.SetResolutionMask(mask); + for (const auto &f : RealisticFrames()) { + runner.Load(f); + RequireIdentical(f.name + " with resolution mask", runner.Cpu(settings), runner.Gpu(settings)); + } + runner.Load(RingFrame("ring_r700_w3_masked", 700, 1, 1)); + RequireIdentical("ring with resolution mask", runner.Cpu(settings), runner.Gpu(settings)); +} + +// The extraction is full of atomics and a lock-free union-find, so the same frame has to come back +// byte-identical however the blocks happen to interleave. Anything less and the spot list - and +// every result downstream of it - would wander between runs of the same data. +TEST_CASE("SpotExtractorGPU_Determinism", "[SpotExtractorGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping SpotExtractorGPU_Determinism"); + return; + } + ParityRunner runner; + const auto settings = ExtractionSettings(); + for (const auto &f : {BlobFrame("blobs_typical", 600, 22), RingFrame("ring_r700_w3", 700, 1, 1), + BandFrame("band_8_rows", 8)}) { + runner.Load(f); + const auto reference = runner.Gpu(settings); + for (int repeat = 0; repeat < 100; repeat++) { + const auto again = runner.Gpu(settings); + INFO(f.name << ", repeat " << repeat); + REQUIRE(again.size() == reference.size()); + REQUIRE(memcmp(again.data(), reference.data(), again.size() * sizeof(DiffractionSpot)) == 0); + } + } +} + +#endif