Jungfraujoch/image_analysis/pixel_refinement/METHODS.md

PixelRefine — methods and improvements

This note documents the changes made to the still-image pixel-refinement integrator (PixelRefine) and, more importantly, why each one was needed. It is written from a methods point of view; the equations are the load-bearing part.

PixelRefine integrates Bragg reflections on a still image by fitting a per-pixel forward model against a known reference intensity set (e.g. F_calc from a deposited model). Unlike a rotation experiment, a still samples only a thin slice of each reflection, so the integrator must (a) model the partiality of that slice, (b) refine the per-image geometry well enough that the high-resolution shoeboxes land on signal, and (c) scale each image onto the reference. Each of those three is where the original code went wrong.

Throughout, a reflection's shoebox is a small box of raw detector pixels I_p with a local flat background B; the (area-normalised) model spot profile at pixel p is P_p, and v_p is the variance used to weight pixel p.

---

0. The forward model

The recorded amplitude of a still reflection is estimated by profile fitting:

J \;=\; \frac{\sum_p w_p\,P_p\,(I_p - B)}{\sum_p w_p\,P_p^{2}},
\qquad w_p = \frac{1}{v_p},
\qquad \operatorname{var}(J) = \frac{1}{\sum_p P_p^{2}/v_p}.

The full (rotation-equivalent) intensity is recovered by dividing out the factors a still does not record,

I \;=\; \frac{J}{p\,B_\mathrm{DW}\,\mathrm{pol}},\qquad
p = \exp\!\left(-\frac{\epsilon_r^{2}}{R_0^{2}}\right),\quad
B_\mathrm{DW}=\exp\!\left(-\frac{B_\mathrm{fac}}{4 d^{2}}\right),

where the partiality p is the fraction of the mosaic block crossing the Ewald sphere, \epsilon_r = 1/\lambda - |S_{hkl}| is the excitation error, R_0 the radial (rocking) width, and \mathrm{pol} the polarisation correction. The per-image model intensity for a pixel is

I_p^\mathrm{model} = G\,I^\mathrm{ref}\,B_\mathrm{DW}\,p\,P^\mathrm{tang}_p\,\mathrm{pol} + B,

with G the per-image scale. The per-image least squares minimises \chi^2 = \sum_p w_p\,(I_p^\mathrm{model}-I_p)^2 over geometry, orientation, G, R.

---

1. De-biased variance (the load-bearing fix)

Symptom. Mean intensities went negative in the high-resolution shells (\langle I/\sigma\rangle down to -12), which a box-sum integrator never does. The per-image scale G also collapsed to 0 on most images, dropping ~80 % of observations.

Cause. Both the extraction and the fit weighted each pixel by its observed count, v_p = I_p. For a background pixel that fluctuated down (I_p < B), v_p is small, so w_p = 1/v_p is large, and its contribution P_p (I_p-B)/v_p < 0 is large in magnitude. Summed over the many (mostly empty) shoebox pixels, this drags J below zero — the classic inverse-observed-count (Poisson-on-data) bias. It bites hardest where the true signal is weakest, i.e. at high resolution. In the fit it manifests differently but identically in origin: the weighted empty pixels make "no signal" (G=0) the cheapest solution, so G\to 0.

Fix. For background-limited (weak) reflections the correct variance is the local background, constant over the shoebox:

v_p = \max(B,\,1)\quad\Longrightarrow\quad
J = \frac{\sum_p P_p\,(I_p-B)}{\sum_p P_p^{2}},

the unbiased uniform-variance estimator. This single change turned \langle I/\sigma\rangle positive at all resolutions and stopped the scale collapse. It is applied to both the extraction weight and the fit weight.

---

2. Prediction band and multiplicity

Symptom. PixelRefine recorded ~4× fewer observations per unique reflection than the classical integrator — completeness was fine, redundancy was not.

Cause. A reflection is given a shoebox only when it lies within a radial band of the Ewald sphere,

\bigl|\,|S_{hkl}| - 1/\lambda\,\bigr| \le \delta .

For randomly oriented stills the number of images on which a given hkl satisfies this is \propto \delta. The default \delta = 5\times10^{-4}\,\text{Å}^{-1} was 4–6× tighter than the classical integrator's \delta = 2\text{–}3\times r_\mathrm{profile}, so each reflection was recorded on 4–6× fewer images.

Fix. Widen to \delta = 2\times10^{-3}\,\text{Å}^{-1}. Multiplicity rose from ~240 k to ~950 k observations and CC$_\mathrm{ref}$ from 49.7 % to 55.9 %. Widening is only safe once the fit is well-behaved (Sections 1, 3, 4); with the original unconstrained fit it caused divergence.

---

3. Regularising the per-image fit

Symptom. Freeing any per-image parameter (orientation, R, even the scalar scale G) collapsed the merged data (CC$_{1/2}$ from 90 % to a few %). The predict↔refine loop with all parameters frozen was, by contrast, byte-identical to extraction-only — proving the fit, not the loop, was at fault.

