Code reference

Parallel beam geometry

radon_parallel(arr, angles[, count, max_count]) Compute the Radon transform (sinogram) of a circular image.
backproject(sinogram, angles[, filtering, …]) 2D backprojection algorithm
fourier_map(sinogram, angles[, intp_method, …]) 2D Fourier mapping with the Fourier slice theorem
integrate(sinogram, angles[, count, max_count]) 2D sum-reconstruction with the Fourier slice theorem
art(sinogram, angles[, initial, iterations, …]) Algebraic Reconstruction Technique
sart(sinogram, angles[, initial, …]) Simultaneous Algebraic Reconstruction Technique

Radon transform

radontea.radon_parallel(arr, angles, count=None, max_count=None)

Compute the Radon transform (sinogram) of a circular image.

The scipy Radon transform performs this operation on the entire image, whereas this implementation requires an input image that has gray-scale values of 0 outside of a circle with diameter equal to the image size.

Parameters:
  • arr (ndarray, shape (N,N)) – the input image.
  • angles (ndarray, length A) – angles or projections in radians
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
Returns:

outarr – Sinogram of the input image. The i’th row contains the projection data of the i’th angle.

Return type:

ndarray of floats, shape (A, N)

See also

scipy.ndimage.interpolation.rotate()
The interpolator used to rotate the image.

Non-iterative reconstruction

Computes the inverse Radon transform with non-iterative techniques. The linear system of equations that describes the forward process can be inverted with several algorithms, most notably the backprojection algorithm radontea.backproject(). The reconstruction is based on the Fourier slice theorem. A Fourier-based interpolation algorithm is implemented in radontea.fourier_map().

Backprojection

radontea.backproject(sinogram, angles, filtering='ramp', weight_angles=True, padding=True, padval=0, count=None, max_count=None, verbose=0)

2D backprojection algorithm

Computes the inverse of the Radon transform using filtered backprojection.

Parameters:
  • sinogram (ndarray, shape (A,N)) – Two-dimensional sinogram of line recordings.
  • angles ((A,) ndarray) – Angular positions of the sinogram in radians.
  • filtering ({'ramp', 'shepp-logan', 'cosine', 'hamming', ) –

    ‘hann’}, optional Specifies the Fourier filter. Either of

    • ”ramp”: mathematically correct reconstruction
    • ”shepp-logan”
    • ”cosine”
    • ”hamming”
    • ”hann”
  • weight_angles (bool) –

    If True, weights each backpropagated projection with a factor proportional to the angular distance between the neighboring projections.

    \[\Delta \phi_0 \longmapsto \Delta \phi_j = \frac{\phi_{j+1} - \phi_{j-1}}{2}\]

    New in version 0.1.9.

  • padding (bool, optional) – Pad the input data to the second next power of 2 before Fourier transforming. This reduces artifacts and speeds up the process for input image sizes that are not powers of 2.
  • padval (float) – The value used for padding. If padval is None, then the edge values are used for padding (see documentation of numpy.pad).
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
  • verbose (int) – Increment to increase verbosity.
Returns:

out – The reconstructed image.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from radontea import backproject
>>> N = 5
>>> A = 7
>>> x = np.linspace(-N/2, N/2, N)
>>> projection = np.exp(-x**2)
>>> sinogram = np.tile(projection, A).reshape(A,N)
>>> angles = np.linspace(0, np.pi, A, endpoint=False)
>>> backproject(sinogram, angles)
array([[ 0.03813283, -0.01972347, -0.02221885,  0.03822303,  0.01903376],
       [ 0.04033526, -0.01591647,  0.02262173,  0.04203002,  0.02123619],
       [ 0.02658797,  0.11375576,  0.41525594,  0.17170226,  0.0074889 ],
       [ 0.04008672,  0.11139209,  0.27650193,  0.16933859,  0.02098764],
       [ 0.02140445,  0.0334597 ,  0.07691547,  0.0914062 ,  0.00230537]])

Fourier mapping

radontea.fourier_map(sinogram, angles, intp_method='cubic', count=None, max_count=None)

