CPU Tomography Module

Tomography tools for nDTomo

@author: Antony Vamvakeros

nDTomo.tomo.conv_tomo.backproject(sinogram, angles, output_size=None, normalize=True)[source]

Performs filtered or unfiltered backprojection of a sinogram.

Parameters:

  • sinogram (ndarray) – 2D array of shape (num_detectors, num_angles).

  • angles (ndarray) – 1D array of projection angles in degrees.

  • output_size (int, optional) – Size (width and height) of the reconstructed image. If None, use num_detectors from sinogram.

Returns:

reconstruction – 2D reconstructed image.

Return type:

ndarray

nDTomo.tomo.conv_tomo.cgls(A, b, x0=None, n_iter=10)[source]

Conjugate Gradient Least Squares (CGLS) algorithm for solving Ax = b.

This method minimizes the least-squares objective:

||Ax - b||^2

It is equivalent to applying the Conjugate Gradient method to the normal equations:

AᵀA x = Aᵀb

Parameters:

  • A (ndarray or sparse matrix, shape (m, n)) – The system matrix. Can be dense or sparse.

  • b (ndarray, shape (m,)) – The right-hand side vector (e.g., sinogram or projection data).

  • x0 (ndarray or None, shape (n,), optional) – Initial guess for the solution. If None, defaults to zeros.

  • n_iter (int, optional) – Number of iterations to perform. Default is 10.

Returns:

x – The reconstructed solution vector.

Return type:

ndarray, shape (n,)

Notes

  • CGLS converges faster than basic iterative methods like ART or SIRT, especially for well-conditioned problems.

  • It avoids explicitly forming AᵀA, which can be ill-conditioned or expensive to compute.

  • This implementation does not support stopping criteria based on tolerance.

nDTomo.tomo.conv_tomo.cgls_functional(sinogram, angles, x0=None, n_iter=10)[source]

CGLS solver using forwardproject and backproject, matching the vectorized logic.

Parameters:

  • sinogram (ndarray) – Input sinogram of shape (num_detectors, num_angles).

  • angles (ndarray) – 1D array of projection angles in degrees.

  • x0 (ndarray or None) – Initial image guess. If None, initialized to zeros.

  • n_iter (int) – Number of CGLS iterations.

Returns:

image – Reconstructed 2D image of shape (num_detectors, num_detectors).

Return type:

ndarray

nDTomo.tomo.conv_tomo.fbpvol(svol, scan=180, theta=None, nt=None)[source]

Reconstructs a stack of images or a single image from sinograms using the filtered backprojection (FBP) algorithm.

The function processes a 2D sinogram or a 3D stack of sinograms, where the third dimension represents either the z-axis or spectral dimension. Reconstruction is performed for each slice using the filtered backprojection method from skimage.transform.iradon.

Parameters:

  • svol (ndarray) – Input sinogram or sinogram volume. Can be either: - A 2D array (n_projections x detector_width) representing a single sinogram. - A 3D array (n_rows x n_projections x n_depths) representing a stack of sinograms.

  • scan (int, optional) – Total range of projection angles (in degrees) over which the sinograms were acquired. Default is 180 degrees.

  • theta (array-like, optional) – Array of projection angles (in degrees). If not provided, the angles are evenly spaced over the range defined by scan.

  • nt (int, optional) – Desired size of the reconstructed images (number of pixels along one dimension). If not provided, it defaults to the number of rows in svol.

Returns:

Reconstructed image or volume: - For a 2D input: A 2D array of shape (nt, nt). - For a 3D input: A 3D array of shape (nt, nt, n_depths).

Return type:

ndarray

Notes

  • The reconstruction is performed using skimage.transform.iradon.

  • The circle parameter in iradon is set to True, assuming the object fits entirely within the reconstruction area.

  • The output dimensions depend on the input: - For 2D sinograms: (nt, nt). - For 3D sinogram stacks: (nt, nt, n_depths).

nDTomo.tomo.conv_tomo.filter_sinogram(sinogram, filter_type='Ramp')[source]

Applies a frequency domain filter to a sinogram.

Parameters:

  • sinogram (ndarray) – 2D array of shape (num_detectors, num_angles).

  • filter_type (str) – Type of filter to apply. Options: - ‘Ramp’ : standard Ram-Lak filter - ‘Shepp-Logan’ : Ramp * sinc

Returns:

filtered_sinogram – The filtered sinogram with the same shape as input.

Return type:

ndarray

nDTomo.tomo.conv_tomo.forwardproject(image, angles, num_detectors=None)[source]

Simulates a forward projection (Radon transform) of a 2D image.

Parameters:

  • image (ndarray) – 2D input image (square or rectangular).

  • angles (ndarray) – 1D array of projection angles in degrees.

  • num_detectors (int, optional) – Number of detectors in the output sinogram. If None, uses image height.

Returns:

