Jungfraujoch/docs/CPU_DATA_ANALYSIS_IMAGE.md

Data analysis: from images to spots (§0–§3)

Part of the CPU/GPU data-analysis reference; the section numbers are continuous across its four parts.

:local:
:depth: 2

0. Getting the image onto the GPU: device-side bitshuffle+LZ4 decoding

Images arrive bitshuffle+LZ4 compressed (HDF5 filter 32008), and everything from §1 onwards runs on the GPU when one is present. Instead of decompressing on the host and uploading the image, the compressed chunk is uploaded — a few MB rather than tens of MB — and decoded on the device. The approach follows Jon Wright (ESRF); the kernels are Jungfraujoch's own.

One kernel does the work: one CUDA block owns one bitshuffle block, from the compressed payload through to finished pixels.

1. LZ4 into shared memory, one warp per bitshuffle block. Blocks are independent, so the parallelism is across them; within the warp every lane runs the same sequence parser over the same bytes, and the literal and match copies are split across the 32 lanes so the stores coalesce. An overlapping match is treated as a pattern of period offset sourced from bytes that already precede the write position, which keeps it parallel rather than a serial byte loop; offset == 1 (a run of one repeated byte, the common case in sparse detector data) and power-of-two offsets avoid the modulo altogether. Because the lanes cooperate on the copies, each one is followed by __syncwarp() — a later match can read bytes another lane wrote, and since Volta that ordering is not implicit.
2. The bitshuffle inverse fused with preprocessing. The whole CUDA block then reads that shared buffer back: one thread owns one group of 8 elements across every byte-plane, so once it has transposed its 8 bytes out of each plane it holds 8 complete elements — and it applies the pixel mask, the error marker and the saturation cap and emits 8 finished int32 pixels directly. Nothing of the block reaches device memory but the pixels — neither the bitshuffled bytes nor the decompressed image is ever materialised. For 8-bit images there is a single plane and the assembly degenerates to a copy.

Decoding into shared memory is worth more than the bandwidth it saves: an LZ4 match reads back bytes written a few sequences earlier, so every copy step is a dependent round trip — tens of cycles in shared memory against hundreds in device memory. It is paid for in residency, because the buffer holds a whole bitshuffle block, and a block larger than 16 kB (larger than either writer this pipeline reads produces) falls back to a pair of kernels instead, the first writing the shuffled image to device memory and the second un-transposing and preprocessing out of it.

The block offsets inside the container can only be discovered by reading the block lengths in order, so that scan stays on the host.

Only BSHUF_LZ4 is decoded on the device. For the zstd variants (BSHUF_ZSTD, BSHUF_ZSTD_RLE, BSHUF_ZSTD_RLE_HUFF), and for uncompressed or float images, BSLZ4DecoderGPU::Supports() returns false and the pipeline decompresses on the host and uploads as before.

The container arrives off the network or off disk and is not trusted. Everything checkable on the host — declared sizes, the block scan, a block size that is not a multiple of 8 elements, a block count the chunk could not hold, trailing bytes — is rejected before any work is queued; the kernel additionally flags a block that did not decode to exactly its declared length, which becomes an exception once the caller has synchronised. That last check matters because the decode buffers are reused frame to frame: a block that stopped early would leave the previous image's most significant byte-plane in place, which reads not as a missing corner but as real pixels several powers of two too bright.

1. Geometry, reciprocal-space mapping, and basic quantities

1.1 Coordinate conventions

For a pixel coordinate (x,y) (in pixels), Jungfraujoch converts to a laboratory direction vector via:

1. shift by the beam-centre pixel (x_\mathrm{beam}, y_\mathrm{beam}) — the PONI, pyFAI's point of normal incidence, which coincides with the direct-beam impact point only for an untilted detector (see Detector geometry),
2. scale by pixel size p (mm),
3. set detector distance D (mm),
4. apply detector orientation rotation R_\mathrm{det} (PyFAI-like parameterization).

The unnormalized detector coordinate (mm) is: $ \mathbf{r}\mathrm{det}(x,y) = \begin{pmatrix} (x-x\mathrm{beam})p\ (y-y_\mathrm{beam})p\ D \end{pmatrix}. $

