Tomographic reconstruction

Last updated
Tomographic reconstruction: Projection, Back projection and Filtered back projection

Tomographic reconstruction is a type of multidimensional inverse problem where the challenge is to yield an estimate of a specific system from a finite number of projections. The mathematical basis for tomographic imaging was laid down by Johann Radon. A notable example of applications is the reconstruction of computed tomography (CT) where cross-sectional images of patients are obtained in non-invasive manner. Recent developments have seen the Radon transform and its inverse used for tasks related to realistic object insertion required for testing and evaluating computed tomography use in airport security. [1]

Contents

This article applies in general to reconstruction methods for all kinds of tomography, but some of the terms and physical descriptions refer directly to the reconstruction of X-ray computed tomography.

Introducing formula

Figure 1: Parallel beam geometry utilized in tomography and tomographic reconstruction. Each projection, resulting from tomography under a specific angle, is made up of the set of line integrals through the object. Tomographic fig1.png
Figure 1: Parallel beam geometry utilized in tomography and tomographic reconstruction. Each projection, resulting from tomography under a specific angle, is made up of the set of line integrals through the object.
Resulting tomographic image from a plastic skull phantom. Projected X-rays are clearly visible on this slice taken with a CT-scan as image artifacts, due to limited amount of projection slices over angles. Ct skull.jpg
Resulting tomographic image from a plastic skull phantom. Projected X-rays are clearly visible on this slice taken with a CT-scan as image artifacts, due to limited amount of projection slices over angles.

The projection of an object, resulting from the tomographic measurement process at a given angle , is made up of a set of line integrals (see Fig. 1). A set of many such projections under different angles organized in 2D is called a sinogram (see Fig. 3). In X-ray CT, the line integral represents the total attenuation of the beam of X-rays as it travels in a straight line through the object. As mentioned above, the resulting image is a 2D (or 3D) model of the attenuation coefficient. That is, we wish to find the image . The simplest and easiest way to visualise the method of scanning is the system of parallel projection, as used in the first scanners. For this discussion we consider the data to be collected as a series of parallel rays, at position , across a projection at angle . This is repeated for various angles. Attenuation occurs exponentially in tissue:

where is the attenuation coefficient as a function of position. Therefore, generally the total attenuation of a ray at position , on the projection at angle , is given by the line integral:

Using the coordinate system of Figure 1, the value of onto which the point will be projected at angle is given by:

So the equation above can be rewritten as

where represents and is the Dirac delta function. This function is known as the Radon transform (or sinogram) of the 2D object.

The Fourier Transform of the projection can be written as

where [2]
represents a slice of the 2D Fourier transform of at angle . Using the inverse Fourier transform, the inverse Radon transform formula can be easily derived.

where is the derivative of the Hilbert transform of

In theory, the inverse Radon transformation would yield the original image. The projection-slice theorem tells us that if we had an infinite number of one-dimensional projections of an object taken at an infinite number of angles, we could perfectly reconstruct the original object, . However, there will only be a finite number of projections available in practice.

Assuming has effective diameter and desired resolution is , a rule of thumb for the number of projections needed for reconstruction is [2]

Reconstruction algorithms

Practical reconstruction algorithms have been developed to implement the process of reconstruction of a three-dimensional object from its projections. [3] [2] These algorithms are designed largely based on the mathematics of the X-ray transform, statistical knowledge of the data acquisition process and geometry of the data imaging system.

Fourier-domain reconstruction algorithm

Reconstruction can be made using interpolation. Assume projections of are generated at equally spaced angles, each sampled at the same rate. The discrete Fourier transform (DFT) on each projection yields sampling in the frequency domain. Combining all the frequency-sampled projections generates a polar raster in the frequency domain. The polar raster is sparse, so interpolation is used to fill the unknown DFT points, and reconstruction can be done through the inverse discrete Fourier transform. [4] Reconstruction performance may improve by designing methods to change the sparsity of the polar raster, facilitating the effectiveness of interpolation.

For instance, a concentric square raster in the frequency domain can be obtained by changing the angle between each projection as follow:

where is highest frequency to be evaluated.

The concentric square raster improves computational efficiency by allowing all the interpolation positions to be on rectangular DFT lattice. Furthermore, it reduces the interpolation error. [4] Yet, the Fourier-Transform algorithm has a disadvantage of producing inherently noisy output.

Back projection algorithm

In practice of tomographic image reconstruction, often a stabilized and discretized version of the inverse Radon transform is used, known as the filtered back projection algorithm. [2]

