Correction steps

PurposeCopy Anchor Link

tomojax.corrections turns detector counts into line integrals and corrects them. tomojax.Frames.corrected(*steps) runs any steps on counts, then (I - D) / (F - D) with each view's flat interpolated between the flat sets taken around it, then steps on transmission, the -log, and steps on line integrals. tomojax.Scan.corrected(*steps) runs line-integral steps on a scan that is already corrected. Views are corrected on the device a batch at a time. Each step is recorded in Scan.corrections and saved with the scan.

HDF5 and TIFF sample and calibration frames stay lazy until correction, so Frames.cropped reads only that detector window of all three. Flat and dark means are accumulated in float64 over bounded batches. A pixel whose flat does not exceed the dark by epsilon and three standard errors of their means (or whose calibration is not finite) is unmeasured. Interior gaps are interpolated in transmission along the same detector row in each view, before phase retrieval or the log; pixels outside the row's measured span take no absorption. The flat_dark record counts unmeasured pixels. Interpolation cannot recover features narrower than the missing span.

Frames.corrected(..., binning=f) averages blocks of the corrected line integrals on the device, so no unbinned stack reaches the host. When only whole-row steps (Stripes) follow the -log and the source reads row windows cheaply (NumPy, uncompressed TIFF, unfiltered HDF5), it corrects slabs of rows read straight from the source; otherwise the steps that need full frames run unbinned first.

tomojax preprocess is a thin layer over these: its options build the same steps.

Public APICopy Anchor Link

Steps, each a frozen dataclass of its settings:

Also exported:

• Step: the base class. A step has a domain ("counts", "transmission" or "line integrals") and an apply on a (views, rows, columns) batch, or, for a step that selects views, a statistic and a select.
• Correction: the record of one step, its settings and what it found.
• correct_frames and correct_projections: the array-level engines behind Frames.corrected and Scan.corrected, returning a Corrected (the line integrals, the records and the indices of the views kept).

See the quickstart for loading frames and choosing steps.

Stripes(..., varying=True) fits gain and offset only at pixels whose view-to-view differences are anomalous. Unselected centred features are preserved, but their constant-offset stripes are left too; sparse views or sharp moving edges can cause false detections. Selected dead pixels receive only the neighbouring mean; their missing angular information is not recovered.

Rings models each detector column as reporting q + a + b (1 - exp(-q)) and fits a and b against a reference built from the neighbouring columns, then inverts the model (two passes). For a stationary parallel full turn, the reference compares equal ranks of the view-sorted row. The reference narrows from 30 columns either side toward 6 where the sample's own sorted values change faster, so sample edges are not taken for stripes; stuck columns are filled from their neighbours first. Given the scan's geometry, it checks for a stationary parallel-beam scan over one evenly sampled full turn: mirror columns about the rotation axis then see the same lines, the odd part of the fitted errors is pure detector error, and large stripes (detected as in Vo, Atwood and Drakopoulos 2018, algorithms 4 and 5) and broad rings are corrected too, each kept only as far as that odd part supports. The record says which axis was used, if any. The CLI binds --rings to the selected and cropped acquisition geometry and applies it before view rejection and detector binning. Otherwise it uses unsorted local quadratics, independent neighbouring-view instruments and spectral regularisation by angular covariance and object power. Built-in angular geometries supply the view ordering and angular interpolation weights within each trajectory; cone segments stay separate. Without such geometry, supply views in continuous acquisition order. At least two views per trajectory are required. This separation is statistical: fine stationary object structure can still resemble detector errors, and broad defects can be left.

Beam hardeningCopy Anchor Link

A polychromatic beam hardens as it crosses material, so thick paths look less attenuating than they are: cupping and dark streaks between dense features. Through one material the measured line integral is q = g(p) = -ln sum_j w_j exp(-k_j p), where p is the material's line integral and the mixture of rates k_j with weights w_j is the detected spectrum seen through the material. Every curve one material can make has this form, and every such curve is increasing, concave and linear at large p. BeamHardening replaces each q by p = g^-1(q), through a tabulated inverse accurate to float32.