The lab-frame vector is: $ \mathbf{r}\mathrm{lab} = R\mathrm{det},\mathbf{r}_\mathrm{det}. $

By this construction (x_\mathrm{beam}, y_\mathrm{beam}) maps to (0,0,D) before the rotation — the point the detector normal through the sample meets — which is what makes it the PONI rather than the direct beam; the two differ by D\tan(\mathrm{tilt}) on a tilted detector. The laboratory frame is fixed the same way everywhere in the system: +z along the beam propagation, +x along increasing pixel column (the fast axis) and +y along increasing pixel row (the slow axis) — a right-handed triple that coincides with XDS's laboratory frame, which is what makes the geometry echo of rugnux directly comparable. The absolute hand of an indexing — and with it the Bijvoet hands of §14.5–§14.6 — follows from this convention.

Let the incident wavevector magnitude be k = 1/\lambda in Å^{-1}, and define: $ \mathbf{S}_0 = (0,0,k). $

The reciprocal-space scattering vector associated with pixel (x,y) is: $ \mathbf{s}(x,y) = k,\frac{\mathbf{r}\mathrm{lab}}{\lVert \mathbf{r}\mathrm{lab}\rVert} - \mathbf{S}_0. $

This \mathbf{s} is the fundamental quantity used for spot finding (resolution filters), indexing, and refinement.

1.2 Two-theta, azimuth, resolution and q

The scattering angle 2\theta is computed from \mathbf{r}_\mathrm{lab} via: $ 2\theta = \mathrm{atan2}!\left(\sqrt{x_\mathrm{lab}^2 + y_\mathrm{lab}^2},; z_\mathrm{lab}\right), $

evaluated as a two-argument arctangent, so the mapping stays correct where a strongly tilted detector's far corner reaches past 2\theta = 90°.

Resolution (Å) at a pixel is: $ d = \frac{\lambda}{2\sin\theta}. $

The magnitude q = 2\pi/d is used for radial binning and ice-ring handling.

1.3 Distance from the Ewald sphere

For a reciprocal lattice point \mathbf{p} (Å^{-1}), define: $ \Delta_\mathrm{Ewald}(\mathbf{p}) = \lVert \mathbf{p} + \mathbf{S}_0\rVert - k. $ 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 §11.1),
• still partiality option in scaling/merging (§10.2).

1.4 Measuring the direct beam before indexing

The beam centre in the file is often a placeholder, and nothing else measured it until post-refinement (§7.5) — by which time a wrong centre has already chosen the lattice. It is now measured on every run and, separately, can be measured from the spots and committed before indexing.

Measured on every run (--beam-center-check, on by default). The isotropy of the scattered background places the centre — the leverage is the curvature of the solvent ring — and the mean image the beam-stop pre-scan already builds is the array the fit runs over, so the measurement costs no frames of its own. The run reports what the file claims, what the background says, how well the fit knows its own answer, and how far this particular geometry can afford to be wrong: the smearing of the accumulated reciprocal-space cloud by a centre error \delta p multiplies the FFT amplitude at an axis of length a by J_0(2\pi\,\delta p\,a/(D\lambda)), and past the first zero the true axis is gone and its harmonic wins. That tolerance runs from well under one pixel to several, so a flat pixel bound means different things on different geometries.

Nothing is committed on the strength of it. The rotation first pass is simply run a second time at the measured centre — the spots are cached and do not depend on the centre, so this costs a fraction of the first pass — and the two lattices are compared. The measured centre is adopted only where the file's centre indexes nothing and the measured one indexes a majority; where both work and disagree, both cells are reported and neither is chosen, because the only arbiter available at that stage is the indexed frame count and it points the wrong way: acceptance is a fractional-Miller test, so a cell twice as long must place every spot twice as accurately to score the same, and a halved axis can index more frames than the true cell. Where the two disagree only in Bravais class at the same primitive volume, that is said separately from a volume ratio that is an axis harmonic, which is the one the centre decides.

Two further things are said before a frame is read, both free: a centre that lands on a masked pixel is wrong about something, and a beam-stop pre-scan that finds exactly zero shadowed pixels is not a beamline without a beam stop but a centre far enough out that every flood seed fell in a gap.