2D Fourier mapping with the Fourier slice theorem

Computes the inverse of the Radon transform using Fourier interpolation. Warning: This is the naive reconstruction that assumes that the image is rotated through the upper left pixel nearest to the actual center of the image. We do not have this problem for odd images, only for even images.

Parameters:
  • sinogram ((A,N) ndarray) – Two-dimensional sinogram of line recordings.
  • angles ((A,) ndarray) – Angular positions of the sinogram in radians equally distributed from zero to PI.
  • intp_method ({'cubic', 'nearest', 'linear'}, optional) –

    Method of interpolation. For more information see scipy.interpolate.griddata. One of

    • ”nearest”: instead of interpolating, use the points closest to the input data
    • ”linear”: bilinear interpolation between data points
    • ”cubic”: interpolate using a two-dimensional poolynimial surface
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
Returns:

out – The reconstructed image.

Return type:

ndarray

See also

scipy.interpolate.griddata()
interpolation method used

Slow integration

radontea.integrate(sinogram, angles, count=None, max_count=None)

2D sum-reconstruction with the Fourier slice theorem

Computes the inverse of the Radon transform by computing the integral in real space.

Parameters:
  • sinogram ((A,N) ndarray) – Two-dimensional sinogram of line recordings.
  • angles ((A,) ndarray) – Angular positions of the sinogram in radians equally distributed from zero to PI.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
Returns:

out – The reconstructed image.

Return type:

ndarray

Iterative reconstruction

Inversion of Radon-based tomography methods using iterative algorithms. The convegence of these algorithms might be slow. The implementation is not optimized.

ART

radontea.art(sinogram, angles, initial=None, iterations=1, count=None, max_count=None)

Algebraic Reconstruction Technique

The Algebraic Reconstruction Technique (ART) iteratively computes the inverse of the Radon transform in two dimensions. The reconstruction technique uses rays of the diameter of one pixel to iteratively solve the system of linear equations that describe the projection process. The binary weighting factors are

  • 1, if the center of the a pixel is within the ray
  • 0, else
Parameters:
  • sinogram (ndarrayy, shape (A,N)) – Two-dimensional sinogram of line recordings.
  • angles (ndarray, length A) – Angular positions of the sinogram in radians. The angles at which the sinogram slices were recorded do not have to be distributed equidistantly as in backproject(). The angles are internaly converted to modulo PI.
  • initial (ndarray, shape (N,N), optional) – The initial guess for the solution.
  • iterations (int) – Number of iterations to perform.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.

Notes

For theoretical backround, see Kak, A. C., & Slaney, M.. Principles of Computerized Tomographic Imaging, SIAM, (2001)

Sec. 7.2: “ART reconstrutions usually suffer from salt and pepper noise, which is caused by the inconsitencies introuced in the set of equations by the approximations commonly used for \(w_{ik}\) ‘s.”

SART

radontea.sart(sinogram, angles, initial=None, iterations=1, count=None, max_count=None)

Simultaneous Algebraic Reconstruction Technique

SART computes an inverse of the Radon transform in two dimensions. The reconstruction technique uses “rays” of the diameter of one pixel to iteratively solve the system of linear equations that describe the image. The weighting factors are bilinear elements. At the beginning and end of each ray, only partial weights are used. The pixel values of the image are updated only after each iteration is complete.

Parameters:
  • sinogram (ndarrayy, shape (A,N)) – Two-dimensional sinogram of line recordings.
  • angles (ndarray, length A) – Angular positions of the sinogram in radians. The angles at which the sinogram slices were recorded do not have to be distributed equidistantly as in backprojection techniques. The angles are internaly converted to modulo PI.
  • initial (ndarray, shape (N,N), optional) – the initial guess for the solution.
  • iterations (integer) – Number of iterations to perform.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.

Notes

Algebraic reconstruction technique (ART) (see art):
Iterations are performed over each ray of each projection. Weighting factors are binary (1 if center of pixel is within ray, 0 else). This leads to salt and pepper noise.
Simultaneous iterative reconstruction technique (SIRT):
Same idea as ART, but for each iteration, the change of the image f is computed for all rays and projections separately and the weights are applied simultaneously after each iteration. The result is a slower convergence but the final image is also less noisy.

