Native SSB: the model-based, tilt-corrected fit

Native SSB: the model-based, tilt-corrected fit#

MetalSSBKernels runs Python’s SSB.fit_aberrations on a prepared native session: the same transfer model, objective, search and branch test, so a native app gives the Python MPS result. The app owns selection and presentation; the shared Swift package owns count extraction, Fourier preparation, the model, the fit and saved numerical state. Python is not needed at runtime.

SSBFitOptions

Python

Method

SSBFitOptions()

fit_aberrations()

model agreement, tilt-corrected (sample tilt and depth spread)

SSBFitOptions(tilt: false)

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

model agreement, tilt and depth held at 0

SSBFitOptions(tilt: false, objective: .phaseVariance)

objective="phase_variance" (retired in Python too)

retired (2026-10-25): saved options decode, validate() refuses a fit

Minimal workflow#

let model = try engine.jointModel(scanDetectorRotationDegrees: rotation)  // the curl estimate
let fit = try SSBModelFit.findAberrations(
  model: model, rotationDegrees: rotation, darkField: darkField,
  scanSampling: (stepRow, stepColumn))                    // both 180-degree branches
let phase = try model.reconstruct(
  fit.fit.point, upsamplingFactor: calibration.ssbUpsamplingFactor())  // phase_of_mean, auto factor

SSBRotationBranches.curlRotationDegrees gives the rotation from exact centre-of-mass moments (dpc.find_optimal_rotation); findAberrations fits it and its opposite branch with the same options, keeps the branch whose dark-field atomic columns have positive phase (else the larger phase skewness), and measures the standard comparison (tiltFitGain) once, for the kept branch. checkRotation: false fits the given rotation alone. Prepare engine with the calibrated geometry and the original unsigned counts, after the loader’s flagged-pixel median fill when the Python session applied one. Preparing the engine again invalidates joint models; retarget moves a model to another rotation without gathering its band evidence again.

Scientific contract#

Quantity

Contract

C10, C12

Angstrom in the engine (SSBModelPoint has nm accessors); f = -C10 at the same reference plane

C10 reference

Model midplane; not a microscope entrance-plane focus measurement

Astigmatism angle

Radians, wrapped to [-pi/2, pi/2)

Tilt

sampleTiltMrad (row, column) in mrad, scan frame

Depth spread

Nonnegative Angstrom; a model depth spread, not a measured thickness

Objective