sinogram – 2D sinogram of shape (num_detectors, num_angles).

Return type:

ndarray

nDTomo.tomo.conv_tomo.paralleltomo(N, theta, p, w)[source]

Creates a 2D tomography test problem using parallel beams.

This function generates a sparse matrix A representing the forward projection operator for a 2D domain discretized into N x N cells. The problem is defined using p parallel rays for each angle in theta.

Parameters:

  • N (int) – Number of discretization intervals in each dimension, resulting in an N x N domain.

  • theta (array-like) – Array of angles in degrees.

  • p (int, optional) – Number of parallel rays for each angle. Default is round(sqrt(2) * N).

  • w (float, optional) – Width of the ray fan (distance from the first ray to the last). Default is sqrt(2) * N.

Returns:

Sparse matrix A with shape (nA * p, N^2), where nA is the number of angles in theta. The rows represent rays, and the columns correspond to cells in the domain.

Return type:

csr_matrix

Notes

  • The generated matrix A can be used in test problems for tomographic reconstruction.

  • This implementation assumes the rays pass through a square domain centered at the origin.

References

      1. Kak and M. Slaney, “Principles of Computerized Tomographic Imaging,” SIAM, Philadelphia, 2001.

  • Original MATLAB code: AIR Tools, DTU Informatics, June 21, 2011.

nDTomo.tomo.conv_tomo.radonvol(vol, scan=180, theta=None)[source]

Computes the Radon transform of a stack of 2D images or a single 2D image.

The Radon transform is calculated along specified projection angles for each 2D slice in the input volume. This function supports both 2D and 3D input data, where the third dimension in the 3D case corresponds to either the z-axis or spectral dimension.

Parameters:

  • vol (ndarray) – Input volume. Can be either: - A 2D array (H x W) representing a single image. - A 3D array (H x W x D) representing a stack of D images.

  • scan (int, optional) – Total range of angles (in degrees) over which projections are computed. Default is 180 degrees.

  • theta (array-like, optional) – Array of projection angles (in degrees). If not provided, the angles are evenly spaced over the range defined by scan.

Returns:

Sinogram or sinogram volume: - For a 2D input: A 2D array of shape (len(theta), W). - For a 3D input: A 3D array of shape (H, len(theta), D).

Return type:

ndarray

Notes

  • The Radon transform is calculated using skimage.transform.radon.

  • The output dimensions depend on the input shape: - For 2D input: (len(theta), width of the input image) - For 3D input: (height of the input image stack, len(theta), depth of the stack)

nDTomo.tomo.conv_tomo.sirt(A, b, x0=None, n_iter=20, relax=1.0)[source]

Simultaneous Iterative Reconstruction Technique (SIRT) with row and column normalization.

This function implements a general version of SIRT using matrix-based normalization inspired by Cimmino/CAV-style scaling. It solves the linear system A @ x ≈ b iteratively.

Each iteration computes the global residual and applies a normalized correction:

x ← x + relax * D * A.T * M * (b - A @ x)

Parameters:

  • A (scipy.sparse matrix (CSR or CSC), shape (m, n)) – The system matrix, typically sparse and non-negative. Rows represent projection rays; columns represent image pixels.

  • b (ndarray, shape (m,)) – The measured projection data (sinogram flattened to 1D).

  • x0 (ndarray or None, shape (n,), optional) – Initial guess for the solution. Defaults to zeros if None.

  • n_iter (int, optional) – Number of SIRT iterations to perform. More iterations yield smoother convergence.

  • relax (float, optional) – Relaxation parameter controlling update magnitude. Values between 0.5–1.0 are typical.

Returns:

x – The reconstructed solution vector (flattened image).

Return type:

ndarray, shape (n,)

Notes

  • The row normalization matrix M scales residual contributions by the inverse of row energy:

    M[i, i] = 1 / ||A[i, :]||²

  • The column normalization matrix D scales the updates to balance pixel sensitivity:

    D[j, j] = 1 / ||A[:, j]||²

  • This global scheme is slower than CGLS but more stable for noisy or incomplete data.

  • Particularly useful in tomographic reconstruction where balancing ray coverage is critical.

nDTomo.tomo.conv_tomo.sirt_functional(sinogram, angles, x0=None, n_iter=20, relax=1.0, epsilon=1e-06)[source]

SIRT solver using forwardproject and backproject, with image size inferred from sinogram.

Parameters:

  • sinogram (ndarray) – Input sinogram of shape (num_detectors, num_angles).

  • angles (ndarray) – 1D array of projection angles in degrees.

  • x0 (ndarray or None) – Initial guess for the image. If None, initialized to zeros.

  • n_iter (int) – Number of SIRT iterations.

  • relax (float) – Relaxation factor.

  • epsilon (float) – Small constant to avoid division by zero.

Returns:

image – Reconstructed 2D image.

Return type:

ndarray