Changelog

All notable changes to this project will be documented in this file.

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

Unreleased

Removed

• anri.refine (and its tutorial), which was never released: refining a map against the pixels will be redone from a design that represents each population’s orientation spread.

Added

• Fly and helical scans, where dty moves while omega turns: anri.io.read_frame_dty(ds) gives each frame’s dty, the DataSet’s (corrected) dty of its row plus the motion within the row from the sparse file’s raw readings (a straight line over the row’s frames). anri.index.dty_offsets turns it into each row’s offset from its nominal dty per omega bin, and anri.index.system (scan["ddty"]) and backproject put each voxel in the rows that really saw it. python -m anri.index does this by default when dty moves (--no-frame-dty to ignore it); step scans are unchanged. On a helical scan (one dty step per turn) the rows’ frames are up to half a step off the row’s mean dty. --dty-source frames bins each frame by its own dty instead of its row’s: a helical rotation then spreads over 2-3 rows, some (row, omega) bins get no frames and some two rotations’ worth (which made the maps’ rings and holes), so dty_offsets(return_exposure=True) also gives each bin’s exposure, which system (scan["exposure"]) and reconstruct scale by. On W2, holes inside the sample: 391 with the DataSet’s rows and the offsets (the default), 813 binned by each frame’s dty, 577 with the exposure.

• python -m anri.index --accel: SQUAREM-accelerated MLEM in the occupancy fit (anri.index.mlem(accel=True); Varadhan & Roland 2008). Thin features converge in about a quarter of the iterations at the same cost per iteration: on the def316l phantom, 50 accelerated iterations match 200 plain ones (twin lamellae, deviance). The deviance never rises above plain MLEM’s.

• python -m anri.index --beam FWHM: each voxel is spread over the dty rows by the beam’s profile across dty (a Gaussian of that FWHM, integrated over the voxel at each omega, as anri.fwd.beam_weight), over as many rows as it reaches (anri.index.beam_rows), instead of linearly over 2 rows. Default 0: the 2-row model, as before. anri.index.system takes n_beam (the fifth element of dims) and the beam in scan (sig_beam, width_beam, voxel); block_voxels takes the entries per prediction. With FWHM = step the two models are close: on the def316l phantom and on a real AP1_1 scan (scored by how well the maps explain the pixels) the maps are as good either way, and the profile costs ~2.4x in the candidate pass. It is for beams wider than the step (overfocusing), where a voxel reaches more rows.

• python -m anri.index --censor C: ImageD11 keeps only pixels above its segmentation cut, so an empty histogram bin had every pixel below the cut. Where the model predicts fewer than C counts in an empty bin, it counts as agreeing with the model, in both MLEMs (anri.index.censored_ratio) and in the deviance and the pruning’s likelihood ratio (anri.index.deviance). Default 0: empty bins are zeros, as before. orientation_mlem and fit_occupancy take censor; orientation_mlem(return_model=True) also returns the fitted histogram.

• python -m anri.index --mask auto|file.npy: fit only the voxels of a sample mask. auto thresholds a quick reconstruction of the histogram as ImageD11’s tomo_2_map makes its whole-sample mask (anri.index.reconstruct: a sinogram of log intensities, Hamming-ramp filtered back-projection; anri.index.threshold_mask: Otsu’s threshold, the largest connected region, its convex hull). A drawn mask: anri.index.draw_mask(image) (a polygon on the image, as ImageD11’s InteractiveMask), saved as a .npy of booleans in reconstruction order. --occupied then works within the mask, and the results hold mask (and recon with auto).

• indexing_parameters.ipynb: the sample mask (threshold slider or drawn) after the rotation axis; the deviance of the orientation fit against the grid step, to choose how fine a grid the fit needs; --beam (with the weights per row it gives) and --censor in the occupancy fit; --occupied within the mask.

• A notebook for choosing the options of python -m anri.index (docs/source/tutorials/indexing_parameters.ipynb), after ImageD11’s S3DXRD notebooks. It runs the indexer’s steps one at a time on an ImageD11 dataset (or renders the 316L phantom), with a diagnostic plot for each option and sliders where a choice is cheap to redo: the 2θ profile and cake with each ring’s window (--rings, --tth-tol); intensity per dty row and per ω (--monitor); y0 from the sinogram’s centre of mass, checked by a filtered back-projection on the indexer’s voxel grid (--y0); lit area against the intensity it holds, and what the η cut removes (--lit, --etacut); chance completeness of every grid step (--grid, --max-chance); completeness, likelihood ratios and kept orientations in Rodrigues space (--prune, --min-comp, --min-lr, --keep); the deviance per iteration, the voxels whose leading orientation still moves, and the occupancy held by the last candidates (--iter, --cand); data, model and residual per row, per ring and as sinograms; and the sample mask and populations (--occupied, --min-frac). It ends by printing the python -m anri.index command with the values chosen, and writes its results as python -m anri.index does (<tag>_params.toml, .npz, _entries.npz, _tmap.h5); on the phantom they match the command line’s to float rounding. Figures are interactive (ipympl, added to the dev extras) and keep their zoom while the sliders move.