Cause. The per-image problem regularised nothing: orientation had no prior, R and G only a lower bound. Three orientation DOF (plus R, G) against a handful of signal pixels per still overfit the noise, and an unconstrained 1/G then scrambled the cross-image merge.

Fix. Anchor each refined parameter to its prior with a data-scaled weight, as ScaleOnTheFly already does for rotation data. The data term has one residual per pixel, so the prior weight must scale with the pixel count:

w_\theta = \sqrt{\frac{N_\mathrm{pix}}{\sigma_\theta^{2}}}\,,
\qquad \chi^2_\mathrm{reg} = w_\theta^{2}\,(\theta-\theta_0)^2 .

Using \sqrt{N_\mathrm{refl}} instead (as in the rotation scaler) is a factor \sqrt{N_\mathrm{pix}/N_\mathrm{refl}}\approx\sqrt{49}\approx 7 too weak and is simply not felt. Applied to:

• Scale: \theta=G,\ \theta_0=1. Prevents 1/G from wandering; restored CC$_{1/2}$ from 4 % back to 87 %.
• Orientation: \theta = Rodrigues vector, \theta_0 = spot-centroid orientation, \sigma_\theta in radians. At \sigma_\theta\!\sim\!1^\circ this gives the best CC$_\mathrm{ref}$; beyond \sim 2^\circ the fit overfits and the merge collapses, so \sigma_\theta is the safety knob.

---

4. Signal-weighting the fit

Cause. Even de-biased, a shoebox is ~80 % empty pixels (≈ 40 of 49 for a radius-3 box). They carry no information on G, R or orientation but add noise and, near the overfitting edge, destabilise the fit.

Fix. Multiply the fit weight by a detector-space Gaussian centred on the predicted spot,

w_p \;\to\; w_p \cdot \exp\!\left(-\frac{r_p^{2}}{2\sigma_s^{2}}\right),
\qquad r_p = \lVert (x_p,y_p) - (x_\mathrm{pred},y_\mathrm{pred})\rVert,

so the signal-bearing core drives the refined parameters (\sigma_s\approx 1.5 px). With the regularised orientation this lifted CC$_\mathrm{ref}$ from 60.5 % to 62.6 %.

---

5. Global orientation + cell-scale sweep

Motivation. The classical refinement (XtalOptimizer) works on spot centroids from spot-finding, which are poor value at high resolution; the local LSQ above is a heavy gradient step that cannot make the ~degree-scale global moves needed to pull a mis-indexed crystal's high-resolution reflections onto their shoeboxes. A small, global sweep that simply asks "at which orientation does the most high-resolution signal appear where the reference says it should?" is structurally better suited.

Score. Pearson CC of the box-summed intensities against the reference, over all matched reflections (not just the strong ones):

\mathrm{CC}_\mathrm{ref}(L) = \operatorname{corr}_{hkl}\bigl(\,I^\mathrm{box}_{hkl}(L),\ I^\mathrm{ref}_{hkl}\,\bigr).

The strong low-resolution reflections anchor the CC (moving them off their boxes collapses it), so the change in CC across the sweep is driven almost entirely by weak high-resolution reflections falling onto — or off — real signal. This is the "appearing out of the void" behaviour.

Geometry-derived bounds (parameter-free). A spot at resolution d sits at detector radius

r(d) = \frac{L}{p}\,\frac{\lambda}{d}\quad[\text{px}],

(L = detector distance, p = pixel size). Both a crystal rotation \delta\theta and a fractional cell-scale \epsilon displace that spot by an amount proportional to its radius:

\Delta_\mathrm{px} = r(d)\,\delta\theta \quad(\text{rotation}),\qquad
\Delta_\mathrm{px} = r(d)\,\epsilon \quad(\text{cell scale}).

So high-resolution spots (large r) move most. The two natural constraints fix the grid:

