Jungfraujoch/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "AdaptiveSpotFinderGPU.h"
#include "AdaptiveThreshold.h"
#include "../../common/JFJochException.h"

namespace {

inline void cuda_err(cudaError_t val) {
    if (val != cudaSuccess)
        throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val));
}

// One thread's merged contribution to a single ring, pushed to the block's shared accumulators.
__device__ __forceinline__ void flush_ring(unsigned long long *s_sum, unsigned long long *s_sum2,
                                           uint32_t *s_count, float *s_sum_corr, float *s_sum2_corr,
                                           bool accumulate_corrected, int b,
                                           unsigned long long r_sum, unsigned long long r_sum2,
                                           uint32_t r_count, float r_sum_corr, float r_sum2_corr) {
    if (r_count == 0) return;   // also covers the initial "no ring yet"
    atomicAdd(&s_sum[b], r_sum);
    atomicAdd(&s_sum2[b], r_sum2);
    atomicAdd(&s_count[b], r_count);
    if (accumulate_corrected) {
        atomicAdd(&s_sum_corr[b], r_sum_corr);
        atomicAdd(&s_sum2_corr[b], r_sum2_corr);
    }
}

// One ring reduction, staging per-ring sums in shared memory (fast path). Shared layout:
//   [ sum(float) | sum2(float) | count(uint32) | sum_corr(float) | sum2_corr(float) ]  x nbins
// The corrected arrays exist only when accumulate_corrected is true (the plain first pass); on the
// sigma-clip passes only the first three are launched/used.
__global__ void reduce_rings_shared(
    const uint16_t *__restrict__ pixel_to_bin,
    const float *__restrict__ corrections,
    const int32_t *__restrict__ image,
    const float *__restrict__ mean,
    const float *__restrict__ sigma,
    float clip_k,
    bool accumulate_corrected,
    unsigned long long *__restrict__ sum, unsigned long long *__restrict__ sum2, uint32_t *__restrict__ count,
    float *__restrict__ sum_corr, float *__restrict__ sum2_corr,
    size_t npix, int nbins) {

    // The raw accumulators are INTEGERS, not floats. A preprocessed pixel is an exact int32 (the
    // masked and saturated sentinels are skipped below), so v and v*v are exact in 64 bits and
    // integer addition is associative - which makes the ring mean and sigma, and therefore the
    // detection threshold, independent of the order the atomics happen to arrive in. With float
    // accumulators the threshold moved in its last bits between runs, and because detection compares
    // an integer pixel value against it, a threshold crossing an integer flipped every pixel of that
    // value in the ring at once.
    // The CORRECTED sums stay float: they are a pixel value times a float correction, so there is no
    // exact integer form. They feed the reported azimuthal profile, not the detection decision.
    extern __shared__ unsigned long long sh[];
    unsigned long long *s_sum = sh;                    // signed value carried as two's complement
    unsigned long long *s_sum2 = &s_sum[nbins];
    uint32_t *s_count = reinterpret_cast<uint32_t *>(&s_sum2[nbins]);
    float *s_sum_corr = reinterpret_cast<float *>(&s_count[nbins]);
    float *s_sum2_corr = &s_sum_corr[nbins];

    for (int i = threadIdx.x; i < nbins; i += blockDim.x) {
        s_sum[i] = 0;
        s_sum2[i] = 0;
        s_count[i] = 0;
        if (accumulate_corrected) {
            s_sum_corr[i] = 0.0f;
            s_sum2_corr[i] = 0.0f;
        }
    }
    __syncthreads();

    // Four pixels per thread, read as one 16-byte and one 8-byte transaction instead of four of
    // each. The image and both tables come straight from cudaMalloc, which aligns to at least the
    // 16 bytes int4/float4 want, and no caller offsets them. npix need not be a multiple of four -
    // the vector loop stops short and the leftovers are done one at a time below, so nothing is read
    // past the end (which would fold uninitialised device memory into the accumulators).
    const size_t stride = static_cast<size_t>(blockDim.x) * gridDim.x;
    const size_t nquad = npix / 4;

    for (size_t q = blockIdx.x * blockDim.x + threadIdx.x; q < nquad; q += stride) {
        const int4 v4 = reinterpret_cast<const int4 *>(image)[q];
        const ushort4 b4 = reinterpret_cast<const ushort4 *>(pixel_to_bin)[q];
        float4 c4 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        if (accumulate_corrected)
            c4 = reinterpret_cast<const float4 *>(corrections)[q];

        const int32_t vq[4] = {v4.x, v4.y, v4.z, v4.w};
        const uint16_t bq[4] = {b4.x, b4.y, b4.z, b4.w};
        const float cq[4] = {c4.x, c4.y, c4.z, c4.w};

        // A ring is several pixels wide, so consecutive pixels along a row usually fall in the same
        // one. Carry a running total for the ring in registers and push it to shared memory only
        // when the ring changes - one set of atomics for the run instead of one per pixel, which is
        // what this kernel is actually limited by. A pixel dropped by one of the tests below does
        // not end a run; it simply contributes nothing.
        int r_b = -1;
        unsigned long long r_sum = 0, r_sum2 = 0;
        uint32_t r_count = 0;
        float r_sum_corr = 0.0f, r_sum2_corr = 0.0f;

        #pragma unroll
        for (int k = 0; k < 4; k++) {
            const int32_t v = vq[k];
            if (v == INT32_MIN || v == INT32_MAX) continue;
            const int b = bq[k];
            if (b >= nbins) continue;
            const float fv = static_cast<float>(v);
            if (clip_k > 0.0f) {
                const float lo = mean[b] - clip_k * sigma[b];
                const float hi = mean[b] + clip_k * sigma[b];
                if (fv < lo || fv > hi) continue;
            }
            if (b != r_b) {
                flush_ring(s_sum, s_sum2, s_count, s_sum_corr, s_sum2_corr, accumulate_corrected,
                           r_b, r_sum, r_sum2, r_count, r_sum_corr, r_sum2_corr);
                r_b = b;
                r_sum = 0; r_sum2 = 0; r_count = 0;
                r_sum_corr = 0.0f; r_sum2_corr = 0.0f;
            }
            // Unsigned, so a run that overflows wraps exactly as the shared accumulator would have:
            // addition mod 2^64 is associative, which is what keeps the regrouping bit-identical.
            r_sum += static_cast<unsigned long long>(static_cast<long long>(v));
            r_sum2 += static_cast<unsigned long long>(static_cast<long long>(v) * v);
            r_count += 1u;
            if (accumulate_corrected) {
                const float cv = fv * cq[k];
                r_sum_corr += cv;
                r_sum2_corr += cv * cv;
            }
        }
        flush_ring(s_sum, s_sum2, s_count, s_sum_corr, s_sum2_corr, accumulate_corrected,
                   r_b, r_sum, r_sum2, r_count, r_sum_corr, r_sum2_corr);
    }

    // The last npix % 4 pixels, one per thread.
    for (size_t idx = 4 * nquad + blockIdx.x * blockDim.x + threadIdx.x; idx < npix; idx += stride) {
        const int32_t v = image[idx];
        if (v == INT32_MIN || v == INT32_MAX) continue;
        const uint16_t b = pixel_to_bin[idx];
        if (b >= nbins) continue;
        const float fv = static_cast<float>(v);
        if (clip_k > 0.0f) {
            const float lo = mean[b] - clip_k * sigma[b];
            const float hi = mean[b] + clip_k * sigma[b];
            if (fv < lo || fv > hi) continue;
        }
        atomicAdd(&s_sum[b], static_cast<unsigned long long>(static_cast<long long>(v)));
        atomicAdd(&s_sum2[b], static_cast<unsigned long long>(static_cast<long long>(v) * v));
        atomicAdd(&s_count[b], 1u);
        if (accumulate_corrected) {
            const float cv = fv * corrections[idx];
            atomicAdd(&s_sum_corr[b], cv);
            atomicAdd(&s_sum2_corr[b], cv * cv);
        }
    }
    __syncthreads();

    for (int i = threadIdx.x; i < nbins; i += blockDim.x) {
        atomicAdd(&sum[i], s_sum[i]);
        atomicAdd(&sum2[i], s_sum2[i]);
        atomicAdd(&count[i], s_count[i]);
        if (accumulate_corrected) {
            atomicAdd(&sum_corr[i], s_sum_corr[i]);
            atomicAdd(&sum2_corr[i], s_sum2_corr[i]);
        }
    }
}

