Tomography

class deepinv.physics.Tomography(angles, img_width, circle=False, parallel_computation=True, normalize=False, fan_beam=False, fan_parameters=None, device=device(type='cpu'), dtype=torch.float32, **kwargs)[source]

Bases: LinearPhysics

(Computed) Tomography operator.

The Radon transform is the integral transform which takes a square image \(x\) defined on the plane to a function \(y=Rx\) defined on the (two-dimensional) space of lines in the plane, whose value at a particular line is equal to the line integral of the function over that line.

Note

The pseudo-inverse is computed using the filtered back-projection algorithm with a Ramp filter. This is not the exact linear pseudo-inverse of the Radon transform, but it is a good approximation which is robust to noise.

Note

The measurements are not normalized by the image size, thus the norm of the operator depends on the image size.

Warning

The adjoint operator has small numerical errors due to interpolation.

Parameters:

angles (int, torch.tensor) – These are the tomography angles. If the type is int, the angles are sampled uniformly between 0 and 360 degrees. If the type is torch.tensor, the angles are the ones provided (e.g., torch.linspace(0, 180, steps=10)).
img_width (int) – width/height of the square image input.
circle (bool) – If True both forward and backward projection will be restricted to pixels inside a circle inscribed in the square image.
parallel_computation (bool) – if True, all projections are performed in parallel. Requires more memory but is faster on GPUs.
normalize (bool) – If True, the outputs are normlized by the image size (i.e. it is assumed that the image lives on [0,1]^2 for the computation of the line integrals). In this case the operator norm is approximately given by \(\|A\|_2^2 \approx \frac{\pi}{2\,\text{angles}}\), If False, then it is assumed that the image lives on [0,im_width]^2 for the computation of the line integrals
fan_beam (bool) – If True, use fan beam geometry, if False use parallel beam
fan_parameters (dict) –
Only used if fan_beam is True. Contains the parameters defining the scanning geometry. The dict should contain the keys:
- ”pixel_spacing” defining the distance between two pixels in the image, default: 0.5 / (in_size)
- ”source_radius” distance between the x-ray source and the rotation axis (middle of the image), default: 57.5
- ”detector_radius” distance between the x-ray detector and the rotation axis (middle of the image), default: 57.5
- ”n_detector_pixels” number of pixels of the detector, default: 258
- ”detector_spacing” distance between two pixels on the detector, default: 0.077
The default values are adapted from the geometry in https://doi.org/10.5281/zenodo.8307932, where pixel spacing, source and detector radius and detector spacing are given in cm. Note that a to small value of n_detector_pixels*detector_spacing can lead to severe circular artifacts in any reconstruction.
device (str) – gpu or cpu.

Examples:

Tomography operator with defined angles for 3x3 image:

>>> from deepinv.physics import Tomography
>>> seed = torch.manual_seed(0)  # Random seed for reproducibility
>>> x = torch.randn(1, 1, 4, 4)  # Define random 4x4 image
>>> angles = torch.linspace(0, 45, steps=3)
>>> physics = Tomography(angles=angles, img_width=4, circle=True)
>>> physics(x)
tensor([[[[ 0.1650,  1.2640,  1.6995],
          [-0.4860,  0.2674,  0.9971],
          [ 0.9002, -0.3856, -0.9360],
          [-2.4882, -2.1068, -2.5720]]]])

Tomography operator with 3 uniformly sampled angles in [0, 360] for 3x3 image:

>>> from deepinv.physics import Tomography
>>> seed = torch.manual_seed(0)  # Random seed for reproducibility
>>> x = torch.randn(1, 1, 4, 4)  # Define random 4x4 image
>>> physics = Tomography(angles=3, img_width=4, circle=True)
>>> physics(x)
tensor([[[[ 0.1650,  1.9493,  1.9897],
          [-0.4860,  0.7137, -1.6536],
          [ 0.9002, -0.8457, -0.1666],
          [-2.4882, -2.7340, -0.9793]]]])

A(x, **kwargs)[source]

Computes forward operator \(y = A(x)\) (without noise and/or sensor non-linearities)

Parameters:: x (torch.Tensor,list[torch.Tensor]) – signal/image
Returns:: (torch.Tensor) clean measurements

A_adjoint(y, **kwargs)[source]

Computes transpose of the forward operator \(\tilde{x} = A^{\top}y\). If \(A\) is linear, it should be the exact transpose of the forward matrix.

Note

If the problem is non-linear, there is not a well-defined transpose operation, but defining one can be useful for some reconstruction networks, such as deepinv.models.ArtifactRemoval.

Parameters:

y (torch.Tensor) – measurements.
params (None, torch.Tensor) – optional additional parameters for the adjoint operator.

Returns:

(torch.Tensor) linear reconstruction \(\tilde{x} = A^{\top}y\).

A_dagger(y, **kwargs)[source]

Computes the solution in \(x\) to \(y = Ax\) using the conjugate gradient method, see deepinv.optim.utils.conjugate_gradient().

If the size of \(y\) is larger than \(x\) (overcomplete problem), it computes \((A^{\top} A)^{-1} A^{\top} y\), otherwise (incomplete problem) it computes \(A^{\top} (A A^{\top})^{-1} y\).

This function can be overwritten by a more efficient pseudoinverse in cases where closed form formulas exist.

Parameters:: y (torch.Tensor) – a measurement \(y\) to reconstruct via the pseudoinverse.
Returns:: (torch.Tensor) The reconstructed image \(x\).

Examples using `Tomography`:

A tour of forward sensing operators

Patch priors for limited-angle computed tomography

Vanilla PnP for computed tomography (CT).

Self-supervised learning with measurement splitting

Learned Primal-Dual algorithm for CT scan.

Tomography

Examples using Tomography:

Examples using `Tomography`: