mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

Add fused GPU adaptive spot finder (azint + spot finding in one pass)

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-25 20:10:45 +02:00
co-authored by Claude Opus 4.8
parent a43f9bda13
commit b65691f312
13 changed files with 812 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CPU_DATA_ANALYSIS.md
+16 -2
View File
@@ -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:
image_analysis/MXAnalysisWithoutFPGA.cpp
+36 -6
View File
@@ -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<BraggIntegrationEngineGPU>(in_experiment, stream);
if (experiment.ROI().size() >= 1)
roi = std::make_unique<ROIIntegrationGPU>(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<AdaptiveSpotFinderGPU>(integration, stream);
fused_adaptive = fused.get();
adaptiveSpotFinder = std::move(fused);
}
}
#endif
adaptiveSpotFinder = std::make_unique<AdaptiveSpotFinderCPU>(integration);
if (!adaptiveSpotFinder)
adaptiveSpotFinder = std::make_unique<AdaptiveSpotFinderCPU>(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<float>(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<float>(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<float>(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<float>(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<Reflection> &predicted, size_t npredicted, int64_t image_number) {
image_analysis/MXAnalysisWithoutFPGA.h
+14 -5
View File
@@ -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<AzIntEngine> azint;
std::unique_ptr<ROIIntegration> roi;
std::unique_ptr<ImageSpotFinder> 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<AdaptiveSpotFinderCPU> 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<ImageSpotFinder> adaptiveSpotFinder;
AdaptiveSpotFinderGPU *fused_adaptive = nullptr;
const bool enable_fused_adaptive_gpu;
IndexAndRefine &indexer;
std::unique_ptr<BraggPrediction> prediction;
std::unique_ptr<BraggIntegrationEngine> bragg_engine;
@@ -56,8 +61,12 @@ class MXAnalysisWithoutFPGA {
std::shared_ptr<CudaStream> 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
image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
+7 -78
View File
@@ -7,64 +7,7 @@
#include <unordered_map>
#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<float>(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<float>(k + 1);
}
} // namespace
#include "AdaptiveThreshold.h"
AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping)
: ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
@@ -142,32 +85,18 @@ std::vector<DiffractionSpot> AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB
const double E = std::max(1.0f, settings.false_pixels_per_frame);
double p = E / static_cast<double>(n_total);
p = std::min(std::max(p, 1e-9), 0.1);
const float z = static_cast<float>(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<double>(p), static_cast<double>(z));
return std::max(gauss, poisson);
};
const float z = static_cast<float>(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<float>(g_mean), static_cast<float>(g_sigma));
const float g_thr = adaptive_threshold::RingThreshold(static_cast<float>(g_mean),
static_cast<float>(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)
image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
+319
View File
@@ -0,0 +1,319 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "AdaptiveSpotFinderGPU.h"
#include "AdaptiveThreshold.h"
#include "../../common/JFJochException.h"
namespace {
inline void cuda_err(cudaError_t val) {
if (val != cudaSuccess)
throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val));
}
// One 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<uint32_t *>(&s_sum2[nbins]);
float *s_sum_corr = reinterpret_cast<float *>(&s_count[nbins]);
float *s_sum2_corr = &s_sum_corr[nbins];
for (int i = threadIdx.x; i < nbins; i += blockDim.x) {
s_sum[i] = 0.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<float>(v);
if (clip_k > 0.0f) {
const float lo = mean[b] - clip_k * sigma[b];
const float hi = mean[b] + clip_k * sigma[b];
if (fv < lo || fv > hi) continue;
}
atomicAdd(&s_sum[b], 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<float>(v);
if (clip_k > 0.0f) {
const float lo = mean[b] - clip_k * sigma[b];
const float hi = mean[b] + clip_k * sigma[b];
if (fv < lo || fv > hi) continue;
}
atomicAdd(&sum[b], 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<float>(v) >= thr[b])
s = true;
}
if (s)
atomicOr(&strong[idx / 32], 1u << (idx % 32));
}
}
} // namespace
AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping,
std::shared_ptr<CudaStream> in_stream)
: ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
static_cast<int32_t>(in_mapping.GetHeight())),
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<size_t>(nbins) * (4 * sizeof(float) + sizeof(uint32_t));
shared_clip = static_cast<size_t>(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<<<reduce_blocks, reduce_threads, shared, *stream>>>(
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<<<reduce_blocks, reduce_threads, 0, *stream>>>(
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<<<blocks, threads, 0, *stream>>>(gpu_sum, gpu_sum2, gpu_count, gpu_mean, gpu_sigma, nbins);
}
// Host reproduction of AdaptiveSpotFinderCPU Stage B, from the clipped raw per-ring stats.
void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &settings) {
int64_t n_total = 0;
double g_sum = 0.0, g_sum2 = 0.0;
for (int b = 0; b < nbins; ++b) {
n_total += host_count[b];
g_sum += 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<double>(n_total);
p = std::min(std::max(p, 1e-9), 0.1);
const float z = static_cast<float>(adaptive_threshold::NormalQuantile(1.0 - p));
const double g_mean = g_sum / n_total;
const double g_sigma = std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * g_mean));
const float g_thr = adaptive_threshold::RingThreshold(static_cast<float>(g_mean),
static_cast<float>(g_sigma), p, z);
host_thr.assign(nbins, 0.0f);
for (int b = 0; b < nbins; ++b) {
if (host_count[b] < adaptive_threshold::MIN_RING_PIXELS) {
host_thr[b] = g_thr;
} else {
const double m = static_cast<double>(host_sum[b]) / host_count[b];
const double var = std::max(0.0, static_cast<double>(host_sum2[b]) / host_count[b] - m * m);
host_thr[b] = adaptive_threshold::RingThreshold(static_cast<float>(m),
static_cast<float>(std::sqrt(var)), p, z);
}
}
}
std::vector<DiffractionSpot> AdaptiveSpotFinderGPU::Run(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings,
const std::vector<bool> &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<<<flag_blocks, flag_threads, 0, *stream>>>(
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);
}
image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
+103
View File
@@ -0,0 +1,103 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// 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 <memory>
#include <vector>
#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<CudaStream> 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<uint16_t> gpu_pixel_to_bin;
CudaDevicePtr<float> gpu_corrections;
// Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold.
CudaDevicePtr<float> gpu_sum;
CudaDevicePtr<float> gpu_sum2;
CudaDevicePtr<uint32_t> gpu_count;
CudaDevicePtr<float> gpu_mean; // per-ring raw mean (clip predicate)
CudaDevicePtr<float> gpu_sigma; // per-ring raw sigma (clip predicate)
// Corrected per-ring accumulators (plain first pass only) -> azimuthal-integration profile.
CudaDevicePtr<float> gpu_sum_corr;
CudaDevicePtr<float> gpu_sum2_corr;
// Per-ring detection threshold (host-computed, uploaded) and the strong-pixel bit buffer.
CudaDevicePtr<float> gpu_thr;
CudaDevicePtr<uint32_t> gpu_strong;
// Host mirrors of the small per-ring transfers.
std::vector<float> host_sum; // clipped raw sum } input to the host threshold computation
std::vector<float> host_sum2; // clipped raw sum^2 }
std::vector<uint32_t> host_count; // clipped raw count }
std::vector<float> host_thr; // per-ring threshold (empty -> frame had no valid pixels)
std::vector<float> prof_sum; // plain corrected sum } azimuthal-integration profile
std::vector<float> prof_sum2; // plain corrected sum^2 }
std::vector<uint32_t> prof_count; // plain pixel count }
CudaRegisteredVector<uint32_t> 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<CudaStream> stream);
~AdaptiveSpotFinderGPU() override = default;
AdaptiveSpotFinderGPU(const AdaptiveSpotFinderGPU &) = delete;
AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete;
std::vector<DiffractionSpot> Run(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings,
const std::vector<bool> &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; }
};
image_analysis/spot_finding/AdaptiveThreshold.h
+91
View File
@@ -0,0 +1,91 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// 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 <algorithm>
#include <cmath>
#include <cstdint>
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<float>(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<float>(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<double>(mean), p, static_cast<double>(z));
return std::max(gauss, poisson);
}
} // namespace adaptive_threshold
image_analysis/spot_finding/CMakeLists.txt
+3 -1
View File
@@ -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()
rugnux/Rugnux.cpp
+6 -3
View File
@@ -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_) {
tests/AdaptiveSpotFinderGPUTest.cpp
+197
View File
@@ -0,0 +1,197 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../common/CUDAWrapper.h"
#ifdef JFJOCH_USE_CUDA
#include <algorithm>
#include <chrono>
#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<int32_t>(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<std::pair<int, int>> SortedCoords(const std::vector<DiffractionSpot> &spots) {
std::vector<std::pair<int, int>> out;
out.reserve(spots.size());
for (const auto &s : spots)
out.emplace_back(static_cast<int>(std::lround(s.RawCoord().y)),
static_cast<int>(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<bool> res_mask(x.GetPixelsNum(), false);
const SpotFindingSettings settings = AdaptiveSettings();
AdaptiveSpotFinderCPU cpu(mapping);
auto stream = std::make_shared<CudaStream>();
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<bool> res_mask(x.GetPixelsNum(), false);
const SpotFindingSettings settings = AdaptiveSettings();
auto stream = std::make_shared<CudaStream>();
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<bool> res_mask(x.GetPixelsNum(), false);
const SpotFindingSettings settings = AdaptiveSettings();
auto stream = std::make_shared<CudaStream>();
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<double, std::milli>(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
tests/CMakeLists.txt
+1
View File
@@ -63,6 +63,7 @@ ADD_EXECUTABLE(jfjoch_test
ResolutionShellsTest.cpp
ImageSpotFinderCPUTest.cpp
ImageSpotFinderGPUTest.cpp
AdaptiveSpotFinderGPUTest.cpp
CalcBraggPredictionTest.cpp
SpotUtilsTest.cpp
LatticeSearchTest.cpp
viewer/JFJochImageReadingWorker.cpp
+1 -1
View File
@@ -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<IndexAndRefine>(curr_experiment, indexing.get(), /*retain_outcomes=*/false);
image_analysis = std::make_unique<MXAnalysisWithoutFPGA>(curr_experiment, *azint_mapping, *dataset->pixel_mask,
*index_and_refine.get());
*index_and_refine.get(), /*enable_fused_adaptive_gpu=*/true);
last_profile_.reset();
}
viewer/widgets/JFJochViewerSettingsDock.cpp
+18
View File
@@ -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<float>(spot_.min_pix_per_spot), 0, "px", page);
auto *maxSpots = new NumberLineEdit(10.0f, 100000.0f, static_cast<float>(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);