LARBED Extraction (pyextal.larbed)

pyextal.larbed is the stage that runs before refinement. It turns a raw EMPAD 4D-STEM acquisition into one 2D rocking curve per reflection — the arrays that LARBEDDiffractionInfo and the rest of the pipeline consume.

It was merged in from the standalone PyLARBED package, so the class and method names are unchanged; only the import path is new:

from pyextal.larbed import LARBEDAnalysis, LARBEDCalibration

---

Where it sits in the pipeline

scan_x256_y256.raw   (EMPAD, 4D: scan_x × scan_y × det_x × det_y)
      │
      ▼
┌──────────────────────────────────────────────┐
│  Stage 0 — Extraction  (pyextal.larbed)      │
│                                              │
│  LARBEDAnalysis                              │
│    → load_data / align_data                  │
│    → assign_gvector + calculate_grid_vectors │
│    → calculate_g_vectors                     │
│    → IntegrateLarbed                         │
│    → to_dinfo(...)  or  save_h5(...)         │
└───────────────┬──────────────────────────────┘
                │  LARBEDDiffractionInfo  /  {sample}.h5  ──▶  /raw
                ▼
        Stage 1 — Pre-refinement   (LARBEDROI, CoarseOptimize, FineOptimize)
                ▼
        Stage 2 — JAX refinement   (jaxEngine)

---

Typical workflow

from pyextal.larbed import LARBEDAnalysis

larbed = LARBEDAnalysis("scan_x256_y256.raw")
larbed.load_data(type=0)
larbed.average_data()
# larbed.align_data()          # optional; only if the direct beam drifts

# 1. Index the reflections. Two basis vectors for a zone-axis pattern …
larbed.assign_gvector((46, 48), (64, 36))          # detector pixel positions
larbed.calculate_grid_vectors(crop=2, order=3, Adjust=True)
larbed.calculate_g_vectors((1, 3, -2), (4, 2, 0))  # matching Miller indices
larbed.plot_grid_vectors(ratio=0.005)              # check the indexing

# 2. Integrate each reflection into a rocking curve
larbed.IntegrateLarbed(ri=3, ro=(6, 8), g=580, m=0.47, Variance=True)
larbed.plot_larbed(ratio=0.5)

assign_gvector takes detector pixel positions; calculate_g_vectors takes the Miller indices of those same reflections, in the same order. Supply two of each for a zone-axis grid, or four to include HOLZ reflections:

larbed.assign_gvector((59.7, 52.3), (52.1, 63.2), (53, 32.2), (62, 62))
larbed.calculate_grid_vectors(order=(8, 8, 0, 0), Adjust=False)
larbed.calculate_g_vectors((-1, 0, 0), (0, 1, 0), (0, 0, 1), (0, 0, 0))

calculate_grid_vectors dispatches on how many basis vectors it was given: four selects the HOLZ grid (and order must then be a 4-tuple), two selects either the peak-snapped grid (Adjust=True) or the plain lattice grid.

Aperture and variance

IntegrateLarbed(ri, ro, g, m, DQE, varB, Variance) sums each reflection over a square aperture of half-width ri. When ro=(inner, outer) is given, the mean of that annulus is subtracted as background.

With Variance=True the counting variance is propagated alongside. The per-pixel model matches pyextal.LucyRichardson.sgnoise(): a signal-dependent term g*m*I/DQE plus a constant readout term g*m*varB, summed over the aperture. varB defaults to 0.

---

Handing off to refinement

Two routes, both replacing the older .npy round-trip:

# (a) straight into the refinement objects
dinfo = larbed.to_dinfo("si110.dat", thickness=1000.0, gl=37.0)
roi = LARBEDROI(dinfo=dinfo, rotation=214, gx=..., gInclude=[...])

# (b) into the same HDF5 file the later stages append to
larbed.save_h5("sample.h5", prefix="YIG_tilt12", crystal_dat="yig.dat")

save_h5 writes /raw/dp, /raw/variance, and /raw/g_vectors — exactly what pyextal.hdf5io.save_raw() defines — so extraction, pre-refinement, and every JAX refinement run share one file per dataset.

gl is the tilt step size, expressed as the number of scan steps spanning the x-axis reflection. Deriving it is covered in Understand how to define/index LARBED dataset; calculate_average_step_length() measures the pixels-per-step half of it.

stack() returns the same three arrays directly if you would rather handle them yourself, and save_larbed() still writes the loose .npy files.

---

MTF calibration and deconvolution

Note

Legacy. Because IntegrateLarbed already integrates over an aperture, deconvolving against the MTF mainly adds noise, and deconv_Larbedstack discards the propagated variance. It is normally skipped.

LARBEDCalibration measures the detector point-spread function from a reference scan taken with no sample and descan off:

calib = LARBEDCalibration("probe_scan.raw")
calib.load_data(type=0)
calib.find_peaks(start=100, end=115)
calib.calculate_average_step_length()
calib.crop_peaks()
calib.interpolate()
calib.calculate_mtf()

larbed.deconv_Larbedstack(mtf_2d=calib.mtf_2d, varB=787, g=580, niter=10)

Two convolution conventions

calculate_mtf returns a Fourier-domain MTF, and pyextal.LucyRichardson.lucy_Richardson() accepts a mode argument selecting how it is applied:

mode	Operator	Use
'fourier' (default)	ifft2(fft2(img) * mtf)	Consistent with a Fourier-domain MTF, and what the CBED path has always used.
'real'	scipy.signal.convolve(img, mtf, mode='same')	What the standalone PyLARBED used. Kept so earlier results reproduce exactly.

The LARBEDAnalysis deconvolution methods pass mode='real', so they reproduce pre-merge PyLARBED output bit for bit. The two operators are genuinely different — a Fourier-domain array fed to the real-space convolution is not the same deconvolution — so switching an existing analysis to 'fourier' will change its numbers.

---

Worked example

example/LARBED/STO_001/STO_001_LARBED.ipynb runs this whole page against a real SrTiO₃ [001] dataset and then continues through pre-refinement, using the indexing parameters validated in the original PyLARBED notebook. It ships with its crystal file, STO_001.dat.

Note the step it takes between extraction and refinement. STO_001.dat leaves both of its beam lists implicit: sele picks the matrix beams geometrically at parse time and aper picks the output reflections by g-length, so both move with the voltage, the tilt, and xaxis. The notebook rewrites it as STO_001_explicit.dat with an explicit hkl block (all 185 beams) and an explicit out block set to exactly gInclude, and refines against that — so the calculation is reproducible beam for beam.

The output beams are written as a prefix of the hkl list, so hkl[:nout] == hklout and the engine’s idptr mapping is the identity. That is worth preserving: bloch_parse resolves each out entry to its index in hkl, and any entry it cannot find is appended to hkl — past the array already allocated from the declared nbeams.

Warning

That notebook’s gl (tilt step size) is a placeholder and its Debye-Waller factors are literature starting values. Set them for your data before trusting a refined number.

---

Module layout

Module	Contents
pyextal.larbed.analysis	LARBEDAnalysis — the extraction workflow
pyextal.larbed.calibration	LARBEDCalibration — MTF and scan-step calibration
pyextal.larbed.empad	Raw I/O, virtual apertures, reciprocal-grid construction
pyextal.larbed.peaks	Least-squares peak models and fitting

Throughout pyextal.larbed.empad, x is the first (slow) detector axis and y the second — not the screen-space x/y of a plotted image. The grid builders take a shape argument used to reject off-detector nodes; it defaults to (124, 124), the cropped EMPAD frame.