From a 4D-STEM scan to an SSB phase

From a 4D-STEM scan to an SSB phase#

Local workflow template. Run this notebook on an Apple silicon Mac (MPS) or an NVIDIA GPU (CUDA), at the microscope or afterwards. SSB needs a GPU, so the documentation site shows the source without running it.

One notebook, one acquisition, seven questions:

  1. What did the detector record, and is it all still on the GPU? (load)

  2. How many mrad is one detector pixel? (detector calibration)

  3. How is the scan turned relative to the detector, including the 180° choice? (scan-detector rotation)

  4. What does the specimen look like in bright field, dark field and centre of mass? (Show4DSTEM)

  5. Which defocus and astigmatism did the probe have, and how is the crystal tilted? (SSB)

  6. Does the phase hold up when you change them by hand? (ShowPtycho)

  7. Can you pick up tomorrow where you stopped today? (save and restore)

The default input is the public gold nanoparticle scan: 512 × 512 probe positions, each recorded on the native 192 × 192 detector. Nothing in this notebook bins, crops or filters the data, and every value is shown as a table or a figure. The assert lines check the data path as the notebook runs.

SSB in quantem: model-based by default#

quantem fits SSB aberrations with a physical model of the measurement, not an image-quality score. For each trial set of aberrations, and the sample’s tilt and thickness, it predicts the double-overlap data the detector should have recorded, then keeps the parameters whose prediction agrees best with the measurement (least squares). This model-based, tilt-corrected fit is the default everywhere, and it is where SSB in quantem is heading.

Call

Method

Use it for

find_aberrations()

model-based, tilt-corrected

the default: every fit

find_aberrations(tilt=False)

model-based, tilt held at 0

thin, untilted samples, or to see what tilt correction changes

find_aberrations(tilt=False, objective="phase_variance")

phase variance, the earlier SSB search

a backup: reproducing older results, cross-checking a fit

Phase variance chooses the aberrations that make the corrected phases of all bright-field pixels agree best: the least variance of the corrected phase across them. It stays as a backup, not an equal alternative: it cannot see sample tilt or thickness, it scores only how consistent the bright-field images are rather than how well a model explains the measurement, and it is about 10x slower. On simulated crystals with known aberrations, the model-based fit lands about 3x closer in C10.

Setup and parameters#

Install once with the line for your GPU, then restart the kernel. Installation lists the options.

# Apple silicon:
# %pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ "quantem.widget[mps]"
# NVIDIA GPU:
# %pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ "quantem.widget[cuda]"

What do you need to know about your acquisition? Three numbers that the diffraction patterns cannot tell you:

  • Beam energy (voltage_kV): it sets the electron wavelength.

  • Convergence semiangle (semiangle_mrad): the edge of the bright-field disk. Section 2 turns it into the detector step.

  • Scan step (scan_sampling_A, in Å): the distance between probe positions.

Two more are optional. detector_sampling_mrad states a calibrated camera length, which then replaces the fit of section 2. scan_detector_rotation_deg states your microscope’s scan rotation; leave it None and section 3 measures it from the data. backend picks the GPU: "auto", "mps" or "cuda".

Point master at your own Arina master file (with its data shards beside it) or a .qem file. Leave it None to use the gold scan.

from pathlib import Path

master = None                       # your *_master.h5 or .qem; None downloads the public gold scan
backend = "auto"                    # "mps" (Apple silicon), "cuda" (NVIDIA) or "auto"
voltage_kV = 300
semiangle_mrad = 30
scan_sampling_A = 0.5               # Å, one number for both scan axes
detector_sampling_mrad = None       # None: fitted from the bright-field disk in section 2
scan_detector_rotation_deg = None   # None: estimated from the centre-of-mass field in section 3

The gold scan is 4.6 GB, downloaded once from the quantem-data dataset at a pinned revision and reused from the local cache afterwards. qw.profile() shows the installed packages and the GPU in use: paste it into any bug report.

import json

import pandas as pd
from huggingface_hub import snapshot_download
from IPython.display import display

import quantem.widget as qw
from quantem.gpu import SSB, dpc, io
from quantem.widget import Show2D, Show4DSTEM

qw.profile()
if master is None:
    folder = snapshot_download(
        "bobleesj/quantem-data", repo_type="dataset",
        revision="5925e4b49b2ca28ff82d6d5a6b9a8bfea9edfb20",
        allow_patterns="4dstem/gold_512/*",
    )
    master = Path(folder) / "4dstem/gold_512/gold_004_master.h5"