• Detector distortion correction, as ImageD11 corrects its peaks: anri.io.read_spatial(ds) reads a DataSet’s distortion maps (a pyFAI detector file, detectorh5, as ImageD11.blobcorrector.get_e2dx_from_h5; else the e2dxfile/e2dyfile EDF maps; spline files are not supported), read_dataset returns their paths, and stream_sparse(..., spatial=(dx, dy)) yields each pixel at (row + dy, col + dx). python -m anri.index and the indexing notebook use the DataSet’s maps whenever it names them, and log which. A pixel or two of distortion is about a ring’s width, so uncorrected rings were wider than the sample makes them. New dependency: fabio.

• anri.io.stream_sparse takes each frame’s omega and dty from the DataSet (dataset_omega, dataset_dty, scans; read_dataset now returns omega too), where the DataSet covers the frame, and the sparse group’s own readings only elsewhere: the positions ImageD11 has binned the scan by. A DataSet scan can be a slice of a sparse group (1.1::[0:1440]), as for a fly scan in one group split into its rotations. The sparse file’s raw readings of a fly scan, where dty drifts through each rotation, aliased between neighbouring rows (gaps and doubled rows in the sinogram). python -m anri.index passes the DataSet’s.

• anri.index.__main__.save_results: the outputs of python -m anri.index (<tag>.npz, _entries.npz and, with ImageD11, _tmap.h5), written by one function that the command line and the notebook share.

• anri.phantom.polycrystal: options for deformed metals, all off by default and drawn from their own random stream, so existing phantoms are unchanged: cells that accumulate from cell to cell like a random walk (cell_walk_deg, from anri.phantom.brownian_field, a 2D Brownian random field of rotation vectors), an intrinsic orientation spread in every voxel (cell_sig_deg, returned as a sig_rot map for the renderer) and bent grains (bend_grains, bend_deg: lattice curvature about one axis). anri.phantom.tensormap turns a phantom into an ImageD11 TensorMap with its truth maps (grain, cell, twin, misorientation from the grain mean, sig_rot), and anri.io.entries_from_tensormap reads an optional sig_rot map. A deformed phantom whose peaks are smeared arcs (bananas), for testing refinement of sub-grain orientations, is in tests/data/phantoms/def316l.

• python -m anri.index writes <tag>_params.toml: the command line, the anri version and git commit (and whether the code had uncommitted changes), every option with its value (defaults included, so a later change of default does not change the record; options left to be worked out from the data are absent there), the values resolved from the data (paths, phase, lattice, y0, voxel size, grid step, ring tolerances) and, rewritten at the end, the results (chance completeness, the completeness cut used, orientations kept, voxels occupied, run time). It is written once the grid is chosen, so an interrupted run leaves it too. New dependency: tomli-w (and tomli for the tests on Python < 3.11).

• anri.fwd.render_peaks(..., origins=...) takes fixed window origins, and anri.fwd._impl.render.window_origins returns the ones it would choose.

• Crystallography as plain functions in anri.crystal, replacing the classes: lattice_parameters and space_group (of a Dans_Diffraction.Crystal, which reads the CIF), B_matrix, symmetry_matrices, reflections (every hkl up to a d*, without systematic absences, sorted by d*), rings (groups sorted d* into rings; a ring is within the tolerance of its first member, so rings never chain) and structure_factors (|F|², with the thermal-factor warning). They return NumPy arrays in float64; B_matrix evaluates the JAX definition (lpars_to_B) in float64 with anri.crystal.float64(), a context manager for any JAX code. The rings match ImageD11’s unitcell.makerings.

• Monitor normalisation: anri.io.stream_sparse(..., monitor="fpico6") multiplies intensities by monitor_ref / monitor frame by frame (default reference: the counter’s mean), as ImageD11’s DataSet.set_monitor; frames without beam are dropped. A counter the sparse file lacks is read from the raw master file’s scan of the same name (anri.io.read_monitor; read_dataset now returns the masterfile). python -m anri.index --monitor fpico6. Flux varying between the dty rows’ scans otherwise makes ring artefacts.