This implementation does NOT use a hamming window to filter the data and to emphasize points at the center of the recon- struction region.

For theoretical backround, see Kak, A. C., & Slaney, M.. Principles of Computerized Tomographic Imaging, SIAM, (2001)

Sec 7.4: “[SART] seems to combine the best of ART and SIRT. […] Here are the main features if SART: First, […] the traditional pixel basis is abandonded in favor of bilinear elements [e.g. interpolation]. Also, for a circular reconstruction region, only partial weights are assigned to the first and last picture elements on the individual rays. To further reduce the noise […], the correction terms are simultaneously applied for all the rays in one projection […].”

Fan beam geometry

fan.radon_fan(arr, det_size[, det_spacing, …]) Compute the Radon transform for a fan beam geometry
fan.get_det_coords(size, spacing) Compute pixel-center positions for 2D detector
fan.get_fan_coords(size, spacing, distance, …) Compute equispaced angular coordinates for 2d detector
fan.lino2sino(linogram, lds[, stepsize, …]) Convert linogram to sinogram for an equispaced detector.
fan.fan_rec(linogram, lds, method[, …]) 2D synthetic aperture reconstruction

Coordinate transforms

radontea.fan.get_det_coords(size, spacing)

Compute pixel-center positions for 2D detector

The centers of the pixels of a detector are usually not aligned to a pixel grid. If we center the detector at the origin, odd images will have a pixel at zero, whereas even images will have two pixels next to the actual center.

Parameters:
  • size (int) – Size of the detector (number of detection points).
  • spacing (float) – Distance between two detection points in pixels.
Returns:

arr – Pixel coordinates.

Return type:

1D ndarray

radontea.fan.get_fan_coords(size, spacing, distance, numang)

Compute equispaced angular coordinates for 2d detector

The centers of the pixels of a detector are usually not aligned to a pixel grid. If we center the detector at the origin, odd images will have a pixel at zero, whereas even images will have two pixels next to the actual center. We want to compute an equispaced array of numang angles that go between the centers of the outmost far pixels of the detector.

Parameters:
  • size (int) – Size of the detector (number of detection points).
  • spacing (float) – Distance between two detection points in pixels.
  • distance (float) – Axial distance from the detector to the angular measurement position in pixels.
  • numang (int) – Number of angles.
Returns:

angles, latpos – Angles and pixel coordinates at the detector.

Return type:

two 1D ndarrays

Notes

Actually one would not need to define spacing and distance, but for convenience, these parameters are separated and an arbitrary uint ‘pixel’ is defined.

radontea.fan.radon_fan(arr, det_size, det_spacing=1, shift_size=1, lS=1, lD=None, return_ang=False, count=None, max_count=None)

Compute the Radon transform for a fan beam geometry

In contrast to radon_parallel(), this function uses (1) a fan-beam geometry (the integral is taken along rays that meet at one point), and (2) translates the object laterally instead of rotating it. The result is sometimes referred to as ‘linogram’.

Problem sketch:

x
^
|
----> z

source       object   detector

          /             . (+det_size/2, lD)
(0,-lS) ./   (0,0)      .
         \              .
          \             . (-det_size/2, lD)

The algorithm computes all angular projections for discrete movements of the object. The position of the object is changed such that its lower boundary starts at (det_size/2, 0) and its upper boundary ends at (-det_size/2, 0) at increments of shift_size.

Parameters:
  • arr (ndarray, shape (N,N)) – the input image.
  • det_size (int) – The total detector size in pixels. The detector centered to the source. The axial position of the detector is the center of the pixels on the far right of the object.
  • det_spacing (float) – Distance between detection points in pixels.
  • shift_size (int) – The amount of pixels that the object is shifted between projections.
  • lS (multiples of 0.5) – Source position relative to the center of arr. lS >= 1.
  • lD (int) – Detector position relative to the center arr. Default is N/2.
  • return_ang (bool) – Also return the angles corresponding to the detector pixels.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