With a sampled discrete system, the inverse Radon transform is

where is the angular spacing between the projections and is a Radon kernel with frequency response .

The name back-projection comes from the fact that a one-dimensional projection needs to be filtered by a one-dimensional Radon kernel (back-projected) in order to obtain a two-dimensional signal. The filter used does not contain DC gain, so adding DC bias may be desirable. Reconstruction using back-projection allows better resolution than interpolation method described above. However, it induces greater noise because the filter is prone to amplify high-frequency content.

Iterative reconstruction algorithm

The iterative algorithm is computationally intensive but it allows the inclusion of a priori information about the system . [2]

Let be the number of projections and be the distortion operator for the th projection taken at an angle . are a set of parameters to optimize the conversion of iterations.

A fan-beam reconstruction of Shepp-Logan Phantom with different sensor spacing. Smaller spacing between the sensors allow finer reconstruction. The figure was generated by using MATLAB. Fan-beam reconstruction of Shepp-Logan Phantom.jpg
A fan-beam reconstruction of Shepp-Logan Phantom with different sensor spacing. Smaller spacing between the sensors allow finer reconstruction. The figure was generated by using MATLAB.

An alternative family of recursive tomographic reconstruction algorithms are the algebraic reconstruction techniques and iterative sparse asymptotic minimum variance.

Fan-beam reconstruction

Use of a noncollimated fan beam is common since a collimated beam of radiation is difficult to obtain. Fan beams will generate series of line integrals, not parallel to each other, as projections. The fan-beam system requires a 360-degree range of angles, which imposes mechanical constraints, but it allows faster signal acquisition time, which may be advantageous in certain settings such as in the field of medicine. Back projection follows a similar two-step procedure that yields reconstruction by computing weighted sum back-projections obtained from filtered projections.

Deep learning reconstruction

The influence of Poisson noise in deep learning reconstruction where Poisson noise causes the U-Net fail to reconstruct an existing high contrast lesion-like object. DeepLearningReconstruction.png
The influence of Poisson noise in deep learning reconstruction where Poisson noise causes the U-Net fail to reconstruct an existing high contrast lesion-like object.

Deep learning methods are widely applied to image reconstruction nowadays and have achieved impressive results in various image reconstruction tasks, including low-dose denoising, sparse-view reconstruction, limited angle tomography and metal artifact reduction. An excellent overview can be found in the special issue [5] of IEEE Transaction on Medical Imaging. One group of deep learning reconstruction algorithms apply post-processing neural networks to achieve image-to-image reconstruction, where input images are reconstructed by conventional reconstruction methods. Artifact reduction using the U-Net in limited angle tomography is such an example application. [6] However, incorrect structures may occur in an image reconstructed by such a completely data-driven method, [7] as displayed in the figure. Therefore, integration of known operators into the architecture design of neural networks appears beneficial, as described in the concept of precision learning. [8] For example, direct image reconstruction from projection data can be learnt from the framework of filtered back-projection. [9] Another example is to build neural networks by unrolling iterative reconstruction algorithms. [10] Except for precision learning, using conventional reconstruction methods with deep learning reconstruction prior [11] is also an alternative approach to improve the image quality of deep learning reconstruction.

Tomographic reconstruction software

Tomographic systems have significant variability in their applications and geometries (locations of sources and detectors). This variability creates the need for very specific, tailored implementations of the processing and reconstruction algorithms. Thus, most CT manufacturers provide their own custom proprietary software. This is done not only to protect intellectual property, but may also be enforced by a government regulatory agency. Regardless, there are a number of general purpose tomographic reconstruction software packages that have been developed over the last couple decades, both commercial and open-source.

Most of the commercial software packages that are available for purchase focus on processing data for benchtop cone-beam CT systems. A few of these software packages include Volume Graphics, InstaRecon, iTomography, Livermore Tomography Tools (LTT), and Cone Beam Software Tools (CST).

Some noteworthy examples of open-source reconstruction software include: Reconstruction Toolkit (RTK), [12] CONRAD, [13] TomoPy, [14] the ASTRA toolbox, [15] [16] PYRO-NN, [17] ODL, [18] TIGRE, [19] and LEAP. [20]

Shown in the gallery is the complete process for a simple object tomography and the following tomographic reconstruction based on ART.

See also

Related Research Articles

<span class="mw-page-title-main">Centripetal force</span> Force directed to the center of rotation

