midas-dct-tt

Differentiable forward + inverse for diffraction contrast tomography (DCT) and topotomography (TT), built on midas-dfxm and midas-diffract.

What it is

DCT and TT are not separate physics. They are the same forward operator as DFXM with the objective removed and the scan convention changed:

So this package is a scan-convention layer plus a projection/reconstruction layer over existing MIDAS primitives — not a new physics stack. The reuse map is in implementation_plan.md §5, and the two reuse claims that matter (na = 0 for the acceptance, M = 1 for the projection) are pinned by tests/test_dfxm_contract.py so an upstream release cannot break them quietly.

The one result to know first

With the tomographic rotation axis parallel to G — the TT alignment — the Bragg condition gives, exactly and for every θ:

angle(k_h, a_hat) = 90° − θ

Parallel-beam tomography needs the projection direction perpendicular to the rotation axis. TT is off by exactly θ, so TT is laminography with a missing cone of half-angle θ, not tomography. Consequences: low-θ (high-energy, low-index) reflections reconstruct better; Friedel-paired settings flip the cone and recover part of it; a reconstructed sphere elongates along the axis.

midas_dct_tt.missing_cone_half_angle_deg carries the full statement; the identity is verified from vectors, independently of the θ that produced them, in tests/test_geometry.py.

Install (editable, from the repo)

pip install -e packages/midas_dct_tt --no-deps    # siblings already editable

Quickstart

import torch
from midas_dct_tt import sphere_grain, attach_uniform_field, tt_alignment

# A planted grain: soft voxel occupancy chi(r) + a perfect-crystal F(r).
grain = attach_uniform_field(sphere_grain(radius_um=5.0, spacing_um=0.5))

# Align its (111) so G lies on the tomographic rotation axis.
G0 = grain.field.reference_G((1, 1, 1))          # sample frame, |G| = 2*pi/d
al = tt_alignment(G0, wavelength_A=0.172979)      # ~71.7 keV, 1-ID/HEXM

print(f"theta        = {float(al.theta_deg):.3f} deg")
print(f"missing cone = {float(al.missing_cone_deg()):.3f} deg half-angle")
print(f"k_h to axis  = {float(al.axis_beam_angle_deg()):.3f} deg  (= 90 - theta)")

# The TT scan: G is invariant, so k_h never moves.
for psi in (0.0, 90.0, 180.0, 270.0):
    lab_positions = grain.positions_lab(al.sample_rotation(psi))

Phase 1 — a topograph, and a full sinogram:

from midas_dct_tt import PlaneDetector, psi_scan, topograph_image, topograph_stack, tt_resolution

det = PlaneDetector(pixel_um=0.5, shape=(96, 96), distance_um=5000.0)

# Ideal grain: the pixel reads path length in micrometers (centre = 2R).
img = topograph_image(grain, al, psi_deg=0.0, detector=det)

# With acceptance and deformation: intensity AND position respond to F(r).
res = tt_resolution(al)                       # DFXM's resolution at na = 0
sino = topograph_stack(grain, al, psi_scan(180), detector=det,
                       hkl=(1, 1, 1), resolution=res)     # (180, 96, 96)

DCT instead — solve the Bragg flashes, then render frames:

from midas_dct_tt import bragg_flashes, dct_frames, dct_omega_scan

flashes = bragg_flashes(G0, wavelength_A=0.172979)   # 2 per 360 deg, or [] if blind
frames = dct_frames([grain], [flashes], wavelength_A=0.172979,
                    detector=PlaneDetector(pixel_um=1.0, shape=(192, 192), distance_um=500.0),
                    omega_centres=dct_omega_scan(720))

Phase 2 — reconstruct the grain shape from the sinogram:

from midas_dct_tt import forward_operator, sirt, reconstruct_differentiable, dice

A = forward_operator(grain.positions, al, psi_scan(36), det,
                     voxel_volume_um3=grain.voxel_volume_um3)
sino = A(grain.occupancy)                       # linear in chi: no sigmoid, no acceptance

chi_sirt = sirt(sino, A, grain.n_voxels, n_iter=60)          # classical baseline
chi_grad, info = reconstruct_differentiable(sino, A, grain.n_voxels, steps=300)
print(dice(chi_sirt, grain.occupancy), dice(chi_grad, grain.occupancy))

Conventions

Inherited unchanged from midas_dfxm so a grain, a field, and a reflection mean the same thing in both packages:

