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.

Parameters

• 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.

Returns

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

Raises

• 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.

Parameters

• 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.

Returns

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

Raises

• 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.

Parameters

• 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.

Parameters

• 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.

Parameters

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

Returns

• 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.

Parameters

• 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.

Parameters

• 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.

Parameters

• 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.

Parameters

• 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.

Parameters

• 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).

Parameters

• 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.

Parameters

• 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.

Parameters

• 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.

Parameters

• 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 64.
• 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.

Parameters

• 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.

Parameters

• 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.

Returns

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

Raises

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

cglsCopy Anchor Link

cgls(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jnp.ndarray | np.ndarray,
    *,
    init_x: jnp.ndarray | np.ndarray | None = None,
    config: CGLSConfig | None = None,
    det_grid: tuple[jnp.ndarray, jnp.ndarray] | None = None,
    out: np.ndarray | None = None,
) -> tuple[jnp.ndarray | np.ndarray, dict[str, object]]

Reconstruct with matched CGLS and report convergence diagnostics.

Inputs and the solve use FP32. Host arrays on CUDA can use the native linear Joseph kernels and return a NumPy volume, avoiding JAX startup and tracing; device arrays, streaming, and other models use the compiled JAX solver. Both paths use the same stationarity and roundoff checks. This public solver returns host diagnostics and is not differentiable.

Cone-beam volumes whose working set exceeds the device, or given out, are solved with their volumes in host memory, streamed through the device in z slabs; the result is then a Fortran-ordered NumPy volume (out itself if given). With a memory-mapped out, the three work volumes are temporary files in its directory.

Parameters

• 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).
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from zero.
• config (CGLSConfig or None, optional): Solver settings. None uses CGLSConfig defaults.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.
• out (numpy.ndarray or None, optional): Writable float32 (nx, ny, nz) host storage for a cone-beam result, including a memmap; it must not overlap the projections.

Returns

• volume (jax.Array or numpy.ndarray): Float32 reconstruction of shape (nx, ny, nz).
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import cgls, CGLSConfig, 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, info = cgls(
    geometry,
    geometry.grid,
    geometry.detector,
    data,
    config=CGLSConfig(iterations=1),
)
volume.shape
# (4, 4, 1)

cgls_multiresCopy Anchor Link

cgls_multires(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jnp.ndarray | np.ndarray,
    *,
    factors: Iterable[int] = (2, 1),
    iterations_per_level: Iterable[int] = (32, 16),
    init_x: jnp.ndarray | None = None,
    config: CGLSConfig | None = None,
    det_grid: tuple[jnp.ndarray, jnp.ndarray] | None = None,
) -> tuple[jnp.ndarray | np.ndarray, dict[str, object]]

Initialize CGLS on coarser grids, then refine against every measured ray.

factors must decrease strictly to 1. Each level gets its explicit positive iteration budget from iterations_per_level; config.iterations is overridden. Other CGLS settings apply unchanged at each level, including damping and streaming. Host arrays remain on the host until each level chooses its placement. Supported host-array CUDA solves and their interpolation avoid JAX startup. Tolerances may stop levels early.

Volume faces, world poses and physical density are preserved. Coarse data select a uniformly strided subset of measured rays without averaging or amplitude scaling. Explicit detector coordinates use the identical subset. A supplied initial volume is at full resolution and is resampled onto the first level. The final solve uses the original grid, detector and all data.

Coarsening is an optional initialization strategy, not a quality guarantee; fine detail and noisy data may require more fine-level iterations. Diagnostics contain every level's actual work and the final level's stopping criterion. Numerical breakdown at any level raises FloatingPointError.

Parameters

• 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).
• factors (iterable of int, optional): Coarsening factors decreasing strictly to 1. Default is (2, 1).
• iterations_per_level (iterable of int, optional): Positive iteration budget at each level. Default is (32, 16); overrides config.iterations.
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from zero.
• config (CGLSConfig or None, optional): Solver settings. None uses CGLSConfig defaults.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.

Returns

• volume (jax.Array or numpy.ndarray): Float32 reconstruction of shape (nx, ny, nz).
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

• ValueError: If shapes, geometry or solver settings are incompatible with this method.
• FloatingPointError: If a level's CGLS solve breaks down numerically.

fbpCopy Anchor Link

fbp(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jnp.ndarray | np.ndarray,
    *,
    config: FBPConfig | None = None,
    det_grid: tuple[jnp.ndarray, jnp.ndarray] | None = None,
) -> jnp.ndarray

Filtered backprojection for parallel rays rotating about a fixed axis.

Cone-beam geometries are reconstructed with tomojax.recon.fdk.fdk using the same filter; backend="cuda" selects its CUDA kernel.

Projections: (n_views, nv, nu) -> attenuation volume (nx, ny, nz). The rotation axis, arc and angular spacing are fitted from the view poses, so tilted (laminography) axes, partial or full turns and irregular angles are weighted exactly on every measured frequency. Repeated directions share their angular weight equally. Frequencies no view measures, such as laminography's missing cone, reconstruct as zero.