A centripetal force is a force that makes a body follow a curved path. The direction of the centripetal force is always orthogonal to the motion of the body and towards the fixed point of the instantaneous center of curvature of the path. Isaac Newton described it as "a force by which bodies are drawn or impelled, or in any way tend, towards a point as to a centre". In Newtonian mechanics, gravity provides the centripetal force causing astronomical orbits.

<span class="mw-page-title-main">Frequency modulation synthesis</span> Form of sound synthesis

Frequency modulation synthesis is a form of sound synthesis whereby the frequency of a waveform is changed by modulating its frequency with a modulator. The (instantaneous) frequency of an oscillator is altered in accordance with the amplitude of a modulating signal.

<span class="mw-page-title-main">Spherical coordinate system</span> Coordinates comprising a distance and two angles

In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a given point in space is specified by three numbers, : the radial distance of the radial liner connecting the point to the fixed point of origin ; the polar angle θ of the radial line r; and the azimuthal angle φ of the radial line r.

<span class="mw-page-title-main">Tautochrone curve</span> Curve for which the time to roll to the end is equal for all starting points

A tautochrone curve or isochrone curve is the curve for which the time taken by an object sliding without friction in uniform gravity to its lowest point is independent of its starting point on the curve. The curve is a cycloid, and the time is equal to π times the square root of the radius over the acceleration of gravity. The tautochrone curve is related to the brachistochrone curve, which is also a cycloid.

Fourier optics is the study of classical optics using Fourier transforms (FTs), in which the waveform being considered is regarded as made up of a combination, or superposition, of plane waves. It has some parallels to the Huygens–Fresnel principle, in which the wavefront is regarded as being made up of a combination of spherical wavefronts whose sum is the wavefront being studied. A key difference is that Fourier optics considers the plane waves to be natural modes of the propagation medium, as opposed to Huygens–Fresnel, where the spherical waves originate in the physical medium.

The Hough transform is a feature extraction technique used in image analysis, computer vision, and digital image processing. The purpose of the technique is to find imperfect instances of objects within a certain class of shapes by a voting procedure. This voting procedure is carried out in a parameter space, from which object candidates are obtained as local maxima in a so-called accumulator space that is explicitly constructed by the algorithm for computing the Hough transform.

<span class="mw-page-title-main">Radon transform</span> Integral transform

In mathematics, the Radon transform is the integral transform which takes a function f defined on the plane to a function Rf 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. The transform was introduced in 1917 by Johann Radon, who also provided a formula for the inverse transform. Radon further included formulas for the transform in three dimensions, in which the integral is taken over planes. It was later generalized to higher-dimensional Euclidean spaces and more broadly in the context of integral geometry. The complex analogue of the Radon transform is known as the Penrose transform. The Radon transform is widely applicable to tomography, the creation of an image from the projection data associated with cross-sectional scans of an object.

In control theory and signal processing, a linear, time-invariant system is said to be minimum-phase if the system and its inverse are causal and stable.

<span class="mw-page-title-main">Phasor</span> Complex number representing a particular sine wave

In physics and engineering, a phasor is a complex number representing a sinusoidal function whose amplitude, and initial phase are time-invariant and whose angular frequency is fixed. It is related to a more general concept called analytic representation, which decomposes a sinusoid into the product of a complex constant and a factor depending on time and frequency. The complex constant, which depends on amplitude and phase, is known as a phasor, or complex amplitude, and sinor or even complexor.

In mathematics, the Hankel transform expresses any given function f(r) as the weighted sum of an infinite number of Bessel functions of the first kind Jν(kr). The Bessel functions in the sum are all of the same order ν, but differ in a scaling factor k along the r axis. The necessary coefficient Fν of each Bessel function in the sum, as a function of the scaling factor k constitutes the transformed function. The Hankel transform is an integral transform and was first developed by the mathematician Hermann Hankel. It is also known as the Fourier–Bessel transform. Just as the Fourier transform for an infinite interval is related to the Fourier series over a finite interval, so the Hankel transform over an infinite interval is related to the Fourier–Bessel series over a finite interval.

In mathematics, a change of variables is a basic technique used to simplify problems in which the original variables are replaced with functions of other variables. The intent is that when expressed in new variables, the problem may become simpler, or equivalent to a better understood problem.

<span class="mw-page-title-main">Sine and cosine</span> Fundamental trigonometric functions

In mathematics, sine and cosine are trigonometric functions of an angle. The sine and cosine of an acute angle are defined in the context of a right triangle: for the specified angle, its sine is the ratio of the length of the side that is opposite that angle to the length of the longest side of the triangle, and the cosine is the ratio of the length of the adjacent leg to that of the hypotenuse. For an angle , the sine and cosine functions are denoted as and .

