{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Radon transform\n\nIn computed tomography, the tomography reconstruction problem is to obtain\na tomographic slice image from a set of projections [1]_. A projection is\nformed by drawing a set of parallel rays through the 2D object of interest,\nassigning the integral of the object's contrast along each ray to a single\npixel in the projection. A single projection of a 2D object is one dimensional.\nTo enable computed tomography reconstruction of the object, several projections\nmust be acquired, each of them corresponding to a different angle between the\nrays with respect to the object. A collection of projections at several angles\nis called a sinogram, which is a linear transform of the original image.\n\nThe inverse Radon transform is used in computed tomography to reconstruct\na 2D image from the measured projections (the sinogram). A practical, exact\nimplementation of the inverse Radon transform does not exist, but there are\nseveral good approximate algorithms available.\n\nAs the inverse Radon transform reconstructs the object from a set of\nprojections, the (forward) Radon transform can be used to simulate a\ntomography experiment.\n\nThis script performs the Radon transform to simulate a tomography experiment\nand reconstructs the input image based on the resulting sinogram formed by\nthe simulation. Two methods for performing the inverse Radon transform\nand reconstructing the original image are compared: The Filtered Back\nProjection (FBP) and the Simultaneous Algebraic Reconstruction\nTechnique (SART).\n\nFor further information on tomographic reconstruction, see:\n\n.. [1] AC Kak, M Slaney, \"Principles of Computerized Tomographic Imaging\",\n IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html\n\n.. [2] Wikipedia, Radon transform,\n https://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform\n\n.. [3] S Kaczmarz, \"Angenaeherte Aufloesung von Systemen linearer\n Gleichungen\", Bulletin International de l'Academie Polonaise\n des Sciences et des Lettres, 35 pp 355--357 (1937)\n\n.. [4] AH Andersen, AC Kak, \"Simultaneous algebraic reconstruction\n technique (SART): a superior implementation of the ART algorithm\",\n Ultrasonic Imaging 6 pp 81--94 (1984)\n\n## The forward transform\n\nAs our original image, we will use the Shepp-Logan phantom. When calculating\nthe Radon transform, we need to decide how many projection angles we wish\nto use. As a rule of thumb, the number of projections should be about the\nsame as the number of pixels there are across the object (to see why this\nis so, consider how many unknown pixel values must be determined in the\nreconstruction process and compare this to the number of measurements\nprovided by the projections), and we follow that rule here. Below is the\noriginal image and its Radon transform, often known as its *sinogram*:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom skimage.data import shepp_logan_phantom\nfrom skimage.transform import radon, rescale\n\nimage = shepp_logan_phantom()\nimage = rescale(image, scale=0.4, mode='reflect', channel_axis=None)\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))\n\nax1.set_title(\"Original\")\nax1.imshow(image, cmap=plt.cm.Greys_r)\n\ntheta = np.linspace(0., 180., max(image.shape), endpoint=False)\nsinogram = radon(image, theta=theta)\ndx, dy = 0.5 * 180.0 / max(image.shape), 0.5 / sinogram.shape[0]\nax2.set_title(\"Radon transform\\n(Sinogram)\")\nax2.set_xlabel(\"Projection angle (deg)\")\nax2.set_ylabel(\"Projection position (pixels)\")\nax2.imshow(sinogram, cmap=plt.cm.Greys_r,\n extent=(-dx, 180.0 + dx, -dy, sinogram.shape[0] + dy),\n aspect='auto')\n\nfig.tight_layout()\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reconstruction with the Filtered Back Projection (FBP)\n\nThe mathematical foundation of the filtered back projection is the Fourier\nslice theorem [2]_. It uses Fourier transform of the projection and\ninterpolation in Fourier space to obtain the 2D Fourier transform of the\nimage, which is then inverted to form the reconstructed image. The filtered\nback projection is among the fastest methods of performing the inverse\nRadon transform. The only tunable parameter for the FBP is the filter,\nwhich is applied to the Fourier transformed projections. It may be used to\nsuppress high frequency noise in the reconstruction. ``skimage`` provides\nthe filters 'ramp', 'shepp-logan', 'cosine', 'hamming', and 'hann':\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\nfrom skimage.transform.radon_transform import _get_fourier_filter\n\nfilters = ['ramp', 'shepp-logan', 'cosine', 'hamming', 'hann']\n\nfor ix, f in enumerate(filters):\n response = _get_fourier_filter(2000, f)\n plt.plot(response, label=f)\n\nplt.xlim([0, 1000])\nplt.xlabel('frequency')\nplt.legend()\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Applying the inverse radon transformation with the 'ramp' filter, we get:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.transform import iradon\n\nreconstruction_fbp = iradon(sinogram, theta=theta, filter_name='ramp')\nerror = reconstruction_fbp - image\nprint(f'FBP rms reconstruction error: {np.sqrt(np.mean(error**2)):.3g}')\n\nimkwargs = dict(vmin=-0.2, vmax=0.2)\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5),\n sharex=True, sharey=True)\nax1.set_title(\"Reconstruction\\nFiltered back projection\")\nax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)\nax2.set_title(\"Reconstruction error\\nFiltered back projection\")\nax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reconstruction with the Simultaneous Algebraic Reconstruction Technique\n\nAlgebraic reconstruction techniques for tomography are based on a\nstraightforward idea: for a pixelated image the value of a single ray in a\nparticular projection is simply a sum of all the pixels the ray passes\nthrough on its way through the object. This is a way of expressing the\nforward Radon transform. The inverse Radon transform can then be formulated\nas a (large) set of linear equations. As each ray passes through a small\nfraction of the pixels in the image, this set of equations is sparse,\nallowing iterative solvers for sparse linear systems to tackle the system\nof equations. One iterative method has been particularly popular, namely\nKaczmarz' method [3]_, which has the property that the solution will\napproach a least-squares solution of the equation set.\n\nThe combination of the formulation of the reconstruction problem as a set\nof linear equations and an iterative solver makes algebraic techniques\nrelatively flexible, hence some forms of prior knowledge can be\nincorporated with relative ease.\n\n``skimage`` provides one of the more popular variations of the algebraic\nreconstruction techniques: the Simultaneous Algebraic Reconstruction\nTechnique (SART) [4]_. It uses Kaczmarz' method as the iterative\nsolver. A good reconstruction is normally obtained in a single iteration,\nmaking the method computationally effective. Running one or more extra\niterations will normally improve the reconstruction of sharp, high\nfrequency features and reduce the mean squared error at the expense of\nincreased high frequency noise (the user will need to decide on what number\nof iterations is best suited to the problem at hand. The implementation in\n``skimage`` allows prior information of the form of a lower and upper\nthreshold on the reconstructed values to be supplied to the reconstruction.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.transform import iradon_sart\n\nreconstruction_sart = iradon_sart(sinogram, theta=theta)\nerror = reconstruction_sart - image\nprint(f'SART (1 iteration) rms reconstruction error: '\n f'{np.sqrt(np.mean(error**2)):.3g}')\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 8.5), sharex=True, sharey=True)\nax = axes.ravel()\n\nax[0].set_title(\"Reconstruction\\nSART\")\nax[0].imshow(reconstruction_sart, cmap=plt.cm.Greys_r)\n\nax[1].set_title(\"Reconstruction error\\nSART\")\nax[1].imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)\n\n# Run a second iteration of SART by supplying the reconstruction\n# from the first iteration as an initial estimate\nreconstruction_sart2 = iradon_sart(sinogram, theta=theta,\n image=reconstruction_sart)\nerror = reconstruction_sart2 - image\nprint(f'SART (2 iterations) rms reconstruction error: '\n f'{np.sqrt(np.mean(error**2)):.3g}')\n\nax[2].set_title(\"Reconstruction\\nSART, 2 iterations\")\nax[2].imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)\n\nax[3].set_title(\"Reconstruction error\\nSART, 2 iterations\")\nax[3].imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)\nplt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 0
}