Backprojection is voxel-driven with bilinear detector interpolation and retains filtered tails through the full output volume, assuming zero raw attenuation beyond the measured detector. This is not a correction for truncated objects. An explicit det_grid instead uses the ray-model adjoint with uniform weights.

Parameters

• 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 (FBPConfig or None, optional): Solver settings. None uses FBPConfig defaults.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.

Returns

• jax.Array: Float32 reconstruction of shape (nx, ny, nz).

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import fbp, FBPConfig, 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 = fbp(geometry, geometry.grid, geometry.detector, np.asarray(data))
volume.shape
# (4, 4, 1)

fbp_hostCopy Anchor Link

fbp_host(
    geometry: ParallelGeometry,
    grid: Grid,
    detector: Detector,
    projections: np.ndarray | jax.Array,
    *,
    config: FBPHostConfig | None = None,
    out: np.ndarray | None = None,
) -> np.ndarray

Reconstruct without storing the full volume or projections on the device.

Accept NumPy arrays, including memmaps, or device arrays and return a float32 NumPy array. Automatically allocated axial output uses Fortran order for contiguous xy-planes; general geometry uses C order for contiguous x-slabs. out may be a writable memmap. Writes complete before returning, but a later failure can leave a partially written output. Input storage must remain unchanged until the call finishes. This host-returning routine is not differentiable. Use fbp for device-resident output.

Every geometry fbp accepts is supported, with the same exact weights: laminography and posed scans reconstruct in x-slabs, built-in parallel scans in axial slabs that read only the detector rows they need. Anisotropic spacing, detector offsets and shifted volume origins are supported. Filter tails use the same zero-extended measurement model as fbp.

Parameters

• 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 (FBPHostConfig or None, optional): Solver settings. None uses FBPHostConfig 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.

Returns

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

Raises

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

fdkCopy Anchor Link

fdk(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jax.Array | np.ndarray,
    *,
    config: FDKConfig | None = None,
) -> jax.Array

Reconstruct a cone-beam scan with FDK; returns an (nx, ny, nz) volume.

projections are (views, nv, nu) line integrals; NumPy and memmap stacks stream to the device views_per_batch views at a time. Full turns on an offset detector (the axis projecting off its centre) get Wang's weights from the conjugate views and reconstruct the field the detector's long side covers. Each view retains its own calibrated source/detector arrangement; angular quadrature covers the whole scan, independently of segment boundaries. Noncircular trajectories retain FDK's approximate character.

Parameters

• 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 (FDKConfig or None, optional): Solver settings. None uses FDKConfig defaults.

Returns

• jax.Array: Float32 reconstruction of shape (nx, ny, nz).

Raises

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

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.

Parameters

• 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.

Returns

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

Raises

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

fista_tvCopy Anchor Link

fista_tv(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jnp.ndarray | np.ndarray,
    *,
    init_x: jnp.ndarray | None = None,
    config: FistaConfig | None = None,
    callback: LossCallback | None = None,
    det_grid: tuple[jnp.ndarray, jnp.ndarray] | None = None,
) -> tuple[jnp.ndarray, dict[str, object]]

Run FISTA with TV regularization using an explicit solver configuration.

If callback is provided, it fires on the first recorded loss sample and on the final recorded loss sample. The callback arguments are (step, loss), where step is the zero-based iteration index that produced loss. When early stopping truncates active iterations, the final callback reports the last active iteration rather than the repeated padded tail entry.

Parameters

• 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).
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from zero.
• config (FistaConfig or None, optional): Solver settings. None uses FistaConfig defaults.
• callback (callable or None, optional): Callback receiving (step, loss) at the first and last logged samples.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.

Returns

• volume (jax.Array): Float32 reconstruction of shape (nx, ny, nz).
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import fista_tv, FistaConfig, 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, info = fista_tv(
    geometry,
    geometry.grid,
    geometry.detector,
    data,
    config=FistaConfig(iterations=1),
)
volume.shape
# (4, 4, 1)

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.

Parameters

• 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.

Returns

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

Raises

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

Examples

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

mlem(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jax.Array | np.ndarray,
    *,
    init_x: jax.Array | np.ndarray | None = None,
    config: MLEMConfig | None = None,
    det_grid: tuple[jax.Array, jax.Array] | None = None,
    out: np.ndarray | None = None,
) -> tuple[jax.Array | np.ndarray, dict[str, object]]

MLEM or OSEM with the exact transpose; see MLEMConfig.

Returns the FP32 volume and an iteration record, without a convergence claim. Wholly unobserved voxels are zero; zero measurements drive observed voxels to zero. Large host measurement stacks stream in bounded view batches. Cone solves whose workspace exceeds the device, or given out, update the volume a z slab at a time, as for cone SIRT; every subset projection sums all slabs before its coupled update. The resident-volume path remains the default whenever its workspace fits.