Two discrete rescues after a first pass that failed. A rotation sweep leaves two things the input often cannot settle, both decidable from the data on the same count the rest of the first pass uses — the right answer indexes and the wrong one does not.

• The rotation-axis sign. A miniCBF header names the axis but gives no direction, and an NXmx vector is only meaningful together with the detector mounting. So after a poor first pass the opposite sign is tried and whichever indexes more validation frames is kept; the decision then holds for the rest of the run, and it flips the axis, not the angles, because prediction reads the axis too. It runs before the long-axis rescue, since with the sign wrong every candidate lattice is wrong.
• The beam centre as a hypothesis (--beam-center-search, on by default, 12 px). After a pass that indexes fewer than half the validation frames, the centre is stepped a pixel at a time along both detector axes and the first rung clearing the same majority is kept. Both directions are searched deliberately: an error across the spindle collapses the indexed fraction and announces itself, while an error along it is the direction the J_0 law calls free — the transform is translation-invariant, so the peaks stay sharp, the run holds nearly every frame indexed, and the lattice fit that follows commits confidently to an axis harmonic. The step is a flat pixel: derived from the J_0 law it would have to use the cell the failed pass returned, which can be a small spurious sub-cell and asks for a step that jumps the lobe being looked for.

The two unknowns are coupled, so where the file's centre indexes nothing and the measured one does not either, the axis sign is asked again at the measured centre: each error otherwise hides the other and the run produces nothing. Where nothing works the file's centre is put back, so a run that fails for some third reason fails at the geometry it was given.

Committed before indexing (--estimate-beam-center). Here the centre is measured from spot positions alone and used in place of the file's. Two exact facts about a rotation sweep supply the two coordinates:

Friedel mates half a turn apart. The Laue condition fixes the component of \mathbf{q} along the beam, q_\parallel = -\lVert\mathbf{q}\rVert^2\lambda/2. Rotating 180° about the spindle \mathbf{m} negates the two components perpendicular to \mathbf{m} and taking -h negates all three, so together they negate only the component along \mathbf{m} and leave q_\parallel untouched. With the spindle perpendicular to the beam, -h therefore diffracts at \varphi+180° exactly where h diffracts at \varphi, and its spot sits at the mirror image of $h$'s along the spindle. This gives the beam coordinate along the spindle. Only the reciprocal lattice's centrosymmetry is needed for the geometry; Friedel's law |F(h)|=|F(-h)| is used separately, to tell a true pairing from an accidental one.

The second crossing. The same reflection meets the Ewald sphere twice, at two angles that are generally not 180° apart, differing only in the sign of the lab component perpendicular to both \mathbf{m} and the beam. This gives the remaining coordinate. The two crossings are separated by a sweep angle fixed by the reflection's own position, which is what identifies genuine pairs.

Neither observable requires a cell or an orientation matrix: each candidate pairing votes for a beam coordinate, and the true value accumulates while wrong pairings scatter. A Friedel pair needs both \varphi and \varphi+180° recorded, so a sweep of S° yields only S-180 degrees' worth of pairs — a sweep of exactly half a turn yields none, and the estimator is refused below a floor on that span.

The mirror is exact in the laboratory frame, so it is sensitive to the spindle direction. A skew of the spindle about the beam spreads the vote rather than shifting it, and is fitted alongside the centre (--no-fit-spindle keeps the axis from the file); a tilt of the spindle towards the beam is measured and reported but not applied, being confounded with the detector tilt until that is fitted too. Nothing inside the fit can tell that the vote settled on the wrong periodic maximum — every frame pair agrees with every other — so the answer is accepted only if it does not move when the search is started from a different position. Frames are sampled away from both ends of the sweep, where shutter synchronisation can spoil an image.

Where the sweep is shorter than half a turn the spot symmetry cannot be formed, and the centre is taken instead from the centroid of the radial background profile, which needs only a few images. Where neither method can measure the centre, the value from the file is kept.

1.5 Finding the beam stop

The beam stop and its holder arm shadow part of the detector. A reflection behind them is attenuated but otherwise ordinary — it integrates low, with a plausible \sigma, and no outlier test catches it — so the shadow is found and masked instead. --detect-beam-stop[=N|off] (on by default, N=60 frames; it reads its own frames, so it is not affected by, and does not affect, the beam-centre pre-scan of §1.4) projects those frames to a per-pixel mean and maximum, and writes the result into the pixel mask as bit 9, from where it excludes those pixels from every later stage.