• python -m anri.index saves the y0 it used (the DataSet’s, or --y0) in its npz, so later steps use the same rotation axis.

• python -m anri.index logs, and saves in its npz, the measured / fitted intensity of each dty row (row_ratio, row_data, row_model). A row that is consistently off makes ring artefacts centred on the rotation axis: steps between rows point to the flux varying between the rows’ scans, a smooth trend with radius to the model. anri.index.fit_occupancy(..., return_model=True) also returns the fitted histogram.

• anri.io.prefetch: read ahead in a background thread, so the next chunk of sparse pixels is read and decompressed while the current one is binned. python -m anri.index uses it.

• anri.index.orientation_mlem: one occupancy per orientation, fitted by MLEM to the histogram with its dty rows summed, and each orientation’s likelihood ratio (how much the fit worsens without it). python -m anri.index now prunes this way by default (--prune likelihood): the grid orientations above chance completeness are fitted, and those with a likelihood ratio above --min-lr (25) are kept. Unlike completeness, this uses the intensities and explains crowded spots jointly. On a crowded phantom with many small grains, 4k orientations kept this way recall 99.6% of the grains (90% of those under 5 µm²), where completeness needs 19k for 98.8% (70%). --prune completeness keeps the old behaviour. With likelihood pruning, --min-comp raises the completeness needed (default: the chance level). That drops the decoys the global fit keeps to soak up what the grid cannot fit, but also small grains.

• Structure factors in indexing: anri.index.ring_table(..., structure) gives each reflection its |F|² from a Dans_Diffraction.Crystal (e.g. read from a CIF), and the fits weight predictions by Lorentz × polarisation × |F|². python -m anri.index --cif.

• anri.index: index scanning-3DXRD data from scratch, giving orientation populations per voxel. The sparse pixels are binned once into a coarse histogram and a row-summed lit map. An orientation grid over the fundamental zone is pruned by completeness, with each prediction’s tolerance computed from the grid spacing, the measured ring widths and the frame step (match_tolerances). Then the kept orientations’ occupancies are fitted per voxel by MLEM, all voxels jointly. Occupancies are sparse (candidates per voxel, from the first MLEM update) and the work runs in blocks of voxels, so memory does not grow with voxels times orientations. Optionally, candidates come from a coarse fit (fit_occupancy(coarse=...)). Each voxel’s occupancy is grouped into populations with a fraction, a mean orientation, a spread and a completeness (populations). python -m anri.index <analysisroot> <sample> <dataset> runs it on an ImageD11 dataset and writes a TensorMap, the populations, and the populations as map entries for the renderer. Tutorials: docs/source/tutorials/phantom.ipynb and indexing.ipynb (stored runs).

• anri.phantom: 2D phantom microstructures (polycrystal): Voronoi grains, cells misoriented a little from their grain, and twin lamellae, on a reconstruction grid. A 316L-like phantom is in tests/data/phantoms/am316l.

• Orientations in anri.crystal: laue_rotations (the proper rotations of a space group’s Laue group, as Cartesian matrices; any crystal system), disorientation, to_fundamental_zone, quaternion and Rodrigues conversions, cubochoric_quaternions (the uniform cubochoric grid of rotations, as orix’s) and orientation_grid. The last covers one fundamental zone: a Rodrigues grid for cubic groups (about half the orientations of a cubochoric grid for the same worst-case spacing), and otherwise the cubochoric grid reduced to the zone plus a thin shell around it, so the zone’s boundary has no gaps. Also allowed_hkls (systematic absences).

• anri.io: read_par, read_dataset (an ImageD11 DataSet’s scan, without ImageD11), read_pars_json, stream_sparse (sparse pixels a chunk at a time, with each frame’s dty row) and tensormap_from_recon (an ImageD11 TensorMap from maps in reconstruction order).

• anri.geom: recon_positions (the sample positions of a reconstruction grid) and sino_shift_and_pad (ImageD11’s padding of a reconstruction).

• Optional per-entry orientation spread "sig_rot" in the renderer’s entries (anri.fwd.render_row, render_peaks): the standard deviation, in radians, of each component of a small isotropic sample-frame rotation of the entry’s lattice. It widens the entry’s peaks in omega and on the detector through the same linearised propagation as the beam’s energy spread and divergence (fine for spreads up to a few degrees), and window size classes (max_frames) account for it. Entries without it render exactly as before, at the same cost.