The units of p are declared: with energy_kev=None the curve's slope at zero is one, so thin paths keep their attenuation and the corrected image is in the unhardened beam's mean attenuation; with energy_kev the rates are relative to the material's attenuation at that energy, and so is the image.

The curve comes from an explicit polynomial or spectrum, each only as good as the tube settings and sample material it is given: the scan itself is not used to check it, so a curve computed for the wrong spectrum or material makes the image worse rather than better.

Nikon coefficients, Zeiss constants and Bruker percentages are retained unapplied in scan.source.geometry_metadata["beam_hardening"], alongside tube settings. Their vendor intensity domains and scaling conventions are not established.

• BeamHardening.from_polynomial(coefficients, upto=q_max) turns an explicit polynomial linearisation p = c0 + c1 q + c2 q^2 + ... into the closest physical curve over 0 <= q <= upto.

• BeamHardening.from_spectrum(energies, weights, material) builds the curve of a sample material (a chemical formula or NIST name) seen by a detected spectrum, with attenuation from xraylib (pip install 'tomojax[spectra]'). tube_spectrum(kv, filters=..., detector=...) models a tungsten-target tube's spectrum behind its filters, as detected by a scintillator:

  Copy Text

  from tomojax.corrections import BeamHardening, tube_spectrum
  
  energies, weights = tube_spectrum(
      90, filters=[("Be", 0.3), ("Al", 0.5), ("Cu", 0.2)], detector=("CsI", 0.6, 4.51)
  )
  step = BeamHardening.from_spectrum(energies, weights, "Al", energy_kev=60)

A second material, scatter and photon starvation break the single-material model; one curve cannot correct them. step.curve(p) and step.inverse(q) evaluate the curve on the host for plotting it against data.

Complex transmission imagesCopy Anchor Link

For ptychography or holotomography, unwrap_phase(complex_images) integrates wrapped phase gradients with nonperiodic boundaries on each image. It returns zero-mean radians on the last two axes, (..., rows, columns). Then remove_phase_ramp(phase, air_mask=air) subtracts a plane fitted to known air pixels; the boolean mask or nonnegative weights may be shared or per view. Crop unreliable probe borders first. Air must exclude the sample, and phase residues produce a least-squares integral rather than an exact unwrap. Both functions return JAX arrays and compose under jax.jit; process large stacks in batches of views.

On a grid measured in detector pixels, convert the corrected phase to line integrals of the refractive-index decrement with -phase * wavelength / (2 * np.pi * pixel_size), then build a tj.Scan. For a grid in physical length units, omit the division by pixel_size. Wavelength, pixel size and grid lengths must use consistent units. For joint alignment, keep -phase in radians on a unit-pixel grid and apply the same conversion factor to the reconstructed volume afterwards, keeping the solver data near unit scale.

DependenciesCopy Anchor Link

Allowed: tomojax._data, tomojax.recon, tomojax.forward, tomojax.backends, tomojax.nuisance, tomojax.motion, tomojax.geometry, tomojax.core, JAX and NumPy.

Forbidden: tomojax.io, tomojax.datasets, tomojax.alignment, the CLI and the workflow modules, which sit above it in the import layers (.importlinter).

InvariantsCopy Anchor Link

• Steps run by domain, counts first and line integrals last, each in the order given within its domain.
• Every step applied is recorded, with its settings, in the scan it produced.
• A step that drops views also drops their geometry; selected views keep their place among the flat fields.

TestsCopy Anchor Link

• tests/test_beam_hardening.py covers the curve and its inverse, each coefficient source and the vendor records.
• tests/test_corrections.py covers loading frames (NXtomo, other HDF5 layouts, APS Data Exchange, Nikon, TIFF), flat interpolation, each step, view selection and cropping, and the saved records.
• tests/test_radial_distortion.py checks the lens map against SciPy's interpolation.
• tests/test_correction_medians.py checks the zinger and stripe medians against SciPy and NumPy references.
• tests/test_correction_batch_memory.py covers automatic batch sizing for phase retrieval.