From 814dff34cbb88cd3b2245523db5652bcb8260a80 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 28 Jun 2026 12:01:56 +0200 Subject: [PATCH] integration: radially-elongated profile for bandwidth-streaked high-res stills 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 --- .../bragg_integration/ProfileIntegrate2D.cpp | 65 +++++++++++++++++-- 1 file changed, 60 insertions(+), 5 deletions(-) diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp index 95cfeed1..39433081 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp @@ -184,23 +184,78 @@ std::vector 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 &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 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 out; out.reserve(npredicted); const auto pol = experiment.GetPolarizationFactor(); + std::vector 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 *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(std::ceil(r2 + 2.0 * std::sqrt(s2r)))); + const int Gf = 2 * Rf + 1; + Pbuf.assign(static_cast(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(xpixel) || y >= static_cast(ypixel)) continue;