master = Path(master)
save_dir = Path("ssb-results") / master.stem   # ShowPtycho's saved calibration for this acquisition

1. Load#

What did the detector record, and is it all still on the GPU? io.load reads the master and its shards straight onto the GPU and keeps every count of every detector pixel there, in lossless ANS-encoded storage: a few GB instead of the tens of GB a dense array would take, with no host copy and no binning. hot_pixel_correction="median" replaces the bad pixels the detector flags with the median of their neighbours, and the load line says how many.

data = io.load(master, backend=backend, hot_pixel_correction="median")
data.calibration.voltage_kV = voltage_kV
data.calibration.semiangle_mrad = semiangle_mrad
data.calibration.scan_sampling_A = scan_sampling_A
if detector_sampling_mrad is not None:
    data.calibration.detector_sampling_mrad = detector_sampling_mrad
assert data.representation.value == "encoded"   # every count stays encoded on the GPU; nothing is decoded whole
data.calibration

What to notice. The table lists each calibration with its source: recorded in the file, or set here. The detector step is still missing; section 2 measures it.

2. Detector calibration#

How many mrad is one detector pixel? The bright-field disk is the image of the probe-forming aperture, so its edge sits at the convergence semiangle: the semiangle divided by the disk’s edge radius in pixels is the detector step in mrad per pixel. fit_detector_sampling() averages every scan position into one mean pattern on the encoded data, fits the disk’s centre and the radius where it falls to half its plateau, and records all three on data.calibration.

fit = data.calibration.fit_detector_sampling()
assert data.calibration.detector_sampling_mrad is not None             # recorded: SSB and Show4DSTEM read it
assert data.calibration.bright_field_center_px == fit.bright_field_center_px
fit.show()

data.calibration

What to notice. The circle sits on the disk’s half-intensity edge, at the semiangle; the axes are in mrad from the disk centre. The table now gives the detector step, the disk centre and its edge radius with the source “Fitted (bright-field disk)”. Every later step reads these values instead of fitting the disk again. A detector step you gave in the parameters cell is kept, and one line says so.

3. Scan-detector rotation#

How is the scan turned relative to the detector? The centre of mass (CoM) of each diffraction pattern shifts when the specimen’s electric field deflects the probe. Rotated into the scan frame by the right angle, that shift field is a gradient field and has no curl; by a wrong angle, swirls appear. SSB(data) estimates the angle that removes the curl and says so in one line, or uses scan_detector_rotation_deg from the parameters cell. The table shows every value the session uses and where it came from.

Predict: rotate every CoM vector by a further 180°, so each one points the opposite way. Does the curl come back?

ssb = SSB(data, backend=backend, scan_detector_rotation_deg=scan_detector_rotation_deg)
ssb.calibration_table()

It does not: reversing every vector leaves a curl-free field curl-free, so the CoM fixes the angle only modulo 180°. The two branches reconstruct the same image with opposite sign and opposite defocus.

Which 180° branch is real? compare_rotation_branches() runs the tilt-corrected fit once on each branch (defocus C10, astigmatism C12 and phi12, with the specimen tilt and a depth spread). Atoms push the electron phase forward, so on the right branch the atomic columns, located on the dark-field image of the same scan, are bright. When the columns do not decide, the branch whose phase has the longer positive tail is kept. The two fits take 10 to 25 s on a 512 × 512 scan; uncomment one fit_options line to constrain both.

Both branches run the default fit, the first row of the table:

Call

Method

Use it for

find_aberrations()

model-based, tilt-corrected

the default: every fit

find_aberrations(tilt=False)

model-based, tilt held at 0

thin, untilted samples, or to see what tilt correction changes

find_aberrations(tilt=False, objective="phase_variance")

phase variance, the earlier SSB search

a backup: reproducing older results, cross-checking a fit

Phase variance chooses the aberrations that make the corrected phases of all bright-field pixels agree best: the least variance of the corrected phase across them. It stays as a backup, not an equal alternative: it cannot see sample tilt or thickness, it scores only how consistent the bright-field images are rather than how well a model explains the measurement, and it is about 10x slower. On simulated crystals with known aberrations, the model-based fit lands about 3x closer in C10.

fit_options = {}
# Limit the defocus search to -50 to 50 nm:
# fit_options = {"search_ranges": {"C10_nm": (-50, 50)}}
# The same, with the astigmatism held at zero:
# fit_options = {"search_ranges": {"C10_nm": (-50, 50)}, "fixed_aberrations": {"C12_nm": 0.0}}
branches = ssb.compare_rotation_branches(**fit_options, verbose=False)
assert ssb.scan_detector_rotation_deg == branches.kept_rotation_deg   # settled: every later step uses this angle
branches.report()

