Radon transform

In computed tomography, the tomography reconstruction problem is to obtain a tomographic slice image from a set of projections [1]. A projection is formed by drawing a set of parallel rays through the 2D object of interest, assigning the integral of the object’s contrast along each ray to a single pixel in the projection. A single projection of a 2D object is one dimensional. To enable computed tomography reconstruction of the object, several projections must be acquired, each of them corresponding to a different angle between the rays with respect to the object. A collection of projections at several angles is called a sinogram, which is a linear transform of the original image.

The inverse Radon transform is used in computed tomography to reconstruct a 2D image from the measured projections (the sinogram). A practical, exact implementation of the inverse Radon transform does not exist, but there are several good approximate algorithms available.

As the inverse Radon transform reconstructs the object from a set of projections, the (forward) Radon transform can be used to simulate a tomography experiment.

This script performs the Radon transform to simulate a tomography experiment and reconstructs the input image based on the resulting sinogram formed by the simulation. Two methods for performing the inverse Radon transform and reconstructing the original image are compared: The Filtered Back Projection (FBP) and the Simultaneous Algebraic Reconstruction Technique (SART).

See also

The forward transform

As our original image, we will use the Shepp-Logan phantom. When calculating the Radon transform, we need to decide how many projection angles we wish to use. As a rule of thumb, the number of projections should be about the same as the number of pixels there are across the object (to see why this is so, consider how many unknown pixel values must be determined in the reconstruction process and compare this to the number of measurements provided by the projections), and we follow that rule here. Below is the original image and its Radon transform, often known as its _sinogram_:

from __future__ import print_function, division

import numpy as np
import matplotlib.pyplot as plt

from skimage.io import imread
from skimage import data_dir
from skimage.transform import radon, rescale

image = imread(data_dir + "/phantom.png", as_grey=True)
image = rescale(image, scale=0.4)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))

ax1.set_title("Original")
ax1.imshow(image, cmap=plt.cm.Greys_r)

theta = np.linspace(0., 180., max(image.shape), endpoint=True)
sinogram = radon(image, theta=theta, circle=True)
ax2.set_title("Radon transform\n(Sinogram)")
ax2.set_xlabel("Projection angle (deg)")
ax2.set_ylabel("Projection position (pixels)")
ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
           extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')

fig.subplots_adjust(hspace=0.4, wspace=0.5)
plt.show()

../_images/plot_radon_transform_1.png

Reconstruction with the Filtered Back Projection (FBP)

The mathematical foundation of the filtered back projection is the Fourier slice theorem [2]. It uses Fourier transform of the projection and interpolation in Fourier space to obtain the 2D Fourier transform of the image, which is then inverted to form the reconstructed image. The filtered back projection is among the fastest methods of performing the inverse Radon transform. The only tunable parameter for the FBP is the filter, which is applied to the Fourier transformed projections. It may be used to suppress high frequency noise in the reconstruction. skimage provides a few different options for the filter.

from skimage.transform import iradon

reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)
error = reconstruction_fbp - image
print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))

imkwargs = dict(vmin=-0.2, vmax=0.2)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
ax1.set_title("Reconstruction\nFiltered back projection")
ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
ax2.set_title("Reconstruction error\nFiltered back projection")
ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()

../_images/plot_radon_transform_2.png

Reconstruction with the Simultaneous Algebraic Reconstruction Technique

Algebraic reconstruction techniques for tomography are based on a straightforward idea: for a pixelated image the value of a single ray in a particular projection is simply a sum of all the pixels the ray passes through on its way through the object. This is a way of expressing the forward Radon transform. The inverse Radon transform can then be formulated as a (large) set of linear equations. As each ray passes through a small fraction of the pixels in the image, this set of equations is sparse, allowing iterative solvers for sparse linear systems to tackle the system of equations. One iterative method has been particularly popular, namely Kaczmarz’ method [3], which has the property that the solution will approach a least-squares solution of the equation set.

The combination of the formulation of the reconstruction problem as a set of linear equations and an iterative solver makes algebraic techniques relatively flexible, hence some forms of prior knowledge can be incorporated with relative ease.

skimage provides one of the more popular variations of the algebraic reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique (SART) [1] [4]. It uses Kaczmarz’ method [3] as the iterative solver. A good reconstruction is normally obtained in a single iteration, making the method computationally effective. Running one or more extra iterations will normally improve the reconstruction of sharp, high frequency features and reduce the mean squared error at the expense of increased high frequency noise (the user will need to decide on what number of iterations is best suited to the problem at hand. The implementation in skimage allows prior information of the form of a lower and upper threshold on the reconstructed values to be supplied to the reconstruction.

from skimage.transform import iradon_sart

reconstruction_sart = iradon_sart(sinogram, theta=theta)
error = reconstruction_sart - image
print('SART (1 iteration) rms reconstruction error: %.3g'
      % np.sqrt(np.mean(error**2)))

fig, ax = plt.subplots(2, 2, figsize=(8, 8.5))
ax1, ax2, ax3, ax4 = ax.ravel()
ax1.set_title("Reconstruction\nSART")
ax1.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
ax2.set_title("Reconstruction error\nSART")
ax2.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)