agreements: per band frequency `

Search

SSBAgreementSearch: the resident multivariate TPE (MT19937, 200 trials, seed 42) of direct_ptycho._resident_search, then a trust-region Newton on the exact gradient and Hessian of the agreement (ssb-low-order-derivatives), in double precision; a joint fit is polished once more from the phi12-mirrored best trial and the best trial in another (phi12, tilt, depth) cluster, kept only when its agreement is higher by more than 1e-6 relative. Only the given rotation is searched: the opposite 180-degree branch gets the exact image of the fit (even-degree chi terms and tilt negated, odd-degree terms and depth kept), which has the same agreement

Fifth order

SSBFifthOrderGradientSearch (higherOrderSearch: .exactGradientCoarseToFine, the default): exact gradients of the 28-coordinate agreement (ssb-fifth-gradient), a trust-region dense BFGS scaled by and started from a Hessian of every 64th BF pixel, the active set held at the box bounds (SSBModelFit.fifthOrderBounds(higherOrderLimit:): the low-order box, then each of the 22 higher-order coordinates within ±SSBFitOptions.higherOrderLimit rad of probe phase at the aperture edge, 20 by default, which no longer binds, or the standard 2), first on every 8th BF pixel and then converged on all of them (a subsample keeps at least 32 pixels, so a small bright field uses all of them); accept/reject and the stop compare Double-summed agreements. The finite-difference Newton (.finiteDifferenceNewton) is retired natively, as in Python: saved options naming it decode, and a fifth-order fit with it is refused. refineToFifthOrder refines a shown low-order fit without its search, again at another limit when asked; onRefinedFit receives the fit before its checks

Half-band check

Default for a ±20 fit (SSBFitOptions.checkHigherOrderConsistency, SSBFit.higherOrderConsistency, SSBHigherOrderConsistency), the same definition on Metal, CUDA and MPS: the fit’s band frequencies are split by ring parity floor(|q| / 0.05 1/A) mod 2 (q and -q in the same half); for each parity the fifth order is refitted at ±20 with the other half’s band weights 0, both from the same full-data low-order start as the ±20 fit; the statistic is the rms over the aperture (|k| <= alpha, 96 x 96 samples of the unit disk), piston removed, of chi_A - chi_B through fifth order (C10 and C12 included, tilt and depth not). Trusted when <= 0.5 rad: the ±20 fit is kept. Otherwise the ±2 fit is returned (options.higherOrderLimit 2) with its higher orders marked not trusted (thick crystals give spurious terms at ±2 too), and the record keeps the checked ±20 fit. Record names as Python: half_band_rms_rad, trusted, shown_limit. Measured data: 0.16-0.35 rad on four frames of a thin two-dimensional crystal, about 1.33 rad on a thick crystal. Cost: two ±20 fits, a third (±2) on fallback; timings higher_order_consistency_s, higher_order_fallback_s

Final phase

phase_of_mean (corrected waves summed, one phase; scan-mean replicas zero). At 1x-4x the interactive redraw kernel (ssb_redraw.metal) serves every term through C56, tilt and depth for cached and streamed spectra: the anti-Hermitian half is evaluated and mirrored, each source frequency visits its output aliases per batch of 32 BF pixels, chi is a per-row Horner polynomial, and SIMD groups skip BF pixels with closed sidebands. It agrees with the earlier per-term kernel within 1e-5 rad (observed about 5e-8 on the synthetic test scene, 4e-6 on the full-size synthetic benchmark). 8x keeps the per-term kernel

Reconstruction methods

SSBReconstructionMethod, ShowPtycho’s three: .ssb (the final phase above), .phase and .complex (regularized direct solves, below). The mean_phase and complex_wave estimators are retired (2026-10-25): their kernels and API are removed; the app opens results saved with them as saved and reconstructs edits with one of the three

Phase

.phase (Python method="phase"): the linear weak-phase model T = exp(i phi) ~ 1 + i phi solved per scan frequency for every supported output alias q (

Complex

.complex (Python method="complex"): the relative log-amplitude eta and the phase phi of T = exp(eta + i phi) solved together (columns t1 + t2 and i (t1 - t2)), amplitude damping on eta; amplitude = exp(eta) (relative; zero-mean eta). CUDA parity: phase as .phase, amplitude below 2e-7 of its mean

Regularization

SSBDirectRegularization: phaseDamping and amplitudeDamping (> 0, default 0.03), smoothnessInverseAngstroms q0 (>= 0, default 1; 0 = frequency-independent prior), detectorWeighting .uniform or .inverseMean (1 / max(mean count, 5 % of the median), normalized to mean 1 over every selected BF pixel, those outside the aperture’s support included, as Python). validate(for:) gives Python’s messages; .ssb ignores it. Results record method and regularization

Caching

The layout (frequencies, aliases, priors, weights, dose) is kept per data, factor, q0 and weighting; the normal equations per optics, rotation and layout, one per method (a complex system also serves .phase). A damping change solves again only (about 1–7 ms at 512 x 512, 1–3x); an optics change prepares again (about 0.17 s phase, 0.17-0.19 s complex at 3x on an Apple M5 Max; complex systems of two or more aliases use one thread per frequency and alias, direct_normal_tiled)

Drag previews

buildDragWindow(along:at:...) builds, one exact preparation per call, a Chebyshev series of the prepared normal equations along one SSBDragControl (C10, C12, phi12, a higher-order magnitude or azimuth, tilt, depth) with every other input held: 29 nodes (50 for an azimuth) over +-h, h from the transfer phase’s spread (C10 +- 6.8 nm, C12 +- 3.4 nm at 300 kV, 30 mrad); up to three tiles per control, a quarter of the device’s working set. dragPreview evaluates it and solves as reconstruct (5-10 ms at 3x); the result’s preview marks it, it agrees with the exact image to about 1e-6 of its std, and it is never saved. Port of Python’s CUDA DragWindow; design note 2026-10-25-ssb-direct-drag-windows.md

Output sampling

Fourier-alias reconstruction at 1, 2, 3, 4 or 8x; not interpolation. SSBOutputSampling.auto (ssbUpsamplingFactor) picks the smallest of 1-4x covering nominal support and caps at 4x, as Python; ssbUpsamplingRequirement gives the uncapped requirement and ssbUpsamplingChoice the factor with the uncapped factor and a capped flag. SSBOutputSampling.choices lists Auto and explicit 1-4x for a viewer; 8x stays an explicit factor

Rotation

Double precision; the angle is converted in double, the BF coordinates rotated in float32 as MPS

Wavelength

Python’s table at 80-400 kV, else the same formula

The TPE port is checked trial by trial: on a synthetic float32 objective that both sides evaluate from the same float32 trial rows, the 200 search trials equal Python’s (SSBAgreementSearchTests, which scripts/ssb_agreement_gate.py also feeds with the current Python records). The Newton after the search has exact derivatives only: on the same objective, given its analytic derivatives, the tilt and tilt=False fits reach its known maximum, and the staged Newton reaches the one-stage Newton’s point with fewer derivative passes. The exact-gradient fifth-order search is checked for a maximum of the Double-summed agreement on the synthetic test scene: a small projected exact gradient, and no rise along any of the 28 coordinates (SSBFifthOrderGradientTests). On three measured 512 x 512 acquisitions it reached a higher Double-summed agreement than the retired finite-difference Newton (+1.2e-7 to +3.8e-7 relative) in about 1 s instead of 65-111 s; the parameters differed along the C10-C30-C50 direction, which the agreement resolves only to about 1e-3 to 1e-2 nm at float32. Measured parity (MetalSSBMeasuredParityTests, given an exported acquisition in SSB_PARITY_FIXTURE) covers the BF disk and fitted detector sampling, the agreement at fixed parameters (1e-7 relative), the default and tilt=False fits on the model’s exact derivatives (tolerances 0.01 nm, 0.2 degrees, 0.02 mrad; the retired stencil Newton measured 1e-5 nm, 1e-4 degrees and 1e-4 mrad, not re-measured since), the curl rotation, the branch decision and the fixed-parameter phases (correlation 0.999999998 or better at 1x and 3x).

Limits#

  • The joint model fits and reconstructs C10 through C56 with tilt and depth spread; an unsupported order or a repeated term is rejected.

  • Supported scans are square 128, 256 and 512. There is no silent binning or padding for this model.

  • Every fit is model-based (the agreement objective), natively and in Python.

  • The native engine has no detector transposition (scan_detector_rotation_flipped is false).

  • No absolute potential, acquisition focus sign or specimen thickness claim follows from contrast or the agreement alone.

Saved results#

MetalSSBJointSavedRun writes a phase-only .jointssb record with source identity, build identity, calibration, ROI, parameters, rotation, output sampling, the fit (SSBBranchedFit, both branches when compared) and a phase SHA-256. Its fit options record higherOrderLimit and checkHigherOrderConsistency, and a fifth-order fit its higherOrderConsistency when computed; records saved before them decode with the standard limit 2, unchecked, and none. Records carrying the retired coherence guard (computeHigherOrderGuard, higherOrderGuard) decode and ignore it. Reopening checks source identity, shape and checksum, then displays the exact stored phase without refitting. Records from the earlier native optimizer still decode, phase-variance fits included (SSBObjective.phaseVariance is kept for decoding). MetalSSBSavedRun (.ssb, the retired non-model engine’s complex object) still opens; MetalSSBEngine.phase(of:) converts its object to the phase.

Reproducible checks#

swift test -c release -Xswiftc -enable-testing --filter MetalSSBKernelsTests
PYTHONPATH=src python scripts/ssb_agreement_gate.py        # the same search trials and exact-Newton fits
SSB_PARITY_FIXTURE=... SSB_PARITY_REFERENCE=... swift test -c release -Xswiftc -enable-testing \
  --filter MetalSSBMeasuredParityTests                     # a measured acquisition, kept outside this repository
SSB_DIRECT_FIXTURE=... [SSB_DIRECT_ALL=1] swift test -c release -Xswiftc -enable-testing \
  --filter SSBDirectMeasuredParityTests                    # ssb, phase and complex against CUDA arrays of the same optics

The joint fixtures (SSB_JOINT_REFERENCE, SSB_JOINT_CALIBRATED_REFERENCE) gate the objective (loss relative error < 5e-5). Their phase images were of the retired mean_phase estimator; they are not compared, and scripts/build_metal_ssb_joint_reference.py no longer writes them (nor needs an upsampling source). Missing external fixtures are reported as skips, never passes. Keep acquisitions, screenshots and raw logs outside this repository.

Retired native engine paths#

The engine’s own non-model SSB and the phase-variance fit were removed on 2026-10-25; native SSB is model-based only. MetalSSBEngine keeps Fourier preparation, the joint-model bridge and the saved-run phase conversion. Removed: MetalSSBEngine.reconstruct(aberrations:scanDetectorRotationDegrees:), phaseVariance, optimize, MetalSSBPhaseVarianceResult, SSBOptimizer with SSBOptimizationProgress, the phase-variance branch of SSBModelFit (and its engine: argument), SSBFit.phaseVarianceLoss, the metal-ssb-benchmark product, the standalone tests/metal SSB harnesses and the frozen phase-variance fit trajectory. Their ssb.metal kernels (ssb_gamma_accumulate*, ssb_finalize_fourier_sum, the chi and cross trigonometric tables, the fused 512 and small-scan phase-loss kernels, the phase-moment and Nyquist-correction kernels and the half-plane transposes) are deleted; ssb.metal keeps the count conversion, row FFTs, transpose, Hermitian extraction and ssb_object_phase. These kernels were never registered (docs/performance/kernels.md), so no registered digest changed.

Retired estimators#

mean_phase (the per-BF-pixel phase average, with its general and fused 1536-point kernels, ssb-mean-phase) and complex_wave (the coherent sum’s amplitude at 1x, C10/C12 only) were removed on 2026-10-25 in favour of the regularized phase and complex reconstructions. docs/performance/kernels.md keeps their last registry rows. Python retired the same two estimators.