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

Bragg integration: drop the 2% sigma floor and carry the background variance

Browse Source 
Two changes to the same variance chain; they are in one commit because the second
exists to remove an assumption the first was breaking, and separating them leaves a
tree that is correct only by luck.

The reported sigma was floored at 2% of the intensity, a per-partial I/sigma cap of
50. It applied only to the box-sum seed, never to the profile fit, so the shipped
default was unaffected - but the combine back-derives each partial's non-signal
variance as sigma^2 - I, and a floored sigma makes that quantity mean nothing. It
then read corr^2 * (0.0004 I^2 - I), which is not a background variance. Measured on
--integrator boxsum: the reported sigma understated the true scatter by up to 16x at
I ~ 21000 counts per partial, and pooled_I amplified a 1 ct/px background drift into
an 11.5% intensity error on the strongest reflections.

What the floor stood in for - that at high intensity the error is systematic rather
than counting - is already carried downstream, twice: the fitted b in
v = a*sigma^2 + (b*I)^2, measured from the data rather than assumed, and
SigmaWithSystematicFloor on the merged sigma. The floor was that idea applied one
level too early with a hardcoded b of 0.02. It arrived without a test or a setter and
was unreachable from the CLI, the API and the config.

The merge now takes the non-signal variance the integrator actually measured instead
of inverting sigma^2 = I + N. That identity is exact for a box sum once the floor is
gone and was never exact for a profile fit, whose sigma^2 = 1/den + (wsum/den)^2 *
bkg_var is formed against a fitted intensity. The value is carried through
BraggFitResult, Reflection and Obs, both engines, both merges, and the process-file
round trip; files written before this change are read with the term absent, which is
what they had.