# Run a second iteration of SART by supplying the reconstruction
# from the first iteration as an initial estimate
reconstruction_sart2 = iradon_sart(sinogram, theta=theta,
                                   image=reconstruction_sart)
error = reconstruction_sart2 - image
print('SART (2 iterations) rms reconstruction error: %.3g'
      % np.sqrt(np.mean(error**2)))

ax3.set_title("Reconstruction\nSART, 2 iterations")
ax3.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
ax4.set_title("Reconstruction error\nSART, 2 iterations")
ax4.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()

../_images/plot_radon_transform_3.png
[1](1, 2) AC Kak, M Slaney, “Principles of Computerized Tomographic Imaging”, IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html
[2]Wikipedia, Radon transform, http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform
[3](1, 2) S Kaczmarz, “Angenaeherte Aufloesung von Systemen linearer Gleichungen”, Bulletin International de l’Academie Polonaise des Sciences et des Lettres 35 pp 355–357 (1937)
[4]AH Andersen, AC Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm”, Ultrasonic Imaging 6 pp 81–94 (1984)

STDOUT


        

STDERR


        

Python source code: download (generated using skimage 0.11dev)

IPython Notebook: download (generated using skimage 0.11dev)

CmZyb20gX19mdXR1cmVfXyBpbXBvcnQgcHJpbnRfZnVuY3Rpb24sIGRpdmlzaW9uCgppbXBvcnQgbnVtcHkgYXMgbnAKaW1wb3J0IG1hdHBsb3RsaWIucHlwbG90IGFzIHBsdAoKZnJvbSBza2ltYWdlLmlvIGltcG9ydCBpbXJlYWQKZnJvbSBza2ltYWdlIGltcG9ydCBkYXRhX2Rpcgpmcm9tIHNraW1hZ2UudHJhbnNmb3JtIGltcG9ydCByYWRvbiwgcmVzY2FsZQoKaW1hZ2UgPSBpbXJlYWQoZGF0YV9kaXIgKyAiL3BoYW50b20ucG5nIiwgYXNfZ3JleT1UcnVlKQppbWFnZSA9IHJlc2NhbGUoaW1hZ2UsIHNjYWxlPTAuNCkKCmZpZywgKGF4MSwgYXgyKSA9IHBsdC5zdWJwbG90cygxLCAyLCBmaWdzaXplPSg4LCA0LjUpKQoKYXgxLnNldF90aXRsZSgiT3JpZ2luYWwiKQpheDEuaW1zaG93KGltYWdlLCBjbWFwPXBsdC5jbS5HcmV5c19yKQoKdGhldGEgPSBucC5saW5zcGFjZSgwLiwgMTgwLiwgbWF4KGltYWdlLnNoYXBlKSwgZW5kcG9pbnQ9VHJ1ZSkKc2lub2dyYW0gPSByYWRvbihpbWFnZSwgdGhldGE9dGhldGEsIGNpcmNsZT1UcnVlKQpheDIuc2V0X3RpdGxlKCJSYWRvbiB0cmFuc2Zvcm1cbihTaW5vZ3JhbSkiKQpheDIuc2V0X3hsYWJlbCgiUHJvamVjdGlvbiBhbmdsZSAoZGVnKSIpCmF4Mi5zZXRfeWxhYmVsKCJQcm9qZWN0aW9uIHBvc2l0aW9uIChwaXhlbHMpIikKYXgyLmltc2hvdyhzaW5vZ3JhbSwgY21hcD1wbHQuY20uR3JleXNfciwKICAgICAgICAgICBleHRlbnQ9KDAsIDE4MCwgMCwgc2lub2dyYW0uc2hhcGVbMF0pLCBhc3BlY3Q9J2F1dG8nKQoKZmlnLnN1YnBsb3RzX2FkanVzdChoc3BhY2U9MC40LCB3c3BhY2U9MC41KQpwbHQuc2hvdygp
CmZyb20gc2tpbWFnZS50cmFuc2Zvcm0gaW1wb3J0IGlyYWRvbgoKcmVjb25zdHJ1Y3Rpb25fZmJwID0gaXJhZG9uKHNpbm9ncmFtLCB0aGV0YT10aGV0YSwgY2lyY2xlPVRydWUpCmVycm9yID0gcmVjb25zdHJ1Y3Rpb25fZmJwIC0gaW1hZ2UKcHJpbnQoJ0ZCUCBybXMgcmVjb25zdHJ1Y3Rpb24gZXJyb3I6ICUuM2cnICUgbnAuc3FydChucC5tZWFuKGVycm9yKioyKSkpCgppbWt3YXJncyA9IGRpY3Qodm1pbj0tMC4yLCB2bWF4PTAuMikKZmlnLCAoYXgxLCBheDIpID0gcGx0LnN1YnBsb3RzKDEsIDIsIGZpZ3NpemU9KDgsIDQuNSkpCmF4MS5zZXRfdGl0bGUoIlJlY29uc3RydWN0aW9uXG5GaWx0ZXJlZCBiYWNrIHByb2plY3Rpb24iKQpheDEuaW1zaG93KHJlY29uc3RydWN0aW9uX2ZicCwgY21hcD1wbHQuY20uR3JleXNfcikKYXgyLnNldF90aXRsZSgiUmVjb25zdHJ1Y3Rpb24gZXJyb3JcbkZpbHRlcmVkIGJhY2sgcHJvamVjdGlvbiIpCmF4Mi5pbXNob3cocmVjb25zdHJ1Y3Rpb25fZmJwIC0gaW1hZ2UsIGNtYXA9cGx0LmNtLkdyZXlzX3IsICoqaW1rd2FyZ3MpCnBsdC5zaG93KCk=
CmZyb20gc2tpbWFnZS50cmFuc2Zvcm0gaW1wb3J0IGlyYWRvbl9zYXJ0CgpyZWNvbnN0cnVjdGlvbl9zYXJ0ID0gaXJhZG9uX3NhcnQoc2lub2dyYW0sIHRoZXRhPXRoZXRhKQplcnJvciA9IHJlY29uc3RydWN0aW9uX3NhcnQgLSBpbWFnZQpwcmludCgnU0FSVCAoMSBpdGVyYXRpb24pIHJtcyByZWNvbnN0cnVjdGlvbiBlcnJvcjogJS4zZycKICAgICAgJSBucC5zcXJ0KG5wLm1lYW4oZXJyb3IqKjIpKSkKCmZpZywgYXggPSBwbHQuc3VicGxvdHMoMiwgMiwgZmlnc2l6ZT0oOCwgOC41KSkKYXgxLCBheDIsIGF4MywgYXg0ID0gYXgucmF2ZWwoKQpheDEuc2V0X3RpdGxlKCJSZWNvbnN0cnVjdGlvblxuU0FSVCIpCmF4MS5pbXNob3cocmVjb25zdHJ1Y3Rpb25fc2FydCwgY21hcD1wbHQuY20uR3JleXNfcikKYXgyLnNldF90aXRsZSgiUmVjb25zdHJ1Y3Rpb24gZXJyb3JcblNBUlQiKQpheDIuaW1zaG93KHJlY29uc3RydWN0aW9uX3NhcnQgLSBpbWFnZSwgY21hcD1wbHQuY20uR3JleXNfciwgKippbWt3YXJncykKCiMgUnVuIGEgc2Vjb25kIGl0ZXJhdGlvbiBvZiBTQVJUIGJ5IHN1cHBseWluZyB0aGUgcmVjb25zdHJ1Y3Rpb24KIyBmcm9tIHRoZSBmaXJzdCBpdGVyYXRpb24gYXMgYW4gaW5pdGlhbCBlc3RpbWF0ZQpyZWNvbnN0cnVjdGlvbl9zYXJ0MiA9IGlyYWRvbl9zYXJ0KHNpbm9ncmFtLCB0aGV0YT10aGV0YSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBpbWFnZT1yZWNvbnN0cnVjdGlvbl9zYXJ0KQplcnJvciA9IHJlY29uc3RydWN0aW9uX3NhcnQyIC0gaW1hZ2UKcHJpbnQoJ1NBUlQgKDIgaXRlcmF0aW9ucykgcm1zIHJlY29uc3RydWN0aW9uIGVycm9yOiAlLjNnJwogICAgICAlIG5wLnNxcnQobnAubWVhbihlcnJvcioqMikpKQoKYXgzLnNldF90aXRsZSgiUmVjb25zdHJ1Y3Rpb25cblNBUlQsIDIgaXRlcmF0aW9ucyIpCmF4My5pbXNob3cocmVjb25zdHJ1Y3Rpb25fc2FydDIsIGNtYXA9cGx0LmNtLkdyZXlzX3IpCmF4NC5zZXRfdGl0bGUoIlJlY29uc3RydWN0aW9uIGVycm9yXG5TQVJULCAyIGl0ZXJhdGlvbnMiKQpheDQuaW1zaG93KHJlY29uc3RydWN0aW9uX3NhcnQyIC0gaW1hZ2UsIGNtYXA9cGx0LmNtLkdyZXlzX3IsICoqaW1rd2FyZ3MpCnBsdC5zaG93KCk=