Tomographic reconstruction

Last updated

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 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 , rule of thumb number of projections needed for reconstruction is [2]

Reconstruction algorithms

Practical reconstruction algorithms have been developed to implement the process of reconstruction of a 3-dimensional object from its projections. [3] [2] These algorithms are designed largely based on the mathematics of the Radon 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 same rate. The Discrete Fourier transform on each projection will yield sampling in the frequency domain. Combining all the frequency-sampled projections would generate a polar raster in the frequency domain. The polar raster will be sparse so interpolation is used to fill the unknown DFT points and reconstruction can be done through 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 radon kernel with frequency response .

The name back-projection comes from the fact that 1D projection needs to be filtered by 1D Radon kernel (back-projected) in order to obtain a 2D signal. The filter used does not contain DC gain, thus 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

Iterative algorithm is computationally intensive but it allows to include a priori information about the system . [2]

Let be the number of projections, be the distortion operator for th projection taken at an angle . are 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 will require 360 degrees range of angles which impose mechanical constraint, however, 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 2 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

For flexible tomographic reconstruction, open source toolboxes are available, such as PYRO-NN, [12] TomoPy, [13] CONRAD, [14] ODL, the ASTRA toolbox, [15] [16] and TIGRE. [17] TomoPy is an open-source Python toolbox to perform tomographic data processing and image reconstruction tasks at the Advanced Photon Source at Argonne National Laboratory. TomoPy toolbox is specifically designed to be easy to use and deploy at a synchrotron facility beamline. It supports reading many common synchrotron data formats from disk through Scientific Data Exchange, [18] and includes several other processing algorithms commonly used for synchrotron data. TomoPy also includes several reconstruction algorithms, which can be run on multi-core workstations and large-scale computing facilities. [19] The ASTRA Toolbox is a MATLAB and Python toolbox of high-performance GPU primitives for 2D and 3D tomography, from 2009–2014 developed by iMinds-Vision Lab, University of Antwerp and since 2014 jointly developed by iMinds-VisionLab (now imec-VisionLab), UAntwerpen and CWI, Amsterdam. The toolbox supports parallel, fan, and cone beam, with highly flexible source/detector positioning. A large number of reconstruction algorithms are available through TomoPy and the ASTRA toolkit, including FBP, Gridrec, ART, SIRT, SART, BART, CGLS, PML, MLEM and OSEM. In 2016, the ASTRA toolbox has been integrated in the TomoPy framework. [20] By integrating the ASTRA toolbox in the TomoPy framework, the optimized GPU-based reconstruction methods become easily available for synchrotron beamline users, and users of the ASTRA toolbox can more easily read data and use TomoPy’s other functionality for data filtering and artifact correction.

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

A centripetal force is a force that makes a body follow a curved path. Its direction 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.

Polar coordinate system Two-dimensional coordinate system where each point is determined by a distance from reference point and an angle from a reference direction

In mathematics, the polar coordinate system is a two-dimensional coordinate system in which each point on a plane is determined by a distance from a reference point and an angle from a reference direction. The reference point is called the pole, and the ray from the pole in the reference direction is the polar axis. The distance from the pole is called the radial coordinate, radial distance or simply radius, and the angle is called the angular coordinate, polar angle, or azimuth. The radial coordinate is often denoted by r or ρ, and the angular coordinate by φ, θ, or t. Angles in polar notation are generally expressed in either degrees or radians.

Spherical coordinate system 3-dimensional coordinate system

In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a point is specified by three numbers: the radial distance of that point from a fixed origin, its polar angle measured from a fixed zenith direction, and the azimuthal angle of its orthogonal projection on a reference plane that passes through the origin and is orthogonal to the zenith, measured from a fixed reference direction on that plane. It can be seen as the three-dimensional version of the polar coordinate system.

Angular displacement

Angular displacement of a body is the angle in radians through which a point revolves around a centre or line has been rotated in a specified sense about a specified axis. When a body rotates about its axis, the motion cannot simply be analyzed as a particle, as in circular motion it undergoes a changing velocity and acceleration at any time (t). When dealing with the rotation of a body, it becomes simpler to consider the body itself rigid. A body is generally considered rigid when the separations between all the particles remains constant throughout the body's motion, so for example parts of its mass are not flying off. In a realistic sense, all things can be deformable, however this impact is minimal and negligible. Thus the rotation of a rigid body over a fixed axis is referred to as rotational motion.

