tomojax.recon

Names resolve on first use, so a Fourier reconstruction does not import JAX and an iterative solve does not import the Fourier/CuPy stack.

projectCopy Anchor Link

project(
    geometry: Geometry,
    volume: jax.Array | np.ndarray,
    *,
    grid: Grid | None = None,
    devices: Device | Sequence[Device] | None = None,
) -> jax.Array

Project a volume through the scan geometry.

volume is (nx, ny, nz) on grid (the geometry's own by default), in the reciprocal of the geometry's length unit. Works for parallel, laminography and cone-beam geometries, with per-view poses and a rolled detector; backproject is its exact transpose. devices (one or several, jax.devices() for all) shares the views among them, each holding the whole volume; the projections come back on the first.

ParametersCopy Anchor Link

geometry : Geometry Parallel, laminography or cone-beam scan geometry, with any per-view poses. volume : array-like Volume of shape (nx, ny, nz) in reciprocal geometry-length units. grid : Grid or None, optional Volume grid. None uses geometry.grid. devices : jax.Device or sequence of jax.Device or None, optional Devices sharing the views. None uses the current default device.

ReturnsCopy Anchor Link

jax.Array Line integrals of shape (views, nv, nu).

RaisesCopy Anchor Link

ValueError If the array shape does not match the grid or scan geometry.

backprojectCopy Anchor Link

backproject(
    geometry: Geometry,
    projections: jax.Array | np.ndarray,
    *,
    grid: Grid | None = None,
    devices: Device | Sequence[Device] | None = None,
) -> jax.Array

Apply the exact transpose of the forward projector.

This is the exact adjoint, not a reconstruction: use tomojax.reconstruct for a volume. devices shares the views among them, as for project, and sums their backprojections.

ParametersCopy Anchor Link

geometry : Geometry Parallel, laminography or cone-beam scan geometry, with any per-view poses. projections : array-like Projection stack of shape (views, nv, nu). grid : Grid or None, optional Volume grid. None uses geometry.grid. devices : jax.Device or sequence of jax.Device or None, optional Devices sharing the views. None uses the current default device.

ReturnsCopy Anchor Link

jax.Array Adjoint volume of shape (nx, ny, nz), not a reconstruction.

RaisesCopy Anchor Link

ValueError If the array shape does not match the grid or scan geometry.

CGLSConfigCopy Anchor Link

class CGLSConfig(
    *,
    iterations: int = 50,
    rtol: float = 1e-06,
    atol: float = 0.0,
    damping: float = 0.0,
    views_per_batch: int = 64,
    backend: Literal['auto', 'jax', 'cuda'] = 'auto',
    projector_model: Literal['auto', 'ray', 'joseph'] = 'auto',
    joseph_interpolation: Literal['linear', 'cubic'] = 'linear',
    gradient_damping: float = 0.0,
    stream_projections: bool | None = None,
    devices: Device | Sequence[Device] | None = None,
)

Control an unconstrained least-squares solve.

Minimize ||A x - y||² + damping² ||x||² + gradient_damping² ||D x||² using the discrete matched adjoint. Stop when the normal-residual norm falls below atol + rtol * initial_normal_residual_norm or after iterations updates. roundoff_limit indicates voxelwise update stagnation or a normal gradient stagnating within an FP32 cancellation estimate, without claiming the requested tolerance. Before reporting convergence or roundoff, recompute the data and normal residuals. A premature convergence indication from recursive residual drift restarts the conjugate direction. Periodic checks near the FP32 noise floor also catch recurrences kept moving by nondeterministic atomic reductions. This criterion measures stationarity, not reconstruction quality. No positivity clipping or TV regularization is applied. auto uses CUDA on CUDA with a canonical detector grid, otherwise JAX. An explicit CUDA request must be supported and never silently falls back.

When both penalties are zero, updates lie in the range of the matched adjoint, up to floating-point error. A converged zero-start solve therefore selects the minimum Euclidean-norm least-squares solution. A supplied initial volume retains its component in the projector's null space: measurements cannot determine that component. A small residual alone does not establish image accuracy, particularly in sparse or tilted scans.

D takes adjacent voxel differences divided by physical voxel spacing, with no edges outside the volume. gradient_damping adds quadratic smoothness; it does not impose positivity or penalize a constant volume. Both penalties default to zero. Their weights are independent of view count; the data term is an unnormalized sum of squared residuals.

projector_model="joseph" samples voxel-centre planes along the dominant voxel direction, interpolates each plane and uses its matched transpose. On CUDA that transpose gathers into voxels without atomic writes. Bilinear interpolation is the default; joseph_interpolation="cubic" selects Keys cubic convolution with a=-1/2 and a 4-by-4 stencil, including negative weights. "ray" is the trilinear ray marcher, the only model supporting explicit detector grids. "auto" selects Joseph unless det_grid is given; both models match analytic line integrals equally well, and Joseph is faster.

stream_projections solves the equivalent normal equations with NumPy or memmap projections read from host memory one view batch at a time, so only volume-sized arrays occupy the device; None reserves space for the resident solver state and verification scratch, streaming when they would exceed free device memory. The normal residual then follows a recurrence between recomputations, each of which streams the data once.

devices (two or more JAX devices) shares the views among them, each holding the whole volume, and sums their backprojections; the projections are then held in device memory, never streamed.

ParametersCopy Anchor Link

iterations : int, optional Maximum number of solver updates. For SIRT and MLEM, an update is a complete subset sweep. Default is 50. rtol : float, optional Nonnegative relative tolerance on the normal-residual norm. Default is 1e-06. atol : float, optional Nonnegative absolute tolerance on the normal-residual norm. Default is 0.0. damping : float, optional Nonnegative weight of the volume-norm penalty; squared in the objective. Default is 0.0. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is 64. backend : Literal['auto', 'jax', 'cuda'], optional Execution backend: "auto", "jax" or "cuda". An explicit accelerator request never silently falls back. Default is 'auto'. projector_model : Literal['auto', 'ray', 'joseph'], optional Projection model: "auto", "ray" or "joseph". Auto selects Joseph unless an explicit detector grid requires rays. Default is 'auto'. joseph_interpolation : Literal['linear', 'cubic'], optional Joseph interpolation: "linear" or Keys "cubic" with a=-1/2. Cubic is not supported by the ray model. Default is 'linear'. gradient_damping : float, optional Nonnegative weight of the physical finite-difference penalty; squared in the objective. Default is 0.0. stream_projections : bool | None, optional Read host projections one view batch at a time. None selects streaming from available device memory; False keeps them resident. devices : jax.Device or sequence of jax.Device or None, optional JAX device or sequence of devices. None uses the current default; multiple devices divide the views and each hold the full volume.

ConeAxisCalibrationCopy Anchor Link

class ConeAxisCalibration(
    axis_offset: float,
    detector_roll_deg: float,
    heights: tuple[float, ...],
    slab_offsets: tuple[float, ...],
)

Calibrated axis offset and detector roll, with the per-slab estimates.

axis_offset and detector_roll_deg are absolute ConeBeam values. slab_offsets[i] is the sharpest axis offset of the slab centred at height heights[i] (physical z), with the calibrated roll applied.

ParametersCopy Anchor Link

axis_offset : float Absolute rotation-axis offset in the geometry length unit. detector_roll_deg : float Absolute detector roll in degrees. heights : tuple[float, ...] Physical z coordinates of the slab centres. slab_offsets : tuple[float, ...] Sharpest axis offset at each slab height with the calibrated roll applied.

ConeAxisCalibration.applyCopy Anchor Link

apply(geometry: ConeGeometry) -> ConeGeometry

Return geometry with the calibrated axis offset and detector roll.

ParametersCopy Anchor Link

geometry : ConeGeometry Geometry to update, without changing its other scan settings.

ReturnsCopy Anchor Link

ConeGeometry Copy with the calibrated absolute axis offset and detector roll.

ConeAxisConfigCopy Anchor Link

class ConeAxisConfig(
    *,
    search_range: float | None = None,
    estimate_roll: bool = True,
    slices: int = 32,
    max_roll_deg: float = 3.0,
    fdk: FDKConfig = (lambda: FDKConfig(filter='hann', views_per_batch=1024))(),
)

Options for calibrate_cone_axis.

search_range bounds the axis offset searched on either side of the geometry's current one, in physical units (default a quarter of the volume width). estimate_roll measures the offset at three heights and fits the detector roll; otherwise only the central slab is used. slices is the thickness of each slab in full-resolution voxels (at most an eighth of the volume), max_roll_deg the largest roll the slabs leave detector rows for, and fdk the filter (Hann by default, so noise does not swamp the sharpness) and backend.

ParametersCopy Anchor Link

search_range : float | None, optional Half-width of the axis-offset search in physical length units. None uses a quarter of the volume width. estimate_roll : bool, optional Estimate roll from three slab heights; False uses only the central slab. Default is True. slices : int, optional Requested slab thickness in full-resolution voxels. Default is 32. max_roll_deg : float, optional Maximum detector roll in degrees allowed for detector-row planning. Default is 3.0. fdk : FDKConfig, optional FDK filter, batching and backend configuration.

FBPConfigCopy Anchor Link

class FBPConfig(
    *,
    filter: str = 'ramp',
    scale: float | None = None,
    views_per_batch: int = 1,
    backend: Literal['auto', 'jax', 'cuda'] = 'auto',
    devices: Device | Sequence[Device] | None = None,
)

Configuration for filtered backprojection.

backend='auto' uses the voxel-driven CUDA kernel on CUDA and JAX otherwise. 'jax' forces the reference; 'cuda' requires CUDA input and never silently changes backend.

scale=None weights each view exactly for the scanned arc (see fbp). A number instead weights every view by that constant and filters along u only, the convention before exact weights; pi / n gives the uniform half-turn result. views_per_batch bounds the views held at once by the host and explicit det_grid paths.

devices (one or several JAX devices) shares the views among them: each filters and backprojects its own, weighted for the whole scan, and the volumes are summed on the first. Not with an explicit det_grid.

ParametersCopy Anchor Link

filter : str, optional Frequency filter: "ramp", "shepp-logan" or "hann". Default is 'ramp'. scale : float | None, optional Constant angular weight per view, or None for geometry-derived weights. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is 1. backend : Literal['auto', 'jax', 'cuda'], optional Backprojection backend: "auto", "jax" or "cuda". Auto uses CUDA acceleration when available. Default is 'auto'. devices : jax.Device or sequence of jax.Device or None, optional JAX device or sequence of devices. None uses the current default; multiple devices divide the views and each hold the full volume.

FBPHostConfigCopy Anchor Link

class FBPHostConfig(
    *,
    slices_per_batch: int | None = None,
    views_per_batch: int = 32,
    filter: str = 'ramp',
    scale: float | None = None,
    backend: Literal['auto', 'jax', 'cuda'] = 'auto',
)

Control slab storage and filter batches for fbp_host.

The output remains in host storage, and so does the input unless it is a device array, which stays resident and is sliced there. For built-in parallel geometry, a slab is slices_per_batch axial (z) slices with only the detector rows they need (balanced slabs of at most 16 when None). For laminography and other geometry, a slab is slices_per_batch x-slices; every view batch passes through the device once per slab. None sizes slabs to about half of the free device memory, so a volume that fits is reconstructed in one pass. Filtering uses views_per_batch views at a time. Axial CUDA slabs with many batches can stage up to 512 MiB of detector rows when those fit in a quarter of free device memory, avoiding per-batch host dispatch. Otherwise rows are read and transferred with at most two input batches on the device. Runtime/compiler caches add memory beyond these arrays. General slabs shrink their depth, then their view batch, after device allocation exhaustion. scale has the same meaning as in FBPConfig.

ParametersCopy Anchor Link

slices_per_batch : int | None, optional Positive slab depth, or None to choose automatically as described above. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is 32. filter : str, optional Frequency filter: "ramp", "shepp-logan" or "hann". Default is 'ramp'. scale : float | None, optional Constant angular weight per view, or None for geometry-derived weights. backend : Literal['auto', 'jax', 'cuda'], optional Backprojection backend: "auto", "jax" or "cuda". Auto uses CUDA acceleration when available. Default is 'auto'.

FDKConfigCopy Anchor Link

class FDKConfig(
    *,
    filter: str = 'ramp',
    backend: str = 'auto',
    views_per_batch: int = 64,
    devices: Device | Sequence[Device] | None = None,
)

FDK options.

filter is ramp, shepp-logan or hann. backend auto uses the CUDA kernel when CuPy and a CUDA device are available. views_per_batch bounds the filtered views held on the device at once. devices (one or several JAX devices) shares the views among them: each filters and backprojects its own, weighted for the whole scan, and the volumes are summed on the first. fdk_host gives each its own slabs.

ParametersCopy Anchor Link

filter : str, optional Frequency filter: "ramp", "shepp-logan" or "hann". Default is 'ramp'. backend : str, optional Execution backend: "auto", "jax" or "cuda". Auto uses CUDA when available. Default is 'auto'. views_per_batch : int, optional Maximum views processed together. Batches also cap FP32 input storage at 64 MiB (or one view) to bound host transfer staging. Default is 64. devices : jax.Device or sequence of jax.Device or None, optional JAX device or sequence of devices. None uses the current default; multiple devices divide the views and each hold the full volume.

FDKHostConfigCopy Anchor Link

class FDKHostConfig(
    *,
    slices_per_batch: int | None = None,
    fdk: FDKConfig = FDKConfig(),
)

Slab options for fdk_host.

slices_per_batch z slices are reconstructed per slab; None sizes slabs to about a third of the free device memory. fdk holds the filter and backend options.

ParametersCopy Anchor Link

slices_per_batch : int | None, optional Positive slab depth, or None to choose automatically as described above. fdk : FDKConfig, optional FDK filter, batching and backend configuration.

FistaConfigCopy Anchor Link

class FistaConfig(
    *,
    iterations: int = 50,
    initial_epochs: int = 0,
    initial_subsets: int = 16,
    tv_weight: float = 0.005,
    regulariser: Regulariser = 'tv',
    huber_delta: float = 0.01,
    lipschitz: float | None = None,
    views_per_batch: int | None = None,
    projector_unroll: int = 1,
    checkpoint_projector: bool = True,
    gather_dtype: str = 'fp32',
    grad_mode: GradMode = 'auto',
    tv_prox_iterations: int = 10,
    recon_rel_tol: float | None = None,
    recon_patience: int = 0,
    power_iterations: int = 3,
    support: jnp.ndarray | None = None,
    nonnegative: bool = False,
    lower_bound: float | None = None,
    upper_bound: float | None = None,
    ray_integrator: str = 'sampled',
    projector_model: ProjectorModel = 'auto',
    backend: ProjectorBackend = 'auto',
    stream_projections: bool | None = None,
    devices: Device | Sequence[Device] | None = None,
)

Configuration for public FISTA/TV reconstruction.

projector_model/backend follow CGLSConfig: "auto" uses batched Joseph (also for ray_integrator="joseph"), with CUDA on CUDA. Explicit grids/exact integration use the reference path, alone honouring gather_dtype, projector_unroll, checkpoint_projector and grad_mode. views_per_batch=None batches 64 views, or one on the reference path.

stream_projections reads NumPy or memmap projections from host memory one view batch at a time, so they never occupy the device whole. None streams when they would take more than 40% of free device memory; True always streams host arrays. Streaming needs the batched operators.

lipschitz is the data term's Lipschitz constant; None estimates it with power_iterations power iterations, started from the backprojected data. Huber-TV's finite-grid gradient bound is added automatically.

initial_epochs optionally initializes nonnegative least squares with accelerated, interleaved ordered subsets before the full-data iterations. Each epoch interleaves nominal angle order into initial_subsets groups; more groups save time but increase gradient error. This needs resident, unmasked parallel or laminographic Joseph data, tv_weight=0, and at least one final iteration. It is not a convergence or image-quality guarantee. Defaults leave it off. Its norm bounds are reused unless lipschitz is supplied. Unmasked resident nonnegative Joseph volumes larger than the data reuse two owned buffers for zero or ordered starts; this changes FP32 summation order.

devices (two or more JAX devices) shares the views among them, each holding the whole volume, and sums their backprojections; one device, or None, runs on the default device. It needs the batched operators and projections held in device memory (no streaming).

ParametersCopy Anchor Link

iterations : int, optional Maximum number of solver updates. For SIRT and MLEM, an update is a complete subset sweep. Default is 50. initial_epochs : int, optional Number of ordered-subset initialization sweeps; zero disables initialization. Default is 0. initial_subsets : int, optional Number of angle-interleaved groups per initialization sweep. Default is 16. tv_weight : float, optional Nonnegative regularization weight. Default is 0.005. regulariser : Regulariser, optional Regularization model: "tv" or "huber_tv". Default is 'tv'. huber_delta : float, optional Positive Huber transition scale when regulariser is "huber_tv". Default is 0.01. lipschitz : float | None, optional Positive data-gradient Lipschitz constant, or None to estimate it by power iteration. views_per_batch : int | None, optional View batch size. None uses 64 for batched operators or one on the reference path. projector_unroll : int, optional Loop unrolling for the explicit-grid or ray-reference projector. Default is 1. checkpoint_projector : bool, optional Recompute intermediate projector values during differentiation on the reference path. Default is True. gather_dtype : str, optional Gather precision on the reference path: "auto", "fp32", "fp16" or "bf16". Default is 'fp32'. grad_mode : GradMode, optional Reference gradient evaluation: "auto", "batched" or "stream". Default is 'auto'. tv_prox_iterations : int, optional Iterations of the TV proximal solve per FISTA update. Default is 10. recon_rel_tol : float | None, optional Relative objective-change tolerance for stopping reconstruction early; None disables it. recon_patience : int, optional Number of successive small-change iterations required for early stopping. Default is 0. power_iterations : int, optional Number of power iterations used to estimate the data-gradient Lipschitz constant. Default is 3. support : array-like | None, optional Volume support mask. Excluded voxels remain zero; None uses the entire grid. nonnegative : bool, optional Constrain the reconstruction to nonnegative values. Default is False. lower_bound : float | None, optional Lower voxel-value bound, or None for no explicit lower bound. upper_bound : float | None, optional Upper voxel-value bound, or None for no upper bound. ray_integrator : str, optional Reference ray integration rule; "sampled" is the default, "exact" uses exact integration, and "joseph" selects Joseph sampling. Default is 'sampled'. projector_model : ProjectorModel, optional Projection model: "auto", "ray" or "joseph". Auto selects Joseph unless an explicit detector grid requires rays. Default is 'auto'. backend : ProjectorBackend, optional Execution backend: "auto", "jax" or "cuda". An explicit accelerator request never silently falls back. Default is 'auto'. stream_projections : bool | None, optional Read host projections one view batch at a time. None selects streaming from available device memory; False keeps them resident. devices : jax.Device or sequence of jax.Device or None, optional JAX device or sequence of devices. None uses the current default; multiple devices divide the views and each hold the full volume.

FourierConfigCopy Anchor Link

class FourierConfig(
    *,
    slices_per_batch: int | None = None,
    backend: Literal['auto', 'numpy', 'cupy'] = 'auto',
)

Choose a portable NumPy reference or CUDA execution with fixed axial slabs.

auto uses CuPy when its current CUDA device matches the selected JAX device and NumPy otherwise, including when CuPy is not installed. Explicit numpy and cupy never switch backends or change the current device. cupy uses single-precision FFTs and an original CUDA interpolation kernel. It requires the optional fourier-cuda12 extra. numpy uses double-precision FFTs and a vectorized interpolation reference. Both return float32 host arrays. slices_per_batch=None chooses up to 16 slices, targeting a 256 MiB padded CUDA FP32 plane stack (at least one slice). NumPy defaults to 16; an explicit integer fixes the requested depth. CUDA allocation exhaustion still halves it down to one slice. FFT plans and allocator caches consume additional memory. Large CUDA calls overlap host row preparation, device transfers, FFT execution and output writes, retaining at most two pending device slabs. A bounded cache retains immutable geometry arrays; each call owns its working buffers. Fractional detector rows are interpolated on CUDA from at most two measured source rows per requested slice; entire projection stacks are not uploaded.

ParametersCopy Anchor Link

slices_per_batch : int | None, optional Positive slab depth, or None to choose automatically as described above. backend : Literal['auto', 'numpy', 'cupy'], optional Execution backend: "auto", "numpy" or "cupy". Explicit choices never change backends. Default is 'auto'.

MLEMConfigCopy Anchor Link

class MLEMConfig(
    *,
    iterations: int = 50,
    subsets: int = 1,
    support: jax.Array | np.ndarray | None = None,
    views_per_batch: int = 32,
    backend: ProjectorBackend = 'auto',
    projector_model: ProjectorModel = 'auto',
)

Reconstruct nonnegative linear Poisson measurements, y ~ Poisson(A x).

subsets=1 is MLEM; larger values use ordered-subsets EM (OSEM). Each iteration visits every view once, in angle-interleaved groups. OSEM accelerates early progress but can cycle; use one subset for convergence. support is a boolean mask broadcastable to the volume shape: voxels outside it stay zero. The default initial volume is one inside the observed support. A supplied zero voxel stays zero under multiplicative updates. tj.reconstruct(..., warm_start=True) instead starts from nonnegative FBP/FDK, including its zero voxels. This can improve early image quality, but fixes the initial zero set for all subsequent EM updates.

These are linear emission-style measurements, not transmission counts with mean I0 * exp(-A x). Applying MLEM to nonnegative log-transformed attenuation data is possible, but does not give it Poisson statistics. The projectors are the same matched pair as CGLSConfig, with nonnegative interpolation weights. A single sensitivity is retained; multiple subset sensitivities are cached within 256 MiB, or in available memory after reserving working buffers when reused across multiple sweeps; otherwise they are recomputed to bound memory use. Exactly independent axial planes share their single-plane sensitivities.

ParametersCopy Anchor Link

iterations : int, optional Maximum number of solver updates. For SIRT and MLEM, an update is a complete subset sweep. Default is 50. subsets : int, optional Number of interleaved view subsets visited per sweep; one uses the simultaneous method. Default is 1. support : array-like | None, optional Boolean support mask broadcastable to the volume shape; excluded voxels remain zero. None uses the entire grid. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is 32. backend : ProjectorBackend, optional Execution backend: "auto", "jax" or "cuda". An explicit accelerator request never silently falls back. Default is 'auto'. projector_model : ProjectorModel, optional Projection model: "auto", "ray" or "joseph". Auto selects Joseph unless an explicit detector grid requires rays. Default is 'auto'.

SIRTConfigCopy Anchor Link

class SIRTConfig(
    *,
    iterations: int = 50,
    subsets: int = 1,
    relaxation: float = 1.0,
    acceleration: bool = False,
    lower_bound: float | None = None,
    upper_bound: float | None = None,
    support: jax.Array | np.ndarray | None = None,
    views_per_batch: int = 64,
    backend: ProjectorBackend = 'auto',
    projector_model: ProjectorModel = 'auto',
)

SIRT, or OS-SART with subsets > 1 (one view per subset is SART).

Each iteration visits all interleaved subsets once, in input view order. A subset updates x += relaxation * C * A.T * R * (data - A*x), where R = 1 / (A*support) and C = 1 / (A.T*ones); zero denominators give zero weights. The forward operator and transpose are exactly matched. relaxation must be in (0, 2); ordered subsets can cycle on inconsistent data, so use fewer subsets or lower relaxation for the final refinement.

acceleration=True applies fixed FISTA momentum to full-data SIRT (subsets=1, relaxation <= 1), reducing the iterations needed for the same weighted least-squares objective without changing the projector.

lower_bound and upper_bound clip each update. support is a volume mask: zero excludes a voxel from reconstruction and fixes it to zero, including in the initial volume; nonzero includes it. Bounds apply only inside that support. Only nonnegative projector weights are used.

ParametersCopy Anchor Link

iterations : int, optional Maximum number of solver updates. For SIRT and MLEM, an update is a complete subset sweep. Default is 50. subsets : int, optional Number of interleaved view subsets visited per sweep; one uses the simultaneous method. Default is 1. relaxation : float, optional Update scale strictly between 0 and 2. Default is 1.0. acceleration : bool, optional Apply FISTA momentum to full-data SIRT; requires subsets=1 and relaxation <=

1. Default is False. lower_bound : float | None, optional Lower voxel-value bound, or None for no explicit lower bound. upper_bound : float | None, optional Upper voxel-value bound, or None for no upper bound. support : array-like | None, optional Volume support mask. Excluded voxels remain zero; None uses the entire grid. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is

backend : ProjectorBackend, optional Execution backend: "auto", "jax" or "cuda". An explicit accelerator request never silently falls back. Default is 'auto'. projector_model : ProjectorModel, optional Projection model: "auto", "ray" or "joseph". Auto selects Joseph unless an explicit detector grid requires rays. Default is 'auto'.

SPDHGConfigCopy Anchor Link

class SPDHGConfig(
    *,
    iterations: int = 400,
    tv_weight: float = 0.005,
    regulariser: Regulariser = 'tv',
    huber_delta: float = 0.01,
    theta: float = 1.0,
    views_per_batch: int = 16,
    seed: int = 0,
    tau: float | None = None,
    sigma_data: float | None = None,
    sigma_tv: float | None = None,
    projector_unroll: int = 1,
    checkpoint_projector: bool = True,
    gather_dtype: str = 'fp32',
    ray_integrator: str = 'sampled',
    projector_model: ProjectorModel = 'auto',
    backend: ProjectorBackend = 'auto',
    nonnegative: bool = False,
    support: jnp.ndarray | None = None,
    log_every: int = 10,
    stream_projections: bool | None = None,
)

Configuration for stochastic primal-dual TV reconstruction.

projector_model and backend choose the projection operator as in FistaConfig: "auto" uses Joseph plane sampling with CUDA on CUDA unless an explicit detector grid or the exact ray integrator requires the ray-model reference path. gather_dtype, projector_unroll and checkpoint_projector apply to the reference path only.

stream_projections keeps NumPy or memmap projections, any weights and the sinogram-sized dual variable in host memory; each iteration reads and writes only its block of views. None streams when the projections would take more than 40% of free device memory. Streaming needs the batched operators.

tau, sigma_data and sigma_tv independently override the primal, data-dual and TV-dual step sizes. None estimates only that missing size; explicit sizes must be finite and positive. Automatic sizes bound every block's squared norm divided by its selection probability. Batched positive-kernel operators use three upper-bound passes per block; reference rays use a one-pass row-sum bound per block. info["A_norm"] is the square root of this probability-scaled bound, also bounding the full-scan operator norm. Supplying overrides makes the caller responsible for the combined step sizes' stability.

ParametersCopy Anchor Link

iterations : int, optional Maximum number of solver updates. For SIRT and MLEM, an update is a complete subset sweep. Default is 400. tv_weight : float, optional Nonnegative regularization weight. Default is 0.005. regulariser : Regulariser, optional Regularization model: "tv" or "huber_tv". Default is 'tv'. huber_delta : float, optional Positive Huber transition scale when regulariser is "huber_tv". Default is 0.01. theta : float, optional Extrapolation factor for the dual gradient. Default is 1.0. views_per_batch : int, optional Number of views processed together to bound temporary storage. Default is 16. seed : int, optional Random seed for stochastic block selection. Default is 0. tau : float | None, optional Positive primal step size, or None to estimate it. sigma_data : float | None, optional Positive data-dual step size, or None to estimate it. sigma_tv : float | None, optional Positive TV-dual step size, or None to estimate it. projector_unroll : int, optional Loop unrolling for the explicit-grid or ray-reference projector. Default is 1. checkpoint_projector : bool, optional Recompute intermediate projector values during differentiation on the reference path. Default is True. gather_dtype : str, optional Gather precision on the reference path: "auto", "fp32", "fp16" or "bf16". Default is 'fp32'. ray_integrator : str, optional Reference ray integration rule; "sampled" is the default, "exact" uses exact integration, and "joseph" selects Joseph sampling. Default is 'sampled'. projector_model : ProjectorModel, optional Projection model: "auto", "ray" or "joseph". Auto selects Joseph unless an explicit detector grid requires rays. Default is 'auto'. backend : ProjectorBackend, optional Execution backend: "auto", "jax" or "cuda". An explicit accelerator request never silently falls back. Default is 'auto'. nonnegative : bool, optional Constrain the reconstruction to nonnegative values. Default is False. support : array-like | None, optional Volume support mask. Excluded voxels remain zero; None uses the entire grid. log_every : int, optional Number of iterations between logged objective samples. Default is 10. stream_projections : bool | None, optional Read host projections one view batch at a time. None selects streaming from available device memory; False keeps them resident.

calibrate_cone_axisCopy Anchor Link

calibrate_cone_axis(
    geometry: ConeGeometry,
    grid: Grid,
    detector: Detector,
    projections: np.ndarray | jax.Array,
    *,
    config: ConeAxisConfig | None = None,
) -> ConeAxisCalibration

Estimate a cone-beam scan's axis offset and detector roll from its projections.

projections are (views, nv, nu) line integrals (NumPy, memmap or JAX). The scan should cover at least 180 degrees plus the fan angle; the sample should have some in-plane structure at the slab heights. Use ConeAxisCalibration.apply to update the geometry.

ParametersCopy Anchor Link

geometry : Geometry Scan geometry, including view angles and any per-view poses. grid : Grid Reconstruction grid in (x, y, z) order, with physical voxel spacing. detector : Detector Detector dimensions and placement in the same length unit as the grid. projections : array-like Measured line integrals with shape (views, detector.nv, detector.nu). config : ConeAxisConfig or None, optional Solver settings. None uses ConeAxisConfig defaults.

ReturnsCopy Anchor Link

ConeAxisCalibration Absolute axis offset and detector roll, with per-slab estimates.

RaisesCopy Anchor Link

ValueError If shapes, geometry or solver settings are incompatible with this method.

cglsCopy Anchor Link

Conjugate-gradient least squares using the exact matched projector pair.

cgls_multiresCopy Anchor Link

Coarse-to-fine CGLS initialization followed by a full-data solve.

fbpCopy Anchor Link

Filtered-backprojection reconstruction routines.

fbp_hostCopy Anchor Link

FBP with host storage and bounded volume slabs on the accelerator.

fdkCopy Anchor Link

Feldkamp-Davis-Kress (FDK) reconstruction for cone-beam scans.

Each projection is cosine-weighted, ramp-filtered along detector rows on the isocentre plane, and backprojected voxel by voxel with the (SOD / U)^2 distance weight, where U is the voxel's distance from the source along the detector normal. Full turns share angular measure across conjugate views; shorter arcs of at least a half turn plus the fan angle use Parker weights with a virtual fan spanning the acquired arc, retaining overscan measurements.

FDK is exact only in the plane through the source orbit; away from it the reconstruction carries the usual cone-beam artefacts, growing with the cone angle. It is a fast first reconstruction and initialiser for the iterative solvers, which model the cone geometry exactly.

fdk_hostCopy Anchor Link

fdk_host(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: np.ndarray,
    *,
    config: FDKHostConfig | None = None,
    out: np.ndarray | None = None,
) -> np.ndarray

FDK with host input and output, reconstructed in z slabs on the device.

Accepts NumPy arrays and memmaps; out may be a writable (nx, ny, nz) FP32 memmap. Each slab filters a bounded band containing the rows it projects onto in any calibrated view, so projections and volume can both exceed device memory. With several config.fdk.devices, each device reconstructs its own slabs, so the volume may exceed every one of them.

Input and output storage must not overlap. Input storage must remain unchanged until the call finishes; sampled projections must be finite in FP32. Completed slabs are written immediately, so a later failure can leave a partially written output. This host-returning routine is not differentiable; use fdk for device-resident output.

ParametersCopy Anchor Link

geometry : Geometry Scan geometry, including view angles and any per-view poses. grid : Grid Reconstruction grid in (x, y, z) order, with physical voxel spacing. detector : Detector Detector dimensions and placement in the same length unit as the grid. projections : array-like Measured line integrals with shape (views, detector.nv, detector.nu). config : FDKHostConfig or None, optional Solver settings. None uses FDKHostConfig defaults. out : numpy.ndarray or None, optional Writable float32 host storage of shape (nx, ny, nz), including a memmap. Must not overlap input storage. None allocates a new Fortran-order output.

ReturnsCopy Anchor Link

numpy.ndarray Float32 reconstruction of shape (nx, ny, nz). Returns out itself when supplied.

RaisesCopy Anchor Link

ValueError If shapes, geometry or solver settings are incompatible with this method.

fista_tvCopy Anchor Link

Public FISTA/TV reconstruction adapter.

fourier_reconstructCopy Anchor Link

fourier_reconstruct(
    geometry: ScanGeometry,
    grid: Grid,
    detector: Detector,
    projections: np.ndarray,
    *,
    config: FourierConfig | None = None,
    out: np.ndarray | None = None,
) -> np.ndarray

Reconstruct a half-turn parallel scan using the Fourier slice theorem.

Projections have shape (views, nv, nu) in physical line-integral units. Angles may be reordered, reversed, offset or replaced by opposing views, with uniform or nonuniform spacing. Redundant half-turn directions are averaged in physical Fourier coordinates after detector phase/orientation. Per-view translations are honored in either detector or object coordinates. Rotational pose corrections and tilted axes are not supported. Detector offsets, anisotropic voxels, shifted/cropped grids and fractional detector-v rows retain their physical coordinates.

This approximate inverse uses six-point Kaiser--Bessel radial interpolation, linear interpolation in measured angle, a detector-Nyquist disk cutoff and a padded Cartesian inverse FFT. DC uses periodic angular trapezoid weights. Missing detector data are zero. Truncation or angular undersampling can cause artifacts; no amplitude fitting, clipping or truncation correction is applied.

Input and output stay in NumPy arrays or memmaps, with float32 output in (nx, ny, nz) order. New outputs are C-contiguous for downstream volume operations; supplied output storage is preserved. Their storage must not overlap. Completed slabs are written immediately, so a later failure may leave partial output. Input storage must remain unchanged during the call. This routine is not differentiable; use the matched projector pair for iterative solves.

ParametersCopy Anchor Link

geometry : Geometry Scan geometry, including view angles and any per-view poses. grid : Grid Reconstruction grid in (x, y, z) order, with physical voxel spacing. detector : Detector Detector dimensions and placement in the same length unit as the grid. projections : array-like Measured line integrals with shape (views, detector.nv, detector.nu). config : FourierConfig or None, optional Solver settings. None uses FourierConfig defaults. out : numpy.ndarray or None, optional Writable float32 host storage of shape (nx, ny, nz), including a memmap. Must not overlap input storage. None allocates a new output.

ReturnsCopy Anchor Link

numpy.ndarray Float32 reconstruction of shape (nx, ny, nz). Returns out itself when supplied.

RaisesCopy Anchor Link

ValueError If shapes, geometry or solver settings are incompatible with this method. ImportError If the requested CuPy backend is not installed.

ExamplesCopy Anchor Link

import numpy as np import tomojax as tj from tomojax.recon import fourier_reconstruct, FourierConfig, project grid = tj.Grid(4, 4, 1, 1.0, 1.0, 1.0) detector = tj.Detector(4, 1, 1.0, 1.0) geometry = tj.ParallelGeometry(grid, detector, np.linspace(0, 180, 8, endpoint=False)) data = project(geometry, np.ones((4, 4, 1), np.float32)) volume = fourier_reconstruct(geometry, geometry.grid, geometry.detector, np.asarray(data)) volume.shape (4, 4, 1)

mlemCopy Anchor Link

Poisson maximum-likelihood reconstruction with matched nonnegative projectors.

sirtCopy Anchor Link

Row/column-normalized simultaneous and ordered-subset algebraic reconstruction.

spdhg_tvCopy Anchor Link

SPDHG/TV reconstruction routine.