From f0cdb027e11f17d426398c154936a6330cb93199 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 6 Aug 2026 19:06:30 +0200 Subject: [PATCH] Ice: default the merge mask off, gate the radial background on smooth ice, and pick detection by geometry Three defaults, each settled by measurement rather than by argument. The arbiter throughout is structure-referenced - anomalous peak height where a crystal can carry it, and otherwise the agreement of the ice bands with a fixed external model against resolution-matched DECOY bands carrying no ice. The band-versus-decoy contrast is used because R-free here tracks completeness, and every one of these switches moves completeness. The damage is real and it localizes: over the rotation battery the ice bands' excess amplitude reaches +9.6% on a smooth-ice crystal and +35% on the worst, while a clean control sits at +0.6% (z +0.45). On the worst crystal, nine of the ten largest excess peaks in a q scan land on hexagonal ring positions. Turning ice handling off leaves the contrast unchanged and forcing it on a clean crystal does not create one, so it is the ice and not the machinery. MERGE-TIME RING MASK -> OFF. It deletes reflections, which no other program does by default - AIMLESS, DIALS, xia2, XDS and CrystFEL all keep ice-band reflections in the merge and exclude them only from the model fit; autoPROC is the sole exception. On the one battery crystal where the mask fires and an anomalous arbiter can score it, dropping the band moved the mean peak height at the known sites by -0.001 +- 0.018 sigma, 2% of the site height, while removing 1149 unique reflections whose mean I/sigma was 3.62 against the dataset's own 3.05 - better than average data - and costing 17 completeness points in that shell. It fires on 5 of 37 crystals, changes no space group, and those 5 disagree in sign: it clearly helps the two most heavily iced, is a wash on two and costs a third. So it stays as a switch, worth setting by hand on a badly iced crystal where it shows in the high shell, but it is not a default. RADIAL BACKGROUND -> AUTO, gated per image. The correction models the background as a function of radius alone, and that is exactly when it works. On a crystal with pure smooth powder ice it removes 43% of the bands' excess amplitude, with the improvement 7x larger inside the bands than outside; on a crystal whose ice is discrete crystallite spots - no smooth ring to model - the excess amplitude GREW by half; on clean data it is inert to four decimals. The two ice channels already separate those morphologies, so --background-radial takes on|off|auto and auto applies it to an image when that image's peak-excluded score reaches --ice-min-score. Auto never engages without such a score, because the plain profile carries the Bragg peaks and cannot support an absolute threshold. Per image rather than per run, and that was tested rather than assumed: the gate fires on 100% and 94% of frames on the two crystals that want it, and on 1.5% of frames - 32 blocks, 23 of them single frames - on the textured-ice crystal. A seam statistic against off + f*(on - off) is null on both mixed runs, every merge statistic is bracketed by the pure arms, and the textured crystal's auto arm lands on `off` rather than on `on`'s harm. A run-level gate would need the score before the pass that integrates, i.e. rotation-only plumbing, and buys nothing measurable. The kernel was already built unconditionally, so flipping the flag per image is free - except on the GPU, where the launches were gated on a construction-time n_rad. That is why the buffers are now allocated whenever the correction could run, and Run() decides per image. DETECTION -> the geometry's default when the file is silent: on for rotation, off for stills, with the command line and then the file taking precedence. A rotation sweep sits on the same rings for the whole run, so ice there is a coherent systematic and the presence gate keeps it inert on a clean crystal; a serial stills run has too few spots per image to spend any on flagging. The master file's key is kept as written rather than collapsed to a bool, so "the file said nothing" is distinguishable from "the file said no" - it used to fall silently to off, taking the exclusion from the scale fit with it. Co-Authored-By: Claude Opus 5 (1M context) --- common/BraggIntegrationSettings.cpp | 4 +- common/BraggIntegrationSettings.h | 19 +++++++-- common/ScalingSettings.h | 22 ++++++++-- docs/CHANGELOG.md | 10 ++++- docs/CPU_DATA_ANALYSIS.md | 34 +++++++++++++-- docs/RUGNUX.md | 5 ++- image_analysis/MXAnalysisWithoutFPGA.cpp | 26 ++++++++++++ .../BraggIntegrationEngine.cpp | 6 ++- .../BraggIntegrationEngine.h | 9 ++++ .../BraggIntegrationEngineGPU.cu | 16 +++++--- reader/HDF5MetadataSource.cpp | 5 ++- reader/JFJochReaderDataset.h | 5 +++ rugnux/rugnux_cli.cpp | 41 ++++++++++++------- 13 files changed, 165 insertions(+), 37 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 7142f13c..e1b92958 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -140,11 +140,11 @@ float BraggIntegrationSettings::GetBackgroundClipNSigma() const { return bkg_clip_nsigma; } -BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(bool input) { +BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(std::optional input) { bkg_radial_correction = input; return *this; } -bool BraggIntegrationSettings::IsBackgroundRadialCorrection() const { +std::optional BraggIntegrationSettings::GetBackgroundRadialCorrection() const { return bkg_radial_correction; } diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index fa24c4ac..24a8f5d2 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -56,7 +56,19 @@ class BraggIntegrationSettings { // mean_annulus(B) - mean_disk(B), evaluated as a fixed kernel over radial offset (O(1), no extra // pixel reads). Measured empty-aperture bias over 9 bands on 3 crystals: 4.33 -> 0.79 counts mean // |bias|, scatter unchanged. - bool bkg_radial_correction = false; + // + // Unset means AUTO, which is the default: apply it per image where that image's peak-excluded ice + // score says a SMOOTH powder ring is present, and not otherwise. The correction models the + // background as a function of radius alone, so it helps exactly where that is true and not + // elsewhere. Measured against a fixed external model, band-versus-decoy-band: on a crystal with + // pure smooth ice it removes 43% of the ice bands' excess amplitude, with the effect 7x stronger + // inside the bands than outside; on a crystal whose ice is discrete crystallite SPOTS - no smooth + // radial ring to model - the excess amplitude instead GREW by half; on clean data it is inert to + // four decimal places. The ice score's two channels separate those two morphologies, so the + // correction is gated on the smooth one. Auto only engages where a peak-excluded score exists + // (adaptive spot finding); the plain profile carries the Bragg peaks and cannot support a + // threshold, so without it auto stays off. + std::optional bkg_radial_correction; // 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 @@ -77,7 +89,7 @@ public: BraggIntegrationSettings& Integrator(IntegratorMode input); BraggIntegrationSettings& BackgroundTrimFraction(float input); BraggIntegrationSettings& BackgroundClipNSigma(float input); - BraggIntegrationSettings& BackgroundRadialCorrection(bool input); + BraggIntegrationSettings& BackgroundRadialCorrection(std::optional input); BraggIntegrationSettings& MaxHKL(std::optional input); @@ -91,6 +103,7 @@ public: [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; [[nodiscard]] float GetBackgroundTrimFraction() const; [[nodiscard]] float GetBackgroundClipNSigma() const; - [[nodiscard]] bool IsBackgroundRadialCorrection() const; + // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction. + [[nodiscard]] std::optional GetBackgroundRadialCorrection() const; [[nodiscard]] std::optional GetMaxHKL() const; }; diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 15ab255f..cf206d60 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -65,9 +65,25 @@ class ScalingSettings { // Merge-time ice-ring mask (FindDecorrelatedIceRings): after a first merge, drop a hexagonal-ice ring // whose merged half-set CC1/2 has collapsed below its resolution shoulders and merge again. Separate // from --detect-ice-rings, which flags ice SPOTS for indexing and keeps ice reflections out of the - // SCALE fit - so this can be turned off on its own to de-confound a merge experiment without changing - // how the data were indexed and scaled. On by default; only consulted when detect_ice_rings is on. - bool ice_ring_merge_mask = true; + // SCALE fit - so this can be turned on or off on its own without changing how the data were indexed + // and scaled. Only consulted when detect_ice_rings is on. + // + // OFF by default. Deleting reflections is not what the field does - AIMLESS, DIALS, xia2, XDS and + // CrystFEL all keep ice-band reflections in the merge, and only autoPROC removes them - and the + // deletion did not pay for itself when it was measured against a structure-referenced metric. On the + // one crystal in the rotation battery where the mask both fires and the anomalous arbiter can score + // it, dropping the ring changed the anomalous peak height by -0.001 +- 0.018 sigma (2% of the mean + // site height), while costing 1149 unique reflections whose mean I/sigma was 3.62 against the + // dataset's own 3.05 - i.e. it deletes better-than-average data. Overall R_meas, CC1/2 and ISa were + // identical to three significant figures either way, and the affected shell went from 82.9% to + // 100.0% complete without it. + // + // It is not useless, which is why the switch stays: over the 37-crystal rotation battery it fires + // on 5, changes no space group, and those 5 disagree in sign - it clearly helps the two most + // heavily iced (one gains 3.8 R_meas and 4.0 CC1/2 points, the other 23 points of high-shell + // CC1/2), is a wash on two and costs a third. It always costs completeness where it fires. So: + // off as a default, worth turning on by hand on a badly iced crystal. + bool ice_ring_merge_mask = false; // Minimum measured ice strength (iceRingScore, 1 = no ice) before any ice-ring handling is applied // at all. The eleven fixed hexagonal bands cover 16-26% of the unique reflections at typical // resolutions REGARDLESS of whether the crystal has ice, so flagging unconditionally taxes clean diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 95eebdfb..91d20c5d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,10 +8,18 @@ This is an UNSTABLE release. It includes many experimental features, as well as * 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. +* Bragg integration: The **uncertainty of the background estimate** is now propagated into `sigma`; both engines omitted it, which understated every reflection's `sigma` by `sqrt(1 + n_signal/n_bkg)` = 1.109 with the shipped stencil. Expect `` to fall by that factor on every dataset. +* Bragg integration: New **radial background correction** (`--background-radial[=on|off|auto]`, default `auto`) for the bias a curved radial background leaves in a flat ring mean - tens of counts on a single reflection sitting on an ice ring. `auto` applies it per image where that image's ice score shows a *smooth* powder ring: measured against a fixed atomic model it removes 43 % of the ice bands' excess amplitude on smooth ice, but on ice made of discrete crystallite spots - where a radius-only background model has nothing to fit - it makes the bias worse, and the gate separates the two. +* rugnux: **Ice-ring handling now runs only when the crystal is measured to have ice** - the eleven fixed bands cost 16-26 % of the unique reflections whether it does or not. Detection uses 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), which catch smooth and textured ice respectively; the previous score was read off the plain azimuthal profile, where strong reflections in a ring's bin score as ice. +* rugnux: The merge-time **ice-ring mask is now off by default** (`--ice-ring-mask[=on|off]`, its own switch rather than sharing `--detect-ice-rings` with indexing). Measured against anomalous peak height it returned -0.001 +- 0.018 sigma while deleting reflections whose mean `I/sigma` was *above* the dataset average, leaving `R_meas`, CC1/2 and ISa unchanged and the affected shell 17 completeness points short. Ice reflections are still kept out of the scale fit and the space-group search, as every established scaling program does. +* rugnux: **Multiplicity was over-reported** whenever the ice-ring mask dropped a band: the per-shell observation counter rode on the R_meas re-walk, which deliberately ignores the mask, so observations the merge had excluded were counted against a unique count that excluded them - and on the GPU path a fully-masked group's resolution is NaN, which `ResolutionShells::GetShell` silently binned as the lowest-resolution shell rather than rejecting. Only counts were affected (intensities, sigmas, R_meas, CC1/2, completeness and ISa were always right), and only on that path. * 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 space-group search no longer **starves on a low-ISa merge**: its fixed `I/sigma >= 3` cut could select nothing at all, since a merged sigma is floored so that no reflection reads above `ISa = 1/b`, leaving every operator correlation undefined and the point group at 1. The cut is now capped at the merge's own quantile, and is unchanged on healthy merges. +* rugnux: The one-line summary **names every space group the data cannot separate** (`I23 (No. 197) or I2(1)3 (No. 199) - indistinguishable from these data`), instead of reporting only the representative as if it had been measured. * rugnux: The per-image mosaicity is fitted from the strongest 250 spots, so the indexing spot budget no longer sets it. -* rugnux: The per-frame crystal orientation is smoothed before scaling and partiality recomputed from it, so per-image refinement noise no longer reaches the merged intensities. +* rugnux: The per-frame crystal orientation is smoothed before scaling and partiality recomputed from it, so per-image refinement noise no longer reaches the merged intensities. Frames that never indexed are now excluded from that smoothing - their all-zero lattice passed the validity check and was averaged into their neighbours', which on a partly-indexing crystal chose a smoothing window nine times too wide. +* rugnux: Rotation post-refinement **reports the goniometer rotation scale** and warns beyond 0.5 %. A stage that turns further than commanded is invisible in the file, since the stored angles are the commanded ones, and reads as crystal drift instead. Report only - nothing is corrected. * rugnux: Stills **partiality post-refinement** added, on by default (`--simple-stills` disables); several non-helping stills scaling/detection knobs removed. * rugnux: **Scaling** hardened against a collapsed per-frame scale on the stills path as well as rotation, and the merged-sigma systematic floor is kept when ISa is too degenerate to report. * rugnux: `--min-image-cc` now works for rotation data (opt-in); new `--search-min-zeta` (rotation default 0.85); reports how close a symmetry axis lies to the spindle. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 4651331a..8db0f1e7 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -262,7 +262,9 @@ A single per-image **ice-ring score** is derived from a radial profile: for each The profile the score is read off is the **peak-excluded** one, not the plain azimuthal integration: where adaptive spot finding runs (§3.2 — the offline and viewer default), the score uses the sigma-clipped per-resolution-ring background that finder already computes for its threshold. This matters more than it sounds. A plain azimuthal profile is a per-ring *mean*, so a few strong low-resolution reflections landing in a ring's bin raise it exactly as ice would; measured over 37 rotation crystals that alone put ice-free crystals at scores of 1.5–4.2, above crystals that really are iced, and the strongest apparent "ice" in the set was a crystal with none. An ice ring is azimuthally smooth and survives the sigma clip, while Bragg peaks do not, so on the clipped profile the same 30 ice-free crystals sit at 1.00–1.22 and the three with confirmed ice at 2.08–2.37. Only where no adaptive finder ran (the FPGA workflow) does the score fall back to the plain profile. -The score is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, and offline it also **gates** ice handling: `--ice-min-score` (default 1.5) is the score a run must reach before ice-ring flagging, the exclusion from the scale fit and the merge-time ice mask are applied at all. The eleven fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. Which rings are then dropped from the merge remains data-driven, from the per-ring merged CC1/2. +The radial profile sees ice only as a **smooth powder ring**. Ice in large crystallites diffracts as discrete spots, leaves the profile flat, and is invisible to the score above, so a second channel is read from the spot list itself: the spots found on the ice bands are counted against the spots found in the ice-free flanks $[w,2w)$ either side of each band, rescaled to the bands' own $q$ width (a flank landing on another ring is dropped with its width). The indicator is the ratio pooled over the run — per image the control is a handful of spots and the ratio means nothing — and it is taken before the spot-count filter, which orders ice spots last and would discard them first. The two channels barely overlap: over the rotation battery the crystals with smooth ice read 2.1–2.4 on the profile and ~1.0 on the spots, the crystals with textured ice ~1.1 and 3.8–17.6, and a clean crystal 1.04 on both. Both counts are stored per image (`spot_count_ice_rings`, `spot_count_ice_control`; HDF5 `/entry/MX/peakCountIceRingRes`, `/entry/MX/peakCountIceRingControl`). + +Both channels are used offline as a **gate** on ice handling: unless the run reaches `--ice-min-score` (default 1.5) on the profile or `--ice-min-spot-ratio` (default 2.0) on the spots — 0 disables a channel — ice-ring flagging, the exclusion from the scale fit and the merge-time ice mask (§10.10) are all skipped. The eleven fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. Which rings are then dropped from the merge remains data-driven, from the per-ring merged CC1/2. A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low. @@ -462,6 +464,8 @@ The refinement above (§7.2) runs per image against that image's spots. For rota The refined pass is written as the canonical `_*` output; the pass-1 (header-geometry) result is kept alongside as `_01_*` for comparison. +**Goniometer rotation scale (report only).** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. Step A already measures it without a new degree of freedom: its residual rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under the same cross-validation that gates the cell move — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it **under-reads** the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error. + --- ## 8. Reflection prediction @@ -554,7 +558,7 @@ $ \hat{b} = \frac{B}{n_B},\qquad \hat{I} = S - n_S \hat{b}, $ -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. +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. @@ -562,6 +566,14 @@ The clip cuts only the high tail, which matters: the **symmetric** trimmed mean 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. +**Radial background correction (opt-in).** A ring mean estimates the background *under* the signal disk correctly only if the background is flat there. The signal disk and the ring are concentric, so for a background that is **linear** in position $\langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk}$ identically — a plane or gradient fit buys exactly nothing. The leading error is the **curvature** of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. Because every reflection uses the same stencil, that error is a fixed kernel over radial offset, +$ +\delta \hat b \;=\; \textstyle\sum_k \kappa_k\, \bar B(r_0+k), +$ +with $\kappa$ obtained once by azimuth-averaging the stencil and $\bar B(r)$ the image's own radial background curve. Applying it costs one short dot product per reflection and no extra pixel reads; correcting the background *scalar* means the box sum, the profile fit and the variance all pick it up. The curve is accumulated from the same annulus pixels the background pass already reads (a pixel's radius is the reflection's radius plus the pixel's projection on the beam→reflection direction, so no per-pixel square root is needed) and specifically from the **clipped** pixels, or it would carry neighbour tails and zingers — which is why the correction is inert under `--integrator boxsum`, that path having no clip pass. + +The model is a function of **radius alone**, so it is applied only where that is true of the background. `--background-radial` takes `on`, `off` or `auto`, and **`auto` is the default**: each image's peak-excluded ice score (§3.3) is taken after spot detection and before integration, and the correction is applied to that image when the score reaches the same `--ice-min-score` gate. Smooth powder ice *is* a radial feature and is corrected; ice made of discrete crystallite spots — which the profile channel is blind to and the spot channel catches — leaves no smooth ring to model, and correcting it makes matters worse. Measured against a fixed atomic model, comparing ice bands with resolution-matched decoy bands carrying no ice: on a crystal with pure smooth ice the correction removes **43 % of the bands' excess amplitude**, and the improvement is **7× larger inside the bands than outside**, which is its stated mechanism; on a crystal whose ice is textured the same correction *increased* the excess amplitude by half; on a clean crystal it is inert to four decimal places. Auto engages only where a peak-excluded score exists (adaptive spot finding, §3.2) — a plain azimuthal profile carries the Bragg peaks and cannot support an absolute threshold, so without one auto leaves the correction off. + ### 9.3 Profile-fitted extraction (default) 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: @@ -634,7 +646,7 @@ Reflections below a minimum partiality can be rejected from merging to avoid uns The per-frame scales $G_i$ are fit by robust (Cauchy) inverse-variance-weighted ratios; there is no explicit $G\approx1$ prior. For rotation datasets, optional smoothing enforces the expectation that scale and mosaicity vary slowly across a sweep: **after** the per-frame fit, $\log G_i$ (and the mosaicity) are replaced by a centred **moving average** over a window spanning a configurable rotation range (XDS DELPHI-like; `--smooth-g`, default 5° for rot3d, off otherwise). It is a post-fit smoothing pass, not a curvature penalty inside the least-squares objective. -The **crystal orientation** is smoothed the same way, and for the same reason. Geometry is re-refined independently on every frame against that frame's spots alone — as few as a dozen on a sparse crystal — so the per-frame orientation carries a real slow drift (crystal slippage, up to ~1.3° across a sweep) on top of fit noise that scales with spots per frame. Before scaling, the per-frame lattices are de-rotated to a common reference, averaged in frame order, rotated back, and every partial's $\Delta\phi$ — hence its partiality — is recomputed from the smoothed lattice. The window is chosen per dataset by leave-one-out cross-validation (does a frame's neighbours predict its orientation?) rather than fixed, because drift and noise both vary by two orders of magnitude between crystals; it is capped, because the per-frame fit also absorbs a real per-frame systematic that smoothing too wide destroys. Refining *less* is not an alternative: with per-image refinement off the space group is lost on several crystals. +The **crystal orientation** is smoothed the same way, and for the same reason. Geometry is re-refined independently on every frame against that frame's spots alone — as few as a dozen on a sparse crystal — so the per-frame orientation carries a real slow drift (crystal slippage, up to ~1.3° across a sweep) on top of fit noise that scales with spots per frame. Before scaling, the per-frame lattices are de-rotated to a common reference, averaged in frame order, rotated back, and every partial's $\Delta\phi$ — hence its partiality — is recomputed from the smoothed lattice. The window is chosen per dataset by leave-one-out cross-validation (does a frame's neighbours predict its orientation?) rather than fixed, because drift and noise both vary by two orders of magnitude between crystals; it is capped, because the per-frame fit also absorbs a real per-frame systematic that smoothing too wide destroys. Only frames that actually indexed take part: a frame that did not carries an all-zero lattice, which is *finite* and so passes a validity check written as a finite test, and would otherwise be both averaged into its neighbours' orientation and scored in the cross-validation that picks the window. Refining *less* is not an alternative: with per-image refinement off the space group is lost on several crystals. A per-frame scale enters every intensity as $1/G$, so a frame whose fit is not determined by its data can amplify it without bound — and $\sigma$ is amplified by the same factor, which makes it invisible to any $\sigma$-based outlier test. A fitted $G$ far below the run's median is therefore treated as *undetermined* rather than as a successful fit, both here and in the separate refit on the combined fulls (§10.6). The bound is a ratio to the run's own median because $G$ is not gauge-fixed: it and the merged means have an exact global multiplicative degeneracy, so no absolute value is meaningful. @@ -733,6 +745,16 @@ A reference dataset (`--reference-mtz`) supplies known intensities for the same **Resolve the indexing (merohedral) ambiguity.** When the lattice symmetry is higher than the crystal's Laue symmetry (e.g. $P3$, $P4$, $P6$, $C2$), more than one indexing of the same lattice is geometrically valid, and the two solutions produce *different* merged intensities that a self-consistent scale cannot tell apart — only an external reference can. The candidate reindexings are the identity together with the twin-law cosets of the metric symmetry (from the unit-cell metric and the Laue group); each is scored by the intensity correlation $\mathrm{CC}_\mathrm{ref}$ of the reindexed merge against the reference, and the data are re-merged in the best-correlating indexing. The reindex is **metric-preserving** — only the $hkl$ labels change, the cell is unchanged — and it is a no-op for a holohedral crystal, which has no twin laws (the lattice and Laue symmetry coincide). For rotation data this is done once, after the space group is determined; the reference is *not* used to scale the rotation merge, which stays self-consistent (its $\mathrm{ISa}$ comes from the data alone). For stills the reference is the per-image scale target of the on-the-fly scaling (§10.2). +### 10.10 Ice rings at the scale and merge stages + +Where the gate of §3.3 has found ice, reflections falling within $\pm w$ in $q$ of a hexagonal-ice band ($w=0.03$ Å$^{-1}$ offline, about the measured ring half-width) are marked. Marked reflections are **excluded where a model is fitted** — the per-frame scale $G$, the per-image correlation, and the $P1$ merge the space-group search runs on — because ice contamination is a *positive bias*, not extra scatter, and a least-squares scale absorbs it into $G$ and into the error-model $b$, where it damages every other reflection on the same frame. They are **kept in the final merge**, which is also what the established scaling programs do by default, so the affected shells keep their completeness. + +Nothing is deleted by default. An **optional** merge-time mask (`--ice-ring-mask[=on|off]`, default **off**, switched separately from ice detection because detection also decides how the data are indexed) will drop a band that can be shown to be dead from the merged data rather than from the presence of ice: after the final merge the half-set $\mathrm{CC}_{1/2}$ is computed inside each ice band and in the shoulders $[w,3w)$ either side of it (shoulder reflections lying on a *neighbouring* band are excluded — the three bands near 1.9 Å are only 0.05–0.06 Å$^{-1}$ apart in $q$, so an unguarded shoulder compares ice against ice), and a band falling more than a margin below its shoulders is dropped and the data re-merged. The margin is 0.10, the 99th percentile of a null measured directly — the identical statistic evaluated at **decoy bands**, $q$ positions carrying no ice ring at all — because at these populations a fixed CC margin is not a fixed significance level, and the parametric (Fisher-$z$) error understates the real scatter of heavy-tailed merged intensities several-fold. + +It is off by default because the deletion did not pay for itself when it was measured against an external reference rather than against the merge's own statistics. On the one rotation-battery crystal where the mask fires *and* an anomalous-peak-height arbiter can score the result, dropping the band changed the mean anomalous density at the known sites by $-0.001\pm0.018\,\sigma$ — about 2 % of the site height — while removing 1149 unique reflections whose mean $I/\sigma$ was 3.62 against the dataset's own 3.05, i.e. better-than-average data. This also matches the field: established scaling programs keep ice-band reflections in the merged output and exclude them only from the model fit. + +It is a switch rather than dead code because it does help some crystals. Over a 37-crystal rotation battery it fires on 5 and changes no space group; of those 5 it clearly helps the two most heavily iced (3.8 points of $R_\mathrm{meas}$ and 4.0 of $\mathrm{CC}_{1/2}$ on one, 23 points of *high-shell* $\mathrm{CC}_{1/2}$ on the other), is a wash on two, and costs a third. It always costs completeness where it fires — 6 to 8 points. Turning it on by hand is reasonable for a badly iced crystal; making it a default is not. + --- ## 11. Mosaicity and “profile radius” monitoring @@ -783,7 +805,7 @@ A **dataset-wide** Wilson $B$ is also estimated over the merged reflections — - **Bragg integration is profile-fitted by default** (per-shell Gaussian profile, Kabsch extraction; §9.3), with plain box summation available as a fallback (`--integrator boxsum`). The profiles are built per frame from that frame's strong spots, which suits fast-feedback and serial/streaming use; a profile shared across many frames (as in full offline workflows) is not currently formed. - **Space-group symmetry** beyond centering absences is not necessarily enforced during prediction/integration unless the space group is supplied and used downstream. -- **Resolution masking and ice rings** are controllable; including ice-ring spots in indexing can improve robustness for some samples but may bias refinement in others. +- **Resolution masking** is controllable, and so is every stage of ice-ring handling (§3.3, §10.10). None of it runs unless the crystal is measured to have ice, because the fixed bands are a fixed cost in unique reflections whether it does or not. - **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted within an excitation-error window and get their partiality from the default-on per-crystal tilt post-refinement (§10.2) — or unit partiality with `--simple-stills`. - **Space-group determination.** When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities. Three tests gate a promotion to higher symmetry, all aimed at the merohedral twin, whose twin law forces non-equivalent reflections together and so mimics symmetry: @@ -791,6 +813,10 @@ A **dataset-wide** Wilson $B$ is also estimated over the merged reflections — 2. **Error-model $b$** (the intensity-proportional systematic). A genuine symmetry step gains multiplicity without inflating $b$; merging a twin law's extra operator inflates it. A $\chi^2$-passing promotion is vetoed when $b$ rises past a bound relative to the confirmed subgroup. 3. **Operator disagreement**, a sigma-free statistic $H=\mathrm{median}\,|I_1-I_2|/(I_1+I_2)$, formed as the ratio of the operators a promotion *adds* to the parent's own, measured on the same reflections. Normalising against the parent divides out the systematic floor that symmetry mates carry on real data, which varies by crystal and by operator; a median is used because a twin perturbs every pair whereas a badly-measured minority perturbs only the tail. Where a candidate has several parents of the same order, it is judged against the worst of them, since a rival subgroup can itself contain the twin laws. + The operator correlations are taken on reflections above an $I/\sigma$ cut, and that cut is **capped at the merge's own $I/\sigma$ quantile** rather than applied as a fixed number. A merged $\sigma$ carries the systematic floor $b|I|$ (§10.4), so no reflection in a merge can read above $\mathrm{ISa}=1/b$; a fixed cut of 3 therefore selects nothing at all on a search merge whose ISa is below 3, leaving every operator correlation undefined and collapsing the point group to 1. The cap keeps at least the strongest quarter and is inert — the cut stays exactly 3.0 — on a healthy merge. + + Several space groups may share an absence pattern exactly. Where they do, the search scores them identically and **all of them are named** in the result rather than one being reported as the answer: some are enantiomorph pairs, which merged intensities cannot distinguish in principle, and others differ only by a screw condition that the centering condition already implies, so the screw has no observable signature at all. The representative reported first is the lowest space-group number, which is a convention and not a measurement. + The Lorentz factor $\zeta$ (§8.3) governs how well a reflection can be measured, so when the spindle lies in a plane of the lattice, an operator permuting the two in-plane axes samples a different mixture of measurement qualities than one that only flips signs. The search is therefore run a second time on a merge of only the well-measured observations (`--search-min-zeta`, rotation default 0.85), and **whichever search found more symmetry is kept** — one-way safe, because discarding observations can starve an operator correlation but never invent one. The filtered merge decides the point group only; absences come from the full merge, since they live in the weak reflections the filter removes. A tie (same order, different symmetry) is reported with both candidates named, for trying in molecular replacement. **Centering** is accepted when the systematically-absent class is weak relative to the present one by *either* of two floor-independent tests: its mean signed $I/\sigma$ well below the present mean, *or* its rate of individually-significant reflections well below the present class's own significant rate. The second test covers weak and low-energy data, where a positive intensity floor (background and profile leakage) lifts the absent class's mean $I/\sigma$ well above zero and, when the present class is itself weak, carries the plain mean ratio past its bound; a false centering fails both tests, its absent class being as strong as the present one. When several centerings pass, they are ranked by their **net** systematic absences (absent minus violating), not the gross absent count, so a super-centering (e.g. $F$ over a true $C$) whose extra, only-half-populated absent class dilutes the strength ratio does not out-rank the correct lower centering. diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 7bec03f2..3df79800 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -211,7 +211,7 @@ Spot finding: | `--spot-low-resolution ` | Low-resolution limit for spot finding, Å (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weak serial data) | | `--min-pix-per-spot ` | Minimum connected strong pixels per spot. **If omitted, min-pix is chosen per image** (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed-spot count × indexed fraction is kept. Give an explicit value to force a fixed min-pix instead. | | `--max-spots ` | Maximum spots kept per image (the strongest ones) and handed to indexing (default: 1000) | -| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling; overrides the dataset/master-file setting (default: use the dataset value). The merge-time mask is separate, see `--ice-ring-mask` | +| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's `detect_ice_rings`, or — where the file carries no such key — **on for rotation and off for stills**. The merge-time mask is separate, see `--ice-ring-mask` | Azimuthal integration (the radial profile behind the per-image ice-ring score): @@ -268,7 +268,7 @@ Scaling and merging: | `--resolution-cc-target ` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) | | `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 10) | | `--min-partiality ` | Minimum partiality to accept a reflection (default: 0.02) | -| `--ice-ring-mask[=on\|off]` | Merge-time ice-ring mask: after a first merge, drop a hexagonal-ice ring whose merged half-set CC1/2 has collapsed below its resolution shoulders and merge again (default: on; only acts when `--detect-ice-rings` is on). `off` disables **only** this mask — ice-spot flagging and the ice exclusion from scaling stay as `--detect-ice-rings` set them | +| `--ice-ring-mask[=on\|off]` | Merge-time ice-ring mask: after a first merge, drop a hexagonal-ice ring whose merged half-set CC1/2 has collapsed below its resolution shoulders and merge again (default: **off**; only acts when `--detect-ice-rings` is on). This is **only** the merge-time mask — ice-spot flagging and the ice exclusion from scaling stay as `--detect-ice-rings` set them. It is off because deleting the band was measured not to pay for itself: it costs completeness in the affected shell and removes reflections whose mean I/σ is *above* the dataset average, while leaving the anomalous signal, R_meas, CC½ and ISa unchanged | | `--ice-min-score ` | Ice-presence gate: the measured per-run ice score (1 = no ice) a dataset must reach before **any** ice handling is applied — the flagging, the exclusion from scaling and the merge-time mask (default: 1.5; 0 = no gate). The eleven fixed hexagonal bands cover 16–26 % of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal only costs completeness | | `--ice-min-spot-ratio ` | The second ice-presence channel: found **spots** on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so `--ice-min-score` alone is blind to it (default: 2.0; 0 disables this channel) | | `--reject-outliers ` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) | @@ -289,6 +289,7 @@ Integration: | `--integration-radius ` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` | | `--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 | +| `--background-radial[=on\|off\|auto]` | Correct the background ring for the **curvature** of the radial background (default `auto`). Disk and 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, since the model is a function of radius alone: on ice made of discrete crystallite spots there is no smooth ring and the correction makes the bias worse. Ignored by `--integrator boxsum` (no clip pass to take the curve from) | | `--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/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 40b81d1b..8445fd4a 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -162,6 +162,30 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); }; + // The radial background correction has to be decided BEFORE this image is integrated, so the + // ice score is taken here rather than with the other per-image quantities at the end of the + // function. It needs the peak-excluded per-ring background, which the adaptive finder has as + // soon as it has detected - so this must be called after detection and before integration. + // Where no such background exists (no adaptive finder), auto leaves the correction off: the + // plain profile carries the Bragg peaks and cannot support an absolute threshold. + const auto decide_radial_background = [this, &spot_finding_settings, &output]() { + if (!bragg_engine->IsBackgroundRadialAuto()) + return; + // Same condition as the score at the end of this function: the adaptive finder holds its + // ring background from whenever it last ran, so requiring that it ran for THIS image is + // what keeps a stale curve out. + if (!spot_finding_settings.adaptive_threshold) + return; + const std::vector &ring_bkg = adaptiveSpotFinder->GetRingBackground(); + if (ring_bkg.empty()) + return; + output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore( + ring_bkg, integration.GetQBinCount(), integration.Settings(), + spot_finding_settings.ice_ring_width_Q_recipA); + bragg_engine->BackgroundRadial(*output.ice_ring_score + >= experiment.GetScalingSettings().GetIceMinScore()); + }; + // A missing min-pix (std::nullopt) means "choose it per image". This applies only to the stills // indexing path (each frame is indexed independently); rotation indexing builds one lattice from // all frames, so it keeps the fixed min-pix and the single-pass finder. @@ -211,6 +235,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // Index and integrate the winning spot list; no spot finding left to do. s.min_pix_per_spot = best_mp; SpotAnalyze(experiment, s, best_spots, output); + decide_radial_background(); indexer.ProcessImage(output, s, *prediction, integrate_fn); indexing_time_s += output.indexing_time_s.value_or(0.0f); } @@ -222,6 +247,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const std::vector spots = finder.Run(*preprocessor_buffer, spot_finding_settings); SpotAnalyze(experiment, spot_finding_settings, spots, output); output.spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - spot_finding_start_time).count(); + decide_radial_background(); if (spot_finding_settings.indexing) indexer.ProcessImage(output, spot_finding_settings, *prediction, integrate_fn); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index b409549b..1a72a439 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -83,7 +83,11 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // exact to the extent the stencil is small against the reflection's radius (r3 = 10 px vs // hundreds). k_diff is the annulus histogram minus the disk histogram, each normalised, so // dot(k_diff, B) is directly mean_annulus(B) - mean_disk(B). - bkg_radial = settings.IsBackgroundRadialCorrection(); + // Unset = auto: start off, and let the analysis raise it per image where the ice score says the + // background really is radial. An engine nobody drives therefore never applies the correction. + const auto radial = settings.GetBackgroundRadialCorrection(); + bkg_radial_auto = !radial.has_value(); + bkg_radial = radial.value_or(false); k_off = static_cast(std::ceil(r3)) + 1; k_diff.assign(2 * k_off + 1, 0.0f); { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 4db28f57..f8cbd9fa 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -108,7 +108,11 @@ protected: // stencil, so mean_annulus(B) - mean_disk(B) of a radial B is a FIXED kernel over radial offset: // bkg_error = sum_k k_diff[k] * B(r0 + k - k_off) // That is one short dot product per reflection and reads no pixels. Built in the constructor. + // The kernel is built whatever the setting says, so bkg_radial can be flipped between images at + // no cost - which is what the auto mode does, applying the correction only to the images whose + // background really is a smooth function of radius (see BackgroundRadial below). bool bkg_radial = false; + bool bkg_radial_auto = false; // settings left it unset: decide per image from the ice score int k_off = 0; // index of offset 0 in k_diff std::vector k_diff; // annulus-minus-disk weight per integer radial offset @@ -129,4 +133,9 @@ public: virtual std::vector Run(const ImagePreprocessorBuffer &image, const std::vector &predicted, size_t npredicted, int64_t image_number) = 0; + + // Turn the radial background correction on or off for the images that follow. The caller owns + // the decision; in the auto mode the analysis sets it per image from that image's ice score. + void BackgroundRadial(bool on) { bkg_radial = on; } + [[nodiscard]] bool IsBackgroundRadialAuto() const { return bkg_radial_auto; } }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 71f2ec5b..9e01139d 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -516,7 +516,9 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment "BraggIntegrationEngineGPU: profile grid exceeds shared memory (r2 too large)"); // Radial background curve: one bin per pixel of distance from the beam, out to the far corner. - if (bkg_radial) { + // Allocated whenever the correction COULD run, so the auto mode can turn it on for an individual + // image; whether it runs for a given image is decided in Run() from the per-image bkg_radial. + if (bkg_radial || bkg_radial_auto) { const double fx = std::max(beam_x, static_cast(xpixel) - beam_x); const double fy = std::max(beam_y, static_cast(ypixel) - beam_y); n_rad = static_cast(std::ceil(std::hypot(fx, fy))) + 2; @@ -598,9 +600,13 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .bkg_trim = bkg_trim, // effective trim fraction (0 for stills), set by the base ctor from settings }; + // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under + // the auto mode they are allocated for every image and used for the ones the ice score selects. + const int rad_n = bkg_radial ? n_rad : 0; + // Pass A: reset accumulators, mask, then box-sum. cuda_err(cudaMemsetAsync(d_mask, 0, npixel, *stream)); - if (n_rad > 0) { + if (rad_n > 0) { cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(float) * n_rad, *stream)); cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); } @@ -610,13 +616,13 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y, d_ok, d_strong, d_has_obs, d_invd2, d_isum, d_ninner, d_rbin, - d_rad_sum, d_rad_cnt, n_rad); + d_rad_sum, d_rad_cnt, rad_n); // Correct the flat annulus background for the curvature of the radial background before anything // downstream (profile fit, variance) reads it. - if (n_rad > 0) + if (rad_n > 0) radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>( - d_rad_sum, d_rad_cnt, n_rad, d_k_diff, static_cast(k_diff.size()), k_off, + d_rad_sum, d_rad_cnt, rad_n, d_k_diff, static_cast(k_diff.size()), k_off, d_isum, d_ninner, d_rbin, d_ok, d_bkg, d_I, n); if (mode != IntegratorMode::BoxSum) { diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index c4205dd4..f119efb7 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -547,8 +547,9 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen auto ring_current_A = master_file->GetOptFloat("/entry/source/current"); if (ring_current_A) dataset->experiment.RingCurrent_mA(ring_current_A.value() * 1000.0); - dataset->experiment.DetectIceRings( - master_file->GetOptBool("/entry/instrument/detector/detectorSpecific/detect_ice_rings").value_or(false)); + dataset->file_detect_ice_rings = + master_file->GetOptBool("/entry/instrument/detector/detectorSpecific/detect_ice_rings"); + dataset->experiment.DetectIceRings(dataset->file_detect_ice_rings.value_or(false)); dataset->experiment.PoniRot1_rad( master_file->GetOptFloat("/entry/instrument/detector/transformations/rot1").value_or(0.0)); dataset->experiment.PoniRot2_rad( diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 3f2f0427..de2ebf49 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -25,6 +25,11 @@ struct JFJochReaderDataset { std::optional error_value; + // The master file's detect_ice_rings key, kept as written rather than collapsed to a bool, so a + // consumer can tell "the file asked for this" from "the file said nothing". Offline processing + // defaults ice handling by geometry when the file is silent (see rugnux_cli). + std::optional file_detect_ice_rings; + std::string jfjoch_release; std::vector az_int_bin_to_q; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index aa115b6e..5abd5ab3 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -78,7 +78,7 @@ void print_usage() { std::cout << " --spot-high-resolution High resolution limit for spot finding. If omitted (or 0), spot finding is not clipped in resolution and extends as far as the detector reaches" << std::endl; std::cout << " --spot-low-resolution Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data)" << std::endl; std::cout << " --max-spots Max spot count per image, the strongest ones, handed to indexing (default: 1000)" << std::endl; - std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling; overrides the dataset/master-file setting (default: use dataset value). The merge-time ice-ring mask is separate, see --ice-ring-mask" << std::endl; + std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's setting, or - where the file says nothing - on for rotation and off for stills. The merge-time ice-ring mask is separate, see --ice-ring-mask" << std::endl; std::cout << std::endl; std::cout << " Indexing" << std::endl; @@ -112,7 +112,7 @@ void print_usage() { std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl; std::cout << " --ice-min-score Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging, the exclusion from scaling and the merge-time mask (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl; std::cout << " --ice-min-spot-ratio Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl; - std::cout << " --ice-ring-mask[=on|off] Drop a hexagonal-ice ring from the merge when its merged half-set CC1/2 has collapsed below its resolution shoulders, then re-merge (default: on; needs --detect-ice-rings on). Off disables ONLY this merge-time mask - ice-spot flagging and the ice exclusion from scaling stay as --detect-ice-rings sets them" << std::endl; + std::cout << " --ice-ring-mask[=on|off] Drop a hexagonal-ice ring from the merge when its merged half-set CC1/2 has collapsed below its resolution shoulders, then re-merge (default: OFF - measured to delete better-than-average reflections for no gain; needs --detect-ice-rings on). This is ONLY the merge-time mask - ice-spot flagging and the ice exclusion from scaling stay as --detect-ice-rings sets them" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; @@ -132,7 +132,7 @@ void print_usage() { 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-radial[=on|off] 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 an ice ring reaches +26 counts on a single reflection. Costs one short dot product per reflection and no extra pixel reads" << std::endl; + std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). 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 << " --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; @@ -536,7 +536,7 @@ static int RunRugnux(int argc, char **argv) { std::optional scale_fulls_arg; // --scale-fulls / --no-scale-fulls; default on for rot3d bool write_process_h5_flag = false; // --write-process-h5; also write _process.h5 when merging std::optional detect_ice_rings; // --detect-ice-rings[=on|off]; unset => use the dataset (file) value - bool ice_ring_mask = true; // --ice-ring-mask[=on|off]; merge-time CC1/2 ice-ring mask + bool ice_ring_mask = false; // --ice-ring-mask[=on|off]; merge-time CC1/2 ice-ring mask std::optional ice_min_score_arg; // --ice-min-score: ice-presence gate on the measured score std::optional ice_min_spot_ratio_arg; // --ice-min-spot-ratio: the same gate on the spot channel std::optional min_q, max_q, q_spacing; // azimuthal integration range / -q spacing (1/A) @@ -572,7 +572,8 @@ static int RunRugnux(int argc, char **argv) { 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 background_radial_arg; // --background-radial: radial curvature correction + bool background_radial_given = false; // --background-radial seen at all (unset => auto) + std::optional background_radial_arg; // when given: set = force on/off, unset = auto 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 @@ -916,12 +917,15 @@ static int RunRugnux(int argc, char **argv) { integration_d_min_arg = parse_double_arg(optarg, "--integration-high-resolution", logger); break; case OPT_BACKGROUND_RADIAL: + background_radial_given = true; if (optarg == nullptr || strcmp(optarg, "on") == 0) background_radial_arg = true; else if (strcmp(optarg, "off") == 0) background_radial_arg = false; + else if (strcmp(optarg, "auto") == 0) + background_radial_arg = std::nullopt; else { - logger.Error("Invalid --background-radial value: {} (expected on|off)", optarg); + logger.Error("Invalid --background-radial value: {} (expected on|off|auto)", optarg); return 1; } break; @@ -1171,14 +1175,17 @@ static int RunRugnux(int argc, char **argv) { experiment.SpaceGroupNumber(space_group_number); if (fixed_reference_unit_cell.has_value()) experiment.SetUnitCell(fixed_reference_unit_cell); - // --detect-ice-rings, applied here as well as on the full path below: this block returns - // before that one runs, so without it the flag is silently ignored by --scale. - if (detect_ice_rings.has_value()) - experiment.DetectIceRings(detect_ice_rings.value()); - // A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling. IndexingSettings indexing_settings; indexing_settings.RotationIndexing(experiment.GetGoniometer().has_value() && !force_still); + + // --detect-ice-rings, applied here as well as on the full path below: this block returns + // before that one runs, so without it the flag is silently ignored by --scale. Same + // precedence as there - command line, then the file, then the geometry's default. + if (detect_ice_rings.has_value()) + experiment.DetectIceRings(detect_ice_rings.value()); + else if (!dataset->file_detect_ice_rings.has_value()) + experiment.DetectIceRings(indexing_settings.GetRotationIndexing()); experiment.ImportIndexingSettings(indexing_settings); // Start from the same defaults the full pipeline uses, so --scale reproduces the merge that wrote @@ -1575,8 +1582,14 @@ static int RunRugnux(int argc, char **argv) { // --detect-ice-rings[=on|off] overrides the value carried in from the dataset (HDF5MetadataSource // sets DetectIceRings from the master file's detect_ice_rings key); with no flag the dataset stands. + // Where the file says nothing at all, the default is the geometry's: on for rotation, off for + // stills. A rotation sweep sits on the same rings for the whole run, so ice there is a coherent + // systematic worth handling, and the ice-presence gate keeps it inert on a clean crystal; a serial + // stills run has too few spots per image to spend any of them on flagging. if (detect_ice_rings.has_value()) experiment.DetectIceRings(detect_ice_rings.value()); + else if (!dataset->file_detect_ice_rings.has_value()) + experiment.DetectIceRings(rotation_indexing); // Scale-fulls refits the per-frame scale on the rotation combined fulls; on by default for rotation // data (where it lifts ISa substantially) and off for stills. --no-scale-fulls overrides. @@ -1696,12 +1709,12 @@ static int RunRugnux(int argc, char **argv) { logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg); } - if (background_radial_arg) { + if (background_radial_given) { BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); - bis.BackgroundRadialCorrection(*background_radial_arg); + bis.BackgroundRadialCorrection(background_radial_arg); experiment.ImportBraggIntegrationSettings(bis); logger.Info("Background ring: radial curvature correction {}", - *background_radial_arg ? "on" : "off"); + background_radial_arg.has_value() ? (*background_radial_arg ? "on" : "off") : "auto"); } SpotFindingSettings spot_settings;