// Same reduction with direct global atomics (used only when nbins is too large to stage in shared
// memory - a rare, high-bin-count configuration).
__global__ void reduce_rings_global(
    const uint16_t *__restrict__ pixel_to_bin,
    const float *__restrict__ corrections,
    const int32_t *__restrict__ image,
    const float *__restrict__ mean,
    const float *__restrict__ sigma,
    float clip_k,
    bool accumulate_corrected,
    unsigned long long *__restrict__ sum, unsigned long long *__restrict__ sum2, uint32_t *__restrict__ count,
    float *__restrict__ sum_corr, float *__restrict__ sum2_corr,
    size_t npix, int nbins) {

    for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x; idx < npix; idx += blockDim.x * gridDim.x) {
        const int32_t v = image[idx];
        if (v == INT32_MIN || v == INT32_MAX) continue;
        const uint16_t b = pixel_to_bin[idx];
        if (b >= nbins) continue;
        const float fv = static_cast<float>(v);
        if (clip_k > 0.0f) {
            const float lo = mean[b] - clip_k * sigma[b];
            const float hi = mean[b] + clip_k * sigma[b];
            if (fv < lo || fv > hi) continue;
        }
        atomicAdd(&sum[b], static_cast<unsigned long long>(static_cast<long long>(v)));
        atomicAdd(&sum2[b], static_cast<unsigned long long>(static_cast<long long>(v) * v));
        atomicAdd(&count[b], 1u);
        if (accumulate_corrected) {
            const float cv = fv * corrections[idx];
            atomicAdd(&sum_corr[b], cv);
            atomicAdd(&sum2_corr[b], cv * cv);
        }
    }
}

