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#

The standard find_aberrations() workflow uses 200 initial TPE trials, then joint exact-gradient refinement through fifth order with tilt and depth. This diagnostic notebook keeps the primary model for its two-rotation comparison and custom search ranges; it reuses that fit for reconstruction without a fifth-order search.

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(search="primary", tilt=False)

model-based, tilt held at 0

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

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 DirectPtycho, 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. DirectPtycho(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?

ptycho = DirectPtycho(data, backend=backend, scan_detector_rotation_deg=scan_detector_rotation_deg)
ptycho.calibration()

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 primary 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 primary fits took 10 to 25 s in earlier 512 × 512 tests; timing depends on the backend and data. The explicit primary-fit ranges below apply to both branches. Edit them to constrain the search.

Both branches use the primary model. The standalone standard API additionally refines through fifth order; it is listed here to distinguish that workflow from this diagnostic comparison:

Call

Method

Use it for

find_aberrations()

model-based, tilt-corrected

the default: every fit

find_aberrations(search="primary", tilt=False)

model-based, tilt held at 0

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

 
fit_options = {
    "search_ranges": {
        "C10_nm": (-30.0, 30.0),
        "C12_nm": (0.0, 20.0),
    },
}
branches = ptycho.compare_rotation_branches(**fit_options, compute_loss=False, verbose=False)
assert ptycho.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=ptycho.scan_detector_rotation_deg,
              scan_detector_rotation_flipped=ptycho.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 primary, model-based fit used here 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 run() applies it without a third search and returns the phase at the native scan sampling, where the fit stays. To look at the phase, run(upsampling_factor="auto") reconstructs the same fit on the finer output grid the scan sampling allows: SSB output sampling, not interpolation.

result = ptycho.run(compute_loss=False)                          # the kept fit, on the native scan grid
shown = ptycho.run(upsampling_factor="auto", compute_loss=False)   # the same fit at the factor the scan sampling allows
result

The result itself displays the fitted parameters. report() returns a full-precision pandas table; report(to_dict=True) returns coefficients under aberrations and the sample geometry separately. Only fitted coefficients appear: this primary-model example has no fitted higher orders. result.aberrations remains the coefficient dictionary. The report includes timing and optimizer diagnostics when recorded. For a fitted result with history, use result.report(show_convergence=True).T to show the curve and a vertical table without refitting.

# Optional: uncomment when you want this report or save/restore step.
# table = result.report()
# parameters = result.report(to_dict=True)

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? ptycho.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.

Additional reconstruction scores are disabled (compute_loss=False). The fitting objective and existing scalar history are retained. Run All writes no calibration or result files. Click Save or uncomment viewer.save_calibration() to save the reviewed calibration.

# SSB is the default on CUDA and MPS. On CUDA, optionally add
# method="phase", phase_regularization=0.03 to compare regularized phase,
# or method="complex", amplitude_regularization=0.03 for amplitude + phase.
viewer = ptycho.show(show_fft=True, save_dir=save_dir, compute_loss=False)
viewer

# Optional: uncomment when you want this report or save/restore step.
# 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. The save and restore cells are commented out so Run All leaves the fitted widget open and writes no result files. Uncomment them after explicitly saving a calibration.

# Optional: uncomment when you want this report or save/restore step.
# viewer.close()
# ptycho.close()
# ptycho = DirectPtycho(data, backend=backend, scan_detector_rotation_deg=scan_detector_rotation_deg)
# # SSB is the default on CUDA and MPS. On CUDA, optionally add
# # method="phase", phase_regularization=0.03 to compare regularized phase,
# # or method="complex", amplitude_regularization=0.03 for amplitude + phase.
# viewer = ptycho.show(show_fft=True, save_dir=save_dir, compute_loss=False)
# viewer

# Optional: uncomment when you want this report or save/restore step.
# 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 ptycho.aberrations == saved["aberrations"]
# {"saved": {key: saved[key] for key in restored_keys},
#  "restored": {key: restored[key] for key in restored_keys}}

What to notice after opting in. The saved and restored columns agree to the last digit, and the session holds the restored model. ptycho.run(compute_loss=False) applies it without fitting again or calculating an additional diagnostic score.

Measured workflow#

Historical primary-model measurements on an Apple silicon Mac (MPS), each notebook run top to bottom with the data already on disk. These timings do not include the current standard API’s fifth-order refinement (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, DirectPtycho(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. run() reuses the fit of section 3 without another search. 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 DirectPtycho(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 primary model with tilt correction; the table at the top gives the explicit primary call for zero tilt. ShowPtycho explains when to trust the fitted tilt.

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