• anri.fwd.render_row(max_frames=...) (and anri.io.simulate_sparse): peaks broad in omega get windows with more frames, in a few size classes (window[0], 2 window[0] + 1, … up to max_frames), so they are not clipped. stats["window_frames"] gives each peak’s window.

• Optional sig_omega in the geometry (and anri.io.geom_from_pars(..., sig_omega=...)): an extra spread of every peak in omega, in degrees.

• A DCT tutorial (docs/source/tutorials/dct.ipynb): a full 360 degree scan of an 80^3 voxel polycrystal with a box beam and a near-field detector. It is stored with its outputs and not re-executed when the docs are built.

• anri.geom.beam_basis: unit vectors along a beam and across it (horizontal and vertical).

• anri.io.beam_from_pars and geom_from_pars take a beam direction, k_in_lab (default lab x).

• anri.fwd.guess_batch_size: the largest batch for render_row that fits in a fraction (default 25%) of the free GPU or host memory, from XLA’s memory analysis of the compiled render step. Adds psutil as a dependency. Raises RuntimeError where XLA gives no usable memory analysis (e.g. jaxlib 0.4.28 on macOS).

• anri.fwd.check_render: checks rendered peaks of your own map and geometry against a Monte Carlo simulation of the beam spreads through the forward model, and reports the worst cell error and window capture for each peak.

• The renderer is public: anri.fwd.render_row, make_row, select_peaks, render_peaks, beam_weight, lorentz and polarisation (previously only importable from anri.fwd._impl.render).

• anri.io.detector_from_pars, gonio_from_pars and beam_from_pars: the detector, goniometer and beam parts of geom_from_pars, usable on their own (e.g. without the renderer’s spreads). detector_from_pars also returns the pixel-to-lab transforms for anri.geom.det_to_lab, and so does geom_from_pars.

• anri.fwd.get_centroid_box_both (and _all_grains_both, _all_both): both Friedel peak centroids of a box-beam forward projection from one call, like get_centroid_scan_both.

Changed

• python -m anri.index: --iter defaults to 50 MLEM iterations (was 10). Features 1-2 voxels thin, such as twin lamellae, are still converging long after the deviance flattens: at 10 iterations a twin lost to its parent in voxels it fills. On a phantom, more iterations keep improving the map (pure voxels wrong: 0.36% at 10, 0.22% at 30, 0.14% at 100).

• anri.index.histogram_pixels takes a list of ((b_e, b_o, n_e, n_o), n_rows) and fills all of them in one pass over the data. python -m anri.index reads the sparse pixels once again (it read them twice, which doubled the slowest step on large scans), and preparing each chunk is cheaper.

• Tests run on CPU by default, as in CI (conftest.py sets JAX_PLATFORMS=cpu unless it is set): on a GPU they mostly waited for XLA:GPU to compile small one-off shapes.

• anri.fwd.render_peaks weights each frame by the beam at that frame’s own omega (the mean omega of the peak’s mass within the frame), not at the peak’s centroid omega. A peak spread in omega (e.g. by "sig_rot") diffracts at different omegas, where a voxel off the rotation axis sits at a different place across the beam, so its intensity spreads over dty rows along the voxel’s sinusoid. For a 0.5 degree spread 50 um off the axis with a 0.25 um beam, the intensity per row was off by ~150% and is now within ~6% of explicitly rotated sub-entries. select_peaks widens its reach across the beam by the voxel’s motion within the omega margin. Sharp peaks are essentially unchanged.

• anri.fwd.render_peaks centres each frame’s pixel window on the peak’s mean position in that frame (it used to centre every frame on the peak’s overall centroid), so a peak that moves across the detector with omega stays inside its window. Peaks within one or two frames are unchanged.

• anri.fwd.render_row is faster on GPUs: duplicate pixels (neighbouring voxels light up the same pixels, ~100x for a grain) are summed on the device, so only unique pixels go to the host, and batches are 1024 peaks per device times a power of 4, so a scan compiles at most a few shapes (XLA:GPU took up to a minute to compile very small batches). A 153-row scan of a 6818-voxel grain renders in 1 minute instead of 3. Output is unchanged up to float rounding; CPUs keep the host merge, as XLA’s CPU sort is slower than NumPy’s.

• For a beam tilted out of the horizontal plane, the scattering origin in the scanning model is where the pencil crosses the voxel’s column: the beam is taken to cross the rotation axis at the voxel’s own height (the layer’s height in a 2D map), so the origin is raised by (k_z / k_x) x. Nothing changes for a horizontal beam.