// Per-ring mean/sigma from the current raw accumulators. Rings with no pixels this pass keep their
// previous value (matches the CPU, which leaves ring_mean/ring_sigma untouched when the count is 0).
__global__ void finalize_rings(const unsigned long long *__restrict__ sum,
                               const unsigned long long *__restrict__ sum2,
                               const uint32_t *__restrict__ count,
                               float *__restrict__ mean, float *__restrict__ sigma, int nbins) {
    for (int b = blockIdx.x * blockDim.x + threadIdx.x; b < nbins; b += blockDim.x * gridDim.x) {
        if (count[b] > 0) {
            // In double, then rounded to float for the clip predicate - the same two steps, in the same
            // order and the same types, as AdaptiveSpotFinderCPU::AccumulateRings.
            const double m = static_cast<double>(static_cast<long long>(sum[b])) / count[b];
            const double var = fmax(0.0, static_cast<double>(sum2[b]) / count[b] - m * m);
            mean[b] = static_cast<float>(m);
            sigma[b] = static_cast<float>(sqrt(var));
        }
    }
}

// Flag strong pixels (value >= ring threshold, or saturated) into the packed bit buffer. Strong
// pixels are sparse, so a plain atomicOr per strong pixel is simpler than warp aggregation and the
// contention is negligible. Mirrors AdaptiveSpotFinderCPU Stage C exactly.
__global__ void flag_strong(const int32_t *__restrict__ image,
                            const uint16_t *__restrict__ pixel_to_bin,
                            const float *__restrict__ thr,
                            uint32_t *__restrict__ strong,
                            size_t npix, int nbins) {
    for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x; idx < npix; idx += blockDim.x * gridDim.x) {
        const int32_t v = image[idx];
        bool s = false;
        if (v == INT32_MAX) {
            s = true;
        } else if (v != INT32_MIN) {
            const uint16_t b = pixel_to_bin[idx];
            if (b < nbins && static_cast<float>(v) >= thr[b])
                s = true;
        }
        if (s)
            atomicOr(&strong[idx / 32], 1u << (idx % 32));
    }
}

} // namespace

AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping,
                                             std::shared_ptr<CudaStream> in_stream)
    : ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
                      static_cast<int32_t>(in_mapping.GetHeight()), false),
      mapping(in_mapping),
      stream(in_stream),
      nbins(in_mapping.GetBinNumber()),
      npix(in_mapping.GetPixelToBin().size()),
      gpu_sum(nbins),
      gpu_sum2(nbins),
      gpu_count(nbins),
      gpu_mean(nbins),
      gpu_sigma(nbins),
      gpu_sum_corr(nbins),
      gpu_sum2_corr(nbins),
      gpu_thr(nbins),
      gpu_strong(OutputSize()),
      host_sum(nbins),
      host_sum2(nbins),
      host_count(nbins),
      host_bkg(nbins, NAN),
      prof_sum(nbins),
      prof_sum2(nbins),
      prof_count(nbins),
      host_sum_reg(host_sum),
      host_sum2_reg(host_sum2),
      host_count_reg(host_count),
      prof_sum_reg(prof_sum),
      prof_sum2_reg(prof_sum2),
      prof_count_reg(prof_count),
      extractor(static_cast<int32_t>(in_mapping.GetWidth()),
                static_cast<int32_t>(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
    // memory and SM count can belong to a different card than the one these kernels launch on.
    int device = 0;
    cuda_err(cudaGetDevice(&device));
    cudaDeviceProp prop{};
    cuda_err(cudaGetDeviceProperties(&prop, device));
    reduce_blocks = 8 * prop.multiProcessorCount;
    // flag_strong stays at four: it is bandwidth-shaped rather than atomic-bound, and eight
    // measured no better (181 vs 175 us/launch).
    flag_blocks = 4 * prop.multiProcessorCount;

    shared_plain = static_cast<size_t>(nbins) * (2 * sizeof(unsigned long long) + 2 * sizeof(float) + sizeof(uint32_t));
    shared_clip = static_cast<size_t>(nbins) * (2 * sizeof(unsigned long long) + sizeof(uint32_t));
    use_shared = (shared_plain < prop.sharedMemPerBlock);

    // Launch as many blocks as the device can actually hold resident, and ask it rather than guess.
    // The two passes have different shared-memory footprints, so they do not fit the same number of
    // blocks: the plain pass carries the corrected rings as well, and a grid sized for the clip pass
    // left it running a second wave at a quarter of the occupancy. The kernel is a grid-stride loop,
    // so any grid is correct - but a block that is not resident is a wave, not parallelism. This also
    // tracks nbins, which is data-driven: a fine q spacing shrinks the number of blocks that fit.
    const auto blocks_per_sm = [&](size_t shared) {
        int bpsm = 0;
        cuda_err(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&bpsm, reduce_rings_shared,
                                                              reduce_threads, shared));
        return std::max(1, bpsm) * prop.multiProcessorCount;
    };
    if (use_shared) {
        reduce_blocks_plain = blocks_per_sm(shared_plain);
        reduce_blocks_clip = blocks_per_sm(shared_clip);
    }

    // Both tables are functions of the detector geometry alone, so they are uploaded once per GPU and
    // shared: the azimuthal-integration engine in the same worker reads the very same two arrays.
    gpu_pixel_to_bin = SharedDeviceTable(mapping.GetPixelToBin().data(), npix,
                                         mapping.GetPixelToBin().data(), *stream);
    gpu_corrections = SharedDeviceTable(mapping.Corrections().data(), npix,
                                        mapping.Corrections().data(), *stream);
}

void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, float clip_k,
                                       bool accumulate_corrected) {
    if (use_shared) {
        const size_t shared = accumulate_corrected ? shared_plain : shared_clip;
        const int blocks = accumulate_corrected ? reduce_blocks_plain : reduce_blocks_clip;
        reduce_rings_shared<<<blocks, reduce_threads, shared, *stream>>>(
            gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma,
            clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr,
            npix, nbins);
    } else {
        reduce_rings_global<<<reduce_blocks, reduce_threads, 0, *stream>>>(
            gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma,
            clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr,
            npix, nbins);
    }
}

void AdaptiveSpotFinderGPU::FinalizeStats() {
    const int threads = 128;
    const int blocks = (nbins + threads - 1) / threads;
    finalize_rings<<<blocks, threads, 0, *stream>>>(gpu_sum, gpu_sum2, gpu_count, gpu_mean, gpu_sigma, nbins);
}