branches.show()

panels = branches.upsampled_phases()   # each branch's own fit at upsampling_factor="auto"
factor = branches.upsampling_factor
Show2D(list(panels),
       labels=[f"{branch.scan_detector_rotation_deg:.1f}°, {factor}x auto" + (" (kept)" if index == branches.kept else "")
               for index, branch in enumerate(branches.branches)],
       sampling=scan_sampling_A / factor, units="Å", ncols=2)

What to notice. In the report, column z is the mean phase at the dark-field columns in standard errors: decisive at ±3, with column effect past ±0.25. The figure’s top row and the viewer show each branch’s own fit reconstructed at upsampling_factor="auto", the smallest factor whose output grid holds the whole SSB bandwidth: 4× for the gold scan’s 0.5 Å step, so 0.125 Å pixels. On the kept branch the atomic columns, circled in the figure, are bright. The tests read the native-grid phases, as do the histograms, where the kept branch is solid. A rotation you gave is kept as given; compare_rotation_branches(apply=True) moves to the selected branch. Bright columns are an assumption: in thick crystals the contrast can invert.

4. Explore with Show4DSTEM#

What does the specimen look like at each scattering angle? Show4DSTEM is the interactive explorer: click a scan position on the right to see its diffraction pattern on the left, and drag the detector on the left to build a virtual image on the right. Its BF, ABF, ADF and HAADF presets are multiples of the edge radius fitted in section 2, around the same centre. Predict: in bright field, are the gold particles brighter or darker than the support film? In dark field?

explorer = Show4DSTEM(data)
explorer

virtual, detectors = {}, {}
for preset in ("bf", "abf", "adf"):
    explorer.apply_preset(preset)
    virtual[preset.upper()] = explorer.virtual_image()
    detectors[preset.upper()] = {"inner radius (px)": explorer.roi_radius_inner, "outer radius (px)": explorer.roi_radius}
explorer.apply_preset("bf")
assert detectors["BF"]["outer radius (px)"] == data.calibration.bright_field_edge_radius_px
display(pd.DataFrame(detectors).T.round(2))
Show2D(list(virtual.values()), labels=list(virtual), link_contrast=False,
       sampling=scan_sampling_A, units="Å", ncols=3)

What to notice. BF sums the inside of the disk, so the particles are dark on a bright background; ABF sums its outer half; ADF sums a ring from the disk edge to twice its radius, so only scattered electrons reach it and the particles are bright. The BF outer radius is the fitted edge radius.

Which way does the specimen push the beam? dpc.run computes the CoM shift and its integral, the iDPC image, turned into the scan frame with the rotation settled in section 3.

com = dpc.run(data, scan_detector_rotation_deg=ssb.scan_detector_rotation_deg,
              scan_detector_rotation_flipped=ssb.scan_detector_rotation_flipped)
Show2D([com.com_row_aligned, com.com_col_aligned, com.phase], labels=["CoM shift, row", "CoM shift, col", "iDPC"],
       link_contrast=False, cmap=["gray", "gray", "inferno"], sampling=scan_sampling_A, units="Å", ncols=3)

What to notice. A thin, weak specimen deflects the beam only where its potential changes, so the CoM outlines the particles; a thick crystal slightly off its zone axis shifts the CoM across its whole area. iDPC turns small offsets into slow, large-scale variations, so read it qualitatively: unlike the SSB phase below, its values are not radians.

5. SSB: probe aberrations and specimen tilt#

Which defocus and astigmatism did the probe have, and how is the crystal tilted? Every bright-field detector pixel sees the specimen through a slightly different probe tilt. With the right aberrations the phase images from all those pixels agree. The default, model-based fit also models a crystal that leans by a sample tilt through a depth spread, the pivot myptycho uses for its multislice reconstruction. Section 3 already kept one branch’s fit, so reconstruct() applies it without a third search and returns the phase at the native scan sampling, where the fit stays. To look at the phase, reconstruct(upsampling_factor="auto") reconstructs the same fit on the finer output grid the scan sampling allows: SSB output sampling, not interpolation.

result = ssb.reconstruct()                          # the kept fit, on the native scan grid
shown = ssb.reconstruct(upsampling_factor="auto")   # the same fit at the factor the scan sampling allows
display(result.report().T)
Show2D(shown.phase, sampling=shown.object_sampling_A[1], units="Å",
       title=f"SSB phase, tilt-corrected, {shown.upsampling_factor}x auto")