Battery, 37 crystals, paired: space groups unchanged, reflection sets unchanged,
median delta zero on R_meas and CC1/2. --integrator boxsum on the reference crystal
goes ISa 8.9 -> 20.2 with a 0.947 -> 1.032.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-08-09 19:10:57 +02:00
co-authored by Claude Opus 5
parent d565b66916
commit f60768d49c
14 changed files with 89 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/BraggIntegrationEngine.cpp
+1 -1
View File
@@ -53,7 +53,6 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe
r2_sq = r2 * r2;
r3 = settings.GetR3();
r3_sq = r3 * r3;
min_sigma_ratio = settings.GetMinimumSigmaInRegardsToI();
R = static_cast<int>(std::ceil(r2));
G = 2 * R + 1;
GG = G * G;
@@ -131,6 +130,7 @@ std::vector<Reflection> BraggIntegrationEngine::Finalize(const std::vector<Refle
refl.I = fr.I;
refl.sigma = fr.sigma;
refl.bkg = fr.bkg;
refl.var_bkg = fr.var_bkg;
if (fr.has_observed) {
refl.observed_x = fr.observed_x;
refl.observed_y = fr.observed_y;
image_analysis/bragg_integration/BraggIntegrationEngine.h
+5 -1
View File
@@ -71,6 +71,11 @@ struct BraggFitResult {
float bkg = 0.0f;
float observed_x = 0.0f; // intensity-weighted centroid (BoxSum mode only)
float observed_y = 0.0f;
// Variance of everything in `sigma` that is NOT the reflection's own Poisson signal (the
// background and the error of its estimate). The merge rebuilds each partial's variance at the
// pooled intensity and needs this term; back-deriving it as sigma^2 - I only works while that
// identity holds exactly, which it does not for a profile fit.
float var_bkg = 0.0f;
bool ok = false;
bool has_observed = false;
};
@@ -86,7 +91,6 @@ protected:
float r1_sq;
float r2, r2_sq;
float r3, r3_sq;
float min_sigma_ratio;
int R, G, GG; // profile-grid half-size, edge (2R+1) and area (G*G)
bool broadband; // a set bandwidth (stills) vs monochromatic (rotation)
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+11 -3
View File
@@ -86,6 +86,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
struct Rough {
double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0;
double bkg_var = 0.0; // variance of the background ESTIMATE itself, bkg / n_bkg
double var_bkg = 0.0; // total non-signal variance carried to the merge
int64_t I_sum = 0; // kept so I can be rebuilt after the radial background correction
int n_inner = 0;
int r_bin = 0; // rounded distance from the beam centre, indexes the radial curve
@@ -194,10 +195,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
// default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
out.bkg_var = out.bkg / n_bkg_used;
out.var_bkg = static_cast<double>(n_inner) * out.bkg
+ static_cast<double>(n_inner) * n_inner * out.bkg_var;
out.I_sum = I_sum;
out.n_inner = static_cast<int>(n_inner);
const double var_bkg_term = static_cast<double>(n_inner) * n_inner * out.bkg_var;
out.sigma = std::max(1.0, out.I * min_sigma_ratio);
out.sigma = 1.0;
if (I_sum > 0) {
out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(I_sum) + var_bkg_term));
out.obs_x = static_cast<double>(I_sum_x) / static_cast<double>(I_sum);
@@ -244,7 +247,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const auto &rh = rough[i];
if (!rh.ok) continue;
results[i] = {static_cast<float>(rh.I), static_cast<float>(rh.sigma), static_cast<float>(rh.bkg),
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y), true, rh.has_obs};
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
static_cast<float>(rh.var_bkg), true, rh.has_obs};
}
return Finalize(predicted, npredicted, results, image_number);
}
@@ -409,15 +413,19 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's
// own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing.
double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var);
double var_bkg = std::max(0.0, 1.0 / den - std::max(0.0, I)
+ (wsum / den) * (wsum / den) * rh.bkg_var);
if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) {
I = rh.I;
sigma = rh.sigma;
var_bkg = rh.var_bkg;
}
// Carry the Pass-A box-sum intensity-weighted centroid (observed spot position) through the
// profile path too - post-refinement uses it as the observed position (beam-centre / distance).
results[i] = {static_cast<float>(I), static_cast<float>(sigma),
static_cast<float>(rh.bkg),
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y), true, rh.has_obs};
static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
static_cast<float>(var_bkg), true, rh.has_obs};
}
return Finalize(predicted, npredicted, results, image_number);
image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
+15 -10
View File
@@ -16,7 +16,6 @@ inline void cuda_err(cudaError_t val) {
struct BraggGpuParams {
int W, H;
float r1_sq, r2, r2_sq, r3, r3_sq;
float min_sigma_ratio;
int R, G, GG;
float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip)
int empirical; // ProfileEmpirical vs ProfileGaussian
@@ -56,7 +55,7 @@ __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, B
__global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
const int32_t *img, const uint8_t *mask, BraggGpuParams p, int n,
int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o,
float *bkgvar_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o,
float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o,
uint8_t *hasobs_o, unsigned long long *invd2mm,
float *isum_o, int *ninner_o, int *rbin_o,
float *rad_sum, int *rad_cnt, int n_rad) {
@@ -235,7 +234,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
// var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree.
const double bkg_var = bkg / (double) n_bkg_used;
const double var_bkg_term = (double) s_ninner * (double) s_ninner * bkg_var;
double sigma = fmax(1.0, I * (double) p.min_sigma_ratio);
double sigma = 1.0;
uint8_t hasobs = 0; double ox = 0.0, oy = 0.0;
if (Isum > 0) {
sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term));
@@ -247,6 +246,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
cy_o[i] = (int) lroundf(cy);
I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg;
bkgvar_o[i] = (float) bkg_var;
varbkg_o[i] = (float) ((double) s_ninner * bkg + var_bkg_term);
isum_o[i] = (float) Isum;
ninner_o[i] = s_ninner;
rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0);
@@ -383,10 +383,10 @@ __global__ void radial_correct(const float *rad_sum, const int *rad_cnt, int n_r
__global__ void fit(const int32_t *img, const float *px_x, const float *px_y,
const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm,
const float *I_seed, const float *sigma_seed, const float *bkg_a,
const float *bkgvar_a, const uint8_t *ok_a,
const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a,
const float *shell_P, const float *global_P,
const float *shell_sigma2, const float *global_sigma2, const int *shell_n,
float *I_o, float *sigma_o, uint8_t *ok_o, BraggGpuParams p, int n) {
float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n) {
const int i = blockIdx.x;
if (i >= n) return;
extern __shared__ float Pbuf[];
@@ -471,13 +471,15 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y,
// estimate's own error (see the CPU engine).
const float wr = s_wsum / s_den;
float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]);
float var_bkg = fmaxf(0.0f, 1.0f / s_den - fmaxf(0.0f, I) + wr * wr * bkgvar_a[i]);
// Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed
// when the profile result diverges from it.
if (fabsf(I - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) {
I = I_seed[i];
sigma = sigma_seed[i];
var_bkg = varbkg_seed[i];
}
I_o[i] = I; sigma_o[i] = sigma; ok_o[i] = 1;
I_o[i] = I; sigma_o[i] = sigma; varbkg_o[i] = var_bkg; ok_o[i] = 1;
} else ok_o[i] = 0;
}
}
@@ -546,6 +548,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
d_sigma = CudaDevicePtr<float>(new_capacity);
d_bkg = CudaDevicePtr<float>(new_capacity);
d_bkg_var = CudaDevicePtr<float>(new_capacity);
d_var_bkg = CudaDevicePtr<float>(new_capacity);
d_isum = CudaDevicePtr<float>(new_capacity);
d_ninner = CudaDevicePtr<int>(new_capacity);
d_rbin = CudaDevicePtr<int>(new_capacity);
@@ -557,6 +560,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
h_px_x.resize(new_capacity); h_px_y.resize(new_capacity); h_d.resize(new_capacity);
h_I.resize(new_capacity); h_sigma.resize(new_capacity); h_bkg.resize(new_capacity);
h_var_bkg.resize(new_capacity);
h_obs_x.resize(new_capacity); h_obs_y.resize(new_capacity);
h_ok.resize(new_capacity); h_has_obs.resize(new_capacity);
capacity = new_capacity;
@@ -588,7 +592,6 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
BraggGpuParams p{
.W = static_cast<int>(xpixel), .H = static_cast<int>(ypixel),
.r1_sq = r1_sq, .r2 = r2, .r2_sq = r2_sq, .r3 = r3, .r3_sq = r3_sq,
.min_sigma_ratio = min_sigma_ratio,
.R = R, .G = G, .GG = GG,
.bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f,
.empirical = empirical ? 1 : 0,
@@ -613,7 +616,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_invd2, GG);
mark_mask<<<n, threads, 0, *stream>>>(d_px_x, d_px_y, d_mask, p, n);
boxsum<<<n, threads, 0, *stream>>>(d_px_x, d_px_y, d_d, img, d_mask, p, n,
d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y,
d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y,
d_ok, d_strong, d_has_obs, d_invd2,
d_isum, d_ninner, d_rbin,
d_rad_sum, d_rad_cnt, rad_n);
@@ -633,14 +636,15 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
d_shell_grid, d_global_grid, d_shell_n, d_global_n,
d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, p);
fit<<<n, threads, fit_shared_bytes, *stream>>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2,
d_I, d_sigma, d_bkg, d_bkg_var, d_ok, d_shell_P, d_global_P,
d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P,
d_shell_sigma2, d_global_sigma2, d_shell_n,
d_I, d_sigma, d_ok, p, n);
d_I, d_sigma, d_var_bkg, d_ok, p, n);
}
cuda_err(cudaMemcpyAsync(h_I.data(), d_I, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream));
cuda_err(cudaMemcpyAsync(h_sigma.data(), d_sigma, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream));
cuda_err(cudaMemcpyAsync(h_bkg.data(), d_bkg, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream));
cuda_err(cudaMemcpyAsync(h_var_bkg.data(), d_var_bkg, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream));
cuda_err(cudaMemcpyAsync(h_ok.data(), d_ok, sizeof(uint8_t) * npredicted, cudaMemcpyDeviceToHost, *stream));
// Pass A always fills the box-sum centroid (observed spot position), so copy it back in all modes -
// post-refinement uses it as the observed position (beam-centre / distance).
@@ -654,6 +658,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
results[i].I = h_I[i];
results[i].sigma = h_sigma[i];
results[i].bkg = h_bkg[i];
results[i].var_bkg = h_var_bkg[i];
results[i].ok = true;
if (h_has_obs[i]) {
results[i].observed_x = h_obs_x[i];
image_analysis/bragg_integration/BraggIntegrationEngineGPU.h
+2 -2
View File
@@ -29,7 +29,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
// --- 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_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, 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;
@@ -51,7 +51,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
// --- 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<float> h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y;
std::vector<uint8_t> h_ok, h_has_obs;
void EnsureCapacity(size_t n);