The shadow is a place where the background is missing, so it is found by comparing each pixel's background against the background at the same radius. The mean is pooled over a small box first — a single pixel of a sparse background carries too few counts to tell a shadow from a Poisson hole, and the stop is much wider than the box — and the comparison is against the median of the pixel's own radius ring, taken over the pixels not already excluded and iterated a few times so the shadow stays out of the baseline it is measured against. A pixel is shadow when the ratio falls below 0.35, and only where the ring has accumulated enough counts for the dip to mean anything. Nothing assumes the stop and the beam are concentric, because only the per-ring comparison is used; a ring lying wholly inside the stop has no unshadowed pixel for its median to find, which is exactly the case where the comparison must fail, and such a ring is shadow in its entirety.

What survives is then shaped into a region: the low pixels connected to the beam centre, bridged across the module gaps the holder arm crosses, grown outward through the partially shadowed penumbra, closed, and with the interior of the disk filled. Last, any pixel that ever recorded a real reflection — a maximum over the frames well above background, in a cluster, so that a single-frame zinger does not count — is given back, because a beam stop cannot have blocked a reflection that was measured.

---

2. Azimuthal integration (radial profiles)

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) 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 and its square: $ 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}}. $

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

Two standard corrections are available:

(i) Solid angle / geometric correction. A flat pixel's solid angle falls off with the incidence angle \alpha between the scattered ray and the detector normal. With the in-plane detector offsets u=(x-x_\mathrm{beam})p and v=(y-y_\mathrm{beam})p — measured from the PONI (§1.1), which is what the tilt-invariance below rests on — and detector distance D, $ \cos\alpha = \frac{D}{\sqrt{u^2+v^2+D^2}},\qquad C_\Omega = \cos^3\alpha, $ applied — like the polarization term below — as a divisor (intensities are scaled by 1/\cos^3\alpha), so pixels at oblique incidence, which subtend a smaller solid angle, are boosted. Because \alpha is evaluated in the detector's own frame it is invariant under detector tilt (\mathrm{rot1}/\mathrm{rot2}/\mathrm{rot3}), matching PyFAI's solidAngleArray and MAX IV azint. It reduces to the commonly quoted \cos^3(2\theta) form only for an untilted detector, where the incidence angle coincides with the scattering angle.

(ii) Polarization correction. With polarization coefficient P (beamline dependent) and azimuth \phi: $ C_\mathrm{pol}(2\theta,\phi) = \frac{1}{2}\left(1+\cos^2(2\theta) - P\cos(2\phi)\left(1-\cos^2(2\theta)\right)\right), $ applied as a divisor to intensities (i.e. scale by 1/C_\mathrm{pol}) when enabled.

2.3 Background estimate for profiles

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).

---

3. Spot finding (strong pixels → Bragg spots)

Spot finding is a two-stage process:

1. Strong-pixel selection using intensity and/or local signal-to-noise criteria.
2. Connected-component labeling (CCL) to group strong pixels into candidate spots, followed by spot-level filtering and feature extraction.

3.1 Strong-pixel detection by local statistics

For each pixel i with value v_i, consider a square window (nominally 31\times 31 pixels) around it. Let the window contain n valid pixels (excluding masked/bad/saturated), and define: $ \Sigma = \sum v,\qquad \Sigma_2 = \sum v^2. $

To avoid biasing the local statistics by the test pixel itself, Jungfraujoch evaluates the pixel against the window with the pixel removed: $ \Sigma' = \Sigma - v_i,\quad \Sigma_2' = \Sigma_2 - v_i^2,\quad n' = n-1. $

A variance-like quantity proportional to n'^2 is formed: $ V = n'\Sigma_2' - (\Sigma')^2, $ and the deviation-from-mean quantity: $ \Delta = v_i n' - \Sigma'. $

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: $ \Delta^2 > V\cdot T^2, $ where T is the configured signal-to-noise threshold.

This is equivalent to a local z-score criterion but implemented in integer arithmetic to be robust and fast.