In geometry, various formalisms exist to express a rotation in three dimensions as a mathematical transformation. In physics, this concept is applied to classical mechanics where rotational kinematics is the science of quantitative description of a purely rotational motion. The orientation of an object at a given instant is described with the same tools, as it is defined as an imaginary rotation from a reference placement in space, rather than an actually observed rotation from a previous placement in space.

<span class="mw-page-title-main">Axis–angle representation</span> Parameterization of a rotation into a unit vector and angle

In mathematics, the axis–angle representation parameterizes a rotation in a three-dimensional Euclidean space by two quantities: a unit vector e indicating the direction (geometry) of an axis of rotation, and an angle of rotation θ describing the magnitude and sense of the rotation about the axis. Only two numbers, not three, are needed to define the direction of a unit vector e rooted at the origin because the magnitude of e is constrained. For example, the elevation and azimuth angles of e suffice to locate it in any particular Cartesian coordinate frame.

The direct-quadrature-zerotransformation or zero-direct-quadraturetransformation is a tensor that rotates the reference frame of a three-element vector or a three-by-three element matrix in an effort to simplify analysis. The DQZ transform is the product of the Clarke transform and the Park transform, first proposed in 1929 by Robert H. Park.

In optics, the Fraunhofer diffraction equation is used to model the diffraction of waves when the diffraction pattern is viewed at a long distance from the diffracting object, and also when it is viewed at the focal plane of an imaging lens.

In physics, and especially scattering theory, the momentum-transfer cross section is an effective scattering cross section useful for describing the average momentum transferred from a particle when it collides with a target. Essentially, it contains all the information about a scattering process necessary for calculating average momentum transfers but ignores other details about the scattering angle.

In mathematical analysis and applications, multidimensional transforms are used to analyze the frequency content of signals in a domain of two or more dimensions.

The problem of reconstructing a multidimensional signal from its projection is uniquely multidimensional, having no 1-D counterpart. It has applications that range from computer-aided tomography to geophysical signal processing. It is a problem which can be explored from several points of view—as a deconvolution problem, a modeling problem, an estimation problem, or an interpolation problem.

<span class="mw-page-title-main">Operation of computed tomography</span>

X-ray computed tomography operates by using an X-ray generator that rotates around the object; X-ray detectors are positioned on the opposite side of the circle from the X-ray source.