Parameters

• 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): Nonnegative linear Poisson measurements (not raw transmission counts), with shape (views, detector.nv, detector.nu).
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from one in observed support.
• config (MLEMConfig or None, optional): Solver settings. None uses MLEMConfig defaults.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.
• out (numpy.ndarray or None, optional): Writable float32 (nx, ny, nz) host storage for a cone-beam result, including a memmap; it must not overlap the projections.

Returns

• volume (jax.Array or numpy.ndarray): Float32 reconstruction of shape (nx, ny, nz), host-owned when the cone solver workspace does not fit the device.
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import mlem, MLEMConfig, 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, info = mlem(
    geometry,
    geometry.grid,
    geometry.detector,
    data,
    config=MLEMConfig(iterations=1),
)
volume.shape
# (4, 4, 1)

sirtCopy Anchor Link

sirt(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jax.Array | np.ndarray,
    *,
    init_x: jax.Array | np.ndarray | None = None,
    config: SIRTConfig | None = None,
    det_grid: tuple[jax.Array, jax.Array] | None = None,
    out: np.ndarray | None = None,
) -> tuple[jax.Array | np.ndarray, dict[str, object]]

Matched FP32 SIRT/OS-SART; each iteration visits every subset once.

Supports parallel, tilted and cone-beam scans on the existing matched projectors. Inputs and updates are FP32. An iteration is a complete sweep, not one subset update. No residual-based stopping or quality claim is made. Large independent axial planes use host slabs and invocation-owned host normalizers. Large host projection stacks stream with exact row recomputation, also allowing reusable column weights to fit in device memory. Cone-beam solves whose working set exceeds the device, or given out, update the volume a z slab at a time; it stays on the device when it fits beside that sweep, and otherwise in host memory as a Fortran-ordered NumPy volume (out if given).

Parameters

• 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).
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from zero.
• config (SIRTConfig or None, optional): Solver settings. None uses SIRTConfig defaults.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.
• out (numpy.ndarray or None, optional): Writable float32 (nx, ny, nz) host storage for a cone-beam result, including a memmap; it must not overlap the projections.

Returns

• volume (jax.Array or numpy.ndarray): Float32 reconstruction of shape (nx, ny, nz).
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import sirt, SIRTConfig, 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, info = sirt(
    geometry,
    geometry.grid,
    geometry.detector,
    data,
    config=SIRTConfig(iterations=1),
)
volume.shape
# (4, 4, 1)

spdhg_tvCopy Anchor Link

spdhg_tv(
    geometry: Geometry,
    grid: Grid,
    detector: Detector,
    projections: jnp.ndarray | np.ndarray,
    *,
    weights: jnp.ndarray | None = None,
    init_x: jnp.ndarray | None = None,
    config: SPDHGConfig | None = None,
    callback: LossCallback | None = None,
    det_grid: tuple[jnp.ndarray, jnp.ndarray] | None = None,
) -> tuple[jnp.ndarray, dict[str, object]]

Reconstruct with stochastic primal-dual weighted L2 and TV updates.

weights must be real, nonnegative and finite when stored in FP32; zero weights exclude unmeasured samples, including nonfinite measurements. Invalid weights fail before solver setup. Boolean/integer masks, standard NumPy floating dtypes and bfloat16 are supported.

The supplied init_x obeys the configured positivity and support map before the first projection, including when iterations=0.

Logged losses estimate the full objective at the primal image used for that iteration's data projection, with inverse block-selection weights. They are samples, not full-scan objectives of the returned primal image.

If callback is provided, it fires on the first logged objective sample and on the final logged objective sample. The callback arguments are (step, loss)`, where `step is the zero-based iteration index that produced loss. Only iterations whose objective was recorded under config.log_every are eligible for callbacks; if a single logged sample exists, the callback fires once.

Parameters

• 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).
• weights (array-like or None, optional): Real, nonnegative weights matching the projection shape; zero excludes unmeasured samples. None weights all samples equally.
• init_x (array-like or None, optional): Initial volume of shape (nx, ny, nz); None starts from zero.
• config (SPDHGConfig or None, optional): Solver settings. None uses SPDHGConfig defaults.
• callback (callable or None, optional): Callback receiving (step, loss) at the first and last logged samples.
• det_grid (tuple of array-like or None, optional): Explicit detector-coordinate arrays (u, v), each of shape (nv, nu). None uses the canonical detector grid.

Returns

• volume (jax.Array): Float32 reconstruction of shape (nx, ny, nz).
• info (dict): Solver diagnostics and effective work; a small residual or iteration count does not establish reconstruction accuracy.

Raises

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

Examples

import numpy as np
import tomojax as tj
from tomojax.recon import spdhg_tv, SPDHGConfig, 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, info = spdhg_tv(
    geometry,
    geometry.grid,
    geometry.detector,
    data,
    config=SPDHGConfig(iterations=1),
)
volume.shape
# (4, 4, 1)