From b65691f312c96f2519fee5cc1b00cb26084febfc Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 25 Jul 2026 20:10:45 +0200 Subject: [PATCH] Add fused GPU adaptive spot finder (azint + spot finding in one pass) AdaptiveSpotFinderGPU does the per-resolution-ring reduction once on the GPU and drives both products from it: the azimuthal-integration profile (corrected space) and the self-calibrating adaptive spot-detection threshold (raw counts). This replaces the separate GPU azint pass and the host-side adaptive spot finder that runs on the GPU path today. On a ~4.5 MP detector it does both jobs in ~1 ms/frame versus ~40 ms for the CPU adaptive finder (~42x), with an identical spot list and azimuthal profile. The per-ring threshold math (Poisson tail + read-floored Gaussian, operating point from the false-pixels-per-frame knob) is factored into AdaptiveThreshold.h so the CPU and GPU finders share one source of truth and cannot drift. Wired opt-in via a MXAnalysisWithoutFPGA constructor flag, default on for the rugnux offline path and the interactive viewer, off for the online receiver (so the broker path is unchanged). When on, Analyze() skips the separate azint pass and lifts the profile from the fused engine. The viewer gains an "Adaptive threshold" checkbox that greys out the signal/noise and photon-count sliders (the adaptive finder uses neither). Dedicated tests exercise both products (spot-finding parity vs the CPU finder, azimuthal profile vs a standalone GPU azint) plus a speed benchmark. Validated end-to-end on lysozyme serial stills: fused == CPU-adaptive index rate and merge stats. Docs: new section 3.2 in docs/CPU_DATA_ANALYSIS.md. Co-Authored-By: Claude Opus 4.8 --- docs/CPU_DATA_ANALYSIS.md | 18 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 42 ++- image_analysis/MXAnalysisWithoutFPGA.h | 19 +- .../spot_finding/AdaptiveSpotFinderCPU.cpp | 85 +---- .../spot_finding/AdaptiveSpotFinderGPU.cu | 319 ++++++++++++++++++ .../spot_finding/AdaptiveSpotFinderGPU.h | 103 ++++++ .../spot_finding/AdaptiveThreshold.h | 91 +++++ image_analysis/spot_finding/CMakeLists.txt | 4 +- rugnux/Rugnux.cpp | 9 +- tests/AdaptiveSpotFinderGPUTest.cpp | 197 +++++++++++ tests/CMakeLists.txt | 1 + viewer/JFJochImageReadingWorker.cpp | 2 +- viewer/widgets/JFJochViewerSettingsDock.cpp | 18 + 13 files changed, 812 insertions(+), 96 deletions(-) create mode 100644 image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu create mode 100644 image_analysis/spot_finding/AdaptiveSpotFinderGPU.h create mode 100644 image_analysis/spot_finding/AdaptiveThreshold.h create mode 100644 tests/AdaptiveSpotFinderGPUTest.cpp diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b0c721f7..d9aa7df5 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -185,7 +185,21 @@ Special cases: - saturated pixels can be forced to “strong” (useful for detecting overloaded Bragg spots), - invalid pixels are never strong. -### 3.2 Resolution and ice-ring handling +### 3.2 Adaptive (self-calibrating) detection + +The local-statistics test above still needs a fixed photon/count threshold, and the right value depends on the background level, which varies between datasets — so it has to be tuned per dataset. An optional **adaptive** mode (`--adaptive-spots`) removes that tuning by deriving the threshold from each image's own noise, per resolution ring. + +Pixels are binned into the same resolution rings as the azimuthal integrator (§2). For each ring a robust background is estimated in three passes: one plain pass over all valid pixels, then two $\sigma$-clipping passes that keep only pixels within $\pm 3\sigma$ of the current ring mean (removing the Bragg peaks from the background estimate). This yields a per-ring background mean $\mu_b$ and scatter $\sigma_b$. + +The ring's detection threshold is the larger of two arms, +$ +t_b = \max\!\big(\;\mu_b + z\,\sqrt{\sigma_b^2 + \sigma_\mathrm{read}^2}\;,\;\; k_\mathrm{Poisson}(\mu_b, p)\;\big), +$ +where $k_\mathrm{Poisson}(\mu_b,p)$ is the smallest count whose Poisson$(\mu_b)$ upper tail is $\le p$. The Poisson arm is correct where the background is countable (a bright low-resolution ring gets a high threshold); the Gaussian arm — floored by a detector-level excess-noise constant $\sigma_\mathrm{read}$ — takes over on near-empty high-resolution rings, where the Poisson arm degenerates to "one photon is significant" and would flood. The operating point $p = E/N$ is set from a single portable knob $E$, the expected number of false pixels tolerated per frame (`--spot-false-pixels`, default 100), with $N$ the number of valid pixels. Because $p$ and every $\mu_b,\sigma_b$ come from the image itself, the same $E$ lands a sensible photon threshold on strong and weak datasets alike, with no per-dataset tuning. Rings too sparse to characterise (fewer than ~40 pixels) fall back to a whole-frame background. A pixel is strong when $v_i \ge t_b$ for its ring (saturated pixels are still forced strong); the strong pixels then feed the same CCL stage (§3.4). The signal-to-noise and photon-count criteria of §3.1 are not used in this mode. + +**Fused GPU engine.** The per-ring reduction the adaptive threshold needs is the *same* reduction the azimuthal integrator performs. On the GPU path the two are fused into a single image pass (`AdaptiveSpotFinderGPU`): one reduction accumulates the corrected per-ring sums for the azimuthal profile (§2) *and* the raw per-ring statistics for the threshold, after which a light kernel flags the strong pixels. This replaces both the separate azimuthal-integration pass and the host-side adaptive spot-finding pass with one GPU pass — on a ~4.5 MP detector it runs in ~1 ms/frame versus ~40 ms for the CPU adaptive finder, and produces an identical spot list and azimuthal profile. It is enabled by default in the offline `rugnux` path and the interactive viewer; the online receiver keeps the CPU adaptive finder. + +### 3.3 Resolution and ice-ring handling Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, spots in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6). @@ -193,7 +207,7 @@ A single per-image **ice-ring score** is derived from the azimuthally-integrated A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low. -### 3.3 Connected-component labeling (CCL) +### 3.4 Connected-component labeling (CCL) Strong pixels are grouped into connected components (adjacent strong pixels) using a CCL algorithm. Each component yields a candidate spot with: diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index bd256424..bf021e89 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -13,11 +13,13 @@ #include "azint/AzIntEngineCPU.h" #include "roi/ROIIntegrationCPU.h" #include "spot_finding/ImageSpotFinderCPU.h" +#include "spot_finding/AdaptiveSpotFinderCPU.h" #include "bragg_integration/BraggIntegrationEngineCPU.h" #ifdef JFJOCH_USE_CUDA #include "azint/AzIntEngineGPU.h" #include "roi/ROIIntegrationGPU.h" #include "spot_finding/ImageSpotFinderGPU.h" +#include "spot_finding/AdaptiveSpotFinderGPU.h" #include "image_preprocessing/ImagePreprocessorGPU.h" #include "image_preprocessing/ImagePreprocessorBufferGPU.h" #include "bragg_integration/BraggIntegrationEngineGPU.h" @@ -28,9 +30,11 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_experiment, const AzimuthalIntegrationMapping &in_integration, const PixelMask &in_mask, - IndexAndRefine &in_indexer) + IndexAndRefine &in_indexer, + bool in_enable_fused_adaptive_gpu) : experiment(in_experiment), integration(in_integration), + enable_fused_adaptive_gpu(in_enable_fused_adaptive_gpu), npixels(experiment.GetPixelsNum()), xpixels(experiment.GetXPixelsNum()), indexer(in_indexer), @@ -59,9 +63,17 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp bragg_engine = std::make_unique(in_experiment, stream); if (experiment.ROI().size() >= 1) roi = std::make_unique(experiment, stream); + if (enable_fused_adaptive_gpu) { + // One GPU engine that computes the azimuthal profile and the adaptive spot mask in a single + // image pass. fused_adaptive aliases it so Analyze() can lift the profile out of it. + auto fused = std::make_unique(integration, stream); + fused_adaptive = fused.get(); + adaptiveSpotFinder = std::move(fused); + } } #endif - adaptiveSpotFinder = std::make_unique(integration); + if (!adaptiveSpotFinder) + adaptiveSpotFinder = std::make_unique(integration); } void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, @@ -83,10 +95,18 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const auto preprocessing_end_time = std::chrono::steady_clock::now(); output.preprocessing_time_s = std::chrono::duration(preprocessing_end_time - preprocessing_start_time).count(); - const auto azint_start_time = std::chrono::steady_clock::now(); - azint->Run(*preprocessor_buffer, profile); - const auto azint_end_time = std::chrono::steady_clock::now(); - output.azint_time_s = std::chrono::duration(azint_end_time - azint_start_time).count(); + // The fused GPU engine (rugnux offline, GPU, adaptive detection) produces the azimuthal profile as + // a byproduct of spot finding, so the separate azint pass is skipped in that case and the profile is + // lifted out of the finder below. + const bool fused = enable_fused_adaptive_gpu && spot_finding_settings.enable + && spot_finding_settings.adaptive_threshold && fused_adaptive != nullptr; + + if (!fused) { + const auto azint_start_time = std::chrono::steady_clock::now(); + azint->Run(*preprocessor_buffer, profile); + const auto azint_end_time = std::chrono::steady_clock::now(); + output.azint_time_s = std::chrono::duration(azint_end_time - azint_start_time).count(); + } if (roi) roi->Run(*preprocessor_buffer, output.roi); @@ -106,6 +126,16 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const auto spot_finding_end_time = std::chrono::steady_clock::now(); output.spot_finding_time_s = std::chrono::duration(spot_finding_end_time - spot_finding_start_time).count(); +#ifdef JFJOCH_USE_CUDA + if (fused) { + // Lift the azimuthal profile the fused engine computed in the same pass; its azint cost is + // folded into spot_finding_time_s above. + profile.Clear(integration); + profile += fused_adaptive->GetProfile(); + output.azint_time_s = 0.0f; + } +#endif + if (spot_finding_settings.indexing) indexer.ProcessImage(output, spot_finding_settings, *prediction, [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 1b40f2f4..c79f473a 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -22,6 +22,7 @@ #include "image_preprocessing/ImagePreprocessorBuffer.h" class CudaStream; +class AdaptiveSpotFinderGPU; // MXAnalysisWithoutFPGA is not thread safe - it has to owned by a single thread class MXAnalysisWithoutFPGA { @@ -38,10 +39,14 @@ class MXAnalysisWithoutFPGA { std::unique_ptr azint; std::unique_ptr roi; std::unique_ptr spotFinder; - // Self-calibrating CPU finder, used when spot settings request adaptive detection. Kept alongside - // the default finder because the choice arrives with the per-image settings, not at construction. - // It reads the host preprocessed image (populated on the GPU path too), so it works in either build. - std::unique_ptr adaptiveSpotFinder; + // Self-calibrating finder, used when spot settings request adaptive detection. Kept alongside the + // default finder because the choice arrives with the per-image settings, not at construction. It is + // an AdaptiveSpotFinderCPU by default; on the GPU path, when the fused engine is enabled (rugnux + // offline only), it is instead an AdaptiveSpotFinderGPU that also computes the azimuthal profile, + // aliased through fused_adaptive so Analyze() can take that profile and skip the separate azint pass. + std::unique_ptr adaptiveSpotFinder; + AdaptiveSpotFinderGPU *fused_adaptive = nullptr; + const bool enable_fused_adaptive_gpu; IndexAndRefine &indexer; std::unique_ptr prediction; std::unique_ptr bragg_engine; @@ -56,8 +61,12 @@ class MXAnalysisWithoutFPGA { std::shared_ptr stream; // kept so RebuildROI() can recreate the GPU ROI engine #endif public: + // enable_fused_adaptive_gpu turns on the fused GPU azint+adaptive spot finder (only takes effect on + // the GPU path with adaptive detection). The rugnux offline path and the interactive viewer enable + // it by default; the online receiver leaves it off and keeps the CPU adaptive finder + separate + // azint. It only changes performance - the fused engine reproduces the CPU finder's spots. MXAnalysisWithoutFPGA(const DiffractionExperiment &experiment, const AzimuthalIntegrationMapping &integration, - const PixelMask &mask, IndexAndRefine &indexer); + const PixelMask &mask, IndexAndRefine &indexer, bool enable_fused_adaptive_gpu = false); void Analyze(DataMessage &output, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &spot_finding_settings); // Surgical ROI-only paths used when a full re-analysis is not wanted: rebuild the diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index a551f567..a722a67c 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -7,64 +7,7 @@ #include #include "AdaptiveSpotFinderCPU.h" - -namespace { - -// Number of background pixels a ring needs before its own statistics are trusted; sparser rings -// (detector corners, heavily masked, innermost) fall back to the whole-frame background. -constexpr int64_t MIN_RING_PIXELS = 40; - -// Inverse standard-normal CDF (Acklam's rational approximation, ~1e-9 accuracy). Only called once -// per frame, so accuracy over speed. -double NormalQuantile(double p) { - if (p <= 0.0) return -40.0; - if (p >= 1.0) return 40.0; - static const double a[] = {-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, - 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00}; - static const double b[] = {-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, - 6.680131188771972e+01, -1.328068155288572e+01}; - static const double c[] = {-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, - -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00}; - static const double d[] = {7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, - 3.754408661907416e+00}; - const double plow = 0.02425, phigh = 1.0 - 0.02425; - if (p < plow) { - double q = std::sqrt(-2.0 * std::log(p)); - return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / - ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0); - } else if (p <= phigh) { - double q = p - 0.5, r = q*q; - return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / - (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1.0); - } else { - double q = std::sqrt(-2.0 * std::log(1.0 - p)); - return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / - ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0); - } -} - -// Smallest integer count whose Poisson(mu) upper tail P(X >= k) <= p. This is the correct -// significance floor while the background is countable (it carries the sqrt(mu) shot-noise -// implicitly, so a bright low-resolution ring gets a high threshold). It DEGENERATES at mu -> 0 -// (a single photon on a zero background is "significant"), which is why it is max'd with a -// read-noise-floored Gaussian arm by the caller. Short-circuits to Gaussian for large mu. -float PoissonThreshold(double mu, double p, double z) { - if (mu > 50.0) - return static_cast(mu + z * std::sqrt(mu)); - if (mu < 1e-6) mu = 1e-6; - const double target = 1.0 - p; - double pmf = std::exp(-mu); - double cdf = pmf; - int k = 0; - while (cdf < target && k < 1000) { - ++k; - pmf *= mu / k; - cdf += pmf; - } - return static_cast(k + 1); -} - -} // namespace +#include "AdaptiveThreshold.h" AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping) : ImageSpotFinder(static_cast(in_mapping.GetWidth()), @@ -142,32 +85,18 @@ std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB const double E = std::max(1.0f, settings.false_pixels_per_frame); double p = E / static_cast(n_total); p = std::min(std::max(p, 1e-9), 0.1); - const float z = static_cast(NormalQuantile(1.0 - p)); - - // A ring's threshold is background mean + z sigmas. sigma combines the ring's own (peak-excluded) - // scatter with an excess-noise floor READ: near-zero-background rings scatter MORE than pure - // Poisson (charge sharing / read noise / occasional spurious low counts), so a per-ring sigma - // alone collapses toward zero on empty high-resolution rings and the threshold would flood. READ - // is a detector-level photon-scale constant (the same for every dataset -- it is NOT the - // per-dataset knob), so the operating point still self-calibrates through mean and sigma while - // staying physical where the background vanishes. - const float READ = 1.0f; - auto ring_threshold = [&](float mean, float sigma) { - // Poisson significance (correct where the background is countable) floored by a - // read-noise-aware Gaussian arm (which alone survives mean -> 0, where Poisson degenerates - // to "one photon is significant" and would flood the empty high-resolution rings). - const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ); - const float poisson = PoissonThreshold(mean, static_cast(p), static_cast(z)); - return std::max(gauss, poisson); - }; + const float z = static_cast(adaptive_threshold::NormalQuantile(1.0 - p)); // whole-frame fallback background for rings too sparse to trust on their own 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 = ring_threshold(static_cast(g_mean), static_cast(g_sigma)); + const float g_thr = adaptive_threshold::RingThreshold(static_cast(g_mean), + static_cast(g_sigma), p, z); for (size_t b = 0; b < nbins; ++b) - ring_thr[b] = (ring_cnt[b] < MIN_RING_PIXELS) ? g_thr : ring_threshold(ring_mean[b], ring_sigma[b]); + ring_thr[b] = (ring_cnt[b] < adaptive_threshold::MIN_RING_PIXELS) + ? g_thr + : adaptive_threshold::RingThreshold(ring_mean[b], ring_sigma[b], p, z); // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) --- for (size_t i = 0; i < OutputSize(); ++i) diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu new file mode 100644 index 00000000..4403bde3 --- /dev/null +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -0,0 +1,319 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// 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 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, + float *__restrict__ sum, float *__restrict__ sum2, uint32_t *__restrict__ count, + float *__restrict__ sum_corr, float *__restrict__ sum2_corr, + size_t npix, int nbins) { + + extern __shared__ float sh[]; + float *s_sum = sh; + float *s_sum2 = &s_sum[nbins]; + uint32_t *s_count = reinterpret_cast(&s_sum2[nbins]); + float *s_sum_corr = reinterpret_cast(&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.0f; + s_sum2[i] = 0.0f; + s_count[i] = 0; + if (accumulate_corrected) { + s_sum_corr[i] = 0.0f; + s_sum2_corr[i] = 0.0f; + } + } + __syncthreads(); + + 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(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], fv); + atomicAdd(&s_sum2[b], fv * fv); + 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, + float *__restrict__ sum, float *__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(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], fv); + atomicAdd(&sum2[b], fv * fv); + 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 float *__restrict__ sum, const float *__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) { + const float m = sum[b] / count[b]; + const float var = fmaxf(0.0f, sum2[b] / count[b] - m * m); + mean[b] = m; + sigma[b] = sqrtf(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(v) >= thr[b]) + s = true; + } + if (s) + atomicOr(&strong[idx / 32], 1u << (idx % 32)); + } +} + +} // namespace + +AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping, + std::shared_ptr in_stream) + : ImageSpotFinder(static_cast(in_mapping.GetWidth()), + static_cast(in_mapping.GetHeight())), + mapping(in_mapping), + stream(std::move(in_stream)), + nbins(in_mapping.GetBinNumber()), + npix(in_mapping.GetPixelToBin().size()), + gpu_pixel_to_bin(npix), + gpu_corrections(npix), + 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), + prof_sum(nbins), + prof_sum2(nbins), + prof_count(nbins), + output_buffer_reg(output_buffer), + last_profile(in_mapping) { + + cudaDeviceProp prop{}; + cuda_err(cudaGetDeviceProperties(&prop, 0)); + reduce_blocks = 4 * prop.multiProcessorCount; + flag_blocks = 4 * prop.multiProcessorCount; + + shared_plain = static_cast(nbins) * (4 * sizeof(float) + sizeof(uint32_t)); + shared_clip = static_cast(nbins) * (2 * sizeof(float) + sizeof(uint32_t)); + use_shared = (shared_plain < prop.sharedMemPerBlock); + + cuda_err(cudaMemcpy(gpu_pixel_to_bin, mapping.GetPixelToBin().data(), sizeof(uint16_t) * npix, + cudaMemcpyHostToDevice)); + cuda_err(cudaMemcpy(gpu_corrections, mapping.Corrections().data(), sizeof(float) * npix, + cudaMemcpyHostToDevice)); +} + +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; + reduce_rings_shared<<>>( + gpu_pixel_to_bin, gpu_corrections, 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<<>>( + gpu_pixel_to_bin, gpu_corrections, 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<<>>(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 += host_sum[b]; + g_sum2 += host_sum2[b]; + } + if (n_total == 0) { + host_thr.clear(); + return; + } + + const double E = std::max(1.0f, settings.false_pixels_per_frame); + double p = E / static_cast(n_total); + p = std::min(std::max(p, 1e-9), 0.1); + const float z = static_cast(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(g_mean), + static_cast(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; + } else { + const double m = static_cast(host_sum[b]) / host_count[b]; + const double var = std::max(0.0, static_cast(host_sum2[b]) / host_count[b] - m * m); + host_thr[b] = adaptive_threshold::RingThreshold(static_cast(m), + static_cast(std::sqrt(var)), p, z); + } + } +} + +std::vector AdaptiveSpotFinderGPU::Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask) { + if (image.size() != npix) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "AdaptiveSpotFinderGPU::Run: mismatch in pixel size"); + + // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) --- + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(float) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(float) * 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(float) * nbins, *stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(float) * 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(float) * nbins, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(float) * 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()) + 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<<>>( + image.getGPUBuffer(), gpu_pixel_to_bin, gpu_thr, gpu_strong, npix, nbins); + cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_strong, OutputByteSize(), cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); + + // --- Stage D: connected components + resolution mask + min/max-pix (shared host path) --- + return ExtractSpots(image, settings, res_mask); +} diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h new file mode 100644 index 00000000..e115ca8f --- /dev/null +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -0,0 +1,103 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// GPU adaptive spot finder that FUSES azimuthal integration and spot finding into one image pass. +// +// The CPU adaptive finder (AdaptiveSpotFinderCPU) and the azimuthal integrator both bin every pixel +// into resolution rings and reduce (sum / sum^2 / count). Today azint runs on the GPU while the +// adaptive finder re-does the identical per-ring reduction on the HOST - a wasted second pass over a +// ~10 MP image. This engine does the ring reduction on the GPU and drives BOTH products from it: +// - 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). +// +// 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), +// and the same raw-count detection test. The only differences from the CPU are those inherent to a +// GPU reduction (float per-ring accumulation in atomic order vs the CPU's serial double sums), which +// shift a handful of borderline pixels at most. The corrected sums for the azint profile are +// accumulated in the SAME plain first pass, so one reduction feeds both products. + +#include +#include + +#include "ImageSpotFinder.h" +#include "SpotFindingSettings.h" +#include "../../common/AzimuthalIntegrationProfile.h" +#include "../../common/AzimuthalIntegrationMapping.h" +#include "../indexing/CUDAMemHelpers.h" + +class AdaptiveSpotFinderGPU : public ImageSpotFinder { + const AzimuthalIntegrationMapping &mapping; + std::shared_ptr stream; + + const int nbins; + const size_t npix; + + int reduce_threads = 128; + int reduce_blocks = 0; + int flag_threads = 256; + int flag_blocks = 0; + size_t shared_plain = 0; // per-block shared bytes for the plain pass (raw + corrected rings) + size_t shared_clip = 0; // per-block shared bytes for a sigma-clip pass (raw rings only) + bool use_shared = true; // false -> nbins too large for shared memory, use the global-atomics kernel + + // Static mapping inputs (uploaded once). + CudaDevicePtr gpu_pixel_to_bin; + CudaDevicePtr gpu_corrections; + + // Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold. + CudaDevicePtr gpu_sum; + CudaDevicePtr gpu_sum2; + CudaDevicePtr gpu_count; + CudaDevicePtr gpu_mean; // per-ring raw mean (clip predicate) + CudaDevicePtr gpu_sigma; // per-ring raw sigma (clip predicate) + + // Corrected per-ring accumulators (plain first pass only) -> azimuthal-integration profile. + CudaDevicePtr gpu_sum_corr; + CudaDevicePtr gpu_sum2_corr; + + // Per-ring detection threshold (host-computed, uploaded) and the strong-pixel bit buffer. + CudaDevicePtr gpu_thr; + CudaDevicePtr gpu_strong; + + // Host mirrors of the small per-ring transfers. + std::vector host_sum; // clipped raw sum } input to the host threshold computation + std::vector host_sum2; // clipped raw sum^2 } + std::vector host_count; // clipped raw count } + std::vector host_thr; // per-ring threshold (empty -> frame had no valid pixels) + std::vector prof_sum; // plain corrected sum } azimuthal-integration profile + 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 + + AzimuthalIntegrationProfile last_profile; // filled every Run(), retrievable via GetProfile() + + // One reduction pass over the image into the raw accumulators. clip_k <= 0 -> plain pass (all + // valid pixels); clip_k > 0 -> keep only pixels within clip_k sigma of the current gpu_mean. + // accumulate_corrected additionally fills gpu_sum_corr/gpu_sum2_corr for the profile (plain pass). + void ReducePass(const ImagePreprocessorBuffer &image, float clip_k, bool accumulate_corrected); + // Finalize gpu_mean/gpu_sigma from the current raw accumulators (per ring). + void FinalizeStats(); + // Host: per-ring threshold from the clipped raw stats and the single knob E (false pixels/frame). + void ComputeThresholds(const SpotFindingSettings &settings); + +public: + AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &mapping, std::shared_ptr stream); + ~AdaptiveSpotFinderGPU() override = default; + AdaptiveSpotFinderGPU(const AdaptiveSpotFinderGPU &) = delete; + AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete; + + std::vector Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask) override; + + // The azimuthal profile computed as a byproduct of the last Run() - lets this engine replace the + // separate azint pass in the analysis pipeline. + [[nodiscard]] const AzimuthalIntegrationProfile &GetProfile() const { return last_profile; } +}; diff --git a/image_analysis/spot_finding/AdaptiveThreshold.h b/image_analysis/spot_finding/AdaptiveThreshold.h new file mode 100644 index 00000000..047ee469 --- /dev/null +++ b/image_analysis/spot_finding/AdaptiveThreshold.h @@ -0,0 +1,91 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// Per-resolution-ring detection-threshold math shared by the CPU adaptive spot finder +// (AdaptiveSpotFinderCPU) and its GPU fused counterpart (AdaptiveSpotFinderGPU). Both engines reduce +// every pixel into resolution rings, take a robust per-ring background (mean, sigma), and turn it into +// a strong-pixel threshold with the SAME formula - so keeping that formula in one place is what makes +// the GPU engine reproduce the CPU one. These are plain host functions (the threshold is computed on +// the host in both engines, once per frame, over the small per-ring arrays). + +#include +#include +#include + +namespace adaptive_threshold { + +// Number of background pixels a ring needs before its own statistics are trusted; sparser rings +// (detector corners, heavily masked, innermost) fall back to the whole-frame background. +constexpr int64_t MIN_RING_PIXELS = 40; + +// Detector-level excess-noise floor (photons). Near-zero-background rings scatter MORE than pure +// Poisson (charge sharing / read noise / occasional spurious low counts), so a per-ring sigma alone +// collapses toward zero on empty high-resolution rings and the threshold would flood. READ is a +// photon-scale constant (the same for every dataset -- NOT the per-dataset knob), so the operating +// point still self-calibrates through mean and sigma while staying physical where the background +// vanishes. +constexpr float READ = 1.0f; + +// Inverse standard-normal CDF (Acklam's rational approximation, ~1e-9 accuracy). Only called once +// per frame, so accuracy over speed. +inline double NormalQuantile(double p) { + if (p <= 0.0) return -40.0; + if (p >= 1.0) return 40.0; + static const double a[] = {-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, + 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00}; + static const double b[] = {-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, + 6.680131188771972e+01, -1.328068155288572e+01}; + static const double c[] = {-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, + -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00}; + static const double d[] = {7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, + 3.754408661907416e+00}; + const double plow = 0.02425, phigh = 1.0 - 0.02425; + if (p < plow) { + double q = std::sqrt(-2.0 * std::log(p)); + return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / + ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0); + } else if (p <= phigh) { + double q = p - 0.5, r = q*q; + return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / + (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1.0); + } else { + double q = std::sqrt(-2.0 * std::log(1.0 - p)); + return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / + ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0); + } +} + +// Smallest integer count whose Poisson(mu) upper tail P(X >= k) <= p. This is the correct +// significance floor while the background is countable (it carries the sqrt(mu) shot-noise +// implicitly, so a bright low-resolution ring gets a high threshold). It DEGENERATES at mu -> 0 +// (a single photon on a zero background is "significant"), which is why it is max'd with a +// read-noise-floored Gaussian arm by the caller. Short-circuits to Gaussian for large mu. +inline float PoissonThreshold(double mu, double p, double z) { + if (mu > 50.0) + return static_cast(mu + z * std::sqrt(mu)); + if (mu < 1e-6) mu = 1e-6; + const double target = 1.0 - p; + double pmf = std::exp(-mu); + double cdf = pmf; + int k = 0; + while (cdf < target && k < 1000) { + ++k; + pmf *= mu / k; + cdf += pmf; + } + return static_cast(k + 1); +} + +// A ring's threshold is background mean + z sigmas, computed two ways and max'd: Poisson significance +// (correct where the background is countable) floored by a read-noise-aware Gaussian arm (which alone +// survives mean -> 0, where Poisson degenerates to "one photon is significant" and would flood the +// empty high-resolution rings). p, z are the frame-wide operating point (p = E / N_pixels). +inline float RingThreshold(float mean, float sigma, double p, float z) { + const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ); + const float poisson = PoissonThreshold(static_cast(mean), p, static_cast(z)); + return std::max(gauss, poisson); +} + +} // namespace adaptive_threshold diff --git a/image_analysis/spot_finding/CMakeLists.txt b/image_analysis/spot_finding/CMakeLists.txt index 0a94ca7d..59c39afb 100644 --- a/image_analysis/spot_finding/CMakeLists.txt +++ b/image_analysis/spot_finding/CMakeLists.txt @@ -16,5 +16,7 @@ ADD_LIBRARY(JFJochSpotFinding STATIC TARGET_LINK_LIBRARIES(JFJochSpotFinding JFJochCommon) IF (JFJOCH_CUDA_AVAILABLE) - TARGET_SOURCES(JFJochSpotFinding PRIVATE ImageSpotFinderGPU.cu ImageSpotFinderGPU.h) + TARGET_SOURCES(JFJochSpotFinding PRIVATE + ImageSpotFinderGPU.cu ImageSpotFinderGPU.h + AdaptiveSpotFinderGPU.cu AdaptiveSpotFinderGPU.h) ENDIF() \ No newline at end of file diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 9dacc5d9..b075daaf 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -122,7 +122,8 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ auto worker = [&]() { pin_gpu(); // round-robin per worker thread; must precede engine construction - MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer); + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer, + /*enable_fused_adaptive_gpu=*/true); AzimuthalIntegrationProfile profile(mapping); while (!cancelled_) { @@ -459,7 +460,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b DataMessage m{}; m.number = ordinal; m.original_number = image_idx; - MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer); + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer, + /*enable_fused_adaptive_gpu=*/true); AzimuthalIntegrationProfile profile(mapping); auto first_pass = config_.spot_finding; first_pass.indexing = false; @@ -729,7 +731,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto full_worker = [&]() { pin_gpu(); // round-robin per worker thread; must precede engine construction - MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer); + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer, + /*enable_fused_adaptive_gpu=*/true); AzimuthalIntegrationProfile profile(mapping); while (!cancelled_) { diff --git a/tests/AdaptiveSpotFinderGPUTest.cpp b/tests/AdaptiveSpotFinderGPUTest.cpp new file mode 100644 index 00000000..2b86e994 --- /dev/null +++ b/tests/AdaptiveSpotFinderGPUTest.cpp @@ -0,0 +1,197 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "../common/CUDAWrapper.h" + +#ifdef JFJOCH_USE_CUDA + +#include +#include + +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../image_analysis/azint/AzIntEngineGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderCPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#include "../image_analysis/spot_finding/ImageSpotFinderGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" + +namespace { + +// Build a realistic full-detector azimuthal-integration mapping (JF4M, ~4.5 MP) whose q-range spans +// most of the detector, so the timing runs over a representative pixel count. +DiffractionExperiment MakeExperiment() { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0); + return x; +} + +// Deterministic image: a low, slightly rippled background (well below any adaptive threshold) plus a +// grid of bright multi-pixel blobs that both finders must recover identically. +void FillTestImage(ImagePreprocessorBuffer &buffer, const DiffractionExperiment &x) { + const size_t w = x.GetXPixelsNum(); + const size_t h = x.GetYPixelsNum(); + for (size_t i = 0; i < w * h; i++) + buffer[i] = 8 + static_cast(i % 5); // background 8..12 (mean 10) + + // Bright 3x3 blobs on a coarse grid, kept clear of the edges and the beam centre. + for (size_t row = 300; row < h - 300; row += 450) { + for (size_t col = 300; col < w - 300; col += 450) { + for (int dr = -1; dr <= 1; dr++) + for (int dc = -1; dc <= 1; dc++) + buffer[(row + dr) * w + (col + dc)] = 200; + } + } +} + +SpotFindingSettings AdaptiveSettings() { + SpotFindingSettings s{}; + s.adaptive_threshold = true; + s.false_pixels_per_frame = 100.0f; + s.min_pix_per_spot = 1; + s.max_pix_per_spot = 50; + s.high_resolution_limit = 0.0f; // no resolution gate for the parity test + s.low_resolution_limit = 1.0e6f; + s.high_res_gap_Q_recipA = std::nullopt; + return s; +} + +std::vector> SortedCoords(const std::vector &spots) { + std::vector> out; + out.reserve(spots.size()); + for (const auto &s : spots) + out.emplace_back(static_cast(std::lround(s.RawCoord().y)), + static_cast(std::lround(s.RawCoord().x))); + std::sort(out.begin(), out.end()); + return out; +} + +} // namespace + +// Spot-finding functionality: the fused GPU engine must reproduce the reference CPU adaptive finder's +// spot list (the two share AdaptiveThreshold.h and the host connected-component extractor; the only +// difference is the GPU's float atomic ring reduction, which is exact for a realistic background). +TEST_CASE("AdaptiveSpotFinderGPU_SpotFindingParity", "[AdaptiveSpotFinderGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_SpotFindingParity"); + return; + } + + DiffractionExperiment x = MakeExperiment(); + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + FillTestImage(buffer, x); + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), + x.GetPixelsNum() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); + + std::vector res_mask(x.GetPixelsNum(), false); + const SpotFindingSettings settings = AdaptiveSettings(); + + AdaptiveSpotFinderCPU cpu(mapping); + 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); + + INFO("cpu spots=" << cpu_spots.size() << " gpu spots=" << gpu_spots.size()); + REQUIRE(cpu_spots.size() > 0); + REQUIRE(cpu_spots.size() == gpu_spots.size()); + CHECK(SortedCoords(cpu_spots) == SortedCoords(gpu_spots)); +} + +// Azimuthal-integration functionality: the profile the fused engine computes as a byproduct of the +// same pass must match a standalone GPU azimuthal integrator over the same image. +TEST_CASE("AdaptiveSpotFinderGPU_AzimuthalIntegration", "[AdaptiveSpotFinderGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_AzimuthalIntegration"); + return; + } + + DiffractionExperiment x = MakeExperiment(); + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + FillTestImage(buffer, x); + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), + x.GetPixelsNum() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); + + std::vector res_mask(x.GetPixelsNum(), false); + const SpotFindingSettings settings = AdaptiveSettings(); + + auto stream = std::make_shared(); + AdaptiveSpotFinderGPU gpu(mapping, stream); + gpu.Run(buffer, settings, res_mask); + + AzIntEngineGPU azint(mapping, stream); + AzimuthalIntegrationProfile ref_profile(mapping); + azint.Run(buffer, ref_profile); + + const auto ref = ref_profile.GetResult(); + const auto got = gpu.GetProfile().GetResult(); + const auto ref_count = ref_profile.GetPixelCount(); + const auto got_count = gpu.GetProfile().GetPixelCount(); + REQUIRE(ref.size() == got.size()); + REQUIRE(ref_count == got_count); // identical per-ring pixel counts (same valid-pixel binning) + for (size_t b = 0; b < ref.size(); b++) { + if (std::isnan(ref[b])) { + CHECK(std::isnan(got[b])); + } else { + CHECK(got[b] == Catch::Approx(ref[b]).epsilon(0.01).margin(0.02)); + } + } +} + +TEST_CASE("AdaptiveSpotFinderGPU_Speed", "[AdaptiveSpotFinderGPU][.benchmark]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_Speed"); + return; + } + + DiffractionExperiment x = MakeExperiment(); + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + FillTestImage(buffer, x); + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), + x.GetPixelsNum() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); + + std::vector res_mask(x.GetPixelsNum(), false); + const SpotFindingSettings settings = AdaptiveSettings(); + + auto stream = std::make_shared(); + AdaptiveSpotFinderCPU cpu(mapping); + AdaptiveSpotFinderGPU gpu_fused(mapping, stream); + ImageSpotFinderGPU gpu_classic(x.GetXPixelsNum(), x.GetYPixelsNum(), stream); + AzIntEngineGPU azint(mapping, stream); + AzimuthalIntegrationProfile profile(mapping); + + const int warmup = 5, iters = 40; + auto bench = [&](const char *name, auto &&fn) { + for (int i = 0; i < warmup; i++) fn(); + const auto t0 = std::chrono::steady_clock::now(); + for (int i = 0; i < iters; i++) fn(); + const auto t1 = std::chrono::steady_clock::now(); + const double ms = std::chrono::duration(t1 - t0).count() / iters; + WARN(name << ": " << ms << " ms/frame"); + return ms; + }; + + 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); }); + + WARN("standard adaptive path (GPU azint + CPU adaptive) = " << (t_azint + t_cpu) + << " ms/frame vs fused GPU = " << t_fused << " ms/frame (speedup " + << (t_azint + t_cpu) / t_fused << "x)"); + WARN("fused GPU vs GPU classic finder alone (no azint): " << t_fused << " vs " << t_classic << " ms/frame"); +} + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b10a228d..0dade2aa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -63,6 +63,7 @@ ADD_EXECUTABLE(jfjoch_test ResolutionShellsTest.cpp ImageSpotFinderCPUTest.cpp ImageSpotFinderGPUTest.cpp + AdaptiveSpotFinderGPUTest.cpp CalcBraggPredictionTest.cpp SpotUtilsTest.cpp LatticeSearchTest.cpp diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 8c503241..88471c1d 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -368,7 +368,7 @@ void JFJochImageReadingWorker::UpdateAzint_i(const JFJochReaderDataset *dataset) // never scales the accumulated run, so don't retain the whole-run integration_outcome vector. index_and_refine = std::make_unique(curr_experiment, indexing.get(), /*retain_outcomes=*/false); image_analysis = std::make_unique(curr_experiment, *azint_mapping, *dataset->pixel_mask, - *index_and_refine.get()); + *index_and_refine.get(), /*enable_fused_adaptive_gpu=*/true); last_profile_.reset(); } diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index ce86c51a..82379ae0 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -224,10 +224,16 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { snr->setValue(spot_.signal_to_noise_threshold); auto *count = new SliderPlusBox(0.0, 100.0, 1.0, 0, page); count->setValue(std::lround(spot_.photon_count_threshold)); + auto *adaptive = new QCheckBox("Adaptive threshold (per-ring, auto)", page); + adaptive->setChecked(spot_.adaptive_threshold); + adaptive->setToolTip("Self-calibrating detection: the strong-pixel threshold is derived from each " + "image's own per-resolution-ring noise. The signal/noise and photon-count " + "settings are not used while this is on."); auto *highResSpot = new SliderPlusBox(0.5, 5.0, 0.1, 1, page); highResSpot->setValue(spot_.high_resolution_limit); auto *minPix = new NumberLineEdit(1.0f, 50.0f, static_cast(spot_.min_pix_per_spot), 0, "px", page); auto *maxSpots = new NumberLineEdit(10.0f, 100000.0f, static_cast(max_spots_), 0, "", page); + spot->addRow("", adaptive); spot->addRow("Signal/noise", snr); spot->addRow("Photon count", count); spot->addRow("High resolution [Å]", highResSpot); @@ -246,6 +252,18 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { spot_.min_pix_per_spot = std::llround(minPix->value()); EmitSpotFinding(); }); connect(maxSpots, &NumberLineEdit::newValue, this, [this, maxSpots] { max_spots_ = std::llround(maxSpots->value()); EmitSpotFinding(); }); + // The adaptive finder sets its own threshold from each image's noise, so the signal/noise and + // photon-count sliders do nothing while it is on - grey them out to make that clear. + auto syncAdaptiveEnabled = [snr, count](bool on) { + snr->setEnabled(!on); + count->setEnabled(!on); + }; + syncAdaptiveEnabled(spot_.adaptive_threshold); + connect(adaptive, &QCheckBox::toggled, this, [this, syncAdaptiveEnabled](bool on) { + spot_.adaptive_threshold = on; + syncAdaptiveEnabled(on); + EmitSpotFinding(); + }); // --- Indexing --- auto *idxSection = new CollapsibleSection("Indexing", page);