• Step = 1 px at the highest resolution: \;\delta\theta_\mathrm{step} = 1/r_\mathrm{max} (finer is below the detector's resolving power).
• Range = the orientation uncertainty \Delta\theta_u (a few px at high res). The number of steps per axis is then

n = \frac{\Delta\theta_u}{\delta\theta_\mathrm{step}} = \Delta\theta_u\, r_\mathrm{max},

a handful of steps, and the lowest-resolution spots move only n\,\delta\theta_\mathrm{step}\,r_\mathrm{min} = \Delta\theta_u\, r_\mathrm{min} \ll 2 px — i.e. the strong anchors stay put by construction.

Lesson learned. Setting the range from "2 px at low resolution" instead gives n = 2\,r_\mathrm{max}/r_\mathrm{min} = 2\,d_\mathrm{low}/d_\mathrm{high} steps — tens of pixels of high-resolution freedom — which lets the per-image CC overfit and degrades the merge. The range must be tied to the (small) orientation uncertainty, not the low-resolution cap.

The sweep is a coordinate descent over the three Rodrigues axes and the cell scale, run before the LSQ. It improves the high-resolution shells (CC$\mathrm{ref}$ +1 to +5 per shell) while preserving CC${1/2}\approx 90$ %. It is an alternative to the loose-\sigma orientation LSQ, not a complement — stacking the two double-moves the orientation and overfits. It is therefore available but off by default.

---

6. Background-estimator bias (the largest σ error; both integrators)

Symptom. Pushed past the true resolution limit (e.g. a 1.8 Å crystal merged to 1.3 Å), the no-signal shells still reported \langle I/\sigma\rangle \approx 4\text{–}6 with \mathrm{CC}_{1/2}\approx 0 — i.e. confident "data" where there is none. Present in both PixelRefine and the classical BraggIntegrate2D.

Cause. Both estimated the local background as the median of the surrounding ring. For a Poisson / right-skewed background the median sits below the mean,

\operatorname{median}(B) < \mathbb{E}[B],

so subtracting it under-subtracts on every pixel. The leftover positive offset is tiny per pixel but coherent, so over an $n_\mathrm{pix}$-pixel peak and a multiplicity-m merge it grows to a fake signal

\langle I\rangle_\mathrm{bias} \;\approx\; n_\mathrm{pix}\,\bigl(\mathbb{E}[B]-\operatorname{median}(B)\bigr),
\qquad
\Bigl(\tfrac{I}{\sigma}\Bigr)_\mathrm{merged,\,bias} \propto \sqrt{m}.

It is worst where the real signal is weakest (high resolution), because there the offset is all that remains — which is exactly the observed signature. Nothing leaks into the high-resolution shells; the background is simply under-estimated.

Fix. Use the mean of the ring (outliers excluded by the existing spot-core mask and saturation sentinels). \langle I/\sigma\rangle then collapses to ~0 wherever \mathrm{CC}\approx0, tracks CC down the shells, and the honest resolution limit becomes visible. This was the single largest contributor to untrustworthy σ — a one-line change in each integrator, not a variance-model problem.

---

7. Error model (global a, b; XDS form)

Counting statistics under-estimate the variance of strong reflections, which carry systematic errors (scaling, partiality, detector) proportional to intensity, not to \sqrt{I}. After merging, this leaves CC$_{1/2}$ and \langle I/\sigma\rangle inconsistent. The standard correction (XDS / DIALS / AIMLESS) inflates the variance with a global two-parameter model:

\sigma'^{\,2} \;=\; a\,\sigma^{2} + \bigl(b\,\langle I\rangle\bigr)^{2},
\qquad
\mathrm{ISa} \;=\; \frac{1}{b}\;=\;\lim_{I\to\infty}\frac{I}{\sigma'} .

Two points matter for an unbiased fit:

• The I^2 term uses the reflection mean \langle I\rangle (constant over its observations), not the per-observation I_i. Using I_i gives a down-fluctuated point a small \sigma' and hence a large 1/\sigma'^2 weight, biasing the merged mean — which collapses CC. (This was the decisive bug in the first attempt.)
• a and b are fit from the spread of symmetry equivalents: for an observation in a group of n, \mathbb{E}[(I_i-\langle I\rangle)^2] = \sigma_i^2(1-h_i) with leverage h_i = w_i/\sum w. Binning by intensity and regressing the bin medians of (I_i-\langle I\rangle)^2/(1-h_i) on (\sigma^2, \langle I\rangle^2) — weighted by 1/\mathrm{dev}^4 so the fit is relative — gives (a, b^2); the relative weight stops the strong bins (which fix b) from swamping the weak bins (which fix a).

It is applied at the merge level (MergeOnTheFly), so both integrators benefit, and jfjoch_process prints the model and ISa. Earlier per-resolution-shell variants were dropped: the standard tools use a single global a, b, and the per-shell version was partly masking the background bias of Section 6.

---

Results (lysozyme jet, 1.8 Å, identical input)

Configuration	N_obs	Compl.	CC$_{1/2}$	CC$_\mathrm{ref}$
Classical integrator (box-sum)	421 k	100 %	81.4 %	60.9 %
PixelRefine — original	61 k	95 %	0.0 %	−0.4 %
PixelRefine — consolidated (this work)	951 k	100 %	80.0 %	62.6 %

The consolidated integrator beats the classical one on the accuracy metric (CC$\mathrm{ref}$) with > 2× the multiplicity, and turns a previously unusable result (CC${1/2}=0$) into a competitive one.

---

Default recipe

Sections 1–5 are PixelRefine-specific; Sections 6–7 act at the integration/merge level and apply to the classical route too.

Field / behaviour	Default	Section
fit/extraction variance	local background B	1
ewald_dist_cutoff	2\times10^{-3}\,\text{Å}^{-1}	2
scale_reg_sigma	2.0	3
orient_reg_sigma_deg	1.0	3
refine_R	false	3
fit_signal_sigma_pix	1.5	4
sweep_orientation	false (available)	5
local background estimator	mean of the ring	6
merge error model	global a,b (ISa printed)	7

orient_reg_sigma_deg (accuracy vs. precision) and ewald_dist_cutoff (multiplicity vs. cost) are the two knobs worth tuning per dataset.