Returns:

outarr – Linogram of the input image. Where N+det_size determines the lateral position of the sample.

Return type:

ndarray of floats, shape (N+det_size,det_size)

See also

scipy.ndimage.interpolation.rotate()
The interpolator used to rotate the image.
radontea.fan.lino2sino(linogram, lds, stepsize=1, det_spacing=1, numang=None, retang=False, count=None, max_count=None)

Convert linogram to sinogram for an equispaced detector.

Parameters:
  • linogram (real 2d ndarray of shape (D, A*)) – Linogram from synthetic aperture measurements.
  • lds (float) – Distance from point source to detector in au.
  • stepsize (float) – Translational increment of object in au (stepsize in D).
  • det_spacing (float) – Distance between detector positions in au.
  • numang (int) – Number of equispaced angles, defaults to linogram.shape[1]
  • retang (bool) – Return the corresponding angles for the sinogram.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
Returns:

sinogram – The distortion-corrected sinogram. If retang is True, then the equispaced angles are returned as well.

Return type:

2d ndarray of shape (D, A)

Notes

This function can be used to convert a linogram obtained with fan-beam tomography to a sinogram, which then can be reconstructed with the backprojection or fourier mapping algorithms.

Reconstruction

radontea.fan.fan_rec(linogram, lds, method, stepsize=1, det_spacing=1, numang=None, count=None, max_count=None, **kwargs)

2D synthetic aperture reconstruction

Computes the inverse of the fan-beam Radon transform using interpolation of the linogram and one of the inverse algorithms for tomography with the Fourier slice theorem.

Parameters:
  • linogram (2d ndarray of shape (D, A)) – Input linogram from the synthetic aperture measurement.
  • lds (float) – Distance in pixels between source and detector.
  • method (callable) – Reconstruction method, e.g. radontea.backproject.
  • numang (int) – Number of angles to be used for the sinogram. A higher number increases quality, but interpolation takes longer. By default numang = linogram.shape[1].
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
  • **kwargs (dict) – Keyword arguments for method.

Volumetric reconstruction

For a slice-wise 3D reconstruction, radontea can use multiprocessing to parallelize the reconstruction process.

Convenience functions

radontea.backproject_3d(sinogram, angles, filtering='ramp', weight_angles=True, padding=True, padval=0, count=None, max_count=None, ncpus=None)

Convenience wrapper for 3D backprojection reconstruction

See backproject() for parameter definitions. The additional parameter ncpus sets the number of CPUs used.

Returns:out – The reconstructed volume.

Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use out.transpose(1, 0, 2).

Return type:ndarray, shape (N,M,N)
radontea.fourier_map_3d(sinogram, angles, intp_method='cubic', count=None, max_count=None, ncpus=None)

Convenience wrapper for 3D Fourier mapping reconstruction

See fourier_map() for parameter definitions. The additional parameter ncpus sets the number of CPUs used.

Returns:out – The reconstructed volume.

Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use out.transpose(1, 0, 2).

Return type:ndarray, shape (N,M,N)

General 3D reconstruction

radontea.volume_recon(func2d, sinogram=None, angles=None, count=None, max_count=None, ncpus=None, **kwargs)

Slice-wise 3D inversion of the Radon transform

Computes the slice-wise 3D inverse of the Radon transform using multiprocessing.

Parameters:
  • func2d (callable) – A method for the slice-wise reconstruction (e.g. backproject()).
  • sinogram (ndarray, shape (A,M,N)) – Three-dimensional sinogram of line recordings. The axis 1 iterates through the M slices. The rotation takes place through axis 1.
  • angles ((A,) ndarray) – Angular positions of the sinogram in radians in the interval [0, PI).
  • max_count (count,) – Can be used to monitor the progress of this function. The value of max_count.value is set initially and the value of count.value is incremented until it reaches the end of the algorithm (max_count.value).
  • **kwargs (dict) – Additional keyword arguments to func2d.
Returns:

out – The reconstructed volume.

Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use out.transpose(1, 0, 2).

Return type:

ndarray, shape (N,M,N)