leonarski_fandClaude Opus 5
0e23fd3ab9
Two independent pieces in the same code path. The background-estimate variance was never propagated. A reflection's background comes from a finite ring of n_b pixels, so subtracting it adds var(B)/n_b per signal pixel - sqrt(1 + n_d/n_b) = 1.109 with the shipped stencil. Both engines omitted it, which is exactly the 1.11-1.19 gap measured between the off-ring scatter and the reported sigma. Three lines each; it affects every dataset, not only iced ones. The radial correction is new and OFF by default (--background-radial). The signal disk and the background ring are concentric, so for any background LINEAR in position <B>_ann == <B>_disk identically and a plane fit buys nothing; the leading error is the CURVATURE of the radial background, which on a sharp ice ring reaches +26 counts on a single reflection. Since every reflection uses the same stencil, that error is a fixed kernel over radial offset - one short dot product per reflection and no extra pixel reads. Validated on empty apertures before any C++: mean |bias| over 9 bands / 3 crystals 4.33 -> 0.79 counts with the scatter unchanged. Three things it cost a battery each to learn, all now in the code: - the radial curve must be accumulated from CLIPPED annulus pixels, inside the clip pass, or it carries neighbour tails and zingers (so it is inert under --integrator boxsum, which has no clip pass); - the GPU version was a 1.8x slowdown from atomicAdd contention on a small radial array - staged in shared memory per block it now costs nothing measurable; - it is battery-NEUTRAL as a default, because the reflections whose bias it fixes are the ones the ice handling already excludes. Hence off by default. CPU/GPU parity extended with two radial sections: 9002 assertions. Also fixes a latent French-Wilson quadrature collapse: j_max = I + 8 sigma on a fixed 400-point grid degenerates to a single cell once sigma >> 50 <I>, giving F = 0.1 sqrt(sigma) with sigmaF -> 0. Harmless today, but any sigma-inflation scheme detonates it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
65 lines
3.2 KiB
C++
65 lines
3.2 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include <cstdint>
|
|
#include <memory>
|
|
#include <vector>
|
|
|
|
#include "BraggIntegrationEngine.h"
|
|
#include "../indexing/CUDAMemHelpers.h"
|
|
|
|
// CUDA engine: reproduces BraggIntegrationEngineCPU up to floating-point precision. Each stage is a
|
|
// kernel with one CUDA block per reflection cooperating over the small window via shared-memory
|
|
// reductions (the natural mapping for thousands of independent, tiny per-spot integrations).
|
|
//
|
|
// Pipeline (profile modes): reset -> mark_mask -> boxsum -> learn_profile -> build_profiles -> fit
|
|
// (the resolution shell is computed inline, so there is no separate shell pass). BoxSum mode stops
|
|
// after boxsum (that pass is the BraggIntegrate2D box integrator and the seed of the profile fit).
|
|
// The preprocessed image already lives on the device (ImagePreprocessorBufferGPU::getGPUBuffer());
|
|
// only the per-frame predicted centres are uploaded.
|
|
class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
|
std::shared_ptr<CudaStream> stream;
|
|
int threads;
|
|
size_t fit_shared_bytes;
|
|
|
|
size_t capacity = 0; // per-reflection device/host arrays hold at least this many reflections
|
|
|
|
// --- per-reflection device arrays (grown by EnsureCapacity) ---
|
|
CudaDevicePtr<float> d_px_x, d_px_y, d_d;
|
|
CudaDevicePtr<int> d_cx, d_cy;
|
|
CudaDevicePtr<float> d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y;
|
|
CudaDevicePtr<float> d_isum; // box-sum raw sum, for the radial correction
|
|
CudaDevicePtr<int> d_ninner, d_rbin;
|
|
CudaDevicePtr<uint8_t> d_ok, d_strong, d_has_obs;
|
|
|
|
// --- radial background curvature correction (see BraggIntegrationEngine) ---
|
|
int n_rad = 0; // radial bins, 0 when the correction is off
|
|
CudaDevicePtr<float> d_rad_sum, d_k_diff;
|
|
CudaDevicePtr<int> d_rad_cnt;
|
|
|
|
// --- fixed-size device arrays ---
|
|
// The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the
|
|
// extraction is Poisson-noise limited, so float reproduces the double CPU path to ~1e-4.
|
|
CudaDevicePtr<uint8_t> d_mask; // per-pixel r2-disk reflection mask
|
|
CudaDevicePtr<float> d_shell_grid, d_global_grid; // learned profile accumulators (N_SHELL*GG, GG)
|
|
CudaDevicePtr<float> d_shell_P, d_global_P; // normalised profiles (empirical mode)
|
|
CudaDevicePtr<float> d_shell_sigma2, d_global_sigma2;
|
|
CudaDevicePtr<int> d_shell_n, d_global_n;
|
|
CudaDevicePtr<unsigned long long> d_invd2; // [min,max] inv-d^2 as monotonic bit patterns
|
|
|
|
// --- host staging (copied back once per frame) ---
|
|
std::vector<float> h_px_x, h_px_y, h_d;
|
|
std::vector<float> h_I, h_sigma, h_bkg, h_obs_x, h_obs_y;
|
|
std::vector<uint8_t> h_ok, h_has_obs;
|
|
|
|
void EnsureCapacity(size_t n);
|
|
|
|
public:
|
|
BraggIntegrationEngineGPU(const DiffractionExperiment &experiment, std::shared_ptr<CudaStream> stream);
|
|
std::vector<Reflection> Run(const ImagePreprocessorBuffer &image,
|
|
const std::vector<Reflection> &predicted, size_t npredicted,
|
|
int64_t image_number) override;
|
|
};
|