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

integration: radially-elongated profile for bandwidth-streaked high-res stills

Browse Source 
With a finite energy bandwidth each reflection is smeared RADIALLY by
sigma_bw = bandwidth_sigma * R_px (R_px = distance from the beam centre, so large
at high resolution): high-resolution spots become radial streaks. The isotropic
per-shell Gaussian both mis-weights them and clips the streak tail on the fixed
profile grid, losing intensity (biased low, noisy).

When a bandwidth is set, fit each reflection with a per-reflection Gaussian
elongated only along its radial direction - sigma^2_radial = sigma^2_intrinsic +
sigma_bw^2, sigma^2_tangential = sigma^2_intrinsic - on a grid grown to hold the
streak. Unlike an isotropic widening this adds no tangential background. It only
engages where the smear exceeds the intrinsic spot (high resolution); low/mid
resolution and monochromatic data (bandwidth 0, e.g. rotation) are untouched.

On the HEWL serial-stills jet (with the background sigma-clip) this lifts the
overall CC-vs-reference 52 -> 55% and the high-resolution I/sig (1.7 A 0.5 -> 1.4),
recovering the 2.0-2.5 A band, with CC1/2 preserved (the per-shot noise the wider
region adds is averaged out by the high serial multiplicity). Rotation ISa 19.1
unchanged.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-28 12:01:56 +02:00
co-authored by Claude Opus 4.8
parent 76e88b0fca
commit 814dff34cb
1 changed files with 60 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/ProfileIntegrate2D.cpp
+60 -5
View File
@@ -184,23 +184,78 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
for (int s = 0; s < N_SHELL; ++s)
shell_P[s] = shell_n[s] >= MIN_STRONG_PER_SHELL ? build_profile(shell_grid[s], shell_n[s]) : global_P;
// Isotropic (intrinsic) width per shell, used to seed the radially-elongated profile below.
auto measure_sigma2 = [&](const std::vector<double> &grid) {
double m2 = 0.0, m2w = 0.0;
for (int dy = -R; dy <= R; ++dy)
for (int dx = -R; dx <= R; ++dx) {
const double g = std::max(0.0, grid[idx(dx, dy)]);
m2 += g * (dx * dx + dy * dy);
m2w += g;
}
return m2w > 0.0 ? std::max(0.25, (m2 / m2w) / 2.0) : 1.0;
};
const double global_sigma2 = global_n > 0 ? measure_sigma2(global_grid) : 1.0;
std::vector<double> shell_sigma2(N_SHELL, global_sigma2);
for (int s = 0; s < N_SHELL; ++s)
if (shell_n[s] >= MIN_STRONG_PER_SHELL)
shell_sigma2[s] = measure_sigma2(shell_grid[s]);
// The energy bandwidth smears a reflection RADIALLY by sigma_bw = bandwidth_sigma * R_px (R_px =
// distance from the beam centre, so large at high resolution). The isotropic profile then both
// mis-weights it and (on the fixed grid) clips its radial tail. In ellipse mode, fit each
// reflection with a per-reflection Gaussian elongated only along the radial direction
// (sigma^2_radial = sigma^2_intrinsic + sigma_bw^2, sigma^2_tangential = sigma^2_intrinsic) on a
// grid grown to hold the streak - no extra tangential background, unlike an isotropic widening.
const double bw_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f;
const float beam_x = geom.GetBeamX_pxl(), beam_y = geom.GetBeamY_pxl();
const bool use_ellipse = !empirical && bw_sigma > 0.0; // only acts when a bandwidth is set
// --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). ---
std::vector<Reflection> out;
out.reserve(npredicted);
const auto pol = experiment.GetPolarizationFactor();
std::vector<double> Pbuf;
for (size_t i = 0; i < npredicted; ++i) {
const auto &rh = rough[i];
if (!rh.ok) continue;
const auto &P = shell_P[rh.shell < 0 ? 0 : rh.shell];
const double B = std::max(rh.bkg, 1.0);
const int sh = rh.shell < 0 ? 0 : rh.shell;
int Rf = R;
const std::vector<double> *Pvec = &shell_P[sh];
if (use_ellipse) {
const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
const double Rpx = std::hypot(rx, ry), sbw = bw_sigma * Rpx;
// Only elongate where the streak is genuinely larger than the intrinsic spot (high
// resolution); at low/mid resolution the smear is sub-pixel and elongating just adds noise.
if (Rpx > 1e-6 && sbw * sbw > shell_sigma2[sh]) {
const double ux = rx / Rpx, uy = ry / Rpx;
const double s2r = shell_sigma2[sh] + sbw * sbw, s2t = shell_sigma2[sh];
Rf = std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(s2r))));
const int Gf = 2 * Rf + 1;
Pbuf.assign(static_cast<size_t>(Gf) * Gf, 0.0);
double gs = 0.0;
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux;
const double g = std::exp(-rad * rad / (2.0 * s2r) - tn * tn / (2.0 * s2t));
Pbuf[(dy + Rf) * Gf + (dx + Rf)] = g;
gs += g;
}
for (double &p : Pbuf) p /= gs;
Pvec = &Pbuf;
}
}
const int Gf = 2 * Rf + 1;
const double B = std::max(rh.bkg, 1.0);
double I = rh.I, den = 0.0;
for (int iter = 0; iter < 4; ++iter) {
double num = 0.0;
den = 0.0;
for (int dy = -R; dy <= R; ++dy)
for (int dx = -R; dx <= R; ++dx) {
const double Pp = P[idx(dx, dy)];
for (int dy = -Rf; dy <= Rf; ++dy)
for (int dx = -Rf; dx <= Rf; ++dx) {
const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
if (Pp <= 0.0) continue;
const int64_t x = rh.cx + dx, y = rh.cy + dy;
if (x < 0 || y < 0 || x >= static_cast<int64_t>(xpixel) || y >= static_cast<int64_t>(ypixel)) continue;