diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 4886611a..1a43be2f 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -44,8 +44,9 @@ class BraggIntegrationSettings { // High-side-only sigma clip: reject ring pixels above mean + this many sqrt(mean). Rejects the same // contamination as the trim - measurably better, in fact - without cutting the low side, so it does // not carry the trim's skew bias. Measured empty-aperture pedestal, counts: plain mean -0.03..-0.20, - // 10% symmetric trim +5.05..+6.34, 4 sigma clip +0.02..+0.54. Applied to monochromatic data - // (rugnux --background-clip); broadband (non-zero bandwidth) data always clip, at their tuned 3 sigma. + // 10% symmetric trim +5.05..+6.34, 4 sigma clip +0.02..+0.54. Whatever is set here is what the + // engine applies (rugnux --background-clip); the front end picks the default, and rugnux lowers it + // to 3 sigma for broadband (non-zero bandwidth) data, where longer spots leak further into the ring. float bkg_clip_nsigma = 4.0f; // Radial background curvature correction. The signal disk and the background annulus are // concentric, so for ANY background linear in position their means are equal - a plane fit buys diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d5d4ee96..7f3ce011 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -13,6 +13,8 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`, `--max-hkl`) instead of a fixed 100, derived per crystal offline and held fixed online. * Bragg integration: The local background ring is made robust with a **high-side sigma clip** (`--background-clip `, default 4) instead of the symmetric trimmed mean; the trim stays reachable with `--background-trim `. * Bragg integration: The **uncertainty of the background estimate** is now propagated into `sigma`, which both engines previously omitted. +* Bragg integration: Profile-fit `sigma` is **no longer inflated on weak reflections** - the fit weights take the signal estimate as it is instead of clamping it at zero, and the per-pixel variance floor is 0.01 counts instead of 1/12. +* Bragg integration: The spot width is now measured **separately along the radial and tangential directions**, and `--bandwidth` no longer changes the background estimator: `--background-clip` and `--background-trim` apply whatever the bandwidth is (the clip still defaults to 3 sigma when a bandwidth is set). * Bragg integration: New **radial background correction** (`--background-radial[=on|off|auto]`, default off); `auto` applies it per image where that image's ice score shows a smooth powder ring. * rugnux: **Ice-ring handling now runs only when the crystal is measured to have ice**, detected on two channels - the spot finder's peak-excluded radial profile (`--ice-min-score`, default 1.5) and found spots on the rings against ice-free flanks (`--ice-min-spot-ratio`, default 2.0). * rugnux: The merge-time **ice-ring mask has been removed**; ice reflections are still kept out of the scale fit and the space-group search, and kept in the final merge. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index ce04b78d..212c653c 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -582,7 +582,7 @@ $ $ with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates $\sigma$ by $\sqrt{1+n_S/n_B}$ — 1.109 with the shipped stencil — uniformly, on every reflection of every dataset. The same term is carried into the profile fit (§9.3), where it adds $(\sum wP/\sum P^2/v)^2\,\mathrm{var}(\hat b)$; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below. -**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ on monochromatic data (`--background-clip`; $n=0$ disables) and $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. +**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), lowered by `rugnux` to $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. That is only a default — the flag sets $n$ whatever the bandwidth is. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. The clip cuts only the high tail, which matters: the **symmetric** trimmed mean it replaced (drop the lowest and highest fraction $f$ of ring pixels, $f=0.10$; still reachable with `--background-trim`, which switches the clip off) is *not* a consistent estimator of the mean of a right-skewed Poisson sample. It sits $\approx0.1$ ct/px **below** the true mean at every level, and with $n_S\approx50$ signal pixels in the $r_1$ disk that under-estimate adds $\approx5$ counts to **every** partial — 0.2 % of the mean partial at 4 Å but $\approx13$ % of the mean and $\approx50$ % of the median partial in the outermost shell. Measured empty-aperture pedestal, in counts: plain mean $-0.03\ldots-0.20$, $10\,\%$ symmetric trim $+5.05\ldots+6.34$, $4\sigma$ clip $+0.02\ldots+0.54$. The same contamination is rejected either way — better, in fact: the trim collapses once contamination exceeds $\approx10\,\%$ of the ring (the same neighbour core shifts it by $+10.1$ ct/px). Note that removing a positive background bias *lowers* $\langle I/\sigma\rangle$ and *raises* edge $R_\text{meas}$, because both are inflated by information-free counts — so neither may be read as evidence against the change. The accuracy gain shows up instead in per-shell agreement with independent processing of the same images, with $CC_{1/2}$ neutral to slightly positive. @@ -601,18 +601,18 @@ The model is a function of **radius alone**, so it is applied only where that is A fixed signal disk captures a *width-dependent* fraction of each spot, which puts a multiplicative floor on the per-observation precision of strong reflections and weights weak reflections poorly. Profile fitting removes this by extracting each intensity against a fitted spot shape, without needing reference intensities. Per frame: 1. **Seed.** Box-sum every reflection (§9.2) to get a rough intensity and observed centroid, and select strong spots (significance $\ge 5$). -2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moment of the strong spots (shell-dependent because spot size grows with resolution); the intrinsic spot is treated as round in the detector plane, the crystal's own anisotropy lying in the rocking direction, which a 2D detector-plane profile does not sample. For `empirical` the profile is instead the averaged, background-subtracted pixel grid of the shell's strong spots, accumulated on their **rounded predicted** positions. For `gaussian` only, the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially: +2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moments of the strong spots (shell-dependent because spot size grows with resolution). The moments are **anisotropic**: each strong spot's pixels are rotated into its *own* radial/tangential frame before being accumulated, giving $\sigma^2_r$ and $\sigma^2_t$ separately. Stacking the spots in the detector frame instead — they sit at every azimuth — averages the two directions away, leaving only $\sigma_r^2+\sigma_t^2$, so radial smearing is read back as a wider *tangential* spot. For `empirical` the profile is instead the averaged, background-subtracted pixel grid of the shell's strong spots, accumulated in the detector frame on their **rounded predicted** positions. For `gaussian` only, the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially: - a finite energy **bandwidth** smears each spot by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ ($R_\mathrm{px}$ = distance from the beam centre, large at high resolution), and - - sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus, on the monochromatic path, a small fixed weak-spot capture term. + - sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus a small fixed weak-spot capture term. - These combine as $\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)$ (tangential unchanged), on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add. The `empirical` profile keeps the fixed per-shell grid and gets none of this. + The two enter as a floor on the measured radial excess: $\sigma^2_\mathrm{radial}=\sigma^2_t+\max\!\left(\sigma^2_r-\sigma^2_t,\ \sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)\right)$, tangential unchanged at $\sigma^2_t$. The measured excess is what the signal disk can resolve; the analytic term takes over for a streak too long to be measured there. The Gaussian is built on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add. The `empirical` profile keeps the fixed per-shell grid and gets none of this. 3. **Fit (Kabsch).** With profile $P$, background $B$ and the shell variance model, the intensity and its uncertainty are $ I = \frac{\sum P\,(c-B)/v}{\sum P^2/v},\qquad \sigma = \sqrt{\frac{1}{\sum P^2/v}},\qquad -v = B + \max(I,0)\,P, +v = \max\!\left(B + I\,P,\ \tfrac{1}{2}B\right), $ -where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the variance floors the background at $1/12$ (the integer-binning pixel-variance floor). The rotation/excitation partiality is carried exactly as in the box-sum path. +where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). The plug-in $I$ enters **as it is**: half-wave rectifying it, $v=B+\max(I,0)P$, lets $v$ — and with it the reported $1/\sum P^2/v$ — respond only to *upward* fluctuations of a noisy estimate, which adds $\approx0.4\,\sigma\sum P^3/(\sum P^2)^2$ to every $\sigma$ whatever the count rate. That offset is invisible on strong reflections and a large fractional inflation on weak ones; the $\tfrac12 B$ clamp keeps $v$ positive without reintroducing it. As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the background term is floored at $0.01$ ct/px — enough to keep $P^2/v$ finite when the ring mean reads exactly zero, which a ring of $n_B$ pixels cannot distinguish from any background below $\approx1/n_B$. The rotation/excitation partiality is carried exactly as in the box-sum path. The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`). diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index c1e3b629..7d239599 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -57,26 +57,25 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe G = 2 * R + 1; GG = G * G; - // A set bandwidth (broadband / stills) vs monochromatic (rotation) splits the treatment: the - // background sigma-clip and radial-elongation terms are path-dependent (see ProfileIntegrate2D). + // The X-ray bandwidth enters ONE place: it smears a reflection radially by bw_sigma * Rpx, which + // the per-reflection Gaussian carries as part of its radial variance. It is not a mode switch - + // the background estimator and the parallax/capture term below are the same whatever the beam is. bw_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f; - broadband = bw_sigma > 0.0; const double c_par = parallax_var_px2(det.GetSensorMaterial(), det.GetSensorThickness_um(), geom.GetWavelength_A(), geom.GetPixelSize_mm() * 1000.0); - c_radial = c_par + (broadband ? 0.0 : bragg_engine::C_CAPTURE); + c_radial = c_par + bragg_engine::C_CAPTURE; F_px = geom.GetDetectorDistance_mm() / std::max(1e-6f, geom.GetPixelSize_mm()); beam_x = geom.GetBeamX_pxl(); beam_y = geom.GetBeamY_pxl(); - use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0); + use_ellipse = !empirical; - // Robust background ring, one estimator or the other (see BraggIntegrationSettings). Broadband - // (non-zero bandwidth: pink-beam / DMM) data keep their tuned 3 sigma high-side clip whatever the - // settings say; monochromatic data - rotation AND stills, the discriminator is the beam, not the - // acquisition mode - take the clip multiplier from settings, and fall back to the symmetric trim - // only when the clip is switched off (rugnux --background-trim). - bkg_clip_nsigma = broadband ? 3.0f : settings.GetBackgroundClipNSigma(); - bkg_trim = (broadband || bkg_clip_nsigma > 0.0f) ? 0.0f : settings.GetBackgroundTrimFraction(); + // Robust background ring, one estimator or the other (see BraggIntegrationSettings): a high-side + // sigma-clip (rugnux --background-clip, the default) or, when the clip is switched off, a + // symmetric trimmed mean (rugnux --background-trim). The caller owns the choice - the engine no + // longer overrides it for broadband data. + bkg_clip_nsigma = settings.GetBackgroundClipNSigma(); + bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); // Radial-offset kernels for the background curvature correction. A stencil pixel at (dx, dy) // sits at radial offset dx*cos(phi) + dy*sin(phi) from the reflection, where phi is the diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 38bf8339..bd5caebf 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -48,7 +48,7 @@ namespace bragg_engine { constexpr int N_SHELL = 6; // resolution shells for per-shell profile learning constexpr double STRONG_I_OVER_SIGMA = 5.0; // strong-spot threshold that seeds the profile constexpr int MIN_STRONG_PER_SHELL = 30; // below this a shell falls back to the global profile -constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (monochromatic only) +constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (coefficient of tan^2(2theta)) // Lower bound on the background term of the Kabsch fit weights (v = max(bkg, floor) + signal). It // guards the background ESTIMATE, not the detector: the r2..r3 ring mean of a high-angle reflection // can come out exactly zero, and v = 0 makes the weights P^2/v diverge. A ring of n pixels cannot @@ -101,7 +101,6 @@ protected: float r3, r3_sq; int R, G, GG; // profile-grid half-size, edge (2R+1) and area (G*G) - bool broadband; // a set bandwidth (stills) vs monochromatic (rotation) double bw_sigma; // bandwidth sigma [dimensionless, * Rpx -> px] float bkg_clip_nsigma; // high-outlier background sigma-clip multiplier (0 = no clip) bool use_ellipse; // radially elongate the per-reflection Gaussian @@ -109,9 +108,9 @@ protected: double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture double F_px; // detector distance expressed in pixels float beam_x, beam_y; - // Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings): the configured - // fraction for monochromatic (rotation) data, forced to 0 for broadband (stills, which keep their - // high-side sigma-clip). 0 = plain ring mean. Read by both the CPU and GPU engines. + // Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings), used only when + // the high-side clip is switched off - the two are alternatives. 0 = plain ring mean. Read by both + // the CPU and GPU engines. float bkg_trim = 0.0f; // --- radial background curvature correction (BraggIntegrationSettings) --- diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index d7ce3f0a..3ee24108 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -262,13 +262,25 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (rough[i].ok) rough[i].shell = shell_of(predicted[i].d); // --- Learn the profile per shell (+ global) from the strong spots. --- + // Two things are learned. The empirical profile is the average grid in the DETECTOR frame, which is + // where it is applied. The width is a pair of second moments taken in each spot's OWN radial / + // tangential frame: a grid stacked in the detector frame is azimuthally averaged, so its is + // sigma_r^2 + sigma_t^2 with no way back, and a radially smeared spot reads as a wide TANGENTIAL + // one. Rotating each contribution into the spot's frame keeps the two apart. + struct Moments { double rad = 0.0, tan = 0.0, w = 0.0; }; + struct Sigma2 { double rad = 1.0, tan = 1.0; }; std::vector> shell_grid(N_SHELL, std::vector(GG, 0.0)); + std::vector shell_mom(N_SHELL); std::vector shell_n(N_SHELL, 0); std::vector global_grid(GG, 0.0); + Moments global_mom; int global_n = 0; for (size_t i = 0; i < npredicted; ++i) { const auto &rh = rough[i]; if (!rh.ok || !rh.strong || rh.I <= 0.0) continue; + const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const double Rpx = std::hypot(rx, ry); + const double ux = Rpx > 1e-6 ? rx / Rpx : 1.0, uy = Rpx > 1e-6 ? ry / Rpx : 0.0; for (int dy = -R; dy <= R; ++dy) for (int dx = -R; dx <= R; ++dx) { const int x = rh.cx + dx, y = rh.cy + dy; @@ -278,65 +290,54 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double v = (static_cast(px) - rh.bkg) / rh.I; shell_grid[rh.shell][grid_idx(dx, dy)] += v; global_grid[grid_idx(dx, dy)] += v; + if (dx * dx + dy * dy >= r1_sq) continue; + const double rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + shell_mom[rh.shell].rad += v * rad * rad; + shell_mom[rh.shell].tan += v * tn * tn; + shell_mom[rh.shell].w += v; + global_mom.rad += v * rad * rad; + global_mom.tan += v * tn * tn; + global_mom.w += v; } ++shell_n[rh.shell]; ++global_n; } - // Isotropic width (2nd moment) of a learned grid: over the r1 disk (monochromatic) or the full - // grid (broadband); = 2 sigma^2 in 2D. The cells are signed: away from the peak a learned - // cell is pure background noise centred on zero, and clamping it at zero turns that noise into a - // positive pedestal spread over the whole domain, which the r^2 weight then reads as extra width. - 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) { - if (!broadband && dx * dx + dy * dy >= r1_sq) continue; - const double g = grid[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; + // Radial and tangential variances from the moments. The domain is the r1 disk, which is + // azimuthally symmetric and so adds no anisotropy of its own. The cells are SIGNED: away from the + // peak a learned cell is background noise centred on zero, and clamping it at zero turns that + // noise into a positive pedestal that the rad^2 / tan^2 weights read as extra width. + auto widths = [](const Moments &m) { + Sigma2 s; + if (m.w > 0.0) { + s.rad = std::max(0.25, m.rad / m.w); + s.tan = std::max(0.25, m.tan / m.w); + } + return s; }; - // Normalised profile (sum = 1): empirical average grid, or an isotropic Gaussian of the measured - // 2nd moment (only used by ProfileEmpirical; ProfileGaussian rebuilds per reflection in Pass B). - auto build_profile = [&](const std::vector &grid, int n) { + // Normalised empirical profile (sum = 1), the average grid over the strong spots of a shell. + // ProfileGaussian does not use it - it rebuilds a per-reflection Gaussian in Pass B. + auto build_profile = [&](const std::vector &grid) { std::vector P(GG, 0.0); - if (n <= 0) return P; double sum = 0.0; for (int k = 0; k < GG; ++k) { - const double g = std::max(0.0, grid[k]); - sum += g; - if (empirical) P[k] = g; + P[k] = std::max(0.0, grid[k]); + sum += P[k]; } - if (sum <= 0.0) return P; - if (empirical) { + if (sum > 0.0) for (double &p : P) p /= sum; - } else { - const double sigma2 = measure_sigma2(grid); - double gsum = 0.0; - for (int dy = -R; dy <= R; ++dy) - for (int dx = -R; dx <= R; ++dx) { - const double g = std::exp(-(dx * dx + dy * dy) / (2.0 * sigma2)); - P[grid_idx(dx, dy)] = g; - gsum += g; - } - for (double &p : P) p /= gsum; - } return P; }; - const std::vector global_P = build_profile(global_grid, global_n); - const double global_sigma2 = global_n > 0 ? measure_sigma2(global_grid) : 1.0; - std::vector> shell_P(N_SHELL); - std::vector shell_sigma2(N_SHELL, global_sigma2); + const std::vector global_P = empirical && global_n > 0 ? build_profile(global_grid) + : std::vector(GG, 0.0); + const Sigma2 global_sigma2 = widths(global_mom); + std::vector> shell_P(N_SHELL, global_P); + 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_P[s] = build_profile(shell_grid[s], shell_n[s]); - shell_sigma2[s] = measure_sigma2(shell_grid[s]); - } else { - shell_P[s] = global_P; - } + if (shell_n[s] < MIN_STRONG_PER_SHELL) continue; + if (empirical) shell_P[s] = build_profile(shell_grid[s]); + shell_sigma2[s] = widths(shell_mom[s]); } // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). --- @@ -352,12 +353,16 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; const double Rpx = std::hypot(rx, ry); const double tan2t = Rpx / F_px; - const double s2t = shell_sigma2[sh]; + const double s2t = shell_sigma2[sh].tan; double s2r = s2t, ux = 1.0, uy = 0.0; bool elong = false; if (use_ellipse) { + // Radial excess over the tangential width: measured where the peak is resolved inside + // the r1 disk, with the analytic bandwidth + parallax/capture term as the floor. The + // analytic term is what carries a streak the disk is too small to measure. const double sbw = bw_sigma * Rpx; - const double radial_extra = sbw * sbw + c_radial * tan2t * tan2t; + const double radial_extra = std::max(shell_sigma2[sh].rad - s2t, + sbw * sbw + c_radial * tan2t * tan2t); if (Rpx > 1e-6 && radial_extra > 0.25) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c132c6b3..63530fa1 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -19,7 +19,6 @@ struct BraggGpuParams { int R, G, GG; float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) int empirical; // ProfileEmpirical vs ProfileGaussian - int broadband; int use_ellipse; float bw_sigma; float c_radial; @@ -33,6 +32,9 @@ __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MA // idea 1: max annulus pixels held in shared memory for the per-block trimmed-mean sort (else flat mean). constexpr int BKG_TRIM_MAX = 512; +// Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL. +constexpr int MOM_STRIDE = 3; + // --- Mark the r2 signal disk of every predicted reflection (race-free: all writes are 1). --- __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, BraggGpuParams p, int n) { const int i = blockIdx.x; @@ -276,12 +278,13 @@ __device__ inline int compute_shell(float d, const unsigned long long *invd2mm) // --- Zero the profile accumulators and seed the inv-d^2 range, in one launch (replaces a handful // of small cudaMemsetAsync calls, which matter when kernel-launch latency is high). --- -__global__ void reset(float *shell_grid, float *global_grid, int *shell_n, int *global_n, +__global__ void reset(float *shell_grid, float *global_grid, float *mom, int *shell_n, int *global_n, unsigned long long *invd2mm, int GG) { for (int k = blockIdx.x * blockDim.x + threadIdx.x; k < N_SHELL * GG; k += blockDim.x * gridDim.x) shell_grid[k] = 0.0f; for (int k = blockIdx.x * blockDim.x + threadIdx.x; k < GG; k += blockDim.x * gridDim.x) global_grid[k] = 0.0f; + if (blockIdx.x == 0 && threadIdx.x < MOM_STRIDE * (N_SHELL + 1)) mom[threadIdx.x] = 0.0f; if (blockIdx.x == 0 && threadIdx.x < N_SHELL) shell_n[threadIdx.x] = 0; if (blockIdx.x == 0 && threadIdx.x == 0) { *global_n = 0; @@ -291,11 +294,16 @@ __global__ void reset(float *shell_grid, float *global_grid, int *shell_n, int * } // --- Learn the profile: each strong spot adds its bkg-subtracted, I-normalised grid to its shell -// (and the global grid). One block per reflection. --- -__global__ void learn_profile(const int32_t *img, const int *cx_a, const int *cy_a, const float *dd, +// (and the global grid), and its second moments to the same shell's moment slot. The grid stays +// in the DETECTOR frame (that is where the empirical profile is applied); the moments are taken in +// the spot's own radial/tangential frame, because a detector-frame stack is azimuthally averaged +// and cannot tell a radially smeared spot from a tangentially wide one (see the CPU engine). +// One block per reflection. --- +__global__ void learn_profile(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_a, const float *bkg_a, const uint8_t *ok_a, const uint8_t *strong_a, - float *shell_grid, float *global_grid, int *shell_n, int *global_n, + float *shell_grid, float *global_grid, float *mom, int *shell_n, int *global_n, BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n || !ok_a[i] || !strong_a[i]) return; @@ -303,57 +311,75 @@ __global__ void learn_profile(const int32_t *img, const int *cx_a, const int *cy if (!(I > 0.0f)) return; const int cx = cx_a[i], cy = cy_a[i], sh = compute_shell(dd[i], invd2mm); const float bkg = bkg_a[i]; + const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; + const float Rpx = sqrtf(rx * rx + ry * ry); + const float ux = Rpx > 1e-6f ? rx / Rpx : 1.0f, uy = Rpx > 1e-6f ? ry / Rpx : 0.0f; float *sg = shell_grid + (size_t) sh * p.GG; + __shared__ float s_rad, s_tan, s_w; + if (threadIdx.x == 0) { s_rad = 0.0f; s_tan = 0.0f; s_w = 0.0f; } + __syncthreads(); + float l_rad = 0.0f, l_tan = 0.0f, l_w = 0.0f; for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { - const int x = cx + (k % p.G - p.R), y = cy + (k / p.G - p.R); + const int dx = k % p.G - p.R, dy = k / p.G - p.R; + const int x = cx + dx, y = cy + dy; if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const float v = ((float) px - bkg) / I; atomicAdd(&sg[k], v); atomicAdd(&global_grid[k], v); + if ((float) (dx * dx + dy * dy) >= p.r1_sq) continue; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + l_rad += v * rad * rad; l_tan += v * tn * tn; l_w += v; + } + // Reduce in shared memory first: every block adds to the same handful of global words, so one + // atomic per block instead of one per thread. + atomicAdd(&s_rad, l_rad); atomicAdd(&s_tan, l_tan); atomicAdd(&s_w, l_w); + __syncthreads(); + if (threadIdx.x == 0) { + atomicAdd(&mom[MOM_STRIDE * sh + 0], s_rad); + atomicAdd(&mom[MOM_STRIDE * sh + 1], s_tan); + atomicAdd(&mom[MOM_STRIDE * sh + 2], s_w); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 0], s_rad); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 1], s_tan); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 2], s_w); + atomicAdd(&shell_n[sh], 1); + atomicAdd(global_n, 1); } - if (threadIdx.x == 0) { atomicAdd(&shell_n[sh], 1); atomicAdd(global_n, 1); } } -// --- Reduce each learned grid to its 2nd-moment width (and, for empirical, a normalised profile). -// One block per grid: blocks [0,N_SHELL) are the shells, block N_SHELL is the global grid. --- -__global__ void build_profiles(const float *shell_grid, const float *global_grid, +// --- Turn the learned moments into radial/tangential variances and, for empirical, normalise the +// learned grid into a profile. One block per grid: blocks [0,N_SHELL) are the shells, block +// N_SHELL is the global one. --- +__global__ void build_profiles(const float *shell_grid, const float *global_grid, const float *mom, const int *shell_n, const int *global_n, - float *shell_P, float *global_P, float *shell_sigma2, float *global_sigma2, + float *shell_P, float *global_P, float *sigma2_r, float *sigma2_t, BraggGpuParams p) { const int b = blockIdx.x; - const float *grid; int nstrong; float *P; float *sig2; - if (b < N_SHELL) { grid = shell_grid + (size_t) b * p.GG; nstrong = shell_n[b]; P = shell_P + (size_t) b * p.GG; sig2 = &shell_sigma2[b]; } - else { grid = global_grid; nstrong = *global_n; P = global_P; sig2 = global_sigma2; } + const float *grid; int nstrong; float *P; + if (b < N_SHELL) { grid = shell_grid + (size_t) b * p.GG; nstrong = shell_n[b]; P = shell_P + (size_t) b * p.GG; } + else { grid = global_grid; nstrong = *global_n; P = global_P; } - __shared__ float s_m2, s_m2w, s_sum; - if (threadIdx.x == 0) { s_m2 = 0.0f; s_m2w = 0.0f; s_sum = 0.0f; } + if (threadIdx.x == 0) { + const float w = mom[MOM_STRIDE * b + 2]; + sigma2_r[b] = w > 0.0f ? fmaxf(0.25f, mom[MOM_STRIDE * b + 0] / w) : 1.0f; + sigma2_t[b] = w > 0.0f ? fmaxf(0.25f, mom[MOM_STRIDE * b + 1] / w) : 1.0f; + } + if (!p.empirical) return; + + __shared__ float s_sum; + if (threadIdx.x == 0) s_sum = 0.0f; __syncthreads(); - - float l_m2 = 0.0f, l_m2w = 0.0f, l_sum = 0.0f; + float l_sum = 0.0f; for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { - const int dx = k % p.G - p.R, dy = k / p.G - p.R; - // The moment reads the SIGNED cell (see the CPU engine): clamping pure-background noise at - // zero leaves a positive pedestal that the r^2 weight reads as extra width. The empirical - // profile still takes the clamped grid - a profile has to be non-negative. - const float raw = grid[k], g = fmaxf(0.0f, raw); - const int r2i = dx * dx + dy * dy; - if (p.broadband || (float) r2i < p.r1_sq) { l_m2 += raw * (float) r2i; l_m2w += raw; } + const float g = fmaxf(0.0f, grid[k]); // a profile has to be non-negative + P[k] = g; l_sum += g; - if (p.empirical) P[k] = g; // pre-store clamped grid for in-place normalisation below } - atomicAdd(&s_m2, l_m2); atomicAdd(&s_m2w, l_m2w); atomicAdd(&s_sum, l_sum); + atomicAdd(&s_sum, l_sum); __syncthreads(); - - if (threadIdx.x == 0) - *sig2 = (nstrong > 0 && s_m2w > 0.0f) ? fmaxf(0.25f, (s_m2 / s_m2w) / 2.0f) : 1.0f; - - if (p.empirical) { - const float sum = s_sum; - const bool normalise = nstrong > 0 && sum > 0.0f; - for (int k = threadIdx.x; k < p.GG; k += blockDim.x) P[k] = normalise ? P[k] / sum : 0.0f; - } + const bool normalise = nstrong > 0 && s_sum > 0.0f; + for (int k = threadIdx.x; k < p.GG; k += blockDim.x) P[k] = normalise ? P[k] / s_sum : 0.0f; } // --- Radial background curvature correction: one thread per reflection, no pixel reads. @@ -388,7 +414,7 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const float *I_seed, const float *sigma_seed, const float *bkg_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, + const float *sigma2_r, const float *sigma2_t, const int *shell_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; @@ -410,15 +436,18 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, for (int k = threadIdx.x; k < p.GG; k += blockDim.x) Pbuf[k] = Psrc[k]; __syncthreads(); } else { - const float s2t = use_shell ? shell_sigma2[sh] : *global_sigma2; + const int si = use_shell ? sh : N_SHELL; // N_SHELL is the global slot + const float s2t = sigma2_t[si]; const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; const float Rpx = sqrtf(rx * rx + ry * ry); const float tan2t = Rpx / p.F_px; float s2r = s2t, ux = 1.0f, uy = 0.0f; bool elong = false; if (p.use_ellipse) { + // Measured radial excess over the tangential width, floored by the analytic bandwidth + + // parallax/capture term (see the CPU engine). const float sbw = p.bw_sigma * Rpx; - const float radial_extra = sbw * sbw + p.c_radial * tan2t * tan2t; + const float radial_extra = fmaxf(sigma2_r[si] - s2t, sbw * sbw + p.c_radial * tan2t * tan2t); if (Rpx > 1e-6f && radial_extra > 0.25f) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; elong = true; } } const int Rf = elong ? min(3 * p.R, (int) ceilf(p.r2 + 2.0f * sqrtf(s2r))) : p.R; @@ -498,8 +527,9 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment d_global_grid(GG), d_shell_P(static_cast(bragg_engine::N_SHELL) * GG), d_global_P(GG), - d_shell_sigma2(bragg_engine::N_SHELL), - d_global_sigma2(1), + d_mom(MOM_STRIDE * (bragg_engine::N_SHELL + 1)), + d_sigma2_r(bragg_engine::N_SHELL + 1), + d_sigma2_t(bragg_engine::N_SHELL + 1), d_shell_n(bragg_engine::N_SHELL), d_global_n(1), d_invd2(2) { @@ -598,7 +628,6 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .R = R, .G = G, .GG = GG, .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f, .empirical = empirical ? 1 : 0, - .broadband = broadband ? 1 : 0, .use_ellipse = use_ellipse ? 1 : 0, .bw_sigma = static_cast(bw_sigma), .c_radial = static_cast(c_radial), .F_px = static_cast(F_px), @@ -616,7 +645,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(float) * n_rad, *stream)); cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); } - reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_invd2, GG); + reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); mark_mask<<>>(d_px_x, d_px_y, d_mask, p, n); boxsum<<>>(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_var_bkg, d_obs_x, d_obs_y, @@ -633,14 +662,15 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (mode != IntegratorMode::BoxSum) { // Pass B: learn (shell computed inline) -> build -> fit. - learn_profile<<>>(img, d_cx, d_cy, d_d, d_invd2, d_I, d_bkg, d_ok, d_strong, - d_shell_grid, d_global_grid, d_shell_n, d_global_n, p, n); + learn_profile<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + d_I, d_bkg, d_ok, d_strong, + d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, p, n); build_profiles<<>>( - 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); + d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, + d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, p); fit<<>>(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_var_bkg, d_ok, d_shell_P, d_global_P, - d_shell_sigma2, d_global_sigma2, d_shell_n, + d_sigma2_r, d_sigma2_t, d_shell_n, d_I, d_sigma, d_var_bkg, d_ok, p, n); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index a9d0e4eb..e698c941 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -45,7 +45,8 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_mask; // per-pixel r2-disk reflection mask CudaDevicePtr d_shell_grid, d_global_grid; // learned profile accumulators (N_SHELL*GG, GG) CudaDevicePtr d_shell_P, d_global_P; // normalised profiles (empirical mode) - CudaDevicePtr d_shell_sigma2, d_global_sigma2; + CudaDevicePtr d_mom; // learned 2nd moments, 3 per shell + global + CudaDevicePtr d_sigma2_r, d_sigma2_t; // radial/tangential widths, N_SHELL + global CudaDevicePtr d_shell_n, d_global_n; CudaDevicePtr d_invd2; // [min,max] inv-d^2 as monotonic bit patterns diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e7a39f82..71aee9ab 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -47,6 +47,10 @@ // in its favour and no measurable cost in wall clock. It is also what jfjoch_viewer already sends. constexpr int64_t RUGNUX_MAX_SPOT_COUNT = 1000; +// Default background-ring high-side clip for broadband (non-zero bandwidth) data, tighter than the +// monochromatic default because a pink-beam reflection is long and its wings reach into the ring. +constexpr float BROADBAND_BACKGROUND_CLIP_NSIGMA = 3.0f; + // Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation. constexpr double SMOOTH_G_DEFAULT_DEG = 5.0; @@ -156,9 +160,9 @@ void print_usage() { std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; std::cout << " --integration-high-resolution High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl; std::cout << " --max-hkl Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl; - std::cout << " --background-clip Monochromatic (rotation + still): high-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3 sigma; ignored by --integrator boxsum" << std::endl; + std::cout << " --background-clip High-side clip of the background ring at mean + n*sqrt(mean) (default 4, or 3 when --bandwidth is set; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl; std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl; - std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl; + std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean. Applies whatever --bandwidth is set to" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; @@ -1847,8 +1851,8 @@ static int RunRugnux(int argc, char **argv) { BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); bis.BackgroundTrimFraction(static_cast(*background_trim_arg)); experiment.ImportBraggIntegrationSettings(bis); - logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip " - "(monochromatic data; broadband always clips)", *background_trim_arg); + logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip", + *background_trim_arg); } if (background_clip_arg) { @@ -1856,6 +1860,15 @@ static int RunRugnux(int argc, char **argv) { bis.BackgroundClipNSigma(static_cast(*background_clip_arg)); experiment.ImportBraggIntegrationSettings(bis); logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg); + } else if (!background_trim_arg && experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f) { + // Broadband (pink-beam / DMM) reflections are long and their wings reach into the background + // ring, so the default clip is tighter there. It is only a default - --background-clip and + // --background-trim override it as they do for any other data. + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.BackgroundClipNSigma(BROADBAND_BACKGROUND_CLIP_NSIGMA); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Background ring: high-side clip at {:.1f} sigma (broadband default)", + BROADBAND_BACKGROUND_CLIP_NSIGMA); } if (background_radial_given) {