• The renderer handles pencil, line and box beams (e.g. DCT) with one model. Every voxel sits at its real lab position for the row’s dty (it used to be moved onto the pencil’s centre line), and anri.fwd.beam_weight (replacing dty_weight) integrates the beam’s profile across it over the voxel: horizontally and vertically, each a flat top blurred by a Gaussian (geom_from_pars(..., sig_beam, width_beam, sig_beam_v, width_beam_v)). Voxels are columns for 2D maps or cubes for 3D maps (voxel_3d). The beam can point anywhere (k_in_lab): it sees the voxel rotated by omega - psi, and travels 1/cos(alpha) further through a column. select_peaks keeps peaks whose voxel is within reach of the beam.

• anri.fwd.polarisation(k_in, k_out, factor) takes the beam direction: horizontal polarisation is across the beam, wherever it points.

Removed

• The crystallography classes anri.crystal.UnitCell, Symmetry, Crystal, Structure and Grain. Use the functions above: e.g. B_matrix(lpars) for Crystal.B, symmetry_matrices(sg) for Symmetry.sym_matrices, reflections + structure_factors + rings for Structure.make_hkls and rings_table. Reflections are no longer filtered by |F|²: keep structure_factors(...) > 0.01 for the old set.

Fixed

• anri.fwd (bin_fractions, truncated_moments, so every rendered pixel): a Gaussian’s mass in a bin far in its upper tail, Phi(b) - Phi(a), cancelled in float32 to about 1e-7 of either sign (2% off at 5 sigma, 0 instead of 5e-10 at 6 sigma). Times a bright peak’s amplitude, and summed over the many voxels sharing a pixel, it made pixels far from peaks wrong, even negative (a log-likelihood of them was NaN). The upper tail now uses Phi(-a) - Phi(-b), and a mass is never negative.

• python -m anri.index: the measured / fitted intensity per dty row (row_ratio) now sums only the eta bins the fit models. It included the data at |sin eta| <= --etacut, which the fit does not predict, so every row read too high.

• anri.io.write_pars writes t_x, t_y, t_z (0) and omegasign (1) when the geometry lacks them: ImageD11 needs them to compute peak geometry (DataSet.update_colfile_pars raised a KeyError on simulated datasets).

• anri.index.inherit_candidates finds each voxel’s coarse neighbours with a KD-tree: sorting every distance took ~2 minutes on a 419 × 419 map, more than the candidate pass that --coarse saves (now ~4 s).

• python -m anri.index writes IPF (ipf_x/y/z) and Euler maps into its TensorMap again, and a ParaView .xdmf beside it, as the sandbox script did. Strain maps are left out: these UBIs are rotations of the nominal lattice, so the strain would be exactly zero.

• anri.fwd.render_row with "sig_rot" failed under newer JAX (shard_map’s check that stacked arrays vary across the same devices): the rotation’s cross-product matrix is now a sum over constant generators.

• Wrong spot positions in float32 on recent NVIDIA GPUs: JAX’s default precision for float32 matrix products there is TF32 (10-bit mantissa), which moved rendered spots ~1000 px from the beam centre by up to 0.3 px on an L40S, a strain error of ~1e-5. anri.utils.setup() now sets jax_default_matmul_precision to "highest" (unless JAX_DEFAULT_MATMUL_PRECISION is set). GPU renders made in float32 without this fix should be redone.

• Beam divergence (ky, kz) was applied along lab y and z, which is only across the beam when it is along lab x: for a tilted beam the vertical divergence shrank by cos(tilt). It is now applied along the beam’s own horizontal and vertical (anri.geom.beam_basis), in every forward model.

• anri.geom.step_grid_from_ybincens always raised: it was jitted, but the size of its grid depends on its inputs.

• render_peaks put up to ~3% of a peak’s intensity in the wrong cells (and lost up to ~7% from broad, truncated peaks): when conditioning fast on slow it held omega at its frame mean, ignoring that slow and omega are correlated within the frame. It now integrates omega out within the frame. Errors against a Monte Carlo of the beam spreads went from 1% (median worst cell) to 0.13%, the Monte Carlo noise.

• Wrong renders on CPU with jaxlib >= 0.11: XLA:CPU’s YNNPACK fusions miscompiled render_peaks for batches of more than a few thousand peaks (in float64, most peaks squeezed into one pixel; in float32, NaNs). anri.utils.setup() now turns them off with --xla_cpu_experimental_ynn_fusion_type=. CPU renders made with jaxlib 0.11.x without this fix should be redone.

• anri.crystal.__all__ and anri.fwd.__all__ listed functions that don’t exist, so from anri.crystal import * failed. anri.fwd now exports propagate_cov_scan_both.