The test is applied in two passes over the image. The first is as described above. The second repeats it with every pixel found strong by the first excluded from the local background — it is treated exactly like a saturated pixel, so it contributes to no window it falls into and stays strong itself. This matters for any spot wide enough to reach into its own background box: on a single pass such a spot inflates the mean and variance it is then tested against, and its outer pixels fail the criterion. Excluding the core recovers them, so the spot is reported with its true extent rather than its brightest few pixels. Both the CPU and GPU implementations run these two passes and return the same spot list for the same frame.

Special cases:

• saturated pixels can be forced to “strong” (useful for detecting overloaded Bragg spots),
• invalid pixels are never strong.

3.2 Adaptive (self-calibrating) detection

The local-statistics test above needs a fixed photon/count threshold whose correct value depends on the background level, which varies between datasets. The adaptive mode (--adaptive-spots; the default in rugnux and in the viewer for both stills and rotation data, --no-adaptive-spots reverts) derives that threshold from each image's own noise, per resolution ring, so no per-dataset value is needed. It admits more spots than the fixed threshold, including genuine reflections that belong to no indexed lattice; these are down-weighted rather than filtered in the per-image geometry fit (§7.4).

Pixels are binned into the same resolution rings as the azimuthal integrator (§2). For each ring a robust background is estimated in three passes: one plain pass over all valid pixels, then two $\sigma$-clipping passes that keep only pixels within \pm 3\sigma of the current ring mean (removing the Bragg peaks from the background estimate). This yields a per-ring background mean \mu_b and scatter \sigma_b.

The ring's detection threshold is the larger of two arms, $ t_b = \max!\big(;\mu_b + z,\sqrt{\sigma_b^2 + \sigma_\mathrm{read}^2};,;; k_\mathrm{Poisson}(\mu_b, p);\big), $ where k_\mathrm{Poisson}(\mu_b,p) is the smallest count whose Poisson$(\mu_b)$ upper tail is \le p. The Poisson arm is correct where the background is countable (a bright low-resolution ring gets a high threshold); the Gaussian arm — floored by a detector-level excess-noise constant \sigma_\mathrm{read} — takes over on near-empty high-resolution rings, where the Poisson arm degenerates to "one photon is significant" and would flood. The operating point p = E/N is set from a single portable knob E, the expected number of false pixels tolerated per frame (--spot-false-pixels, default 100), with N the number of valid pixels. Because p and every \mu_b,\sigma_b come from the image itself, the same E lands a sensible photon threshold on strong and weak datasets alike, with no per-dataset tuning. Rings too sparse to characterise (fewer than ~40 pixels) fall back to a whole-frame background. A pixel is strong when v_i \ge t_b for its ring (saturated pixels are still forced strong); the strong pixels then feed the same CCL stage (§3.4). The signal-to-noise and photon-count criteria of §3.1 are not used in this mode.

Because detection reads the pixel's ring, a pixel that falls outside the azimuthal-integration q range has no ring and can never be strong: the integration range bounds what adaptive detection can see. Both upper limits are therefore optional and default to the detector itself — the azimuthal integration runs to the highest q any pixel of the detector reaches (--azim-max-q unset), and spot finding is not clipped in resolution (--spot-high-resolution unset), for rotation data as well as stills. Setting either one narrows detection accordingly — appropriate for weak, high-background data, where the spots admitted at the detector edge are dominated by noise.

Fused GPU engine. The per-ring reduction the adaptive threshold needs is the same reduction the azimuthal integrator performs. On the GPU path the two are fused into a single image pass (AdaptiveSpotFinderGPU): one reduction accumulates the corrected per-ring sums for the azimuthal profile (§2) and the raw per-ring statistics for the threshold, after which a light kernel flags the strong pixels. One GPU pass therefore replaces both the separate azimuthal-integration pass and the host-side adaptive spot-finding pass, at a small fraction of the CPU finder's cost per frame and producing the same spot list and azimuthal profile. It is enabled by default in the offline rugnux path, the interactive viewer and the online receiver.