• lab x along the incident beam, z vertical up, y = z × x;
• k_in = (2π/λ) x̂, and reciprocal vectors carry the 2π convention (|G| = 2π/d);
• v_lab = R @ v_sample;
• micrometers, degrees, Ångström (wavelength and lattice only).

Two conventions are hazards rather than choices, and both are named constants with loud docstrings rather than literals:

• DCT ω sign. The 1-ID aero stage is clockwise, so recorded ω must be negated (DCT_OMEGA_SIGN_AERO). A wrong sign here does not degrade any residual — it mirrors the reconstruction, undetectably.
• Soft-boundary volume bias. Occupancy is sigmoid(−sdf/w), which over-counts volume on a convex surface by A·(2H)·w²·π²/6. Measured, ∝ w², and it inverts below w ≈ 0.25 × spacing where the grid stops resolving the ramp. Do not quote a grain volume from a soft occupancy without correcting or stating w. See logits_from_signed_distance.

Status

Phases 0-3 complete — conventions/geometry/containers; the forward model (TT topographs, DCT frames, acceptance, extinction, projector); shape reconstruction (SIRT + differentiable, Friedel pairing, missing-cone study); and the deformation inverse with a pre-registered identifiability study. 343 tests, all three device backends verified: CPU + MPS on macOS, CUDA on an H200 (sentosa, torch 2.11/cu128, PYTHONPATH overlay on the shared MIDAS env).

Validated against two independent implementations:

• skimage.transform.radon for the projector — 4e-12 relative at 0°/90°, ~0.2% RMS elsewhere (its interpolation blur, which we don't incur).
• pymicro for the DCT scan geometry — our Bragg-flash angles match Orientation.dct_omega_angles to 1.1e-13 degrees across 3 orientations × 3 reflections, with no fitted constant.

plus an internal NumPy oracle to 1e-13. A sphere projects to its chord, 2R, and projected mass is conserved exactly.

Phase 3 refuted this package's original novelty claim. Exact finite-strain inversion is correct — its difference from the linearised model is exactly the omitted O(|H|²) term, verified to better than 1% — but the difference reaches only 2.2% at |H| = 5% and needs |H| ≈ 24% to reach 10%, below the measurement floor throughout. Separately, a fixed TT setting goes dark at |H| ≈ 2e-3 (rotation), 13× before the linearisation matters. Pre-registered in dev/paper/PREREGISTER.md; full result in dev/paper/RESULTS_phase3.md.

What survives: a validated differentiable forward model, a clean identifiability result (one reflection constrains 3 of 9 components of F; three non-coplanar give all 9; collinear add nothing), and the finding that TT's spatial resolution resolves a strain sign that intensity alone cannot.

Three results that change how you should use this package:

• A Friedel partner is a lab-frame mirror, not an inversion. The two lab-frame G vectors share (x, y) and have opposite z (cos = +0.79). Pair in the sample frame — a lab-frame antiparallel test returns zero pairs and fails silently. BraggFlash carries G_sample for exactly this.
• Do not measure the missing cone by reconstructing a sphere. It detects nothing (Dice 1.0000, zero elongation, at any θ) because same-grid noiseless data is an inverse crime. Measure the operator's Fourier response instead — tests/test_missing_cone.py shows ~10× suppression inside the cone with the width tracking θ.

And two from Phase 1 about reflection selection:

• Low θ trades tomographic coverage against strain contrast. Both scale as cot θ: fcc (111) at 71.7 keV gives a 2.4° missing cone but σ_par/|Q| = 6.4e-3; at 17 keV, 10.0° and 1.5e-3. Neither choice dominates, so reflection selection wants Fisher-information design rather than a rule of thumb.
• TT strain contrast is divergence-limited, not bandwidth-limited at HEXM energies (the cot²θ·div² term beats 4ε² by ~2000×). A narrower monochromator will not sharpen it; collimation will.

Tests

KMP_DUPLICATE_LIB_OK=TRUE python -m pytest tests/ -q

Markers: unit, contract (sibling API surfaces), device (CPU/CUDA/MPS), autograd (gradcheck), slow (deselected by default; run with -m slow).

tests/test_radon_crosscheck.py needs scikit-image and skips without it. tests/test_pymicro_crosscheck.py uses a frozen reference table rather than importing pymicro, which pins numpy<2; regenerate with examples/pymicro_reference.py.