What to notice. The report lists C10 (defocus), C12 and phi12 (twofold astigmatism and its direction), the sample tilt along the scan rows and columns, the depth spread, and the rotation SSB kept with the test that chose it. A fitted depth spread is a model parameter, not a measured thickness. The phase is that fit at 4× auto on the gold scan: the title names the factor, the scale bar uses the scan step divided by it, and the particles are bright.

6. Review in ShowPtycho#

Does the phase hold up when you change the aberrations by hand? ssb.show() opens ShowPtycho on the fitted session. Drag C10, C12, phi12 or the rotation and the phase and its FFT update live. The Sample panel holds the fitted tilt and depth spread, and the launch table names each value’s source. Save writes calibration.json to save_dir; save_calibration() does the same from Python, so a top-to-bottom run also leaves a saved calibration.

viewer = ssb.show(show_fft=True, save_dir=save_dir)
viewer

viewer.save_calibration()
calibration = viewer.calibration()
{key: calibration[key] for key in ("aberrations", "scan_detector_rotation_deg", "scan_detector_rotation_flipped",
                                   "sample_tilt_mrad", "propagator_tilt_mrad", "depth_spread_nm",
                                   "detector_sampling_mrad", "scan_sampling_A", "upsampling_factor", "loss")}

What to notice. In the FFT the lattice spots should be sharp in every direction: a defocus error blurs them all, astigmatism one direction more than the other. calibration() returns the model on screen, slider changes included, without refitting or saving; its detector step is the one fitted in section 2.

7. Save and restore#

Can you pick up tomorrow where you stopped today? Close the session, open a new one on the same acquisition, and show it with the same save_dir. ShowPtycho finds calibration.json, checks that it belongs to this acquisition and launch calibration, and restores it before the first preview, with one line saying what it restored. No aberration search runs, so this takes seconds.

viewer.close()
ssb.close()
ssb = SSB(data, backend=backend, scan_detector_rotation_deg=scan_detector_rotation_deg)
viewer = ssb.show(show_fft=True, save_dir=save_dir)
viewer

saved = json.loads((save_dir / "calibration.json").read_text())
restored = viewer.calibration()
restored_keys = ("aberrations", "scan_detector_rotation_deg", "sample_tilt_mrad", "depth_spread_nm")
assert all(restored[key] == saved[key] for key in restored_keys)   # restored exactly, without a new fit
assert ssb.aberrations == saved["aberrations"]
{"saved": {key: saved[key] for key in restored_keys},
 "restored": {key: restored[key] for key in restored_keys}}

What to notice. The saved and restored columns agree to the last digit, and the session holds the restored model, so ssb.reconstruct() applies it without fitting. Restarting the kernel and running the notebook again ends the same way.

Measured workflow#

Measured on an Apple silicon Mac (MPS), each notebook run top to bottom with the data already on disk (the gold download, 4.6 GB, comes once before this):

Section

Gold scan

A 512 × 512 scan with a 192 × 192 detector

1. Load

4.6 s

1.8 s

2. Detector calibration

1.8 s

0.2 s

3. Rotation estimate, SSB(data)

6.3 s

0.6 s

3. Two-branch fits

22.1 s

8.4 s

4. Show4DSTEM, presets and DPC

11.4 s

0.8 s

6. ShowPtycho and save

1.0 s

0.8 s

7. Reopen and restore

8.2 s

2.1 s

Whole notebook

58.5 s

17.8 s

The detector calibration is the mean pattern of every scan position, computed on the encoded data, and a fit on that one pattern. reconstruct() takes no time here: it applies the fit of section 3. The table predates the auto-factor displays: on the gold scan the two 4× branch panels (2048 × 2048 each) add 0.4 s and the 4× phase of section 5 about 0.2 s. Reopening repeats the rotation estimate of SSB(data) before the saved model replaces it.

Next steps#

  • Share a viewer. viewer.export("ssb-results/showptycho") writes a folder that opens in Chrome without Python: double-click ShowPtycho.command. See Export and run ShowPtycho.

  • Iterative ptychography. viewer.calibration() is the starting model for multislice ptychography; the SSB to iterative regions notebook hands each fitted region to its own run.

  • Thick, tilted crystals. The fit above is the model-based, tilt-corrected default; the table at the top lists tilt=False and the phase-variance backup. ShowPtycho explains when to trust the fitted tilt.

  • Close the data when you are done, to free the GPU: ssb.close(), then data.close().