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);