docs: update analysis references to match current code
Reconcile docs/CPU_DATA_ANALYSIS.md and docs/RUGNUX.md with the current image_analysis/ and rugnux CLI (verified section-by-section against the code): unified profile-fit Bragg integration engine, multi-lattice indexing, azimuthal phi binning, radial parallax/bandwidth profile with sub-pixel centring, rot3d capture-fraction handling, automatic CC1/2 resolution cutoff, and the new rugnux options; fix the section numbering and cross-references; remove the never-implemented French-Wilson and still-partiality descriptions. Delete the stale in-source design notes (ICE_RING_DETECTION, BRAGG_INTEGRATION_ENGINE, NEXTGEN_INTEGRATOR) and fix the two code comments that pointed at them; the BraggIntegrationEngine header no longer claims it is 'not yet wired'. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 4.8
@@ -9,7 +9,7 @@
|
||||
// profile-fits with a measured-width Gaussian (Kabsch-style) - more accurate intensities than the
|
||||
// classical uniform BoxSum; validated on HEWL anomalous (stronger S/Cl peaks vs box-sum). BoxSum is
|
||||
// the simpler, faster fallback. ProfileEmpirical learns the profile per resolution shell from strong
|
||||
// spots - see bragg_integration/BRAGG_INTEGRATION_ENGINE.md.
|
||||
// spots - see docs/CPU_DATA_ANALYSIS.md (Bragg integration).
|
||||
enum class IntegratorMode { BoxSum, ProfileGaussian, ProfileEmpirical };
|
||||
|
||||
class BraggIntegrationSettings {
|
||||
|
||||
@@ -6,6 +6,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
|
||||
* jfjoch_viewer: Open datasets on the WSL2/UNC filesystem (paths starting `\\`); write processing outputs next to the input file, with a Browse button and independent `_process.h5` / merged `.mtz`/`.cif` toggles; and show the determined space group in the merge-statistics window.
|
||||
* rugnux: Accept an absolute `-o` output prefix in offline processing.
|
||||
* Packaging: The self-contained Linux viewer `.tgz` now bundles cuFFT, so it runs without a system CUDA toolkit (`.deb`/`.rpm` are unchanged, distro-managed).
|
||||
* Docs: Bring the analysis references up to date with the code. `docs/CPU_DATA_ANALYSIS.md` now reflects the unified profile-fit Bragg integration engine, multi-lattice indexing, azimuthal phi binning, the radial parallax/bandwidth profile with sub-pixel centring, the rot3d capture-fraction handling and the automatic CC1/2 resolution cutoff, and drops the descriptions of features that were never implemented (French-Wilson amplitudes, the still excitation-error partiality model); `docs/RUGNUX.md` documents the new `--resolution-cutoff`/`--resolution-cc-target`/`--resolution-shells`, `--min-captured-fraction`, `--mosaicity`, `--reference-column`, the azimuthal correction toggles and the geometry-override options, and corrects the `-N` default. The outdated in-source design notes (ICE_RING_DETECTION, BRAGG_INTEGRATION_ENGINE, NEXTGEN_INTEGRATOR) are removed.
|
||||
|
||||
### 1.0.0-rc.157
|
||||
This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use.
|
||||
|
||||
+65
-82
@@ -11,7 +11,7 @@ This document describes the crystallographic algorithms implemented in Jungfrauj
|
||||
5. Bravais lattice / centering inference,
|
||||
6. geometry and lattice refinement,
|
||||
7. reflection prediction (still and rotation),
|
||||
8. Bragg integration by either 2D summation or reference-driven profile fitting,
|
||||
8. Bragg integration by either 2D box summation or profile fitting (Kabsch, reference-free),
|
||||
9. scaling and merging,
|
||||
10. merge-level error modelling and outlier rejection,
|
||||
11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref).
|
||||
@@ -87,29 +87,30 @@ $
|
||||
$
|
||||
Jungfraujoch uses $|\Delta_\mathrm{Ewald}|$ as an operational proxy for excitation error. This appears in:
|
||||
- still prediction (accept if $|\Delta_\mathrm{Ewald}|\le \Delta_\mathrm{cut}$),
|
||||
- profile radius estimation (see §7.1),
|
||||
- still partiality option in scaling/merging (§9.3).
|
||||
- profile radius estimation (see §11.1),
|
||||
- still partiality option in scaling/merging (§10.2).
|
||||
|
||||
---
|
||||
|
||||
## 2. Azimuthal integration (radial profiles)
|
||||
|
||||
Azimuthal integration produces a 1D radial profile $I(q)$ or $I(d)$ by histogramming pixels into radial bins. Pixels are **not split** across bins; each pixel contributes wholly to a single bin.
|
||||
Azimuthal integration produces a radial profile $I(q)$ or $I(d)$ by histogramming pixels into radial bins. Pixels are **not split** across bins; each pixel contributes wholly to a single bin. By default the profile is purely radial (a single azimuthal bin), but the azimuth can optionally be split into up to 512 $\phi$ sectors (`azim_bins`, `--azim-phi-bins`), giving a **2D $q\times\phi$ profile** that exposes azimuthal anisotropy such as detector shadowing or sample texture.
|
||||
|
||||
### 2.1 Histogram estimator
|
||||
|
||||
Let bin index $b(x,y)\in\{0,\dots,B-1\}$ be precomputed from $q(x,y)$ (or equivalently from $d(x,y)$). For each bin $b$:
|
||||
Let bin index $b(x,y)$ be precomputed from $q(x,y)$ (or equivalently from $d(x,y)$) and, when $\phi$ sectors are enabled, the azimuth $\phi(x,y)$ — so $b = b_q + b_\phi B_q$. For each bin $b$:
|
||||
|
||||
- accumulate corrected intensity:
|
||||
- accumulate corrected intensity and its square:
|
||||
$
|
||||
S_b = \sum_{(x,y):\,b(x,y)=b} I(x,y)\,C(x,y),
|
||||
S_b = \sum_{(x,y):\,b(x,y)=b} I(x,y)\,C(x,y),\qquad
|
||||
S^{(2)}_b = \sum I(x,y)^2\,C(x,y)^2,
|
||||
$
|
||||
- and count:
|
||||
$
|
||||
N_b = \#\{(x,y):\,b(x,y)=b \text{ and pixel is valid}\}.
|
||||
$
|
||||
|
||||
A simple mean profile is then $ \bar{I}_b = S_b / N_b$ (when $N_b>0$). Invalid pixels (masked, saturated, detector error codes) are excluded.
|
||||
The profile reports both the mean $\bar{I}_b = S_b / N_b$ (when $N_b>0$) and a per-bin sample standard deviation $\sigma_b = \sqrt{(S^{(2)}_b - S_b^2/N_b)/(N_b-1)}$ (a spread/error estimate for each radial point). Invalid pixels (masked, saturated, detector error codes) are excluded.
|
||||
|
||||
### 2.2 Corrections applied
|
||||
|
||||
@@ -117,8 +118,9 @@ Two standard corrections are available:
|
||||
|
||||
**(i) Solid angle / geometric correction.** A commonly used approximation for flat detectors gives a $\cos^3(2\theta)$ factor:
|
||||
$
|
||||
C_\Omega(2\theta) = \cos^3(2\theta).
|
||||
C_\Omega(2\theta) = \cos^3(2\theta),
|
||||
$
|
||||
applied — like the polarization term below — as a **divisor** (intensities are scaled by $1/\cos^3(2\theta)$), so pixels at higher $2\theta$, which subtend a smaller solid angle, are boosted.
|
||||
|
||||
**(ii) Polarization correction.** With polarization coefficient $P$ (beamline dependent) and azimuth $\phi$:
|
||||
$
|
||||
@@ -129,7 +131,7 @@ applied as a divisor to intensities (i.e. scale by $1/C_\mathrm{pol}$) when enab
|
||||
|
||||
### 2.3 Background estimate for profiles
|
||||
|
||||
A background estimate is derived from the integrated profile using the azimuthal integration settings (details depend on the configured estimator). This background is used for monitoring and diagnostics; it is **not** the same as local Bragg-spot background used in summation integration (§8).
|
||||
A background estimate is derived from the profile as its mean intensity over a fixed low-to-mid $Q$ window (default $2\pi/5$ to $2\pi/3$ Å$^{-1}$). This background is used for monitoring and diagnostics; it is **not** the same as the local Bragg-spot background used in summation integration (§9.2).
|
||||
|
||||
---
|
||||
|
||||
@@ -163,6 +165,7 @@ $
|
||||
|
||||
A pixel is considered strong if:
|
||||
- it is above a photon/count threshold, and
|
||||
- its window contains enough valid neighbours (more than 100), so the local statistics are meaningful, and
|
||||
- $\Delta>0$, and
|
||||
- the squared deviation exceeds a scaled variance:
|
||||
$
|
||||
@@ -178,7 +181,7 @@ Special cases:
|
||||
|
||||
### 3.2 Resolution and ice-ring handling
|
||||
|
||||
Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, pixels in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6).
|
||||
Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, spots in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6).
|
||||
|
||||
A single per-image **ice-ring score** is derived from the azimuthally-integrated radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). (A significance/z-score was considered but is uninformative here: with many photons any real ice ring is highly significant, so the discriminating quantity is the ice *magnitude*, i.e. this ratio.) It is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, distinct from the merge-time ice masking, which is data-driven from the per-ring merged CC1/2.
|
||||
|
||||
@@ -211,7 +214,7 @@ Jungfraujoch supports two complementary indexing strategies:
|
||||
1. **FFT-based indexing** (Rossmann-type): does not require an a priori unit cell; suitable for unknown samples.
|
||||
2. **Fast-feedback indexing** (TORO-like): requires an approximate unit cell; optimized for speed and feedback.
|
||||
|
||||
Both feed into a common robust refinement/selection stage which maximizes the number of inliers under an indexing tolerance.
|
||||
Both feed into a common robust refinement/selection stage which maximizes the number of inliers under an indexing tolerance, and which can return **more than one lattice** per image (multi-lattice indexing; see §5.4).
|
||||
|
||||
### 4.1 Indexed-spot decision (inlier test)
|
||||
|
||||
@@ -225,7 +228,7 @@ Let $(h,k,l)=(\mathrm{round}(h_f),\mathrm{round}(k_f),\mathrm{round}(l_f))$ and
|
||||
$
|
||||
\delta^2 = (h_f-h)^2 + (k_f-k)^2 + (l_f-l)^2.
|
||||
$
|
||||
A spot is indexed if $\delta^2 \le \tau^2$, where $\tau$ is the configured tolerance.
|
||||
A spot is indexed if $\delta^2 < \tau^2$, where $\tau$ is the configured tolerance.
|
||||
|
||||
For indexed spots, the reciprocal lattice point $\mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^*$ is used to compute $\Delta_\mathrm{Ewald}(\mathbf{p})$ (stored as a diagnostic and later used in profile-radius estimation).
|
||||
|
||||
@@ -250,7 +253,7 @@ where $L_\mathrm{max}$ is the maximum expected real-space unit-cell edge (Å). T
|
||||
|
||||
### 5.2 FFT peak picking and candidate vectors
|
||||
|
||||
For each direction, the FFT magnitude spectrum is computed; peaks correspond to periodicities along $\mathbf{u}_d$. Each direction yields a candidate real-space length $L$ with maximum spectral magnitude (subject to $L\ge L_\mathrm{min}$).
|
||||
For each direction, the FFT magnitude spectrum is computed; peaks correspond to periodicities along $\mathbf{u}_d$. Each direction yields a candidate real-space length $L$ chosen **not** by raw magnitude but by **maximum prominence above a running-mean local background** (subtracting the broad low-frequency envelope that otherwise dominates on weak or pink-beam frames), subject to $L\ge L_\mathrm{min}$.
|
||||
|
||||
Candidate vectors are $\mathbf{v}_d = L_d\,\mathbf{u}_d$.
|
||||
|
||||
@@ -258,22 +261,16 @@ A collinearity filter removes nearly parallel vectors (e.g. within 5°) and atte
|
||||
|
||||
### 5.3 Lattice reduction and cell candidates
|
||||
|
||||
Triples of candidate vectors are combined to form candidate bases $(\mathbf{A},\mathbf{B},\mathbf{C})$. A simple reduction is applied:
|
||||
$
|
||||
\mathbf{B} \leftarrow \mathbf{B} - \mathrm{round}\!\left(\frac{\mathbf{B}\cdot\mathbf{A}}{\mathbf{A}\cdot\mathbf{A}}\right)\mathbf{A},
|
||||
$
|
||||
$
|
||||
\mathbf{C} \leftarrow \mathbf{C} - \mathrm{round}\!\left(\frac{\mathbf{C}\cdot\mathbf{A}}{\mathbf{A}\cdot\mathbf{A}}\right)\mathbf{A}
|
||||
- \mathrm{round}\!\left(\frac{\mathbf{C}\cdot\mathbf{B}}{\mathbf{B}\cdot\mathbf{B}}\right)\mathbf{B}.
|
||||
$
|
||||
|
||||
Candidates are filtered by allowed length and angle ranges.
|
||||
Triples of candidate vectors are combined to form candidate bases $(\mathbf{A},\mathbf{B},\mathbf{C})$, each reduced to its **Niggli-reduced cell** (Gruber-vector reduction) before comparison, and filtered by allowed length and angle ranges. Two passes are run: a standard pass forms shortest-vector triples from the ~30 strongest filtered directions; if the best cell then indexes fewer than half the spots, a **widened fallback** anchors the two shortest axes and lets the third range over up to ~60 candidate vectors (deduplicated by Niggli cell), catching large, elongated or superstructure cells the first pass misses.
|
||||
|
||||
### 5.4 Robust refinement and best-cell selection
|
||||
|
||||
Candidate bases are refined against observed spots using an iterative inlier‑focused least‑squares procedure (trimmed/contracting threshold). The output cell is chosen to:
|
||||
1. maximize the number of indexed spots under the tolerance $\tau$, and
|
||||
2. break ties by a refined score (smaller residual threshold/score is preferred).
|
||||
Candidate bases are refined against observed spots using an iterative inlier‑focused least‑squares procedure (trimmed/contracting threshold). Candidates are then ranked:
|
||||
1. more indexed spots wins — **unless** two candidates index within ~10 % of each other, in which case
|
||||
2. the **smaller-volume** cell is preferred (when the volumes differ by more than ~5 %), avoiding a doubled supercell, then
|
||||
3. the smaller refinement score, then the spot count again.
|
||||
|
||||
Selection is **not limited to a single lattice**: after the best cell is accepted, further lattices are added as separate crystals provided fewer than ~40 % of their indexed spots overlap an already-accepted lattice (up to two extra by default), so split or multi-lattice crystals are indexed rather than discarded.
|
||||
|
||||
An optional reference unit cell (if supplied) restricts acceptance to cells within a relative distance tolerance in edge lengths (permutation-invariant).
|
||||
|
||||
@@ -292,9 +289,9 @@ If not, Jungfraujoch attempts to infer the most plausible Bravais lattice type f
|
||||
The output includes:
|
||||
- a conventional cell,
|
||||
- crystal system (triclinic, monoclinic, …),
|
||||
- centering symbol $P, A, B, C, I, F, R$.
|
||||
- centering symbol (one of $P, C, I, F, R$; the $A/B$ variants are not emitted here — they are handled only later as prediction absences, §8.4).
|
||||
|
||||
This stage provides centering information used for systematic absences in prediction (§7.3) and for reporting.
|
||||
This stage provides centering information used for systematic absences in prediction (§8.4) and for reporting.
|
||||
|
||||
**Note.** In ambiguous or special cases, forcing space group to $P1$ (no symmetry assumptions) is recommended.
|
||||
|
||||
@@ -315,6 +312,8 @@ The refinement jointly optimizes, depending on mode and constraints:
|
||||
- crystal orientation (a global rotation),
|
||||
- unit-cell parameters, with constraints determined by inferred crystal system.
|
||||
|
||||
By default only the beam center, unit cell and crystal orientation are refined; the detector distance, tilt angles and rotation-axis direction are held fixed unless explicitly enabled. A lighter **orientation-only** mode refines just the crystal orientation (with a weak small-rotation prior on the poorly-determined out-of-plane component), for stills whose geometry is already trusted.
|
||||
|
||||
For higher symmetries, constraints are enforced, e.g.
|
||||
- cubic: $a=b=c,\ \alpha=\beta=\gamma=90^\circ$,
|
||||
- tetragonal: $a=b$,
|
||||
@@ -341,15 +340,15 @@ For oscillation/rotation data, each image corresponds to a rotation angle $\phi$
|
||||
$
|
||||
\mathbf{s}_\mathrm{obs,ref} = R(\phi)\,\mathbf{s}_\mathrm{obs},
|
||||
$
|
||||
with $R(\phi)$ constructed from the axis-angle representation of the goniometer model.
|
||||
with $R(\phi)$ constructed from the axis-angle representation of the goniometer model. The angle $\phi$ is taken at the centre of each frame's oscillation (the frame angle plus half the oscillation width).
|
||||
|
||||
### 7.4 Multi-stage tightening of inlier tolerance
|
||||
|
||||
Refinement is performed in stages with decreasing acceptance tolerance for including reflections (e.g. from coarse to fine), which stabilizes convergence when starting from imperfect indexing and approximate geometry.
|
||||
Refinement is performed in stages with decreasing acceptance tolerance for including reflections (three stages, indexing tolerance $0.3\to0.2\to0.1$), which stabilizes convergence when starting from imperfect indexing and approximate geometry.
|
||||
|
||||
---
|
||||
|
||||
## 9. Bragg integration
|
||||
## 8. Reflection prediction
|
||||
|
||||
Jungfraujoch predicts reflection positions for integration by enumerating Miller indices within a resolution cutoff and accepting those that satisfy a diffraction condition model.
|
||||
|
||||
@@ -368,6 +367,8 @@ $
|
||||
$
|
||||
Accepted reflections are projected to the detector by intersecting the diffracted direction $\mathbf{S}=\mathbf{S}_0+\mathbf{p}$ with the detector plane, using the current geometry.
|
||||
|
||||
When the beam has a finite energy bandwidth, this window is **broadened radially per reflection**: the cutoff is combined in quadrature with a bandwidth smear, $\sqrt{\Delta_\mathrm{cut}^2 + (3\,\sigma_\mathrm{bw})^2}$, where $\sigma_\mathrm{bw}\propto|p_z|$ (the reciprocal-space depth along the beam, growing as $\sim 1/d^2$). This keeps high-resolution reflections — smeared by the bandwidth into radial streaks — from being clipped. The same $\sigma_\mathrm{bw}$ is deconvolved from the measured profile radius (§11.1), so it is not double-counted.
|
||||
|
||||
### 8.3 Rotation prediction (Laue equation + partiality model)
|
||||
|
||||
For rotation/oscillation datasets, Jungfraujoch solves for rotation angles $\phi$ where the rotated reciprocal lattice point satisfies the Ewald-sphere condition. In an XDS-like notation, define:
|
||||
@@ -422,7 +423,7 @@ Pixels are classified by their squared distance $r^2=(x-x_p)^2+(y-y_p)^2$:
|
||||
- **signal region:** $r^2 < r_1^2$,
|
||||
- **background annulus:** $r_2^2 \le r^2 < r_3^2$.
|
||||
|
||||
Invalid pixels (masked/bad/saturated) are excluded from both sums.
|
||||
Invalid pixels (masked/bad/saturated) are excluded from both sums. In addition, pixels lying inside the signal disk ($r<r_2$) of any *other* predicted reflection are removed from this reflection's background annulus, so a neighbouring spot cannot leak into the background estimate.
|
||||
|
||||
### 9.2 Box summation (seed and fallback)
|
||||
|
||||
@@ -437,23 +438,27 @@ $
|
||||
\hat{b} = \frac{B}{n_B},\qquad
|
||||
\hat{I} = S - n_S \hat{b},
|
||||
$
|
||||
with a Poisson-like uncertainty $\sigma(\hat{I})\approx\sqrt{S}$ (floored at 1). 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}\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.
|
||||
|
||||
The background mean is computed with a single **high-outlier reject** (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, 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.
|
||||
For the **profile-fit path on broadband (still) data**, the background mean is additionally 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. The reject is **not** applied to plain box summation (`--integrator boxsum`) or to monochromatic/rotation data.
|
||||
|
||||
### 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:
|
||||
|
||||
1. **Seed.** Box-sum every reflection (§9.2) to get a rough intensity and observed centroid, and select strong spots (significance $\ge 5$).
|
||||
2. **Build the profile.** From the strong spots, form a profile **per resolution shell**: an isotropic Gaussian of the measured second moment (the default), or an empirical averaged grid. The width is shell-dependent because spot size grows with resolution; the intrinsic spot is essentially round in the detector plane (per-detector-region and crystal-anisotropy profiles were evaluated and add nothing — the real crystal anisotropy lives in the discarded rocking direction). When a finite energy bandwidth is set, however, it smears each spot **radially** by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ (distance from the beam centre, large at high resolution), turning high-resolution spots into radial streaks. There the profile is elongated *only along the radial direction* per reflection ($\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2$, tangential unchanged) on a grid grown to hold the streak, capturing it without the tangential background an isotropic widening would add.
|
||||
2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moment of the strong spots (shell-dependent because spot size grows with resolution); the intrinsic spot is essentially round in the detector plane (per-detector-region and crystal-anisotropy profiles were evaluated and add nothing — the real crystal anisotropy lives in the discarded rocking direction). For `empirical` the profile is instead the averaged, centroid-aligned, background-subtracted pixel grid of the shell's strong spots. Either way the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially:
|
||||
- a finite energy **bandwidth** smears each spot by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ ($R_\mathrm{px}$ = distance from the beam centre, large at high resolution), and
|
||||
- sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus, on the monochromatic path, a small fixed weak-spot capture term.
|
||||
|
||||
These combine as $\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)$ (tangential unchanged), on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add.
|
||||
3. **Fit (Kabsch).** With profile $P$, background $B$ and the shell variance model, the intensity and its uncertainty are
|
||||
$
|
||||
I = \frac{\sum P\,(c-B)/v}{\sum P^2/v},\qquad
|
||||
\sigma = \sqrt{\frac{1}{\sum P^2/v}},\qquad
|
||||
v = B + \max(I,0)\,P,
|
||||
$
|
||||
where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated. The rotation/excitation partiality is carried exactly as in the box-sum path.
|
||||
where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the variance floors the background at $1/12$ (the integer-binning pixel-variance floor). The rotation/excitation partiality is carried exactly as in the box-sum path.
|
||||
|
||||
The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`).
|
||||
|
||||
@@ -465,7 +470,7 @@ For integrated reflections, polarization correction can be applied as a multipli
|
||||
|
||||
## 10. Scaling and merging
|
||||
|
||||
After per-image integration, Jungfraujoch scales observations and merges them into unique reflections. The design is intentionally compatible with XDS/XSCALE concepts, while supporting both still and rotation partiality models.
|
||||
After per-image integration, Jungfraujoch scales observations and merges them into unique reflections. The design is intentionally compatible with XDS/XSCALE concepts, and handles both still and rotation data.
|
||||
|
||||
### 10.1 Observation model
|
||||
|
||||
@@ -484,40 +489,30 @@ A least-squares objective is minimized:
|
||||
$
|
||||
\sum_{ij} \left(\frac{I_{ij}^{\mathrm{pred}} - I_{ij}^{\mathrm{obs}}}{\sigma_{ij}}\right)^2
|
||||
$
|
||||
with regularization on $G_i$ and optional smoothness constraints (particularly meaningful for rotation series).
|
||||
solved by robust (Cauchy) weighted least squares, with optional post-fit smoothing of the per-frame scales for rotation series (§10.3).
|
||||
|
||||
### 10.2 Partiality models available
|
||||
### 10.2 Partiality models
|
||||
|
||||
Jungfraujoch supports several partiality choices:
|
||||
The partiality applied is fixed by the data type and scaling stage, not chosen from a user menu:
|
||||
|
||||
1. **Rotation partiality** (XDS-like; see §8.3):
|
||||
1. **Rotation partiality** (XDS-like; see §8.3), used for the per-frame scaling of rotation partials:
|
||||
$
|
||||
P_{ij} = \frac{1}{2}\left[
|
||||
\mathrm{erf}\!\left(\frac{\Delta\phi_{ij}+\Delta\phi/2}{\sqrt{2}\,\sigma_{M,i}/\zeta_{ij}}\right) -
|
||||
\mathrm{erf}\!\left(\frac{\Delta\phi_{ij}-\Delta\phi/2}{\sqrt{2}\,\sigma_{M,i}/\zeta_{ij}}\right)
|
||||
\right].
|
||||
$
|
||||
Mosaicity $\sigma_{M,i}$ can be refined per image group with bounds.
|
||||
The mosaicity $\sigma_{M,i}$ is **measured once per image at indexing** (MLE, §11.2) and held fixed during scaling — only smoothed in frame order (§10.3), never re-refined (it is degenerate with the scale $G$; §11.2).
|
||||
|
||||
2. **Still partiality** (excitation-error proxy):
|
||||
$
|
||||
P_{ij} = \exp\!\left(-\frac{\Delta_\mathrm{Ewald}^2}{R_i^2}\right),
|
||||
$
|
||||
where $R_i^2$ is a refined width parameter (bounded).
|
||||
2. **Unity** ($P_{ij}=1$): used for the scale-on-fulls refit (§10.6), where each observation is already a complete reflection.
|
||||
|
||||
3. **Unity**: $P_{ij}=1$.
|
||||
|
||||
4. **Fixed**: use the per-reflection partiality carried from prediction.
|
||||
3. **Fixed**: use the per-reflection partiality carried from prediction. Still/serial images are predicted with $P=1$, so their scaling is effectively unity/fixed — there is no excitation-error still-partiality model.
|
||||
|
||||
Reflections below a minimum partiality can be rejected from merging to avoid unstable corrections.
|
||||
|
||||
### 10.3 Regularization and smoothness
|
||||
### 10.3 Smoothing of per-frame scales
|
||||
|
||||
To stabilize scale determination, a weak prior $G_i\approx 1$ is used. For rotation datasets, optional smoothness encourages slowly varying scales and mosaicity:
|
||||
$
|
||||
\log G_{i-1} - 2\log G_i + \log G_{i+1} \approx 0,
|
||||
$
|
||||
(and similarly for mosaicity), reflecting the expectation of gradual changes during a rotation scan.
|
||||
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.
|
||||
|
||||
### 10.4 Merging estimator
|
||||
|
||||
@@ -538,12 +533,13 @@ An internal-consistency term can inflate uncertainties when multiple observation
|
||||
### 10.5 Merging statistics
|
||||
|
||||
Per-shell and overall merging statistics are computed on corrected intensities, including:
|
||||
- number of observations,
|
||||
- number of unique reflections,
|
||||
- number of observations and of unique reflections, and multiplicity,
|
||||
- mean $I/\sigma(I)$,
|
||||
- an R$_\mathrm{meas}$-like quantity derived from within‑HKL deviations (shell-binned).
|
||||
- $R_\mathrm{meas}$ (the redundancy-independent Diederichs–Karplus form) from within‑HKL deviations,
|
||||
- $\mathrm{CC}_{1/2}$ (half-set correlation) and, when a reference dataset is supplied, $\mathrm{CC}_\mathrm{ref}$,
|
||||
- completeness against the enumerated reflections for the cell and symmetry.
|
||||
|
||||
Completeness requires enumeration of possible reflections given a unit cell and symmetry; where this is not fully available, completeness may be reported as 0 or omitted.
|
||||
The error model is refined as $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2$ with a systematic floor $\sigma\ge b|I|$; the asymptotic signal-to-noise $\mathrm{ISa}=1/b$ is reported and written to the output files.
|
||||
|
||||
### 10.6 Rotation datasets: combining partials into fulls (3D integration)
|
||||
|
||||
@@ -552,7 +548,7 @@ In a rotation scan a reflection is recorded as a series of *partials* spread acr
|
||||
The combine groups each reflection's partials into rocking events (contiguous runs of frames) and reduces each event to one full:
|
||||
|
||||
- **De-biased weighted sum.** Partials are combined by inverse-variance weighting, where each partial's variance is its background-noise component plus the *model* signal shared across the event (Kabsch profile-fit form). Using the shared model signal rather than the individual down-fluctuating intensity stops weak partials from being over-weighted, which would otherwise inflate the merged error model. The weights depend on the full, so the estimate is iterated.
|
||||
- **Captured fraction.** The partiality summed over the event, $f=\min(1,\sum_j p_j)$, measures how completely the rocking curve was sampled; it replaces a per-partial minimum-partiality cut, because an event seen over only a few percent of its curve is unreliable however many frames it spans.
|
||||
- **Captured fraction.** The partiality summed over the event, $f=\min(1,\sum_j p_j)$, measures how completely the rocking curve was sampled. A full whose curve was captured below a threshold (`--min-captured-fraction`, default 0.7 for rotation) is dropped — an event seen over only a small fraction of its curve is unreliable however many frames it spans. (The per-partial minimum-partiality cut of §10.2 still applies upstream, in the per-frame scaling.)
|
||||
- **Capture-aware uncertainty.** A full captured incompletely ($f<1$) is extrapolated and biased high. The unobserved fraction is charged as an extra systematic uncertainty, $\sigma^2 \leftarrow \sigma^2 + \big(c\,(1-f)\,I\big)^2$, so the merge down-weights these extrapolated fulls and the error model treats their scatter as expected. It is enabled by default for the rotation path.
|
||||
|
||||
The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4). Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$.
|
||||
@@ -579,7 +575,7 @@ The estimated mosaicity feeds the rotation prediction (how many frames each refl
|
||||
|
||||
---
|
||||
|
||||
## 12. Wilson statistics and French–Wilson treatment
|
||||
## 12. Auxiliary statistics: ⟨I/σ(I)⟩ and Wilson plot
|
||||
|
||||
### 12.1 Per-shell ⟨I/σ(I)⟩
|
||||
|
||||
@@ -597,24 +593,6 @@ $
|
||||
$
|
||||
A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of $B$, subject to basic quality checks (e.g. $R^2$ threshold).
|
||||
|
||||
### 12.3 French–Wilson (posterior expectation of I and |F|)
|
||||
|
||||
To mitigate negative intensities and obtain physically meaningful amplitudes, Jungfraujoch implements a French–Wilson style Bayesian treatment using per-shell mean intensity as a prior scale.
|
||||
|
||||
For each merged observation $I_\mathrm{obs}$ with uncertainty $\sigma$, the posterior over true intensity $I\ge 0$ is:
|
||||
$
|
||||
p(I\mid I_\mathrm{obs}) \propto p(I)\,\exp\!\left(-\frac{(I_\mathrm{obs}-I)^2}{2\sigma^2}\right),
|
||||
$
|
||||
with priors differing between acentric and centric cases (standard Wilson distributions).
|
||||
|
||||
Numerical quadrature over a scaled intensity variable is used to compute posterior moments:
|
||||
- $\langle I\rangle$,
|
||||
- $\langle |F|\rangle = \langle \sqrt{I}\rangle$,
|
||||
and an amplitude uncertainty estimate via:
|
||||
$
|
||||
\sigma_F \approx \sqrt{\langle I\rangle - \langle |F|\rangle^2}.
|
||||
$
|
||||
|
||||
---
|
||||
|
||||
## 13. Practical notes and limitations
|
||||
@@ -622,4 +600,9 @@ Numerical quadrature over a scaled intensity variable is used to compute posteri
|
||||
- **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.
|
||||
- **Rotation vs still modes** differ substantially in prediction and scaling because partiality is angle-driven in rotation data and excitation-error-driven in still data.
|
||||
- **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 scaled with unit partiality.
|
||||
- **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.
|
||||
- **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ are computed and written to the merged mmCIF as a twinning diagnostic.
|
||||
- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default).
|
||||
- **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it.
|
||||
- **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream.
|
||||
|
||||
+26
-4
@@ -51,7 +51,9 @@ the first pass.
|
||||
lattice symmetry). No-reference scaling additionally emits per-iteration `<prefix>_iterN_scale.dat`.
|
||||
|
||||
Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the
|
||||
console.
|
||||
console. By default the written resolution is trimmed automatically where CC1/2 falls off
|
||||
(`--resolution-cutoff cc-logistic`, CC1/2 target 0.30); set `--scaling-high-resolution` to fix the
|
||||
limit by hand, or `--resolution-cutoff off` to keep the full range.
|
||||
|
||||
## Re-scaling and re-merging (`rugnux --scale`)
|
||||
|
||||
@@ -114,7 +116,7 @@ General:
|
||||
| Option | Description |
|
||||
| --- | --- |
|
||||
| `-o, --output-prefix <txt>` | Output file prefix (default: `output`) |
|
||||
| `-N, --threads <num>` | Number of worker threads (default: 1) |
|
||||
| `-N, --threads <num>` | Number of worker threads (default: all hardware threads) |
|
||||
| `-s, --start-image <num>` | First image to process (default: 0) |
|
||||
| `-e, --end-image <num>` | Last image to process (default: all) |
|
||||
| `-t, --stride <num>` | Process every *n*-th image (default: 1) |
|
||||
@@ -135,7 +137,7 @@ Spot finding:
|
||||
| `--spot-threshold <num>` | Photon-count threshold for spot finding (default: 10) |
|
||||
| `--spot-high-resolution <num>` | High-resolution limit for spot finding, Å (default: 1.5) |
|
||||
| `--max-spots <num>` | Maximum spot count (default: 250) |
|
||||
| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots and exclude ice-ring reflections from scaling/merging; overrides the dataset setting (default: use the dataset value) |
|
||||
| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use the dataset value) |
|
||||
|
||||
Azimuthal integration (the radial profile behind the per-image ice-ring score):
|
||||
|
||||
@@ -145,6 +147,8 @@ Azimuthal integration (the radial profile behind the per-image ice-ring score):
|
||||
| `--azim-min-q <num>` | Minimum Q, 1/Å |
|
||||
| `--azim-max-q <num>` | Maximum Q, 1/Å |
|
||||
| `--azim-phi-bins <num>` | Number of azimuthal (phi) bins (default: 1) |
|
||||
| `--polarization-correction <on\|off>` | Enable/disable the azimuthal polarization correction |
|
||||
| `--solid-angle-correction <on\|off>` | Enable/disable the azimuthal solid-angle correction |
|
||||
|
||||
Indexing:
|
||||
|
||||
@@ -179,14 +183,20 @@ Scaling and merging:
|
||||
| `--scale-fulls` / `--no-scale-fulls` | rot3d: refit a per-frame scale on the combined fulls (XDS order, Unity model); on by default for rotation data, off for stills |
|
||||
| `--smooth-g[=deg]` | rot3d: smooth the per-frame scale *G* over a degree range before the 3D combine (XDS DELPHI-like; default 5° for rotation, 0 = off) |
|
||||
| `--capture-uncertainty <num>` | rot3d: systematic sigma on under-captured fulls, ~num·(1−captured_fraction)·I (default: 1.0 for rotation, 0 otherwise) |
|
||||
| `--scaling-high-resolution <num>` | High-resolution limit for scaling, Å (default: no limit) |
|
||||
| `--min-captured-fraction <num>` | 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) |
|
||||
| `--scaling-high-resolution <num>` | High-resolution limit for scaling, Å — manual override (default: no limit; disables the automatic cutoff below) |
|
||||
| `--resolution-cutoff <txt>` | Automatic high-resolution cutoff for the written reflections and reported shells: `cc-logistic` \| `off` (default: `cc-logistic`; ignored when `--scaling-high-resolution` is set) |
|
||||
| `--resolution-cc-target <num>` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) |
|
||||
| `--resolution-shells <num>` | Number of resolution shells in the reported statistics table (default: 10) |
|
||||
| `--min-partiality <num>` | Minimum partiality to accept a reflection (default: 0.02) |
|
||||
| `--reject-outliers <num>` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) |
|
||||
| `--reject-delta-cchalf <num>` | Drop images with ΔCC1/2 below mean − N·stddev (default: off) |
|
||||
| `--min-image-cc <num>` | Per-image CC limit, percent (default: no limit) |
|
||||
| `--mosaicity <num>` | Diagnostic: fix the scaling mosaicity (°) instead of using the per-image seed |
|
||||
| `--scaling-iterations <num>` | Scaling iterations with no reference data (default: 3) |
|
||||
| `--scaling-output <txt>` | Reflection output format: `cif` (mmCIF, default) \| `mtz` \| `txt` |
|
||||
| `-z, --reference-mtz <file>` | Reference MTZ (enables reference-driven scaling) |
|
||||
| `--reference-column <label>` | Reference MTZ column to use (default: auto — F-model, else IMEAN/I/…) |
|
||||
| `--write-process-h5` | Also write the (large) `_process.h5` when merging (default: only `.mtz`/`.cif`) |
|
||||
|
||||
Integration:
|
||||
@@ -196,3 +206,15 @@ Integration:
|
||||
| `--integrator <txt>` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) |
|
||||
| `--integration-radius <r>` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` |
|
||||
| `--bandwidth <num>` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) |
|
||||
|
||||
Geometry overrides (defaults are taken from the input file; override them to reprocess with a corrected geometry):
|
||||
|
||||
| Option | Description |
|
||||
| --- | --- |
|
||||
| `--beam-x <num>` | Beam centre X (pixel) |
|
||||
| `--beam-y <num>` | Beam centre Y (pixel) |
|
||||
| `--detector-distance <num>` | Detector distance (mm) |
|
||||
| `--wavelength <num>` | Wavelength (Å) |
|
||||
| `--rot1 <num>` | PONI detector rotation 1 (rad) |
|
||||
| `--rot2 <num>` | PONI detector rotation 2 (rad) |
|
||||
| `--polarization <num>` | Polarization factor |
|
||||
|
||||
@@ -1,71 +0,0 @@
|
||||
# Auto-gated ice-ring detection & estimation (design)
|
||||
|
||||
_Status: designed and Python-prototyped (2026-06-30), not yet implemented in C++. The current
|
||||
`--detect-ice-rings` flag uses the fixed tabulated `ICE_RING_RES_A` with a 0.03 q half-width and must
|
||||
be enabled by hand; this document describes making it automatic and data-driven._
|
||||
|
||||
## Motivation
|
||||
|
||||
`--detect-ice-rings` today is a manual on/off flag over a fixed ring list. Two improvements:
|
||||
1. **Auto-gate** it — turn the handling on only when ice is actually present, so it can default to
|
||||
`auto` without risking dropping good reflections on ice-free data.
|
||||
2. **Estimate** the real ring positions/widths from the data instead of the table (rings drift slightly
|
||||
with detector geometry, λ and ice texture/temperature).
|
||||
|
||||
Both fall out of the azimuthally-averaged radial **background** profile, which the `azint` engine already
|
||||
produces (`AzimuthalIntegrationProfile`) — ice rings are sharp peaks in it, distinct from Bragg spots
|
||||
(which are localised, not azimuthally complete).
|
||||
|
||||
## Method (validated in Python on raw NXmx frames)
|
||||
|
||||
1. Accumulate a **clipped-mean radial background profile** over ~30 frames spread across the sweep
|
||||
(clip at ~20 counts to kill Bragg/hot pixels; mask the beamstop and bad pixels). This is exactly
|
||||
`AzimuthalIntegrationProfile::GetResult()`.
|
||||
2. Smooth-baseline it (rolling median, ±25 bins) and take `excess = profile − baseline`;
|
||||
`rel = excess / max(baseline, ε)`.
|
||||
3. Map radius → d → q (`θ = 0.5·atan(R·pix/dist)`, `d = λ/2sinθ`, `q = 2π/d`).
|
||||
|
||||
### Detection statistic (self-calibrating z-score)
|
||||
|
||||
For each tabulated ice line, `z = (max rel within ±0.03 q of the line − μ_ctrl) / σ_ctrl`, where the
|
||||
control population is all valid q at least 0.06 away from every ice line (in-image noise — so the
|
||||
threshold needs no per-detector tuning).
|
||||
|
||||
| dataset | max line z | #lines z>5 | verdict |
|
||||
|---|---:|---:|---|
|
||||
| EP0117 (diffuse ice) | 11.8 | 6 | ICE |
|
||||
| EP0124 (sharp 3 rings) | 17.0 | 3 | ICE |
|
||||
| EP0210 (control) | 1.6 | 0 | no ice |
|
||||
| MyoB (control) | 0.8 | 0 | no ice |
|
||||
|
||||
**Gate: enable ice handling iff ≥2 ice lines have z > 5.** Large, clean margin.
|
||||
|
||||
### Position & width estimation
|
||||
|
||||
- **Positions:** `ICE_RING_RES_A` is accurate for the JUNGFRAU/λ=0.953 Å — measured centroid offset
|
||||
Δq = +0.004 ± 0.011 q (~0.25 %). No recalibration needed; the table stands.
|
||||
- **Widths:** measured ring FWHM is **0.055–0.068 q** ⇒ half-width ~0.027–0.034. The fixed 0.02 q
|
||||
under-covered the strong low-res rings (already bumped to **0.03** as the fallback). Best is per-ring
|
||||
`half_width = 0.5·FWHM` clamped to [0.012, 0.040], estimated from the peak shape.
|
||||
|
||||
## Proposed wiring (small, no new data path)
|
||||
|
||||
- New `image_analysis/azint/IceRingDetect.{h,cpp}`:
|
||||
`IceRingModel DetectIceRings(const AzimuthalIntegrationProfile&)` → `{bool present; vector<IceRing>
|
||||
rings;}` with per-ring `{centre_q, half_width_q}`.
|
||||
- Add an `IsOnIceRing(d, const vector<IceRing>&)` overload in `common/Definitions.h` (per-ring centre
|
||||
+ width), alongside the existing fixed-table one.
|
||||
- `tools/rugnux_cli.cpp`: change `--detect-ice-rings` to take `auto|on|off` (default `auto`); in
|
||||
`auto`, run `DetectIceRings` on the accumulated profile and set
|
||||
`experiment.DetectIceRings(model.present)` + feed `model.rings` to the spot-finder
|
||||
(`MarkIceRings`) and the scaling exclusion path.
|
||||
- Online (broker): accumulate the profile over the first N frames, freeze the model, apply thereafter.
|
||||
|
||||
Prototype scripts + the `ice_ring_detect.diff` sketch are in the session scratchpad
|
||||
(`scratchpad/ice_detect/`).
|
||||
|
||||
## Note
|
||||
|
||||
For the EP0117 test crystal specifically, ice is a **secondary** issue — its dominant problem is
|
||||
radiation damage across the sweep, addressed by per-observation outlier rejection (`--reject-outliers`),
|
||||
not ice handling. Auto-gated ice detection matters for genuinely ice-dominated, otherwise-good crystals.
|
||||
@@ -1,107 +0,0 @@
|
||||
# BraggIntegrationEngine — GPU box + profile-fit integrator
|
||||
|
||||
Status: **standalone, not yet wired into the pipeline.** This is a GPU (CUDA) reimplementation of
|
||||
`ProfileIntegrate2D.cpp` (the default Kabsch profile-fit integrator, ~142 ms/frame on CPU) and
|
||||
`BraggIntegrate2D.cpp` (the box-sum integrator) under one roof, following the
|
||||
`AzIntEngine` / `ROIIntegration` base + CPU + GPU pattern. It reproduces the CPU result up to
|
||||
floating-point precision and runs in **< 2 ms/frame** for realistic reflection counts.
|
||||
|
||||
## Files (all in `image_analysis/bragg_integration/`)
|
||||
|
||||
| file | role |
|
||||
|------|------|
|
||||
| `BraggIntegrationEngine.{h,cpp}` | base class: extracts the fixed per-experiment config (r-radii, ellipse coefficients, geometry) and owns `Finalize()` (polarization + `image_scale_corr` + `vector<Reflection>` assembly) |
|
||||
| `BraggIntegrationEngineCPU.{h,cpp}` | plain-C++ engine: the fallback **and** the numeric oracle for the GPU |
|
||||
| `BraggIntegrationEngineGPU.{h,cu}` | the CUDA engine |
|
||||
| `tests/BraggIntegrationEngineGPUTest.cpp` | GPU-vs-CPU equivalence test + hidden `[bragg_bench]` benchmark |
|
||||
|
||||
`CMakeLists.txt` builds the base + CPU always; the `.cu` + GPU header are added only under
|
||||
`JFJOCH_CUDA_AVAILABLE` (mirrors `azint`/`roi`). No change needed to `image_analysis/CMakeLists.txt`
|
||||
(the aggregate `JFJochImageAnalysis` already links `JFJochBraggIntegration` and
|
||||
`JFJochImagePreprocessing`, so the `ImagePreprocessorBuffer` symbols resolve at the final link).
|
||||
|
||||
## Interface
|
||||
|
||||
```cpp
|
||||
class BraggIntegrationEngine {
|
||||
BraggIntegrationEngine(const DiffractionExperiment&);
|
||||
virtual std::vector<Reflection> Run(const ImagePreprocessorBuffer& image, // preprocessed int32
|
||||
const std::vector<Reflection>& predicted, size_t npredicted,
|
||||
int64_t image_number) = 0;
|
||||
};
|
||||
// CPU: BraggIntegrationEngineCPU(experiment)
|
||||
// GPU: BraggIntegrationEngineGPU(experiment, std::shared_ptr<CudaStream>)
|
||||
```
|
||||
|
||||
The integrator mode comes from `experiment.GetBraggIntegrationSettings().GetIntegrator()`:
|
||||
`BoxSum` (= BraggIntegrate2D), `ProfileGaussian` (default), `ProfileEmpirical`. The box sum is also
|
||||
Pass A / the seed of the two profile modes, so it always runs.
|
||||
|
||||
## Key design points
|
||||
|
||||
- **Input is the preprocessed int32 buffer**, NOT the raw `CompressedImage`. Masked/bad pixels are
|
||||
`INT32_MIN`, saturated are `INT32_MAX` (a pixel is valid iff `v != INT32_MIN && v != INT32_MAX`).
|
||||
This is the same buffer `AzIntEngineGPU`/`ROIIntegrationGPU` consume. **Consequence:** the `±1`
|
||||
special/saturation band that `ProfileIntegrate2D` rejects on the lossy `CompressedImage` is handled
|
||||
upstream by the preprocessor instead — correct for the online path; sanity-check it for offline
|
||||
`rugnux` if that ever uses this engine on a lossy-compressed stored file.
|
||||
- Config uses raw detector dims — `xpixel = GetXPixelsNum()`, `npixel = GetPixelsNum()` — matching how
|
||||
`MXAnalysisWithoutFPGA` sizes the `ImagePreprocessorBuffer` and the frame the predicted `predicted_x`
|
||||
live in.
|
||||
- **GPU: one CUDA block per reflection**, shared-memory reductions. Kernels (all in the `.cu`):
|
||||
`reset → mark_mask → boxsum → learn_profile → build_profiles → fit` (Kabsch, 4 iters). The
|
||||
resolution shell is computed inline (`compute_shell`), so there is no separate shell pass. `BoxSum`
|
||||
mode stops after `boxsum`.
|
||||
- **Precision:** the hot path (profile learning + Gaussian fit) is `float` (FP64 is throttled on
|
||||
consumer GPUs; the extraction is Poisson-noise limited, so float matches the double CPU to ~1e-4).
|
||||
`double` is kept only for the box-sum background mean / centroid, where it drives weak intensities.
|
||||
- Per frame only the predicted centres (`predicted_x/y`, `d`) are uploaded; the image is already on
|
||||
the device via `ImagePreprocessorBufferGPU::getGPUBuffer()`. Device buffers grow-only
|
||||
(`EnsureCapacity`). Everything runs on the passed-in `CudaStream`, one sync at the end.
|
||||
|
||||
## Build & test
|
||||
|
||||
```
|
||||
cd build && cmake . # pick up the new CMake sources
|
||||
make -j$(nproc) jfjoch_test
|
||||
cd tests
|
||||
./jfjoch_test "BraggIntegrationEngineGPU_MatchesCPU" # GPU==CPU, 4 modes (skips if no GPU)
|
||||
./jfjoch_test "[bragg_bench]" # perf sweep vs reflection count
|
||||
```
|
||||
|
||||
The equivalence test compares I / sigma / bkg within a small tolerance (float vs double + unordered
|
||||
atomic summation of the learned profile → ~1e-4, not bit-exact — unlike the integer ROI test).
|
||||
|
||||
## Performance (measured on a throttled RTX PRO 2000 Blackwell **laptop**, WSL2 — worst case)
|
||||
|
||||
~0.5 ms @528 refl, 1.1 ms @2600, 2.3 ms @5325 (stress). ~0.45 ms of that is fixed WSL2 kernel-launch
|
||||
latency; on a native-Linux datacenter GPU it is sub-millisecond throughout. 15–20× vs the CPU engine.
|
||||
(First cut was 10.8 ms with FP64 → 3.8 ms with float → the above after collapsing ~19 CUDA API
|
||||
calls/frame into one `reset` kernel + inline shell.)
|
||||
|
||||
## To integrate (the binding — do this on the target machine)
|
||||
|
||||
The engine is deliberately *not* called anywhere yet. Wiring it in:
|
||||
|
||||
1. **Provide the preprocessed buffer + stream to the integration call.** Today
|
||||
`MXAnalysisWithoutFPGA::Analyze` (`image_analysis/MXAnalysisWithoutFPGA.cpp`, ~line 101-104) wraps
|
||||
`preprocessor_buffer->getBuffer()` back into a `CompressedImage` and hands it to
|
||||
`indexer.ProcessImage(...)`. To use this engine you instead want to pass the
|
||||
`ImagePreprocessorBuffer` (the GPU subclass carries the device image) and the `CudaStream` through
|
||||
to where integration happens.
|
||||
2. **Swap the integrator call.** `IndexAndRefine.cpp` (~line 337-339) currently does:
|
||||
```cpp
|
||||
i_outcome.reflections = mode == IntegratorMode::BoxSum
|
||||
? BraggIntegrate2D(...) // CompressedImage
|
||||
: ProfileIntegrate2D(...);
|
||||
```
|
||||
Replace with a `BraggIntegrationEngine::Run(buffer, prediction.GetReflections(), nrefl, msg.number)`
|
||||
on a member engine constructed once (GPU when `get_gpu_count() > 0`, else CPU — mirror how
|
||||
`MXAnalysisWithoutFPGA` chooses `AzIntEngineGPU` vs `AzIntEngineCPU`). Construct with the same
|
||||
`stream` the other GPU engines share.
|
||||
3. **Keep the old functions** (`ProfileIntegrate2D`/`BraggIntegrate2D`) until the engine is validated
|
||||
end-to-end against them on real data (`rugnux --integrator ... --dump-observations`, A/B vs
|
||||
`XDS_ASCII.HKL`), since the `±1`-band difference above is the one behavioural change.
|
||||
|
||||
The output `vector<Reflection>` is identical in shape (I, sigma, bkg, partiality, d, ...), so
|
||||
`ScaleImage`/`CalcISigma`/`Combine3D`/merge downstream consume it unchanged.
|
||||
@@ -25,9 +25,10 @@
|
||||
// ProfileEmpirical-> per-shell learned empirical profile fit
|
||||
// The box sum is also the seed pass (Pass A) of the two profile modes, so it always runs.
|
||||
//
|
||||
// This class is intentionally standalone: it is NOT yet wired into IndexAndRefine. It takes a
|
||||
// preprocessed image + the predicted reflections and returns the same vector<Reflection> shape
|
||||
// (I, sigma, bkg, partiality, ...) that the downstream scaling/merge consumes unchanged.
|
||||
// This is the Bragg integrator used by the pipeline (bound in MXAnalysisWithoutFPGA: the GPU
|
||||
// engine when a device is present, otherwise the CPU engine). It takes a preprocessed image +
|
||||
// the predicted reflections and returns the vector<Reflection> (I, sigma, bkg, partiality, ...)
|
||||
// that the downstream scaling/merge consumes unchanged.
|
||||
// =============================================================================
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@@ -1,391 +0,0 @@
|
||||
# Next-generation integrator & the XDS-quality investigation
|
||||
|
||||
_Last updated 2026-06-25. Supersedes the 2026-06-24 design draft — `ProfileIntegrate2D` is now
|
||||
built, wired, and the default integrator; this document records the design **and** everything the
|
||||
data-quality investigation established (especially the dead ends, so they are not re-tried)._
|
||||
|
||||
---
|
||||
|
||||
## TL;DR
|
||||
|
||||
- **The question that started this:** jfjoch's merged data was ~3–4× worse than XDS on per-observation
|
||||
precision (R-meas, ISa, anomalous peak height) on the HEWL test crystals. Where is the gap, and can
|
||||
we close it **without** XDS-style 3D shoeboxes (Filip's philosophy: integrate per-frame, combine at
|
||||
the end, keep ~300 Hz throughput)?
|
||||
- **Diagnosis (decisive):** our **mean intensities are accurate everywhere** (`CC(ours, XDS) = 0.994`,
|
||||
≥0.986 in every resolution shell) and **unbiased** (ANODE: our anomalous peaks sit on the same real
|
||||
atoms as XDS). The gap is **per-observation noise**, localized to the **2D box-sum integration** of
|
||||
**strong** reflections — a ~18% multiplicative floor that a fixed disk imposes by capturing a
|
||||
width-dependent fraction of each spot.
|
||||
- **What moved the needle:** swapping box-sum for a **profile-fitting integrator** (`ProfileIntegrate2D`,
|
||||
Gaussian profile). That is the single lever that ever improved intensity quality. Everything on the
|
||||
*scaling* side was neutral or negative.
|
||||
- **Current best on the test crystal** (HEWL `fixed_master.h5`, anomalous S/Cl peak height as a fraction
|
||||
of XDS, `1.0` = match): **0.61×**, ISa **11.0**, via Gaussian profile-fit + `-R100` + `rot3d` +
|
||||
`--scale-fulls`. Box-sum baseline was **0.46×**, ISa 8.4. XDS is 1.0× / ISa ~28.
|
||||
- **The remaining gap is XDS's full XSCALE/CORRECT scaling machinery**, whose individual pieces are each
|
||||
neutral on this well-behaved crystal — not a quick win, and not a bias to fix.
|
||||
|
||||
---
|
||||
|
||||
## Philosophy & constraints (Filip)
|
||||
|
||||
1. **Per-frame integrate, combine at the end.** Jungfraujoch integrates each frame as a 2D image and
|
||||
combines the per-frame partials of a rotation reflection into one "full" afterwards. We deliberately
|
||||
do **not** batch frames into XDS-style 3D shoeboxes in the online path. Push per-frame + combine as
|
||||
far as it goes.
|
||||
2. **Throughput matters.** ~300 Hz end-to-end was demonstrated; that is an order of magnitude over XDS.
|
||||
Judge every improvement by its **performance penalty**, not just quality (this is why global scaling
|
||||
was removed — see below).
|
||||
3. **Modular, not bolted-on.** The integrator is swappable (`BoxSum` ↔ `ProfileGaussian` ↔
|
||||
`ProfileEmpirical`) behind a single setting; scaling, the rot3d combine, and the merge consume a
|
||||
plain `Reflection{I, σ, partiality, d}` and are untouched by the integrator choice. Both **stills and
|
||||
rotation** benefit.
|
||||
4. **Simple, readable code.** No speculative/defensive machinery; neutral experiments get reverted.
|
||||
|
||||
---
|
||||
|
||||
## The pipeline today
|
||||
|
||||
```
|
||||
per-frame 2D integration (BraggIntegrate2D | ProfileIntegrate2D) <- the lever
|
||||
-> per-image scaling (ScaleOnTheFly: Ceres fit G, B_resolution, mosaicity, wedge)
|
||||
-> rot3d combine (Combine3D: weight-sum per-frame partials -> one full per event)
|
||||
-> scale fulls (optional) (ScaleOnTheFly again, Unity model, on the combined fulls)
|
||||
-> merge + error model + stats (MergeOnTheFly: RefineErrorModel, CC1/2, R-meas, CCref, ISa)
|
||||
```
|
||||
|
||||
Relevant knobs (`rugnux`, mirrored in the viewer's Processing-settings tabs):
|
||||
- `--integrator boxsum|gaussian|empirical` — default **gaussian** (`IntegratorMode::ProfileGaussian`).
|
||||
- `-P fixed|rot|rot3d` — partiality model. `rot3d` = `PartialityModel::Rotation` **+** the orthogonal
|
||||
`ScalingSettings::combine_3d` bool (see "Why rot3d is a bool" below).
|
||||
- `--scale-fulls` — second scaling pass on the combined fulls (XDS order). Implies `-M`, needs `rot3d`.
|
||||
- `-R<N>` — first-pass rotation-indexing images; default **100**.
|
||||
- `-M` — run scaling; `-A` — anomalous (keep Friedel mates separate).
|
||||
|
||||
### Current integrator: `ProfileIntegrate2D`
|
||||
|
||||
Drop-in for `BraggIntegrate2D`, same `vector<Reflection>` output, **no reference intensities**.
|
||||
Selected at `IndexAndRefine.cpp:347`. Per frame:
|
||||
|
||||
- **Pass A — rough + pick.** Box-sum every predicted reflection (rough `I`, background, centroid);
|
||||
select strong spots (significance ≥ 5).
|
||||
- **Pass B — build the profile**, per resolution shell, from the strong spots:
|
||||
- `ProfileGaussian` (**default, the keeper**): a Gaussian of the measured second moment.
|
||||
- `ProfileEmpirical`: the averaged, centroid-aligned, background-subtracted pixel grid.
|
||||
- **Pass C — profile-fit extract (Kabsch).** `I = Σ P(c−B)/v ÷ Σ P²/v`, de-biased variance
|
||||
`v = B + max(I,0)·P` (iterated), `σ = √(1/Σ P²/v)`. Rotation `partiality` carried exactly as
|
||||
`BraggIntegrate2D` (so `rot3d` works unchanged).
|
||||
|
||||
**Why it removes the floor:** profile weighting integrates out to the spot tail (the width-varying part
|
||||
a fixed disk clips) while down-weighting the noisy edge (avoiding the wide-box background penalty). The
|
||||
floor is multiplicative and lives in strong reflections, so matched weighting + de-biased variance is
|
||||
the right lever — **not** a tighter aperture (a tight profile *lost* to a generous box for stills, the
|
||||
0.4-px centroid-undersampling floor).
|
||||
|
||||
**Measured, do not re-try:**
|
||||
- **Anisotropy adds nothing — SUPERSEDED 2026-06-30 (round 3 below).** The original test was confounded
|
||||
by an over-wide (contaminated) width; once the width is de-contaminated, a **radial** elongation does
|
||||
help at high resolution. Per-detector-region profiles still add nothing.
|
||||
- **The empirical profile under-performs the Gaussian** as built (per-frame, integer-pixel-binned,
|
||||
centred on `round(predicted)`): sub-pixel smear + too few strong spots per frame/shell. It would need
|
||||
sub-pixel alignment + cross-frame accumulation to compete; not currently worth it.
|
||||
- _(Note: the header comment in `ProfileIntegrate2D.h` still says "DRAFT / not yet wired in" — that is
|
||||
stale; it is wired and default.)_
|
||||
|
||||
### The rot3d combine (`Combine3D`)
|
||||
|
||||
Buckets per-frame partials by raw `(h,k,l)`, splits into contiguous-frame rocking events (frame gap ≤2),
|
||||
inverse-variance weight-sums each event's corrected partials into **one full**
|
||||
(`F = Σ wⱼ·I_corr,j / Σ wⱼ`, `σ = 1/√Σ wⱼ`; partiality = 1, `image_scale_corr` = 1), dropping events
|
||||
with `Σ pⱼ < min_partiality`.
|
||||
|
||||
- **De-biased Poisson variance (kept):** weight each frame by `var = max(0, corr²(σ²−I)) + corr·max(0,F)`
|
||||
(background noise + **model** signal, iterated 3×) instead of `σ_corr²`. Using the model signal (not
|
||||
the down-fluctuating observed `I`) stops weak partials being over-weighted. ISa 5.3→6.4 (ref),
|
||||
6.7→8.4 (no-ref).
|
||||
- **Validation that the grouping is structurally right:** jfjoch produces ~the **same observation count
|
||||
as XDS** (936k vs 961k, multiplicity 25.1 vs 25.9). The residual gap is pure per-observation noise,
|
||||
not bad event grouping.
|
||||
|
||||
**Why `rot3d` is a `ScalingSettings::combine_3d` bool, not a 4th `PartialityModel` enum.** Many sites
|
||||
gate on `PartialityModel::Rotation` *exactly* (e.g. `GetRotationWedgeForScaling`, rotation prediction).
|
||||
A 4th enum value silently broke per-image rotation handling → NaN scales, all partials rejected. So
|
||||
`-P rot3d` sets `PartialityModel::Rotation` **and** `combine_3d = true`; the pipeline only ever sees
|
||||
`Rotation`, and the post-pass keys off `GetScalingSettings().GetCombine3D()`. **Do not reintroduce a
|
||||
rot3d enum value.**
|
||||
|
||||
### Scale-on-fulls (`--scale-fulls`) — XDS order
|
||||
|
||||
XDS/DIALS scale on the 3D-integrated **fulls**; jfjoch scales per-frame **partials** then combines, so
|
||||
the per-frame scale `G` is fit via `G·partiality·B·lp·Itrue − Iobs` — rocking-curve/partiality **model
|
||||
error leaks into G**, and the combined fulls inherit that partial scale (never scaled *as* fulls). Fix:
|
||||
a second `ScaleOnTheFly` pass on the combined fulls with the **Unity** model (`G·Itrue − I_full`, no
|
||||
partiality), between combine and merge. **This is the only post-integration win:** anomalous 0.53→0.57×,
|
||||
ISa 9.4→10.5 — and it improves **both** precision and accuracy (unlike outlier rejection, which trades
|
||||
accuracy for CC½). Kept.
|
||||
|
||||
---
|
||||
|
||||
## The diagnosis — where the gap is (and is not)
|
||||
|
||||
The investigation localized the gap with three independent instruments (see Tooling) and a per-reflection
|
||||
comparison against `XDS_ASCII.HKL` (same 4/mmm ASU). Findings:
|
||||
|
||||
**Pipeline-wide systematic-error map (ours vs XDS, per reflection):**
|
||||
- `CC(ours, XDS) = 0.994` overall, ≥0.986 in **every** resolution shell → **mean intensities are accurate
|
||||
everywhere.** The deficiencies are systematic *scale trends*, not random error.
|
||||
- Ratio ours/XDS by resolution: 3.0 Å **0.94 (strong, −6%)** → 1.8 Å 1.02 → 1.45 Å 1.19 → 1.2 Å
|
||||
**1.27 (weak, +27%)**. Two distinct problems:
|
||||
1. **Strong low-res:** mean ~right but ~6% low (mild saturation/non-linearity at high counts) **and**
|
||||
per-observation scatter R-meas **13% vs XDS 2.9%** — averages out (CC fine) but kills **anomalous**
|
||||
(Bijvoet differences need <1% agreement). This is the integration-quality floor.
|
||||
2. **Weak high-res:** **+20–27% high = background under-subtraction.** Code-confirmed: the integrator
|
||||
subtracts a **flat ring mean** (`ProfileIntegrate2D.cpp` `out.bkg = bkg_sum/n_bkg`), so an unmodeled
|
||||
background **gradient** under the peak biases weak reflections high, worst at high res.
|
||||
|
||||
**Anomalous accuracy (ANODE, external benchmark):** averaged anomalous peak height at the model atoms,
|
||||
ours vs XDS — SD_MET 6.05 vs 13.23, SG_CYS 5.56 vs 11.79, CL_CL 4.06 vs 9.58 (XDS ~2.2× stronger). **But
|
||||
both put their strongest peaks on the same real atoms** (<0.2 Å off; non-anomalous atoms in the noise).
|
||||
→ our data is **accurate / unbiased**; the gap is **precision (noise)**, ~2.2× on the anomalous ΔF,
|
||||
consistent with R-meas (~3×) and ISa (~3.5×). We are not systematically wrong, just ~2.2× noisier.
|
||||
|
||||
**Where the strong-reflection floor lives:** it is in **strong** reflections (jf/XDS R-meas ratio 2.1×
|
||||
weak → 6.5× top-5%), a ~18% **multiplicative** floor, **flat across resolution**, **per-observation**
|
||||
(per-image/scaling is only 6.6% of 20.7%), **immune to box radius** (tight/default/wide all ≈ 18–24%),
|
||||
and present even in fully-captured spots (so not partiality/capture). It survives every box size →
|
||||
it is the **box-sum method** (uniform pixel weighting), confirmed by the profile-fit win.
|
||||
|
||||
---
|
||||
|
||||
## Experiment ledger — what was tried (the dead ends are the value)
|
||||
|
||||
Scored on the HEWL test crystal, anomalous S/Cl peak height as ×XDS (`1.0` = XDS), unless noted.
|
||||
Baseline progression that **worked**:
|
||||
|
||||
| change | effect | verdict |
|
||||
|---|---|---|
|
||||
| box-sum → **2D Gaussian profile-fit** | 0.46 → **0.53×**, ISa 8.4 → 9.4, R-meas 18.9 → 14.6% | ✅ **the lever** (committed default) |
|
||||
| **`rot3d` combine** (+ de-biased Poisson var) | ISa 1.7 → 5.3 → 6.4; R-meas (partials 67% → fulls 18.9%) | ✅ keep |
|
||||
| **`--scale-fulls`** (XDS order) | 0.53 → **0.57×**, ISa 9.4 → 10.5 (precision **and** accuracy) | ✅ keep |
|
||||
| **`-R100`** (more first-pass indexing images) | 0.57 → **0.61×**, ISa 10.5 → 11.0 (~free; sweet spot ~100) | ✅ keep (committed default) |
|
||||
| masked-sentinel fix (`±1` band) | neutral on lyso; correctness (gap-clipping reflections elsewhere) | ✅ keep (committed) |
|
||||
|
||||
Things that **did not work** (do not re-prototype without new evidence):
|
||||
|
||||
| change | result | why it failed |
|
||||
|---|---|---|
|
||||
| 3D box-sum / single 3D background | R-meas identical, ISa bump is a σ-model artifact | the 3D box-sum is reconstructible from the 2D partials; recovers nothing |
|
||||
| **3D profile-fit** (`Integrate3D`, voxel Kabsch) | 0.53 → **0.56× (+6%)**, diminishing (box→2D +15%, 2D→3D +6%) | integration cannot bridge the rest; offline-only, slower, dropped |
|
||||
| **global (joint) scaling** | 0.54× (= 0.53 alternating), **+120 s** | block-coordinate loop already reaches the joint optimum; **removed (perf)** |
|
||||
| **absorption surface** (per-image / global quadratic in detector xy) | 0.54×, ISa 9.4 → 7.3 (worse) | radially degenerate with per-HKL resolution structure → "intensity sponge"; we're unbiased, a position-scale can't help |
|
||||
| **learned rocking-curve profile** (on partials) | ISa 9.4 → 6.7; with `--scale-fulls` 0.37× | partials are already divided by the model partiality; the learned `c(t)` spikes at curve edges and up-weights noise |
|
||||
| **per-batch geometry refinement** (30° batches) | 0.61 → 0.56–0.57× | batch deltas confirm real smooth drift, **but** per-image refine wants the most *stable* start; a local-but-noisy per-batch start hurts |
|
||||
| **sinusoidal geom-drift model** (beam/distance vs φ) | ~0 drift found; 0.57–0.58× | this crystal is well-centred (beam sub-pixel, distance 0.04%); physically sound, no drifty data to demonstrate |
|
||||
| background **plane fit** (a + gₓ·dx + g_y·dy over the ring) | anomalous identical; +27% weak bias persists | the weak bias is **not** a linear gradient (it's radial / profile-vs-summation); no superfluous code |
|
||||
| per-frame **B-factor** decay (`--scale-fulls -B`) | 0.57 → 0.56× | no significant radiation damage; the scalar `G` already absorbs it |
|
||||
| outlier rejection (`--reject-outliers`, `--reject-delta-cchalf`) | CC½ ↑ **but anomalous ↓** (0.53 → 0.49×) | **the accuracy trap** — improves precision, degrades accuracy; not Bijvoet-safe. Do not enable for anomalous |
|
||||
| raise `--min-partiality` (0.02 → 0.6) | 0.53 → 0.35× | low-partiality partials *feed* the scale fit; partiality model is good, post-refining it doesn't help |
|
||||
| cut high resolution | hurts (1.2 Å > 2.0 Å) | high-res helps; not over-weighted noise |
|
||||
| more scaling iterations; error-model tuning | flat | scaling converged; the error model is a **thermometer** (fits the scatter, doesn't reduce it) |
|
||||
| beam-stop masking | anomalous unchanged | only ~5 low-res reflections affected on this well-behaved crystal |
|
||||
| classical post-refinement | already done | cell (a=78.23, c=37.78 ≈ model), orientation (XtalOptimizer), and **mosaicity** (Ceres in `ScaleOnTheFly`) are already refined |
|
||||
|
||||
**Radiation damage ruled out:** anomalous rises monotonically with dose (first 300 frames 0.24× → full
|
||||
1800 0.61×) → redundancy-dominated, not damage-limited; per-frame B fits flat.
|
||||
|
||||
**`-R100` is refinement, not indexing:** FFT saturates by ~33 first-pass images; the 10.5→11.0 gain is
|
||||
refining the **global** geometry against more spots (30→50 populated images). **Per-image refinement is
|
||||
essential** and fits real per-frame signal: `-r none` (global only) collapses to ISa 5.1 / 0.31×;
|
||||
`-r orientation` 0.56×; default `BeamCenter` (per-image orient + beam) 11.0 / 0.61×. Finer + more DOF
|
||||
wins — the per-image beam refit is really a per-frame position correction (crystal motion /
|
||||
goniometer sphere-of-confusion).
|
||||
|
||||
---
|
||||
|
||||
## Tooling (reusable instruments)
|
||||
|
||||
- **`--dump-observations <file>`** (`rugnux`): `Combine3D` writes the unmerged fulls
|
||||
(`h k l I σ d n_frames captured_fraction peak_frame`) for direct A/B vs `XDS_ASCII.HKL`.
|
||||
- **`lyso_test/anomalous_scoreboard.sh <jfjoch_anomalous.hkl> [label]`**: one command → SHELXC + ANODE
|
||||
under `qemu-x86_64-static` (the SHELX binaries hit the legacy vsyscall page; this WSL2 kernel is
|
||||
`vsyscall=none`, so run them under qemu) → prints `⟨d″/σ⟩` per shell + averaged anomalous peak height
|
||||
at SD_MET/SG_CYS/CL_CL vs a cached XDS reference. ~30 s/rerun. **This is the accuracy benchmark** —
|
||||
anomalous ΔF is tiny, so it is exquisitely sensitive to systematic error, and it caught the outlier-
|
||||
rejection accuracy trap that CC½/R-meas hid.
|
||||
- **Per-reflection compare** (`ccp4-python` + gemmi, P43212 ASU): `CC(ours, XDS)` per shell and the
|
||||
ours/XDS ratio-vs-resolution that produced the systematic-error map.
|
||||
- **`[timing]` log lines** in the scaling/merge post-pass (per sub-step, plus the serial `MergeAll` vs
|
||||
parallel `ScaleAllImages` split) — to find performance bottlenecks.
|
||||
|
||||
---
|
||||
|
||||
## Performance notes
|
||||
|
||||
- **Scaling is already multi-threaded.** `ScaleOnTheFly` parallelizes across images (all cores); each
|
||||
per-image Ceres solve is single-threaded (correct grain).
|
||||
- **Global (joint) scaling was removed** (2026-06-25): 122 s vs 3.5 s for the alternating loop (~35×) for
|
||||
**no** quality gain (0.54 ≈ 0.53). Deleted everywhere (`GlobalScale.{h,cpp}`, the CLI flag, the
|
||||
`ScalingSettings` field, the viewer checkbox).
|
||||
- **The serial post-pass bottlenecks** (the real MT targets, if pursued): `Combine3D` (~1.8 s),
|
||||
`MergeAll` (~1.4 s × scaling iterations), the final merge (~0.5 s).
|
||||
- **Latent perf trap (open):** stats computation blows up to ~110 s when **no space group is set**
|
||||
(huge unmerged set → `SearchSpaceGroup` auto-detection over all of it; with a SG, 0.14 s). Needs a
|
||||
guard/cap — not yet touched.
|
||||
|
||||
---
|
||||
|
||||
## Where the remaining gap is — honest assessment
|
||||
|
||||
Box-sum → 2D profile (+15%) → 3D profile (+6%) shows integration improvements with **diminishing
|
||||
returns**; integration alone cannot reach 1.0×. ANODE says we are unbiased, just ~2× noisier. Every
|
||||
isolated *scaling* improvement (global fit, absorption surface, per-frame B, more iterations, error-model
|
||||
tuning) was neutral or negative on this well-behaved crystal. The conclusion: **the residual ~1.8× is
|
||||
XDS's full XSCALE/CORRECT machinery acting together** (multi-dimensional smooth scale + per-batch
|
||||
corrections + tighter error model, fit jointly), whose individual pieces are each neutral here — a large
|
||||
reimplementation with low odds of a quick win.
|
||||
|
||||
**Candidate directions, with current expectation:**
|
||||
1. **Weak high-res background modelling** (the +20–27% bias) — radial / profile-vs-summation, **not** a
|
||||
linear gradient (the plane fit was neutral). A genuine target, but the simple fix doesn't work.
|
||||
2. **Strong-reflection per-observation profile scatter** — the integration-quality floor; 2D profile-fit
|
||||
already took the easy part, 3D profile gives a diminishing +6% and needs a 3D pixel shoebox (departs
|
||||
from the per-frame philosophy; reader-blocked online).
|
||||
3. **A learnable model.** The systematic-error map is, in effect, the "what to do" answer for an ML
|
||||
integrator: the weak background bias, the strong-reflection saturation/response, and the per-observation
|
||||
profile scatter are all **model-learnable** corrections — none is a simple hand-coded fix.
|
||||
|
||||
**What is settled and should not be re-litigated:** we are unbiased (ANODE); 3D summation/background is
|
||||
null; global scaling, absorption surface, per-batch geometry, learned rocking profile, per-frame B,
|
||||
outlier rejection, partiality cuts, and beam-stop masking are all dead on this data. The lever that ever
|
||||
worked is **2D profile-fit integration**; the keepers are **Gaussian profile-fit + rot3d + scale-fulls +
|
||||
-R100** (≈ 0.61× anomalous, ISa 11.0).
|
||||
|
||||
---
|
||||
|
||||
## Round 3 (2026-06-30): the profile width was wrong — de-contaminate + radial parallax ellipse
|
||||
|
||||
Re-targeting **weak high-resolution** reflections (ligand-screening; metric = high-res CC1/2 / CCref, the
|
||||
only metric that can exceed XDS), a profile dump (per strong spot: the 13×13 background-subtracted,
|
||||
intensity-normalised grid + per-shell learned width, in the integrator's **own** units) showed the fitted
|
||||
Gaussian was **3–8× wider than the true spot**. The width came from the full-13×13 second moment, which is
|
||||
(a) **neighbour-contaminated** — the learning grid had no neighbour mask, catastrophic at low res where
|
||||
spots crowd near the beam (32× at shell 0) — and (b) inflated by **rectified background noise** in the far
|
||||
corners (lever arm `dx²+dy²` up to 72).
|
||||
|
||||
Splitting the spot 2nd-moment into **radial vs tangential** (about the beam) separated two effects cleanly:
|
||||
the **tangential (isotropic)** width grows with resolution and saturates (mosaicity/divergence ∝ R_px),
|
||||
while the **radial excess is pure sensor parallax**, ∝ `tan²(2θ)`, giving σ_z ≈ 106 µm for 450 µm Si
|
||||
(physical: below the `t/√12 = 130 µm` uniform limit). The contaminated isotropic width was accidentally OK
|
||||
at high res (wide) but far too wide at low res — so simply narrowing it helped mid-res but *truncated* the
|
||||
radial parallax streak and hurt the weakest high-res shell.
|
||||
|
||||
**The fix (shipped):**
|
||||
1. **De-contaminate the width** — measure the 2nd-moment over the signal disk **`r < r1`** (the same disk
|
||||
the intensity is summed over), which inherently excludes neighbours and caps the radial tail. Gated to
|
||||
the **monochromatic path**; the broadband (stills) path keeps the generous full-grid width (sparse
|
||||
spots, the centroid-undersampling floor wants it). A neighbour-mask on the learning grid or an
|
||||
iterative `r<kσ` window were considered but the disk cut is simpler and the iterative one *diverges*
|
||||
at high res (the isotropic window chases the radial streak).
|
||||
2. **Radial parallax ellipse, always-on** — extend the (was stills-only bandwidth) ellipse with an
|
||||
**analytic, material-aware** parallax term `c_par·tan²(2θ)`, `c_par = Var(z)/pixel²` from sensor
|
||||
thickness + material + energy (`parallax_var_px2`; Si 450 µm → 0.66, JUNGFRAU Si 320 µm → 1.42, CdTe
|
||||
tiny). Plus a fixed **weak-spot capture** term (`C_CAPTURE = 2.5`, monochromatic only) absorbing the
|
||||
per-frame radial position scatter of weak high-res spots; the metric is a broad plateau (parc 2–4 ≈
|
||||
equal) so it needs no tuning. On the broadband path the bandwidth streak supplies this, so it is omitted.
|
||||
|
||||
**Result (HEWL rotation `fixed_master.h5` @ 1.0 Å, vs the round-2 baseline):** ISa **13.5 → 15.7**, CC1/2
|
||||
1.12 Å **91.3 → 95.9**, 1.05 Å **83.0 → 85.2**; external **CCref** band 1.00–1.19 Å **88.1 → 89.9**, CCxds
|
||||
**93.4 → 94.8** (per-shell CCref 1.00 Å 60→66, 1.05 Å 82→85.5); low/mid res flat, overall R-meas 9.4 → 8.7.
|
||||
Generalises: sharp serial stills (Jet8, JUNGFRAU 320 µm) get a small clean gain from the parallax term;
|
||||
broadband stills (Jet5, 1% DMM) are neutral (bandwidth dominates). The de-contamination is rotation-only
|
||||
because it *hurts* stills (their ellipse is on; tightening the seed starves it).
|
||||
|
||||
This is the biggest single integration lever since box→Gaussian, and it confirms parallax is real and
|
||||
analytic.
|
||||
|
||||
**Sub-pixel centring (added).** The profile is now built per reflection and centred on the **predicted
|
||||
sub-pixel position** (the learning/fit grid is otherwise binned about `round(predicted)`, mis-placing a
|
||||
shared profile by up to 0.5 px; the predicted position is noise-free geometry, unlike the observed
|
||||
centroid, which hurt). Small clean win on every metric — rotation ISa 15.7→16.2, CCref band 89.9→90.0,
|
||||
CCxds 94.8→95.0 (high-res 1.00 Å CCref 66.0→66.9), sharp serial stills 1.68 Å CC1/2 61.6→62.5, and the
|
||||
**anomalous S peak 0.92×→0.92× XDS** (SG_CYS/CL_CL nudge up) — i.e. it does not trade accuracy for it.
|
||||
The companion idea of **de-broadening** the learned width by the known 1/12 px²/axis integer-binning
|
||||
floor was tested and **rejected** (neutral-to-slightly-negative: the learned width was not over-broad
|
||||
enough to matter, so subtracting 1/12 just over-narrows).
|
||||
|
||||
For the record, the combined round-3 lift on the **anomalous** accuracy benchmark (which the re-target had
|
||||
de-prioritised): S peak height **0.85× → 0.92× XDS** (pre-fix → profile-width fix → +sub-pixel), the
|
||||
closest to XDS yet — confirming the weak-high-res work raised genuine accuracy, not just precision.
|
||||
|
||||
Open follow-up (deferred, untested): a neighbour mask on the learning accumulation as the root-cause
|
||||
alternative to the `r1` window (keeps this spot's genuine [r1,r2] tail; needs a multiplicity mask to
|
||||
separate self from neighbour in the overlap annulus).
|
||||
|
||||
### Reciprocal-space profile width — `--reciprocal-profile` (tried and REMOVED 2026-07-02)
|
||||
|
||||
**The idea.** The per-shell width is learned in **pixels**, so it changes ~4× with resolution (the
|
||||
geometric projection of a roughly-constant reciprocal-space relrod) and must be binned per shell, which
|
||||
**starves** at high resolution / on sparse data. Mapping the tangential 2nd-moment into **reciprocal
|
||||
space** (Jacobian `g_tan² = cos²2θ`) should remove the projection and give a transferable mosaicity law
|
||||
`σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`), fit once per frame — like XDS/DIALS
|
||||
integrating in a reciprocal profile coordinate. It was shipped as an off-by-default dial after being
|
||||
tuned on the single sharp HEWL crystal, where it was metric-neutral.
|
||||
|
||||
**Why it was removed — an 11-crystal mosaicity-stratified re-test (`/data/rotation_test`, off vs on vs a
|
||||
de-contaminated variant, plus a per-frame dump of the fitted widths) showed it is net-negative and cannot
|
||||
work per-frame.** Three linked, proven facts:
|
||||
1. **The `C|q|²` mosaicity term — the whole point — is unfittable per-frame.** The fitted curvature `a2`
|
||||
came out ≈ 0 (often slightly *negative*) on **every** crystal, with **zero correlation** to the XDS
|
||||
mosaicity (0.09→0.42°). Strong spots live at low `q` (q̄ ≈ 0.24–0.47) where `η²q²` is invisible; the
|
||||
curvature only appears at high `q`, where there are ~0 strong spots to constrain it. The law degenerates
|
||||
to a straight line.
|
||||
2. With `a2 ≈ 0`, the high-res width becomes a **blind `1/cos²2θ` extrapolation → 2–4× wider than
|
||||
per-shell** (e.g. px²: pding4_001 1.07 flat → 2.0–3.0; cytC 1.10 → 2.7; MyoB 0.89 → 3.4). Meanwhile the
|
||||
per-shell path *starves* at high res (n_strong<30 ⇒ flat global fallback that even **shrinks** vs
|
||||
mid-res) — which turns out to be **accidentally correct**: weak, crowded high-res spots want a *narrow*
|
||||
aperture, not the true (wide) spot shape.
|
||||
3. The over-wide profile pulls background into weak spots ⇒ **per-observation scatter (R-meas) rises**,
|
||||
CC1/2 drops in reliable, high-multiplicity shells (pding4_001 mult 25, pding4_003 mult 12, EcwtCQ066,
|
||||
MyoB); the mean intensities survive (CC-overall flat). A confirmatory cap (clamp the reciprocal width at
|
||||
the widest well-sampled per-shell value) **recovers the regression** — direct proof the over-widening is
|
||||
the harm. The only "improvements" seen (EcwtAL500, EP_cs_01-17 overall CC) were entirely in
|
||||
noise-dominated shells (mult 2–3, CC<20 %, R-meas>100 %). **No crystal reliably wins.**
|
||||
|
||||
The per-frame paradigm structurally cannot fit mosaicity, and un-starving the high-res width is net harm,
|
||||
so the dial (CLI flag, `BraggIntegrationSettings::reciprocal_profile`, and the per-frame fit block) was
|
||||
deleted. A real win would need cross-frame / rolling-window pooling to constrain `a2`, or a two-pass
|
||||
mosaicity seed from scaling — both against the per-frame / online + stills philosophy. Also dropped along
|
||||
the way: a separate **radial** relrod fit (radial is divergence + rotation-smear, not mosaicity → no gain).
|
||||
|
||||
## Lineage: PixelRefine (removed 2026-06-25)
|
||||
|
||||
`ProfileIntegrate2D` is the surviving half of an earlier experimental still-integrator, **PixelRefine**,
|
||||
which did three things bundled together: a reference-free **profile + extraction** (`J = Σ P(I−B)/v ÷
|
||||
Σ P²/v`, de-biased variance, plus a radial partiality) and a reference-based **per-image scaling**
|
||||
(`J ≈ G·B_DW·p·pol·I_ref`). The decomposition observation was that the first two are independent of the
|
||||
reference and identical in spirit to what a profile-fitting integrator does — so they were lifted out as
|
||||
`ProfileIntegrate2D`, while scaling stayed in the shared `ScaleOnTheFly` (which already does reference
|
||||
scaling when a reference is supplied: `IndexAndRefine::ReferenceIntensities` builds a
|
||||
`ScaleOnTheFly(experiment, reference)` and the per-image pass calls it on whatever the integrator
|
||||
produced — **integrator-agnostic**).
|
||||
|
||||
**The decomposition won, so PixelRefine was deleted.** On the LysozymeJet5 serial stills, the default
|
||||
Gaussian path (`ProfileIntegrate2D` + the same reference scaling) **matched or beat** whole-PixelRefine
|
||||
in every per-shell CC½ (overall 95.7% vs 91.9%), ISa (1.6 vs 1.2) and R-meas (98.5% vs 175%), with
|
||||
CCref a tie — same indexing, same uniques. PixelRefine collected ~1.8× more observations (a wider Ewald
|
||||
band) but they were noisier and didn't help. Whole-PixelRefine was also incompatible with the rot3d
|
||||
combine (it emitted `image_scale_corr = 1/G` with no rotation partiality → garbage σ), which the
|
||||
decomposed integrator avoids by emitting `partiality` like `BraggIntegrate2D`.
|
||||
|
||||
**Lessons inherited from PixelRefine (kept here so they survive the deletion):**
|
||||
- **Mean (not median) background** in both integrators — the median of a skewed background
|
||||
under-subtracts and biases weak intensities positive (fake ⟨I/σ⟩ at CC≈0). Single biggest
|
||||
"untrustworthy σ" cause.
|
||||
- **De-biased variance** (weight by background, not observed counts) — fixed negative high-res ⟨I/σ⟩
|
||||
and scale collapse.
|
||||
- **A tight profile loses to a generous box for stills** — the ~0.4-px centroid-undersampling floor of
|
||||
~2×2 spots is a real sampling floor, not a fixable misprediction; recentring/adaptive-width add noise.
|
||||
The lever is matched weighting + de-biased variance, never a tight aperture.
|
||||
- **Per-image R (profile-width) refinement is futile** — G↔R is intrinsically degenerate; refine G/B
|
||||
freely but never R per image (it slides the per-image scale and wrecks the merge).
|
||||
Reference in New Issue
Block a user