Online. spot_finding_settings in the REST API carries adaptive_threshold and false_pixels_per_frame, so the mode is reachable from the broker and from the web frontend as well as from rugnux and the viewer. It defaults to off online, unlike rugnux and the viewer, because the broker serves both workflows and only one of them can run it: spots are found in software only on the DECTRIS/SIMPLON path, while the JUNGFRAU and EIGER workflows find them on the FPGA at its own fixed threshold. Setting adaptive_threshold on those is refused with an error rather than accepted and ignored, so a detection setting that had no effect cannot be mistaken for one that did.

3.3 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, 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 a radial profile: for each hexagonal-ice powder ring (see Where the ring positions come from below), 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 Å⁻¹ spacing — in q = 2\pi/d, like every q in this document (§1.2) — --azim-q-spacing, so the rings are well resolved). The reported quantity is the ice magnitude rather than a significance: with many photons any real ice ring is statistically significant, so significance does not discriminate.

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 a rotation battery, that alone ranked ice-free crystals above crystals that really are iced. An ice ring is azimuthally smooth and survives the sigma clip, while Bragg peaks do not, so on the clipped profile ice-free crystals sit near 1 and crystals with confirmed ice above 2. Only where no adaptive finder ran (the FPGA workflow) does the score fall back to the plain profile.

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: smooth ice reads high on the profile and ~1 on the spots, textured ice the reverse, and a clean crystal ~1 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 and the exclusion from the scale fit (§10.10) are skipped. The fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, and more than that on a detector reaching past 1.5 Å, so handling ice on a clean crystal is a pure loss.

Where the ring positions come from. The eleven bands from 3.895 to 1.522 Å are the measured positions of Moreau, Atakisi & Thorne (Acta Cryst D77, 2021, 540–554). That list ends at 1.522 Å by its own scope — the paper states that hexagonal ice "has 11 diffraction rings between 4 and 1.5 Å resolution", and its subject was detecting ice in deposited data rather than masking it — not because ice stops there. On a detector that reaches further, the rings it does not list are the ones left in the data.

The eight bands below it are calculated, since past that paper there is nothing measured to copy. Enumerating hkl from the ice Ih cell is not enough: ice Ih is P6_3/mmc with oxygen on 4f, and most of what enumeration emits is extinguished by the oxygen sublattice rather than by the space group — which is why (004) at 1.830 Å and (104) at 1.657 Å are missing from the measured list even though they lie inside its range and its reflection conditions allow them (for even l the 4f structure factor carries a factor \cos 2\pi l z at every hk, and z \approx 1/16 kills l=4 — (104) along with (004)). Structure factors are computed instead — oxygen only, the hydrogens being half-occupancy disordered and weak to X-rays — and the lines kept are those reaching 3 % of the strongest. That rule reproduces the measured eleven exactly, and every line it drops inside their range computes to zero, which is what makes it trustworthy below 1.522 Å. The cell is that of Röttger et al. (Acta Cryst B50, 1994, 644–648). The list stops at 1.170 Å because below it the calculated real lines fall to 2–3 % while the extinct ones rise to about 1 %, and an oxygen-only calculation cannot separate them any further.

One consequence is worth stating: the profile score is the strongest ring's ratio, a maximum over the bands, so a longer list can only raise it. The gate at 1.5 is therefore read against a list of this length, and lengthening it again would need the gate re-checked.

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.

3.4 Connected-component labeling (CCL)

Strong pixels are grouped into connected components (adjacent strong pixels) using a CCL algorithm. Each component yields a candidate spot with:

• centroid (x,y) (often intensity-weighted),
• pixel count (spot size),
• integrated spot intensity proxy (sum of pixel values),
• resolution d at the centroid (or mean over pixels),
• and quality flags (e.g. ice-ring classification).

Spot-level filters include minimum/maximum pixel count and resolution limits.