Synchrotron radiation

Synchrotron radiation is the electromagnetic radiation emitted when charged particles are accelerated radially, e.g., when they are subject to an acceleration perpendicular to their velocity. It is produced, for example, in synchrotrons using bending magnets, undulators and/or wigglers. If the particle is non-relativistic, the emission is called cyclotron emission. If the particles are relativistic, sometimes referred to as ultrarelativistic, the emission is called synchrotron emission. Synchrotron radiation may be achieved artificially in synchrotrons or storage rings, or naturally by fast electrons moving through magnetic fields. The radiation produced in this way has a characteristic polarization and the frequencies generated can range over the entire electromagnetic spectrum, which is also called continuum radiation.

Tautochrone curve

A tautochrone 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.

Inverted pendulum

An inverted pendulum is a pendulum that has its center of mass above its pivot point. It is unstable and without additional help will fall over. It can be suspended stably in this inverted position by using a control system to monitor the angle of the pole and move the pivot point horizontally back under the center of mass when it starts to fall over, keeping it balanced. The inverted pendulum is a classic problem in dynamics and control theory and is used as a benchmark for testing control strategies. It is often implemented with the pivot point mounted on a cart that can move horizontally under control of an electronic servo system as shown in the photo; this is called a cart and pole apparatus. Most applications limit the pendulum to 1 degree of freedom by affixing the pole to an axis of rotation. Whereas a normal pendulum is stable when hanging downwards, an inverted pendulum is inherently unstable, and must be actively balanced in order to remain upright; this can be done either by applying a torque at the pivot point, by moving the pivot point horizontally as part of a feedback system, changing the rate of rotation of a mass mounted on the pendulum on an axis parallel to the pivot axis and thereby generating a net torque on the pendulum, or by oscillating the pivot point vertically. A simple demonstration of moving the pivot point in a feedback system is achieved by balancing an upturned broomstick on the end of one's finger.

Tomography Imaging by sections or sectioning using a penetrative wave

Tomography is imaging by sections or sectioning through the use of any kind of penetrating wave. The method is used in radiology, archaeology, biology, atmospheric science, geophysics, oceanography, plasma physics, materials science, astrophysics, quantum information, and other areas of science. The word tomography is derived from Ancient Greek τόμος tomos, "slice, section" and γράφω graphō, "to write" or, in this context as well, " to describe." A device used in tomography is called a tomograph, while the image produced is a tomogram.

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.

Radon 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 linear algebra, linear transformations can be represented by matrices. If is a linear transformation mapping to and is a column vector with entries, then

Phasor

In physics and engineering, a phasor, is a complex number representing a sinusoidal function whose amplitude (A), angular frequency (ω), and initial phase (θ) are time-invariant. 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.

Etendue or étendue is a property of light in an optical system, which characterizes how "spread out" the light is in area and angle. It corresponds to the beam parameter product (BPP) in Gaussian beam optics.

Sinusoidal plane-wave solutions are particular solutions to the electromagnetic wave equation.

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.

Axis–angle representation

In mathematics, the axis–angle representation of a rotation parameterizes a rotation in a three-dimensional Euclidean space by two quantities: a unit vector e indicating the direction of an axis of rotation, and an angle θ describing the magnitude 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.

In algebra, casus irreducibilis is one of the cases that may arise in attempting to solve polynomials of degree 3 or higher with integer coefficients, to obtain roots that are expressed with radicals. It shows that many algebraic numbers are real-valued but cannot be expressed in radicals without introducing complex numbers. The most notable occurrence of casus irreducibilis is in the case of cubic polynomials that are irreducible over the rational numbers and have three real roots, which was proven by Pierre Wantzel in 1843. One can decide whether a given irreducible cubic polynomial is in casus irreducibilis using the discriminant Δ, via Cardano's formula. Let the cubic equation be given by

TomoPy

TomoPy is an open-sourced Python toolbox to perform tomographic data processing and image reconstruction.

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.

