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

integration: sigma-clip the background ring to de-bias high-resolution stills

Browse Source 
Weak high-resolution still reflections were systematically over-subtracted: a
bandwidth-streaked high-res spot (or a neighbour) leaks into the r2-r3 background
annulus and biases its mean high, so the subtracted background is too large and
the merged high-resolution intensities go negative (seen as reproducibly negative
<I/sig> at 100% completeness and high multiplicity past ~1.9 A).

Add one high-outlier sigma-clip pass to the box-sum background (reject ring pixels
above mean + 3*sqrt(mean), recompute) so the contamination no longer inflates it.
A clean Poisson background is essentially unchanged (~0.1% exceed the cut). On the
HEWL serial-stills jet this de-biases the high-res band - <I/sig> 2.03 A 0.9 -> 1.6,
1.79 A -0.1 -> +0.7 - extends the usable resolution ~2.2 -> ~2.0 A and improves
overall R-meas 130 -> 124%, with CC1/2 and CC-vs-reference neutral. The rotation
crystal is unchanged (ISa 19.1), its clean backgrounds being barely clipped.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-28 11:50:31 +02:00
co-authored by Claude Opus 4.8
parent 7bb25979e6
commit 76e88b0fca
1 changed files with 14 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/ProfileIntegrate2D.cpp
+14
View File
@@ -39,6 +39,8 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6
int64_t I_sum = 0, n_inner = 0, n_inner_valid = 0; int64_t I_sum = 0, n_inner = 0, n_inner_valid = 0;
double bkg_sum = 0.0; double bkg_sum = 0.0;
int n_bkg = 0; int n_bkg = 0;
std::vector<double> bkg_vals;
bkg_vals.reserve(128);
for (int64_t y = y0; y <= y1; ++y) for (int64_t y = y0; y <= y1; ++y)
for (int64_t x = x0; x <= x1; ++x) { for (int64_t x = x0; x <= x1; ++x) {
const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y); const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
@@ -52,11 +54,23 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6
if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue; if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
bkg_sum += static_cast<double>(px); bkg_sum += static_cast<double>(px);
++n_bkg; ++n_bkg;
bkg_vals.push_back(static_cast<double>(px));
} }
} }
if (n_inner_valid != n_inner || n_bkg <= 5) if (n_inner_valid != n_inner || n_bkg <= 5)
return out; return out;
out.bkg = bkg_sum / n_bkg; out.bkg = bkg_sum / n_bkg;
// One high-outlier sigma-clip pass on the background ring. A bandwidth-streaked high-resolution spot
// (or a neighbour) leaks into the annulus and biases the mean high, which over-subtracts and drives
// the merged high-resolution intensities negative; rejecting pixels above mean + 3*sqrt(mean) removes
// that contamination (and barely touches a clean Poisson background, where ~0.1% exceed the cut).
{
const double thr = out.bkg + 3.0 * std::sqrt(std::max(out.bkg, 1.0));
double s = 0.0;
int n = 0;
for (const double v : bkg_vals) if (v <= thr) { s += v; ++n; }
if (n > 5) out.bkg = s / n;
}
out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner) * out.bkg; out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner) * out.bkg;
out.sigma = std::max(1.0, out.I * min_sigma_ratio); out.sigma = std::max(1.0, out.I * min_sigma_ratio);
if (I_sum > 0) if (I_sum > 0)