The upper bound is on how large a round spot may be, not on how bright one may be. Under the self-calibrating threshold (§3.2) a component's area above the contour grows as \sigma^2\ln(A/T) — without bound in the peak amplitude A — so an upper bound on pixel count alone becomes an intensity ceiling: measured on a strongly diffracting crystal, footprints run 3 px at 30–100 counts to 50 px above 10 000, four times the slope the fixed local-box test gives, and the old bound of 50 discarded the brightest reflections on every image. The bound is therefore 200 (the same as CrystFEL peakfinder8's --max-pix-count; XDS has no such parameter at all and guards on shape instead, with SPOT_MAXIMUM-CENTROID), and a component above 50 pixels must in addition fill at least a fifth of the square its bounding box fits inside. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, which is what an upper bound was ever protecting against. Below 50 pixels no shape is asked for, so everything accepted before still is. The shape test is inert on every dataset it has been measured on — it exists to bound the shape of what the larger size bound now admits, on data carrying arcs or tracks, not because the crystals measured here needed it. The test is integer arithmetic, so the host and the GPU extractor agree by construction.

The host implementation (StrongPixelSet::sparseccl) is the SparseCCL of the ACTS/traccc project: it runs over the strong pixels sorted row-major, uses a sliding window over the previous line and a union-find whose root is each component's lowest index. On the GPU the same labelling runs on the device (SpotExtractorGPU): the packed strong-pixel bitmask is compacted into that same sorted list without leaving the card, each pixel finds its at most four earlier 8-neighbours by binary search, and a lock-free union-find with path halving labels them. Only the finished spot list — a few hundred entries — comes back to the host, instead of the whole bitmask (2.26 MB per frame at 18 MP). The two implementations produce the same components, in the same order, with the same pixel counts and intensities; tests/SpotExtractorGPUParityTest.cpp holds them to it. The device version is also insensitive to frame content: the host sliding window becomes quadratic when many pixels light up in one detector line — a hot module, or a diffraction ring where it runs tangent to a row — which costs tens to hundreds of milliseconds on such a frame, while the device version stays under a millisecond.

3.5 Adaptive per-image minimum spot size

The minimum-pixels-per-spot filter (§3.4) trades sensitivity against noise: a small value keeps faint one- or two-pixel spots — real signal on strong data, but detector noise on high-background frames — while a larger value keeps only well-formed spots. The best value is dataset-dependent, so for serial-stills indexing it can be chosen per image rather than fixed. The frame is indexed three times, at min-pix 3, 2 and 1, and the setting that maximises

\frac{n_\mathrm{indexed}^2}{n_\mathrm{total}} \quad\text{(indexed-spot count weighted by indexed fraction)}

is kept; the frame is then integrated once at that min-pix. The fraction factor discounts the extra spots a smaller min-pix admits unless the lattice actually explains them, so strong frames keep their real weak spots (extending resolution) while noise-flooded frames stay strict. Because min-pix filters the connected components after detection, strong-pixel detection AND the connected-component labelling both run once per frame, and the three attempts only repeat the spot-level filter; the azimuthal profile is the one that single detection pass computed. The winning attempt's spot list is kept rather than re-extracted, so the frame that is integrated is exactly the frame that was scored. This is a stills-only, indexing-path option — rotation indexing builds one global lattice from all frames and keeps a fixed min-pix. In rugnux it is the default; giving an explicit --min-pix-per-spot pins a fixed value instead.

3.6 Predicting the resolution the merged data will reach

A per-image resolution estimate is read off the finished spot list. It predicts how far the merged data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by \sqrt{I} — the intensity is a summed photon count, so \sqrt{I} is its Poisson significance — the 1/d^2 is found beyond which a fraction f=0.30 of that weight lies, and the estimate is that resolution taken 2.25\times further in 1/d. It is deliberately not limited to what the detector records: the quantile sits in the middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal where the detector stops before the diffraction does, and on such a run it reads finer than the detector corner. The dataset value is the median over images.

Both constants carry a mechanism. A quantile from the middle of the distribution measures the shape of the fall-off, which is the crystal's own \exp(-B/2d^{2}), where the extreme end of it measures where detection stops — a threshold that moves with the exposure and with how many reflections the unit cell puts on a frame. And merging averages many observations of each reflection, so intensities go on being measurable a fixed factor in 1/d past the point at which one image's spot finder still detects them; that factor is the 2.25. Both are calibrated on rotation data against the resolution at which per-shell CC1/2 falls through 0.30, and the estimate is good to about 0.2 Å there. It is a prediction and not a measurement of what a run achieved: nothing downstream is cut on it, and it is reported alone (rugnux SPOT_RESOLUTION_ESTIMATE, and per image in the stream, the plots and HDF5).