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

Bragg integration: propagate the background-estimate uncertainty, add an opt-in radial background correction

Browse Source 
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>
This commit is contained in:
leonarski_f
2026-08-06 15:44:13 +02:00
co-authored by Claude Opus 5
parent 52f0e58cae
commit 0e23fd3ab9
10 changed files with 309 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+72 -5
View File
@@ -85,6 +85,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). ---
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
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
int cx = 0, cy = 0, shell = -1;
bool ok = false, strong = false, has_obs = false;
};
@@ -92,6 +96,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
std::vector<int32_t> bkg_vals; // reused per reflection for the trimmed-mean background (idea 1)
// Radial background curve, accumulated from the annulus pixels this pass already reads. A pixel's
// radius is the reflection's radius plus the pixel's projection on the beam->reflection direction,
// so no per-pixel sqrt is needed.
const int n_rad = static_cast<int>(std::ceil(std::hypot(std::max<double>(beam_x, W - beam_x),
std::max<double>(beam_y, H - beam_y)))) + 2;
std::vector<double> rad_sum;
std::vector<int> rad_cnt;
if (bkg_radial) {
rad_sum.assign(n_rad, 0.0);
rad_cnt.assign(n_rad, 0);
}
for (size_t i = 0; i < npredicted; ++i) {
const auto &r = predicted[i];
Rough out;
@@ -100,6 +116,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - r3 - 1.0)));
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + r3 + 1.0)));
// Unit vector beam -> reflection: a stencil pixel's radial offset is its projection on it.
const double rx = r.predicted_x - beam_x, ry = r.predicted_y - beam_y;
const double r0 = std::hypot(rx, ry);
const double ux = r0 > 1e-6 ? rx / r0 : 1.0, uy = r0 > 1e-6 ? ry / r0 : 0.0;
out.r_bin = std::clamp(static_cast<int>(std::lround(r0)), 0, n_rad - 1);
int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0;
double bkg_sum = 0.0;
int n_bkg = 0;
@@ -124,6 +146,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
}
}
int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it)
if (n_inner_valid == n_inner && n_bkg > 5) {
out.bkg = bkg_sum / n_bkg;
if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5) {
@@ -137,6 +160,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
double s = 0.0;
for (size_t t = lo; t < hi; ++t) s += bkg_vals[t];
out.bkg = s / static_cast<double>(hi - lo);
n_bkg_used = static_cast<int>(hi - lo);
}
} else if (do_clip) {
// One high-outlier sigma-clip pass on the background ring: reject pixels above
@@ -151,14 +175,31 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
if (refl_mask[y * W + x]) continue;
const int32_t px = img[y * W + x];
if (!valid(px)) continue;
if (static_cast<double>(px) <= thr) { s += px; ++n; }
if (static_cast<double>(px) <= thr) {
s += px;
++n;
if (bkg_radial) {
const double off = (x - r.predicted_x) * ux + (y - r.predicted_y) * uy;
const int b = std::clamp(static_cast<int>(std::lround(r0 + off)), 0, n_rad - 1);
rad_sum[b] += static_cast<double>(px);
++rad_cnt[b];
}
}
}
if (n > 5) out.bkg = s / n;
if (n > 5) { out.bkg = s / n; n_bkg_used = n; }
}
out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner) * out.bkg;
// I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its
// error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
// 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.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);
if (I_sum > 0) {
out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(I_sum)));
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);
out.obs_y = static_cast<double>(I_sum_y) / static_cast<double>(I_sum);
out.has_obs = true;
@@ -176,6 +217,27 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
rough[i] = out;
}
// --- Radial background curvature correction. The annulus mean is blind to the curvature of the
// radial background (a linear background cancels between the concentric disk and annulus), so
// correct it by mean_annulus(B) - mean_disk(B) taken from the curve just accumulated. Reads no
// pixels: one short dot product per reflection. ---
if (bkg_radial) {
for (size_t i = 0; i < npredicted; ++i) {
auto &rh = rough[i];
if (!rh.ok) continue;
// An empty bin contributes the reflection's own background, so a fully empty
// neighbourhood gives corr == 0 exactly (the kernel weights sum to zero).
double corr = 0.0;
for (int k = 0; k < static_cast<int>(k_diff.size()); ++k) {
const int b = std::clamp(rh.r_bin + k - k_off, 0, n_rad - 1);
const double v = rad_cnt[b] > 0 ? rad_sum[b] / rad_cnt[b] : rh.bkg;
corr += static_cast<double>(k_diff[k]) * v;
}
rh.bkg -= corr; // annulus mean -> mean over the signal disk
rh.I = static_cast<double>(rh.I_sum) - static_cast<double>(rh.n_inner) * rh.bkg;
}
}
// --- BoxSum mode is BraggIntegrate2D: emit the rough result directly. ---
if (mode == IntegratorMode::BoxSum) {
for (size_t i = 0; i < npredicted; ++i) {
@@ -317,10 +379,11 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const int Gf = 2 * Rf + 1;
const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR);
double I = rh.I, den = 0.0;
double I = rh.I, den = 0.0, wsum = 0.0;
for (int iter = 0; iter < 4; ++iter) {
double num = 0.0;
den = 0.0;
wsum = 0.0;
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
@@ -332,6 +395,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
const double v = B + std::max(0.0, I) * Pp;
num += Pp * (static_cast<double>(px) - rh.bkg) / v;
den += Pp * Pp / v;
wsum += Pp / v;
}
if (den > 0.0) I = num / den; else break;
}
@@ -341,7 +405,10 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// iteration has no real peak to lock onto and can manufacture intensity the box sum never sees.
// Keep the profile intensity only if it agrees with the summation seed within the margin;
// otherwise fall back to the summation, which is robust there.
double sigma = std::sqrt(1.0 / den);
// 1/den is the profile-fit variance with the background taken as exact. The fit is
// 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);
if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) {
I = rh.I;
sigma = rh.sigma;