// Host reproduction of AdaptiveSpotFinderCPU Stage B, from the clipped raw per-ring stats.
void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &settings) {
    int64_t n_total = 0;
    double g_sum = 0.0, g_sum2 = 0.0;
    for (int b = 0; b < nbins; ++b) {
        n_total += host_count[b];
        g_sum += static_cast<double>(static_cast<int64_t>(host_sum[b]));
        g_sum2 += static_cast<double>(host_sum2[b]);
    }
    if (n_total == 0) {
        host_thr.clear();
        std::fill(host_bkg.begin(), host_bkg.end(), NAN);
        return;
    }

    const double E = std::max(1.0f, settings.false_pixels_per_frame);
    double p = E / static_cast<double>(n_total);
    p = std::min(std::max(p, 1e-9), 0.1);
    const float z = static_cast<float>(adaptive_threshold::NormalQuantile(1.0 - p));

    const double g_mean = g_sum / n_total;
    const double g_sigma = std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * g_mean));
    const float g_thr = adaptive_threshold::RingThreshold(static_cast<float>(g_mean),
                                                          static_cast<float>(g_sigma), p, z);

    host_thr.assign(nbins, 0.0f);
    for (int b = 0; b < nbins; ++b) {
        if (host_count[b] < adaptive_threshold::MIN_RING_PIXELS) {
            host_thr[b] = g_thr;
            host_bkg[b] = NAN;
        } else {
            const double m = static_cast<double>(static_cast<int64_t>(host_sum[b])) / host_count[b];
            const double var = std::max(0.0, static_cast<double>(host_sum2[b]) / host_count[b] - m * m);
            host_thr[b] = adaptive_threshold::RingThreshold(static_cast<float>(m),
                                                           static_cast<float>(std::sqrt(var)), p, z);
            host_bkg[b] = static_cast<float>(m);
        }
    }
}

void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,
                                   const SpotFindingSettings &settings) {
    if (image.size() != npix)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "AdaptiveSpotFinderGPU::Detect: mismatch in pixel size");

    // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) ---
    cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_mean, 0, sizeof(float) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_sigma, 0, sizeof(float) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_sum_corr, 0, sizeof(float) * nbins, *stream));
    cuda_err(cudaMemsetAsync(gpu_sum2_corr, 0, sizeof(float) * nbins, *stream));

    ReducePass(image, 0.0f, true);   // plain pass also fills the corrected profile accumulators
    FinalizeStats();

    // Snapshot the plain corrected profile (and its pixel count) before the raw accumulators are
    // re-zeroed for the sigma-clip passes.
    cuda_err(cudaMemcpyAsync(prof_sum.data(), gpu_sum_corr, sizeof(float) * nbins, cudaMemcpyDeviceToHost, *stream));
    cuda_err(cudaMemcpyAsync(prof_sum2.data(), gpu_sum2_corr, sizeof(float) * nbins, cudaMemcpyDeviceToHost, *stream));
    cuda_err(cudaMemcpyAsync(prof_count.data(), gpu_count, sizeof(uint32_t) * nbins, cudaMemcpyDeviceToHost, *stream));

    for (int pass = 0; pass < 2; ++pass) {
        cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * nbins, *stream));
        cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * nbins, *stream));
        cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * nbins, *stream));
        ReducePass(image, 3.0f, false);
        FinalizeStats();
    }

    // Snapshot the clipped raw stats that drive the threshold.
    cuda_err(cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(unsigned long long) * nbins, cudaMemcpyDeviceToHost, *stream));
    cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(unsigned long long) * nbins, cudaMemcpyDeviceToHost, *stream));
    cuda_err(cudaMemcpyAsync(host_count.data(), gpu_count, sizeof(uint32_t) * nbins, cudaMemcpyDeviceToHost, *stream));
    cuda_err(cudaStreamSynchronize(*stream));

    // --- Stage B: per-ring threshold on the host (shared with the CPU finder) ---
    ComputeThresholds(settings);

    // The profile is a byproduct even when the frame has no valid pixels for detection.
    last_profile.Clear(mapping);
    last_profile.Add(prof_sum, prof_sum2, prof_count);

    if (host_thr.empty()) {
        // 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;
    }

    // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) ---
    cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream));
    cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream));
    flag_strong<<<flag_blocks, flag_threads, 0, *stream>>>(
        image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins);
    // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream,
    // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host.
}

void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector<bool> &mask) {
    ImageSpotFinder::SetResolutionMask(mask);
    extractor.SetResolutionMask(res_mask_bits);
}

const std::vector<DiffractionSpot> &AdaptiveSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image,
                                                                             const SpotFindingSettings &settings) {
    extractor.Extract(gpu_strong, image.getGPUBuffer(), settings, components);
    return components;
}