Operation of computed tomography

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. Megherbi, N., Breckon, T.P., Flitton, G.T., Mouton, A. (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. 8901. SPIE. pp. 1–7. doi:10.1117/12.2028506 . Retrieved 5 November 2013.CS1 maint: multiple names: authors list (link)
  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 and Ye, Jong Chu and Mueller, Klaus and 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.CS1 maint: multiple names: authors list (link)
  6. Gu, Jawook and Ye, Jong Chul (2017). Multi-scale wavelet domain residual learning for limited-angle CT reconstruction. Fully3D. pp. 443–447.CS1 maint: multiple names: authors list (link)
  7. Huang Y., Würfl T., Breininger K., Liu L., Lauritsch G., Maier A. (2018). Some Investigations on Robustness of Deep Learning in Limited Angle Tomography. MICCAI. doi:10.1007/978-3-030-00928-1_17.CS1 maint: multiple names: authors list (link)
  8. Maier, Andreas K and Syben, Christopher and Stimpel, Bernhard and Wuerfl, Tobias and Hoffmann, Mathis and Schebesch, Frank and Fu, Weilin and Mill, Leonid and Kling, Lasse and Christiansen, Silke (2019). "Learning with known operators reduces maximum error bounds". Nature Machine Intelligence. 1 (8): 373–380. doi:10.1038/s42256-019-0077-5. PMC   6690833 . PMID   31406960.CS1 maint: multiple names: authors list (link)
  9. Tobias Wuerfl and Mathis Hoffmann and Vincent Christlein and Katharina Breininger and Yixing Huang and Mathias Unberath and 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.
  10. J. Adler and 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.
  11. Huang Y., Preuhs A., Lauritsch G., Manhart M., Huang X., Maier A. (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.CS1 maint: multiple names: authors list (link)
  12. 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:2019arXiv190413342S. doi:10.1002/mp.13753. PMC   6899669 . PMID   31389023.CS1 maint: multiple names: authors list (link)
  13. Gursoy D, De Carlo F, Xiao X, and Jacobsen C (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.CS1 maint: multiple names: authors list (link)
  14. A. Maier, H. G. Hofmann, M. Berger, P. Fischer, C. Schwemmer, H. Wu, K. Mueller, J. Hornegger, J. Choi, C. Riess, A. Keil, A. Farhig (2013). "CONRAD - A software framework for cone-beam imaging in radiology". Medical Physics. 40 (11): 111914. doi:10.1118/1.4824926. PMC   3820625 . PMID   24320447.CS1 maint: multiple names: authors list (link)
  15. Van Aarle, W., Palenstijn, W.J., De Beenhouwer, J., Altantzis T., Bals S., Batenburg K. J., and J. Sijbers (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. PMID   26057688.CS1 maint: multiple names: authors list (link)
  16. W. Van Aarle, W J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K. J. Batenburg, and J. Sijbers (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 . PMID   27828452.CS1 maint: multiple names: authors list (link)
  17. 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.
  18. De Carlo F, Gursoy D, Marone F, Rivers M, Parkinson YD, Khan F, Schwarz N, Vine DJ, Vogt S, Gleber SC, Narayanan S, Newville M, Lanzirotti T, Sun Y, Hong YP, Jacobsen C (2014). "Scientific Data Exchange: a schema for HDF5-based storage of raw and analyzed data". Journal of Synchrotron Radiation. 22 (6): 35–47. doi:10.1107/S160057751401604X. PMID   25343788.CS1 maint: multiple names: authors list (link)
  19. Bicer T, Gursoy D, Kettimuthu R, De Carlo F, and Foster I (2016). "Optimization of tomographic reconstruction workflows on geographically distributed resources". Journal of Synchrotron Radiation. 23 (4): 997–1005. doi:10.1107/S1600577516007980. PMC   5315096 . PMID   27359149.CS1 maint: multiple names: authors list (link)
  20. Pelt DM, Gursoy D, Batenburg KJ, De Carlo F, Palenstijna WJ, and Sijbers J (2016). "Integration of TomoPy and the ASTRA toolbox for advanced processing and reconstruction of tomographic synchrotron data". Journal of Synchrotron Radiation. 23 (3): 842–849. doi:10.1107/S1600577516005658. PMC   5315009 . PMID   27140167.CS1 maint: multiple names: authors list (link)

Further reading