diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 00f0fb0d..e1fad22f 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -103,6 +103,8 @@ BraggIntegrationSettings &BraggIntegrationSettings::BackgroundTrimFraction(float check_min("Background trim fraction", input, 0.0); check_max("Background trim fraction", input, 0.49); // must leave a central majority after trimming bkg_trim_fraction = input; + if (input > 0.0f) + bkg_clip_nsigma = 0.0f; // the two ring estimators are alternatives, not a stack return *this; } @@ -124,3 +126,16 @@ BraggIntegrationSettings &BraggIntegrationSettings::MaxHKL(std::optional in std::optional BraggIntegrationSettings::GetMaxHKL() const { return max_hkl; } + +BraggIntegrationSettings &BraggIntegrationSettings::BackgroundClipNSigma(float input) { + check_finite("Background clip nsigma", input); + check_min("Background clip nsigma", input, 0.0); + bkg_clip_nsigma = input; + if (input > 0.0f) + bkg_trim_fraction = 0.0f; // the two ring estimators are alternatives, not a stack + return *this; +} + +float BraggIntegrationSettings::GetBackgroundClipNSigma() const { + return bkg_clip_nsigma; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 4f917155..074168f2 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -29,13 +29,24 @@ class BraggIntegrationSettings { std::optional d_min_limit_A; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; - // Symmetric trimmed-mean fraction for the r2..r3 background ring: drop the lowest and highest this - // fraction of ring pixels before averaging. Resists the high-side contamination (neighbour-spot - // wings, tails, zingers) that biases the plain ring mean up and makes it over-subtract weak - // high-angle reflections. Applied to monochromatic data (rotation and still); the integration engine - // keeps the tuned high-side sigma-clip for broadband (non-zero bandwidth) data instead. 0 = plain ring - // mean (rugnux --background-trim). - float bkg_trim_fraction = 0.10f; + // The r2..r3 background ring is estimated with ONE of two robust means, never both: a high-side + // sigma-clip (bkg_clip_nsigma, the default) or a symmetric trimmed mean (bkg_trim_fraction). Setting + // either through its setter clears the other, so whichever was asked for last is the one in force; + // with both at 0 the ring is a plain mean. + // + // Symmetric trimmed-mean fraction: drop the lowest and highest this fraction of ring pixels before + // averaging. Robust to the high-side contamination (neighbour-spot wings, tails, zingers) that + // biases a plain ring mean up, but a symmetric trim is NOT a consistent estimator of the mean of a + // right-skewed (Poisson) sample - it sits ~0.1 ct/px low at every level, which with ~50 ring pixels + // adds ~5 counts to every partial. Kept reachable (rugnux --background-trim) for back compatibility; + // 0.10 was the shipped value. + float bkg_trim_fraction = 0.0f; + // 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. + float bkg_clip_nsigma = 4.0f; // Half-width of the hkl cube the predictor walks: every reflection with |h|,|k|,|l| <= this is // tested against the Ewald sphere, and nothing outside it can ever be predicted. An axis is // truncated once a/d_min exceeds this, and the GPU cost is the cube (2n+1)^3 of candidates, so @@ -55,6 +66,7 @@ public: BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); BraggIntegrationSettings& BackgroundTrimFraction(float input); + BraggIntegrationSettings& BackgroundClipNSigma(float input); BraggIntegrationSettings& MaxHKL(std::optional input); @@ -67,5 +79,6 @@ public: [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; [[nodiscard]] float GetBackgroundTrimFraction() const; + [[nodiscard]] float GetBackgroundClipNSigma() const; [[nodiscard]] std::optional GetMaxHKL() const; }; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a8d4bc9b..95eebdfb 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -7,6 +7,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Spot finding: Serial stills pick `--min-pix-per-spot` per image; connected components run on the GPU; detection is configurable over the API, where a `high_resolution_limit` of 0 means "no limit" instead of throwing. * Resolution limits: Bragg integration, azimuthal integration and spot finding all default to **as far as the detector reaches**, replacing a fixed 1.0 Å integration limit and a 1.5 Å rotation spot-finding limit that silently discarded everything beyond them; `--integration-high-resolution` and `--spot-high-resolution` still set one by hand. * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`) instead of a fixed 100, derived per crystal offline from the refined cell (`--max-hkl` overrides); the broker keeps a fixed bootstrap so a live acquisition has a predictable per-image cost. +* Bragg integration: The local background ring is now made robust with a **high-side sigma clip** (`--background-clip `, default 4) instead of the symmetric trimmed mean, which is biased low on Poisson data and added ~5 counts to every partial. The trim stays reachable with `--background-trim `. Expect `` to fall and edge `R_meas` to rise - that is the removed bias, not a regression; per-shell agreement with independent processing improves. * rugnux: De-novo **space-group search** substantially more robust - centering ranked by net absences and judged on absent-class strength, merohedral-twin over-promotion vetoed, and genuine high-symmetry groups recovered on weak data. * rugnux: The space-group search takes systematic absences from the merge of all observations, needs at least three control reflections on an axial row to claim a **screw axis**, and no longer alters the production merge. * rugnux: The per-image mosaicity is fitted from the strongest 250 spots, so the indexing spot budget no longer sets it. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 6b56b11c..e564c9df 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -552,9 +552,11 @@ $ $ with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. 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. -**Trimmed-mean background (monochromatic, default on).** On monochromatic data — rotation *and* still (the discriminator is the beam, not the acquisition mode) — the ring mean $\hat{b}$ is by default replaced with a **symmetric trimmed mean**: the ring pixels are sorted, the lowest and highest fraction $f$ are dropped, and the central $(1-2f)$ are averaged ($f=0.10$ by default, `--background-trim`; $f=0$ restores the plain mean). 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. The plain mean reads high there because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants; dropping the extreme ring pixels removes that bias while a clean Poisson ring is essentially unchanged. The trimmed estimate has a slightly higher variance than the plain mean, so the gain is at the resolution edge and the cost, where any, is in already-clean shells. It is applied to the shared background used by both the box sum and the profile fit. +**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. -For the **profile-fit path on broadband (non-zero bandwidth: pink-beam / DMM) data**, the trimmed mean is *not* used; instead the background mean is computed with a single high-outlier reject (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, then recompute): a bandwidth-streaked high-resolution spot or a close neighbour can leak into the ring and bias the mean high, over-subtracting and driving weak high-resolution intensities negative. A clean Poisson background is essentially unchanged by the cut. This $\sigma$-clip is the one robustification that is *not* applied to plain box summation (`--integrator boxsum`); the trimmed mean above is computed in the shared background pass and so applies to box summation too. +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. + +Both estimators are computed in the shared background pass, but only the trim reaches plain box summation: the high-side clip is skipped for `--integrator boxsum`, which therefore uses the plain ring mean unless `--background-trim` is given. ### 9.3 Profile-fitted extraction (default) diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index c8baafa1..a0a25dc0 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -284,7 +284,8 @@ Integration: | --- | --- | | `--integrator ` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) | | `--integration-radius ` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` | -| `--background-trim ` | Monochromatic (rotation + still): symmetric trimmed-mean fraction for the background ring, 0≤f<0.5 (default 0.10; 0 = plain mean) — removes the high-side bias that over-subtracts weak high-angle spots | +| `--background-clip ` | Monochromatic (rotation + still): high-side clip of the background ring at `mean + n·√mean` (default 4; 0 = off). The default background estimator — it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3σ; ignored by `--integrator boxsum` | +| `--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 | | `--integration-high-resolution ` | High-resolution limit for prediction and integration. Omitted (or 0) means integration extends as far as the detector reaches — which is what the predictor can place on the detector anyway, since it rejects reflections that miss it. Set a value to integrate less than the detector offers | | `--max-hkl ` | Predict reflections with \|h\|,\|k\|,\|l\| ≤ `n` (max 511). By default this is derived per crystal from the refined cell as `ceil(max(a,b,c)/d_min) + 1`, which is the exact bound: the predictor keeps only \|q\| ≤ 1/d_min and `h = a·q`, so no reflection can lie outside it and no candidate inside it is wasted on a shorter axis. Set it only to override that | | `--bandwidth ` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) | diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index f89337bb..58a4e2fb 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -58,7 +58,6 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // background sigma-clip and radial-elongation terms are path-dependent (see ProfileIntegrate2D). bw_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f; broadband = bw_sigma > 0.0; - apply_bkg_clip = broadband; const double c_par = parallax_var_px2(det.GetSensorMaterial(), det.GetSensorThickness_um(), geom.GetWavelength_A(), geom.GetPixelSize_mm() * 1000.0); @@ -68,10 +67,13 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe beam_y = geom.GetBeamY_pxl(); use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0); - // Trimmed-mean background applies to monochromatic data - both rotation AND stills (the discriminator - // is the beam, not the acquisition mode). Broadband (non-zero bandwidth: pink-beam / DMM) data keep the - // tuned high-side sigma-clip instead, so the trim is forced off there. The fraction comes from settings. - bkg_trim = broadband ? 0.0f : settings.GetBackgroundTrimFraction(); + // 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(); polarization = experiment.GetPolarizationFactor(); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index bcac7217..a2c39c9a 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -91,7 +91,7 @@ protected: bool broadband; // a set bandwidth (stills) vs monochromatic (rotation) double bw_sigma; // bandwidth sigma [dimensionless, * Rpx -> px] - bool apply_bkg_clip; // stills-only high-outlier background sigma-clip + float bkg_clip_nsigma; // high-outlier background sigma-clip multiplier (0 = no clip) bool use_ellipse; // radially elongate the per-reflection Gaussian double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index d4371d71..2b2c9b14 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -55,12 +55,13 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, return Finalize(predicted, npredicted, results, image_number); const int W = static_cast(xpixel), H = static_cast(ypixel); - const bool do_clip = apply_bkg_clip && mode != IntegratorMode::BoxSum; + const bool do_clip = bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum; // Symmetric trimmed-mean background fraction (BraggIntegrationSettings, rugnux --background-trim): // replaces the r2..r3 ring MEAN with an f-trimmed mean, robust to the high-side contamination // (neighbour wings, tails) that biases the plain mean up and makes it over-subtract weak high-angle - // reflections. Already forced to 0 for broadband/stills by the base ctor (they keep the clip below). + // reflections. The base ctor has already forced it to 0 whenever the clip below is in force, which + // is the default - the two are alternatives. const double bkg_trim_frac = bkg_trim; auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); }; @@ -138,9 +139,9 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, out.bkg = s / static_cast(hi - lo); } } else if (do_clip) { - // One high-outlier sigma-clip pass on the background ring (stills-only): reject pixels - // above mean + 3*sqrt(mean) to strip a bandwidth-streaked neighbour that biases the mean. - const double thr = out.bkg + 3.0 * std::sqrt(std::max(out.bkg, 1.0)); + // One high-outlier sigma-clip pass on the background ring: reject pixels above + // mean + n*sqrt(mean) to strip a neighbour core or a zinger that biases the mean. + const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(std::max(out.bkg, 1.0)); double s = 0.0; int n = 0; for (int y = y0; y <= y1; ++y) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index de4401de..42d5131e 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -18,7 +18,7 @@ struct BraggGpuParams { float r1_sq, r2, r2_sq, r3, r3_sq; float min_sigma_ratio; int R, G, GG; - int do_clip; // background sigma-clip (stills, profile modes) + float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) int empirical; // ProfileEmpirical vs ProfileGaussian int broadband; int use_ellipse; @@ -110,7 +110,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (threadIdx.x == 0) { s_accept = (s_ninner_valid == s_ninner && s_nbkg > 5) ? 1 : 0; s_bkg = s_accept ? (s_bkgsum / (double) s_nbkg) : 0.0; - s_thr = s_bkg + 3.0 * sqrt(fmax(s_bkg, 1.0)); + s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0)); s_clipsum = 0.0; s_clipn = 0; } __syncthreads(); @@ -164,8 +164,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, __syncthreads(); } - // Second ring pass for the stills sigma-clip (re-reads the annulus; avoids storing bkg values). - if (s_accept && p.do_clip && !do_trim) { + // Second ring pass for the high-side sigma-clip (re-reads the annulus; avoids storing bkg values). + if (s_accept && p.bkg_clip_nsigma > 0.0f && !do_trim) { double c_l = 0.0; int cn_l = 0; for (int t = threadIdx.x; t < area; t += blockDim.x) { const int x = x0 + t % bw, y = y0 + t / bw; @@ -185,7 +185,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (!s_accept) { ok_o[i] = 0; strong_o[i] = 0; hasobs_o[i] = 0; return; } double bkg = s_bkg; - if (p.do_clip && s_clipn > 5) bkg = s_clipsum / (double) s_clipn; + if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) bkg = s_clipsum / (double) s_clipn; const long long Isum = (long long) s_Isum; const double I = (double) Isum - (double) s_ninner * bkg; double sigma = fmax(1.0, I * (double) p.min_sigma_ratio); @@ -491,7 +491,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .r1_sq = r1_sq, .r2 = r2, .r2_sq = r2_sq, .r3 = r3, .r3_sq = r3_sq, .min_sigma_ratio = min_sigma_ratio, .R = R, .G = G, .GG = GG, - .do_clip = (apply_bkg_clip && mode != IntegratorMode::BoxSum) ? 1 : 0, + .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, diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 0e73f3bd..91d5769d 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -120,10 +120,13 @@ std::string RugnuxCommandLine(const ProcessConfig &config, std::ostringstream radii; radii << bragg.GetR1() << "," << bragg.GetR2() << "," << bragg.GetR3(); add("--integration-radius", radii.str()); - // Background trim defaults to 0.10 in the CLI, so emit it whenever the GUI value differs (a - // custom fraction, or 0 when the box is unchecked) to reproduce the GUI's choice faithfully. - if (bragg.GetBackgroundTrimFraction() != 0.10f) + // Background ring: the CLI defaults to the 4 sigma high-side clip, so emit a flag only when the + // GUI chose otherwise. The trim is the alternative estimator (it clears the clip), and a clip of + // 0 with no trim is the plain ring mean, which needs the flag to be reproduced. + if (bragg.GetBackgroundTrimFraction() > 0.0f) add("--background-trim", num(bragg.GetBackgroundTrimFraction())); + else if (bragg.GetBackgroundClipNSigma() != BraggIntegrationSettings().GetBackgroundClipNSigma()) + add("--background-clip", num(bragg.GetBackgroundClipNSigma())); if (const auto max_hkl = bragg.GetMaxHKL()) add("--max-hkl", std::to_string(*max_hkl)); // Unset means "to the detector edge"; emitting the resolved number would pin it to this run's diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 65417482..c27f4386 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -128,7 +128,8 @@ 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-trim Monochromatic (rotation + still): symmetric trimmed-mean fraction for the background ring (0<=f<0.5, default 0.10; 0 = plain mean). Removes the high-side bias that over-subtracts weak high-angle spots (broadband data keep the sigma-clip instead)" << 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-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 << " --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; @@ -171,6 +172,7 @@ enum { OPT_REDO_ROTATION_SPOTS, OPT_FORCE_ROTATION_LATTICE, OPT_ROTATION_NO_POSTREFINE, + OPT_BACKGROUND_CLIP, OPT_REFINE_GEOMETRY, OPT_BANDWIDTH, OPT_INTEGRATION_RADIUS, @@ -278,6 +280,7 @@ static option long_options[] = { {"search-min-zeta", required_argument, nullptr, OPT_SEARCH_MIN_ZETA}, {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, + {"background-clip", required_argument, nullptr, OPT_BACKGROUND_CLIP}, {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET}, {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS}, @@ -553,6 +556,7 @@ static int RunRugnux(int argc, char **argv) { std::optional search_min_zeta_arg; // --search-min-zeta; rotation default below int64_t scaling_iter = 3; std::optional forced_rotation_lattice; + std::optional background_clip_arg; // --background-clip: background-ring high-side sigma clip std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on @@ -878,6 +882,9 @@ static int RunRugnux(int argc, char **argv) { case OPT_INTEGRATION_HIGH_RES: integration_d_min_arg = parse_double_arg(optarg, "--integration-high-resolution", logger); break; + case OPT_BACKGROUND_CLIP: + background_clip_arg = parse_double_arg(optarg, "--background-clip", logger); + break; case OPT_INTEGRATOR: if (strcmp(optarg, "boxsum") == 0) integrator_mode = IntegratorMode::BoxSum; else if (strcmp(optarg, "gaussian") == 0) integrator_mode = IntegratorMode::ProfileGaussian; @@ -1599,8 +1606,15 @@ 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 trimmed-mean fraction set to {:.2f} (rotation; stills keep the sigma-clip)", - *background_trim_arg); + logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip " + "(monochromatic data; broadband always clips)", *background_trim_arg); + } + + if (background_clip_arg) { + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.BackgroundClipNSigma(static_cast(*background_clip_arg)); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg); } SpotFindingSettings spot_settings; diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index f38bf30d..78257b87 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -78,7 +78,10 @@ Scene BuildScene(size_t width, size_t height, int spacing = 60) { return s; } +// clip_nsigma 0 selects the OTHER background-ring estimator, the symmetric trim, so the two branches +// the CPU and GPU each implement separately are both covered. DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional bandwidth_fwhm, + float clip_nsigma = 4.0f, const DetectorSetup &det = DetJF(2)) { DiffractionExperiment experiment(det); // DetJF(2) (small) keeps the correctness test fast experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV) @@ -86,12 +89,17 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b experiment.BandwidthFWHM(bandwidth_fwhm); BraggIntegrationSettings settings; settings.Integrator(mode); + if (clip_nsigma > 0.0f) + settings.BackgroundClipNSigma(clip_nsigma); + else + settings.BackgroundTrimFraction(0.10f); experiment.ImportBraggIntegrationSettings(settings); return experiment; } -void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm) { - const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm); +void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, + float clip_nsigma = 4.0f) { + const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); @@ -145,6 +153,7 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { SECTION("ProfileGaussian mono") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt); } SECTION("ProfileGaussian broadband") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.03f); } SECTION("ProfileEmpirical") { CompareCpuVsGpu(IntegratorMode::ProfileEmpirical, std::nullopt); } + SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); } } // Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU @@ -154,7 +163,8 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { WARN("No CUDA GPU present. Skipping benchmark"); return; } - const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::ProfileGaussian, std::nullopt, DetJF4M()); + const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::ProfileGaussian, std::nullopt, + 4.0f, DetJF4M()); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index 10454670..aab113cc 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -520,16 +520,16 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() { auto *r3 = new NumberLineEdit(1.0f, 30.0f, bragg_.GetR3(), 1, "px", this); auto *radii = new QHBoxLayout(); radii->addWidget(r1); radii->addWidget(r2); radii->addWidget(r3); - // Background trim: replace the r2..r3 ring mean with a symmetric trimmed mean (drop the lowest and - // highest fraction of ring pixels), which removes the high-side bias that over-subtracts weak - // high-angle reflections. Checkbox + fraction; 0 = plain mean. + // Background trim: use the old symmetric trimmed mean for the r2..r3 ring (drop the lowest and + // highest fraction of ring pixels) in place of the default high-side sigma clip. Setting the trim + // clears the clip, so unchecked simply leaves the default estimator in place. const float trim = bragg_.GetBackgroundTrimFraction(); auto *bkgTrim = new QCheckBox("Background trim", this); bkgTrim->setChecked(trim > 0.0f); bkgTrim->setToolTip("Estimate the local Bragg background with a symmetric trimmed mean of the " - "background ring (drop the lowest and highest fraction of pixels) instead of the " - "plain mean, removing the high-side bias that over-subtracts weak high-angle " - "reflections. 0.10 recommended; unchecked = plain mean."); + "background ring instead of the default high-side sigma clip. A symmetric trim " + "is biased low on Poisson data and adds a few counts to every partial, so this " + "is for back compatibility; unchecked = the 4 sigma clip."); auto *bkgTrimFrac = new NumberLineEdit(0.01f, 0.49f, trim > 0.0f ? trim : 0.10f, 2, "", this); bkgTrimFrac->setEnabled(bkgTrim->isChecked()); auto *trimRow = new QHBoxLayout(); @@ -545,7 +545,11 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() { bragg_.Integrator(gaussian->isChecked() ? IntegratorMode::ProfileGaussian : IntegratorMode::BoxSum); bragg_.R1(static_cast(r1->value())).R2(static_cast(r2->value())) .R3(static_cast(r3->value())); - bragg_.BackgroundTrimFraction(bkgTrim->isChecked() ? static_cast(bkgTrimFrac->value()) : 0.0f); + // Either estimator clears the other, so unchecking has to put the default clip back explicitly. + if (bkgTrim->isChecked()) + bragg_.BackgroundTrimFraction(static_cast(bkgTrimFrac->value())); + else + bragg_.BackgroundClipNSigma(BraggIntegrationSettings().GetBackgroundClipNSigma()); emit braggChanged(bragg_); }; connect(gaussian, &QCheckBox::toggled, this, [emitBragg] { emitBragg(); });