Jungfraujoch/image_analysis/spot_finding/SpotExtractorGPU.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// 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 integers, term for
//     term as DiffractionSpot::AddPixel does them, so there is no rounding for the two compilers to
//     disagree about.
// 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 <cstdint>
#include <memory>
#include <vector>

#include "../../common/DiffractionSpot.h"
#include "../indexing/CUDAMemHelpers.h"
#include "SpotFindingSettings.h"
#include "../image_preprocessing/PreprocessedPixel.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 {
    int64_t x;
    int64_t y;
    int64_t photons;
    int64_t max_photons;
    int32_t pixel_count;
    int32_t padding;
};

class SpotExtractorGPU {
    std::shared_ptr<CudaStream> 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<uint32_t> gpu_res_mask;      // packed, bit set = pixel excluded
    CudaDevicePtr<uint32_t> gpu_block_count;
    CudaDevicePtr<uint32_t> gpu_block_offset;
    CudaDevicePtr<uint32_t> gpu_nstrong;
    CudaDevicePtr<uint32_t> gpu_index;         // strong pixels, sorted by flat index
    CudaDevicePtr<int32_t>  gpu_value;
    CudaDevicePtr<uint32_t> gpu_parent;        // union-find parent
    CudaDevicePtr<uint32_t> gpu_root;
    CudaDevicePtr<uint32_t> gpu_label;         // compact label, indexed by root
    CudaDevicePtr<int32_t>  gpu_count;         // pixels per component
    CudaDevicePtr<SpotExtractorGPUSpot> gpu_spot;
    CudaDevicePtr<SpotExtractorGPUSpot> gpu_spot_out;
    CudaDevicePtr<uint32_t> gpu_nspot;

    CudaHostPtr<uint32_t> host_nspot;
    CudaHostPtr<SpotExtractorGPUSpot> host_spot;      // SPOT_PREFIX entries, pinned
    std::vector<SpotExtractorGPUSpot> overflow_spot;  // only for a frame with more spots than that

public:
    SpotExtractorGPU(int32_t width, int32_t height, std::shared_ptr<CudaStream> stream);

    void SetResolutionMask(const std::vector<uint32_t> &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, PixelView gpu_image,
                 const SpotFindingSettings &settings, std::vector<DiffractionSpot> &spots);
};