References

  1. Najla Megherbi; Toby P. Breckon; Greg T. Flitton; Andre Mouton (October 2013). "Radon Transform based Metal Artefacts Generation in 3D Threat Image Projection" (PDF). Proc. SPIE Optics and Photonics for Counterterrorism, Crime Fighting and Defence. Vol. 8901. SPIE. pp. 1–7. doi:10.1117/12.2028506. S2CID   14001672 . Retrieved 5 November 2013.
  2. 1 2 3 4 5 Dudgeon and Mersereau (1984). Multidimensional digital signal processing. Prentice-Hall.
  3. Herman, G. T., Fundamentals of computerized tomography: Image reconstruction from projection, 2nd edition, Springer, 2009
  4. 1 2 R. Mersereau, A. Oppenheim (1974). "Digital reconstruction of multidimensional signals from their projections". Proceedings of the IEEE. 62 (10): 1319–1338. doi:10.1109/proc.1974.9625. hdl: 1721.1/13788 .
  5. Wang, Ge; Ye, Jong Chu; Mueller, Klaus; Fessler, Jeffrey A (2018). "Image reconstruction is a new frontier of machine learning". IEEE Transactions on Medical Imaging. 37 (6): 1289–1296. doi:10.1109/TMI.2018.2833635. PMID   29870359. S2CID   46931303.
  6. Gu, Jawook; Ye, Jong Chul (2017). Multi-scale wavelet domain residual learning for limited-angle CT reconstruction. Fully3D. pp. 443–447.
  7. Yixing Huang; Tobias Würfl; Katharina Breininger; Ling Liu; Günter Lauritsch; Andreas Maier (2018). Some Investigations on Robustness of Deep Learning in Limited Angle Tomography. MICCAI. doi:10.1007/978-3-030-00928-1_17.
  8. Maier, Andreas K; Syben, Christopher; Stimpel, Bernhard; Wuerfl, Tobias; Hoffmann, Mathis; Schebesch, Frank; Fu, Weilin; Mill, Leonid; Kling, Lasse; Christiansen, Silke (2019). "Learning with known operators reduces maximum error bounds". Nature Machine Intelligence. 1 (8): 373–380. arXiv: 1907.01992 . doi:10.1038/s42256-019-0077-5. PMC   6690833 . PMID   31406960.
  9. Tobias Wuerfl; Mathis Hoffmann; Vincent Christlein; Katharina Breininger; Yixing Huang; Mathias Unberath; Andreas Maier (2018). "Deep Learning Computed Tomography: Learning Projection-Domain Weights from Image Domain in Limited Angle Problems". IEEE Transactions on Medical Imaging. 37 (6): 1454–1463. doi:10.1109/TMI.2018.2833499. PMID   29870373. S2CID   46935914.
  10. J. Adler; O. Öktem (2018). "Learned Primal-Dual Reconstruction". IEEE Transactions on Medical Imaging. 37 (6): 1322–1332. arXiv: 1707.06474 . doi:10.1109/TMI.2018.2799231. PMID   29870362. S2CID   26897002.
  11. Yixing Huang; Alexander Preuhs; Guenter Lauritsch; Michael Manhart; Xiaolin Huang; Andreas Maier (2019). Data Consistent Artifact Reduction for Limited Angle Tomography with Deep Learning Prior. Machine Learning for Medical Image Reconstruction. arXiv: 1908.06792 . doi:10.1007/978-3-030-33843-5_10.
  12. Reconstruction Toolkit (RTK)
  13. Maier, Andreas; Hofmann, Hannes G.; Berger, Martin; Fischer, Peter; Schwemmer, Chris; Wu, Haibo; Müller, Kerstin; Hornegger, Joachim; Choi, Jang-Hwan; Riess, Christian; Keil, Andreas; Fahrig, Rebecca (2013). "CONRAD - A software framework for cone-beam imaging in radiology". Medical Physics. 40 (11): 111914. Bibcode:2013MedPh..40k1914M. doi:10.1118/1.4824926. PMC   3820625 . PMID   24320447.
  14. Gürsoy, Doǧa; De Carlo, Francesco; Xiao, Xianghui; Jacobsen, Chris (2014). "TomoPy: A framework for the analysis of synchrotron tomographic data". Journal of Synchrotron Radiation. 22 (5): 1188–1193. Bibcode:2014SPIE.9212E..0NG. doi:10.1107/S1600577514013939. PMC   4181643 . PMID   25178011.
  15. van Aarle, Wim; Palenstijn, Willem Jan; De Beenhouwer, Jan; Altantzis, Thomas; Bals, Sara; Batenburg, K. Joost; Sijbers, Jan (October 2015). "The ASTRA Toolbox: a platform for advanced algorithm development in electron tomography". Ultramicroscopy. 157: 35–47. doi:10.1016/j.ultramic.2015.05.002. hdl: 10067/1278340151162165141 . PMID   26057688.
  16. van Aarle, Wim; Palenstijn, Willem Jan; Cant, Jeroen; Janssens, Eline; Bleichrodt, Folkert; Dabravolski, Andrei; De Beenhouwer, Jan; Joost Batenburg, K.; Sijbers, Jan (2016). "Fast and flexible X-ray tomography using the ASTRA toolbox". Optics Express. 24 (22): 35–47. Bibcode:2016OExpr..2425129V. doi: 10.1364/OE.24.025129 . hdl: 10067/1392160151162165141 . PMID   27828452.
  17. Syben, Christopher; Michen, Markus; Stimpel, Bernhard; Seitz, Stephan; Ploner, Stefan; Maier, Andreas (2019). "PYRO-NN: Python Reconstruction Operators in Neural Networks". Medical Physics. 46 (11): 5110–5115. arXiv: 1904.13342 . Bibcode:2019MedPh..46.5110S. doi:10.1002/mp.13753. PMC   6899669 . PMID   31389023.
  18. "Odlgroup/Odl". GitHub .
  19. Released by the University of Bath and CERN.
    Biguri, Ander; Dosanjh, Manjit; Hancock, Steven; Soleimani, Manuchehr (2016-09-08). "TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction". Biomedical Physics & Engineering Express. 2 (5): 055010. doi: 10.1088/2057-1976/2/5/055010 . ISSN   2057-1976.
  20. Kim, Hyojin; Champley, Kyle (2023). "Differentiable Forward Projector for X-ray Computed Tomography". ICML. arXiv: 2307.05801 .

Further reading