Abstract
Computed tomography (CT) has been widely applied in medical diagnosis, nondestructive evaluation, homeland security, and other science and engineering applications. Image reconstruction is one of the core CT imaging technologies. In this review paper, we systematically reviewed the currently publicly available CT image reconstruction open source toolkits in the aspects of their environments, object models, imaging geometries, and algorithms. In addition to analytic and iterative algorithms, deep learning reconstruction networks and open codes are also reviewed as the third category of reconstruction algorithms. This systematic summary of the publicly available software platforms will help facilitate CT research and development.
Introduction
Computed tomography (CT) can be used to visualize and assess objects without physical damage, which offers many benefits in related fields. In medical diagnosis, CT plays an important role in clinical routine examinations. Biomedical research often relies on micro-CT scanning devices to study the effects of substances on small animals [1]. In industries, CT can be used for nondestructive testing: cavities or cracks can be detected, particles in materials can be analyzed, and the internal and external geometric shapes of complex parts can also be reconstructed [2].
CT system mainly includes scanning devices and a data processing system. The data processing system processes the data obtained by scanning devices. Image reconstruction is the core part of data processing. Figure 1 shows the number of times that the phrase “CT image reconstruction” has been searched on the Web of Science since 2001. Each bar represents a span of every three years. In 2001–2003, the results show that “CT image reconstruction” has been searched 1938 times. By 2016–2018, the results show a much higher amount, a total of 7835 searches. It can be seen from the diagram that there have been more and more studies on CT image reconstruction in the past 18 years, and it has been increasing every year that goes by.
Retrieval results with the phrase “CT Image Reconstruction” on the Web of Science from 2001 to 2018.
Commercial CT image reconstruction packages are reliable and usually have user-friendly operation interfaces. These performances are important to many CT users. However, those commercial packages are not flexible for researchers and high-end users. Moreover, the prices of commercial packages are too high. Typically, a reconstruction software comes bundled with a scanner, which usually encrypts the collected and reconstructed data. Besides, commercial software usually cannot support user-defined scanning trajectories and data truncations.
Many research teams have developed their open source image reconstruction toolkits. These toolkits provide an expressway for a beginner to start the research. In this paper, open source reconstruction toolkits, CT system models, some reconstruction algorithms, and deep learning issues are briefly reviewed and discussed. In section 2, some important concepts and CT imaging system models are introduced, including projection geometries, object models and (back-) projection models. Then section 3 discusses the reconstruction algorithms included in these toolkits. Section 4 introduces each toolkit from the aspects of CT models and algorithms. Section 5 discusses deep learning techniques that apply to image reconstruction. Finally, section 6 concludes this work.
A CT system is usually determined by X-ray source, detector, and scanning trajectory. In CT image reconstruction, a mathematical model is necessary for describing a CT system. CT system model is also important to the simulation study. In this section, some CT system model factors will be introduced.
X-ray source and detector
Parallel beam geometry is an architecture for the first generation CT. Because only one X-ray focal spot and one detector element are used, translation and rotation are needed to get complete projection data set, and the scanning speed is very slow. The second generation CT consists of an X-ray source and a short detector. Translation is still needed while the rotation times are reduced. To avoid translation, wide aperture fan beam CT that can cover the whole object was developed. Curved-detector with equal angle intervals is usually used in medical CT scanners while line-detector with equal distance intervals is common in industrial CT systems. To further reduce the scanning time and increase the resolution along z direction, cone beam CT with multi-row detectors or planar detectors were developed. In some CT systems, more than one source-detector pairs are used for reducing scanning time or realizing multi-energy imaging [3, 4].
Scanning trajectory
Scanning trajectory generally refers to the trajectory of the ray source relative to the object. In commercial CT scanners, circular and helical scanning trajectories are commonly used. In cone beam CT, circular scanning cannot get complete projection data of a 3D object. The slice far from the trajectory plane suffers from cone beam artifacts. To solve this problem, saddle trajectory [5] and compositing-circling trajectory [6] were developed. Helical scanning trajectory that is common in medical CT scanners is designed for long objects [7]. In industrial field, to scan the objects with special shapes or satisfy on-site conditions, some special scanning trajectories were developed [8].
Object modeling
Object modeling is needed in the process of reconstruction or projection simulation. Generally, the models are divided into two categories: analytical model and discrete model. Analytical model utilizes two-dimensional or three-dimensional analytical functions to describe objects, such as Shepp-Logan phantom [9]. In the discrete model [10], the object is discretized by different basis functions like pixel, voxel, and blob [11].
Projection/backprojection model
Projection and backprojection are the cores of CT data simulation and image reconstruction. Projection is an operator that computes the projection of an object model. The projection of a ray is the sum of the points in the object along this ray. Backprojection is the adjoint of projection. It assigns the projection of a ray to each point in this ray.
In analytical object model, for each ray, the projection can be simulated by the following steps: (1) According to the expressions of the ray and each regular region with constant value, calculating the intersection length of the ray and each region. (2) Getting the projection of the ray passing through each region by multiplying the length by region value. (3) Summing these projections to get the projection of this ray.
When an object is expressed with discrete model, for instance, pixel model, the following models are usually used to describe the projection/backprojection: pixel driven, ray driven, distance driven, area model, etc. For more details, please refer to [12]. For more realistic models, noise and scatter should be considered [13–15].
Reconstruction algorithms
Image reconstruction algorithms include analytic reconstruction algorithms and iterative reconstruction algorithms [11]. In analytic algorithms, filters and back-projectors are needed. The iterative algorithms include forward projection, error correction, back-projection, and updating the image data.
Analytic algorithms
FBP algorithm is typically applied to 2D parallel-beam and fan-beam CT, based on 2D Fourier central slice theorem. The direct Fourier inverse transform is also given in some reconstruction toolkits. The Rho-filtered layergram reconstruction method is for image deblurring that obtained by back-projection [16]. The linogram method is a faster reconstruction algorithm than FBP [17, 18].
For circular cone-beam geometry, the Feldkamp-Davis-Kress (FDK) algorithm is commonly used. The FDK algorithm converts the 3D reconstruction into the 2D reconstruction (a fan-beam reconstructions on a tilted plane determined by the source and a filtering line) [19]. For central slice, FDK is the same as FBP. For spiral/helical cone-beam scanning, the first approximate algorithm was proposed in 1991 [20]. From 2002 to 2004 two exact reconstruction algorithms based on PI lines were proposed [21–23]. A detailed review of approximate and exact cone-beam reconstruction with standard and non-standard spiral scanning can be found in [24].
Iterative algorithms
In general, iterative reconstruction scheme consists of a projection model and a backprojection of the error in projection domain [25]. Iterative reconstruction algorithms have advantages in reducing image noise and various artifacts. With the use of prior information, image quality can be greatly improved, especially for sparse or incomplete data. Since the computational hardware has been greatly developed, it is possible for iterative algorithms to be popularly used.
Algebraic reconstruction technique (ART) [26] updates the reconstructed results ray by ray while simultaneous iterative reconstruction technique (SIRT) and simultaneous ART (SART) update the reconstructed results by the average error of all rays. Ordered subset (OS) -based methods update the reconstructed results one subset by one subset. It converges faster as the number of subsets increases. When the idea of ordered subset is applied to SIRT and SART, it leads to OS-SIRT and OS-SART. Each of these ART-like iterative reconstruction algorithms can be seen as a specific case of general Landweber scheme [27].
Based on the parameter estimation theory, statistical iterative reconstruction algorithms have a good ability in denoising. Some statistical methods, like maximum likelihood expectation-maximization (MLEM) algorithm, are based on the maximum likelihood principle [28]. Expectation maximization (EM) method is used for parameter estimation. Ordered subset EM (OSEM) algorithms are also proposed to increase the rate of convergence. EMAP, a maximum a posteriori probability (MAP) algorithm based on modified EM algorithm, has a faster convergence rate and smoother image quality [29].
The iterative image reconstruction algorithm is generally based on an optimization model. Some statistical methods, for example, weighted least-squares (WLS) [30] algorithm that considers second-order statistical properties, are based on least squares principle. In the family of least squares algorithms, methods based on iterative coordinate descent (ICD) [31] are successfully applied to image reconstruction. Iterated conditional modes (ICM) method is essentially an ICD algorithm. It has some important advantages, for example, ICM can add robust results to OSL (one-step-late, based on Gibbs prior, is unstable when the smoothing parameter getting larger) algorithm [32]. In the optimization model, the objective functions can be solved by quadratic optimization techniques [33]. Conjugate gradient least squares (CGLS) algorithm is one of the quadratic optimization techniques. It has advantages in convergence rate, simplicity and potential for parallelization, compared to general gradient-based methods [34]. Steepest descent methods also belong to quadratic optimization techniques. The projection onto steepest descent (POSD) algorithm is a steepest descent algorithm that is projection controlled and the extensive POSD algorithms have better robustness [35]. Gradient ascent (GA) is a gradient-based optimization algorithm that can be used in imaging processing. Fast iterative shrinkage-thresholding algorithm (FISTA) is for inverse problem that is attractive due to computational simplicity and a global rate of convergence [36].
Regularization methods are to find solutions to the optimization models for ill-posed problems. Total variation (TV) minimization is a regularization method that has advantages of preserving the sharp edges and denoising [37]. Soft threshold filtering algorithm is also applied to limited angle reconstruction [38]. The projection onto convex sets (POCS) is to find an intersection of several well-defined closed convex sets as a solution [39]. Alternating direction method of multipliers (ADMM) can be applied to the distributed convex optimization of large-scale problems [40].
Image quality metrics
Mean squared error (MSE) and mean square deviation (MSD) are used for describing the difference between the estimator and what is estimated. Similarly, root-MSE (RMSE) and root-MSD (RMSD) can be also used for evaluating the quality of reconstruction images. Signal-to-noise ratio (SNR) is used as a physical measure of the sensitivity of an imaging system. The signal amplitude of each pixel is the amount in which patch of the image is elevated, relative to the mean background signal. SNR represents the integrated signal over a region of interest (ROI). Structural similarity (SSIM) [41] is an index for measuring the similarity between two images. Mean structural similarity (MSSIM) is an average of the SSIM used to measure the similarity between two images in terms of brightness, contrast, and structure. Pearson correlation coefficient is to measure the linear correlation between two images. Universal quality index (UQI) [42] is used to evaluate the similarity between the reconstructed image of the ROI and the phantom image. The value of UQI ranges from 0 to 1. If the UQI value is closer to 1, the reconstruction image is closer to the real image.
CT image reconstruction open source toolkits
CT image reconstruction open source toolkits are investigated from the introduction websites, corresponding publications and documents in software packages. The survey results are briefly summarized in Table 1, carried out in alphabetical order. The research content of each toolkit is mainly carried out from the following items corresponding to the following notes:
Survey results of 15 software from item (1) to (11)
Survey results of 15 software from item (1) to (11)
Note: “*” indicates that the survey results of corresponding toolkits are unknow. “/” indicates that the toolkit does not support corresponding function. “divonx;X” indicates that the toolkit relies on the X toolkit to support related function.
Object model: (A) Analytical model, (B) Discrete model: (B1) Pixel, (B2) Voxel, (B3) Blob; Projection geometry: (A) Circle trajectory, (B) Helical trajectory, (C) Arbitrary trajectory; (D) Curved-detector, (E) Line-detector, (F) Multi-row curved-detector, (G) Planar-detector; Projection model: (A) Analytical method, (B) Pixel driven, (C) Ray driven, (D) Strip model; (E) GPU support; (F) Noise simulation, (G) Scattering simulation, (H) Artifacts simulation; Projection data pre-processing model: (A) Denoise, (B) Dead pixel correction, (C) Detector consistency correction, (D) Artifact correction, (E) Beam hardening correction, (F) FOV correction, (G) Phase retrieval, (H) Flat field correction; Backprojection model: (A) Pixel driven, (B) Ray driven, (C) Strip model, (D) GPU support; Post-processing: (A) Denoise, (B) Threshold segmentation, (C) 3D cutting, (D) CT number, (E) Quantitative analysis, (F) Artifacts removal; Special function: (A) Quarter pixel offset, (B) Simulating polychromatic X-ray, (C) Support arbitrary scanning axis, (D) Scattering correction, (E) Respiratory motion correction, (F) Leverage nonlinear interval data acquisition, (G) Limited angle problem; Features other than CT: (A) PET reconstruction, (B) SPECT reconstruction, (C) MRI reconstruction; Image quality metrics: (A) Standardized root mean square distance measurement, (B) Normalized average absolute distance measure, (C) the worst case distance measure over a 2×2 pixel area, (D) Structure accuracy, (E) Pointwise accuracy, (F) Hit-ratio, (G) Pearson correlation coefficient, (H) RMSE, (I) MSSIM, (J) UQI. Reconstruction algorithms: Analytical algorithms (FBP, FDK, etc.), Iterative algorithms, etc.; Software environment: Operation system (Windows, Linux), Compiling languages (C, C++, MATLAB, Python, etc.), GUI (Graphical User Interface), software structure. Year: The year for the latest version of each toolkit.
All scale tomographic reconstruction Antwerp (ASTRA) toolbox [43] provides an efficient and flexible open source codes for CT image reconstruction. From 2010 to 2016, the toolbox was developed and maintained by the iMinds-Vision Lab in Antwerp University. In 2014, the CWI institute in Netherlands joined in the ASTRA toolbox developing group.
ASTRA toolbox supports Windows and Linux operating systems, and compiling languages include C++, MATLAB and Python. Three main layers are involved in the toolbox [44]: (1) Efficient algorithm building module is provided in the bottom layer, for instance, projectors and back-projectors with GPU hardware accelerator; (2) In the layer 2, there are reconstruction algorithms written in C++; (3) MATLAB interface and Python interface are provided in the top layer for end-users.
ASTRA toolbox is used for three-dimensional reconstruction, supporting parallel beam geometry and cone beam geometry. It also supports fan beam geometry with line-detector. One of the features of ASTRA toolbox is that it supports arbitrary trajectory for projection geometry [45]. In ASTRA toolbox, projection geometry is defined by the coordinate of the source, the detector center, and the principal axis of the detector plane (typically horizontal and vertical, specified as two 3D vectors). The magnitude of each vector corresponds to the pixel value of every detector.
The ASTRA toolbox supports ray-driven and strip model for computing projection data, including the following three projection operators: “line”: Calculate the length of the ray through a pixel as the contribution; “linear”: Linear interpolation is performed between the two nearest voxels of the ray and the column/row intersection as the contribution of the column/row to the ray; “strip”: The area surrounded by two adjacent rays passing through a pixel is calculated as the contribution of the pixel.
In ASTRA toolbox, parallel beam supports the above three projectors, fan beam supports “line” and “strip” projectors, and cone beam only supports “line” projector. The toolbox supports Poisson noise simulation. ASTRA toolbox also includes the following algorithms: Two-dimensional reconstruction algorithms include FBP, SIRT, ART, SART, CGLS, EM; Three-dimensional reconstruction algorithms include FDK, SIRT, CGLS. FBP, ART, and SART support parallel beam and fan beam geometries. SIRT and CGLS support parallel, fan beam, and all 3D geometries.
Filters in analytical reconstruction algorithms include ram-lak, shepp-logan, cosine, hamming, hann, tukey, lanczos, triangular, gaussian, barlett-hann, blackman, nuttall, blackman-harris, blackman-nuttall, flat-top, kaiser, parzen.
The toolbox includes a discrete algebraic reconstruction technique (DART) that can segment the reconstructed image by setting a threshold. DART can be applied to dense nanoparticles segmentation problem.
CCTT
Copenhagen computed tomography toolbox (CCTT) is a CT image reconstruction toolbox based on GPU acceleration that contains reconstruction application programs [46]. CCTT runs on Windows or Linux. The toolkit was written in Python and implemented by Numpy, PyCUDA and PyOpenCL. In advanced versions, all intensive computing is done in the PyCUDA or PyOpenCL kernel [47]. Thus, the core Python and NumPy functions are only used for convenient scripting and basic wrapping. Researchers aim to offer an open CT reconstruction toolbox that can be applied to all possible CT system, including education, research, and industries. The toolbox should be easy, flexible, and fast to achieve this goal [48].
Data pre-processing and post-processing can be done efficiently by the relevant GPU plugins in the toolbox. The toolbox provides a flexible framework and a plugin-in structure to make it easy to extend the toolbox with new reconstruction algorithms.
CCTT supports helical cone beam geometry and circle trajectory volume geometry.
CCTT supports detector consistency correction, dead pixel correction and flat field correction. CCTT contains the following reconstruction algorithms: Center FDK and FDK: Center FDK utilizes FDK algorithm for cone beam reconstruction in central fan beam projecting. Filters include “hamming”, “ram-lak”, “shepp-logan”, “cosine”, “hann”, and “skip”. Katsevich: Katsevich algorithm is a cone beam reconstruction algorithm applied to helical CT.
CIL
The core imaging library (CIL) [49] is a set of modules in the CT datasets analysis workflow. The library is to provide imaging tools for CT imaging community. Those tools in CIL can be integrated into existing workflows, such as SAVU (a Python package assisting with the processing and reconstruction of parallel-beam tomography data [50]). CIL was written in C++, Python, and MATLAB. CIL includes algorithms for pre-processing module, reconstruction module, data analysis module, and visualization module. The algorithms are provided by CT imaging community. The codes are redesigned by the core CCPi [51] staff.
The pre-processing module [52] provides beam hardening removal algorithms for X-ray tomography [53]. This module can also be applied to limited angle CT problem. In the reconstruction module, CIL provides framework that aims to give a unified platform to use, to distribute, and to document. The framework provides basic converters to read data that is generated from different scanners. Currently, it supports NeXus format and Nikon format. With the framework and some plugins [54], such as CCPi plugins and ASTRA plugins, each module is independent that the developers can use these modules independently to analyze their data.
The reconstruction module also provides some iterative algorithms [55]: CGLS, ML-EM, SIRT. Besides, three extensive algorithms of CGLS algorithms are included: CGLS based on convolution, CGLS based on Tikhonov regularization, and CGLS based on TV method. The reconstruction module includes a CCPi-Regularisation toolkit (CCPi-RGL) [56, 57] for image denoising and inpainting problems. CCPi-RGL supports GPU hardware for acceleration. To guarantee a better performance in SNR and resolution for the ill-posed problems, CCPi-RGL provides 2D/3D regularization strategies that can be used with proximal splitting methods [58], such as the primal-dual hybrid gradient (PDHG) algorithm, Douglas-Rachford splitting algorithm, ADMM, FISTA.
CIL provides algorithms for suppression of ring artifacts. The data analysis module provides quantification algorithms for calculating characteristics and accessible volume in a range of sphere sizes (3D image). CIL also provides a segmentation algorithm based on contour tree [59].
CTSim
CTSim [60] simulates the process of the X-ray going through virtual object and utilizes all kinds of algorithms to reconstruct the original image from the projection. It has a wide array of image analysis and image processing functions. CTSim can be tested in Windows and Linux operating systems. The toolkit was written in C++ and equipped with GUI by wxWindows graphical toolbox. CT Sim also provides analytical models that simulate the scanned object model in combination with some simple geometric objects. A phantom is usually made up of multiple elements: rectangle, triangle, ellipse, sector and segment. Besides, Herman phantom, Shepp-Logan phantom, and Unit Pulse phantom are also given by CTSim.
CTSim can simulate single circle scanning trajectory. The software supports line-detector and curved-detector fan beam geometries. CTSim calculates the projection and backprojection by linear interpolation. Operations of reconstructing include Fourier inverse transformation, filtering backprojection. As for the filters, Shepp-Logan, sinc, Hamming, Hanning, Cosine, Triangle and Bandlimit are included in CTSim. In addition, the software provided three image quality metrics for analyzing: Standardized root mean square distance measurement, normalized average absolute distance measure, and the worst case distance measure over a 2×2 pixel area.
MIRT
The Michigan image reconstruction toolbox (MIRT) is a collection of open source algorithms for image reconstruction (and related imaging problems) [61]. The toolbox was written in Mathwork’s MATLAB. It can be tested in Windows and Linux operating systems. Besides Shepp-Logan phantom, MIRT contains voxel discrete model and blob discrete model for simulation. MIRT can be used for parallel beam geometry and cone beam geometry. MIRT supports line-detector or curved-detector for data generation. The toolbox can be used for circular trajectory reconstruction and helical CT. MIRT also can build models for planar detector and multi-row curved-detector cone beam image reconstruction. MIRT provides linear projectors for calculating projection data, as well as noise simulation and scattering simulation.
The toolbox supports quarter-detector shift resulting in a duplication of the spatial resolution in a 360 degree rotation of the sampling unit. Nonuniform sampling can also be carried out by setting parameters. Moreover, for the data processing of non-uniform sampling, MIRT software gives a nonuniform fast Fourier transform (NUFFT) toolbox [62]. Algorithms in MIRT are mainly applied to polychromatic CT, including FBP and WLS. The toolbox provides some filters: “hann”, “ramp” and “hamming”. The toolbox also includes the following imaging functions other than CT: Iterative and non-iterative algorithms for tomographic imaging (PET, SPECT, X-ray CT); Methods for magnetic resonance (MR) image reconstruction, including compensation for off-resonance effects (field inhomogeneity, susceptibility, etc.); Methods for MR RF pulse design, including MRI spectral-spatial pulse design for phase-precompensatory slice selection. Iterative image restoration tools; Methods for B-spline based image registration with regularization to encourage the deformation to be invertible (diffeomorphic).
OpenRecon
OpenRecon is an open source toolkit for image reconstruction [63]. The toolkit was developed for sharing a standard platform equipped with key analytic reconstruction methods for different scanning modes and imaging geometries. It is implemented by C++ and MATLAB. While the intensive computational parts are implemented by C++ dynamic link libraries, MATLAB serves as the basic interface and calls the dynamic link functions directly.
OpenRecon provides analytical object models. The input of the reconstruction algorithms can be either from numerical simulations or real CT scanners. It supports helical/spiral multi-slice/cone-beam scanning modes. The toolkit can be used for flat-panel detector and multi-row curved-detector imaging geometries. It also provides pixel-driven in the backprojection. OpenRecon contains analytical reconstruction algorithms (FDK) and numerical projection generations.
OSCaR
Open source cone-beam reconstructions (OSCaR) [64] is a cone beam image reconstruction package developed by American Association for Physicists in Medicine (AAPM). The package was implemented in MATLAB and the GUI was written by MATLAB. Researchers are aimed to provide and implement an intuitive, open source GUI image reconstruction software with flexibility, universality, and ease of use.
OSCaR is used for cone beam geometry following circle scanning trajectory with the detector panel perpendicular to the trajectory plane. OSCaR software supports detector consistency correction by adjusting the offset of the center of the detector panel, the pitch angle and rotation angle of the detector panel. In summary, OSCaR is mainly used for FDK cone beam reconstruction and the filters include Ram-Lak, Shepp-Logan, Cosine, Hamming and Hann [65].
RTK
Reconstruction toolkit (RTK) [66] is an open source reconstruction software based on insight toolkit (ITK). The RTK Toolbox supports GPU acceleration and also provides a set of command line applications defined by Gengotop, a tool that generates command line option resolution codes for C programs [67]. RTK uses CUDA to implement several time-consuming filters. Part of the codes are taken from Plastimatch (an open source toolkit for image computation [68]) and NiftyRec (an open source toolkit for emission tomography [69]). There is also a filter implemented by OpenCL in RTK.
Currently, RTK supports only cone beam geometry with circle scanning trajectory. Deviations from this strict geometry can be described for each projection image by means of nine degrees of freedom in total: 3 coordinates for the position of the punctual source, 3 coordinates for the flat panel position and 3 angles for the flat panel orientation. The software provides Gaussian noise simulation. RTK software provides many filters for denoising, including conditional median image filter, average out of ROI image filter, TV denoising image filter, conjugate gradient image filters, Laplace filters and additive Gaussian noise filter. In addition, RTK software contains FDK and SART reconstruction methods. RTK software provides tools for respiratory motion correction and preprocessing of raw data for scatter correction.
SNARK14
SNARK14 is a programming system for the reconstruction of 2D image from 1D projections. SNARK14 is the updated version of SNARK09, including many features of SNARK09. SNARK14 toolkit was developed to provide a unified framework to implement the image reconstruction algorithms and evaluate the performance for researchers. SNARK14 was written by C++ with GUI and can be tested in Linux operation system. The developers provide SNARK14 of VirtualBox version and the installation tutorial in the website [70].
In SNARK14, the object models include analytical models combined by several fundamental elements, including rectangles, ellipses, triangles, segments and sectors. Each element is defined by its center coordinate, lengths of two perpendicular directions, orientation angle and density (in case of CT, density represents the linear attenuation coefficient). The densities of overlapping areas are added together. Object models also include pixel model and blob model. SNARK14 supports parallel beam and divergent beam scanning geometry of circle trajectory for 2D reconstruction. The toolkit uses ray-driven based on the length of intersection of rays and pixels, strip model and analytical method for calculating the projection data. The software provides noise simulation and scattering simulation. SNARK14 also supports the beam hardening correction.
SNARK14 also has built-in image reconstruction algorithms and it provides users with custom options to create their own reconstruction algorithms and tests. Reconstruction methods are usually divided into two categories: transformation methods and series expansion methods. Transformation methods include FBP, Rho-filtered layergram and other various deblurring methods that can be used with this algorithm, Fourier methods and Linogram based on projection theorem; Series expansion methods include ART, SIRT, Quadratic optimization techniques, EMAP.
The series expansion methods can be used in both pixel model and blob model, the transformation methods can be only used in pixel model reconstruction. Besides built-in algorithms, users can add in 10 user-defined algorithms in SNARK09. The major advances in SNARK14 are the capability of optimizing iterative algorithms to improve performance, the capability of applying SART to projection data, added several new options to the statistical evaluation and added several new stopping criteria.
For single reconstructions, SNARK09 [71] provides means for the evaluation of some quantitative measures of the overall difference between a digitized test phantom and its reconstruction. Such an evaluation can be performed either over the entire region of the image or over selected areas, and can be also restricted to pixels whose densities fall within a user-selected range. On this basis, the updated version of SNARK14 adds several statistical evaluation indicators and introduces them in the SNARK14 usage manual [72]. The evaluation indicators include structural accuracy, hit-ratio, and pointwise accuracy.
SNARK14 can customize the number of rays while scanning and the spacing of the detectors can be defined as uniform or non-uniform. It can also simulate polychromatic X-rays for PET simulation.
SPRAY
Spectral X-ray image reconstruction (SPRAY) toolbox contains many functions written by MATLAB. SPRAY is developed for spectral X-ray nonlinear image reconstruction [73]. SPRAY uses the data measured by different energy levels to reconstruct the three-dimensional entity. It mainly solves two problems: material decomposition problem and CT image reconstruction problem. SPRAY implements the regularized weighted least squares Gauss-Newton algorithm (RWLS-GN) [74] and provides several different regularization methods to process. The algorithm is less sensitive to noise and improved contrast-to-noise ratio of the gadolinium image, it efficiently solves the nonlinear decomposition problem for spectral CT, which produces new possibilities. In addition, SPRAY also has a built-in thorax model.
TIGRE
Tomographic iterative GPU-based reconstruction (TIGRE) [75, 76] is a MATLAB and Python/CUDA toolbox for fast and accurate 3D tomographic reconstruction. TIGRE toolbox operates in Window or Linux system. The main construction module of the toolbox is projection and back-projection operators at the bottom layer, which are optimized for GPU acceleration; the algorithm module is on the top layer written by MATLAB; in order to communicate between lower layer, hardware-oriented CUDA and higher layer, design-oriented MATLAB, a set of the so-called MEX functions are needed. Moreover, Python version is also provided in the TIGRE toolbox which is still under developed.
TIGRE toolbox contains analytical models. Analytical models include Shepp-Logan phantom and modified Shepp-Logan phantoms. There are two modified versions: One is modified in the contrast, the other is combined with “Katsevich-Type algorithm of variable radius spiral cone beam CT”. The toolkit contains a built-in thorax model.
TIGRE toolbox supports circle and helical trajectory. TIGRE toolbox provides three-dimensional projection data calculation methods, i.e., ray-driven method: one is to calculate the length of the intersection of rays and voxels, the other is to perform trilinear interpolation of sampling points on rays. These two calculations have been optimized for GPU computing using CUDA. The toolbox also has Poisson noise simulation and Gaussian noise simulation.
TIGRE toolbox provides detector consistency correction. Moreover, system detectors rotate slightly due to mechanical inaccuracies, in order to address this problem, TIGRE toolbox uses three rotating angles (spinning rotation angle, self-offset rotation angle and self-offset pitch angle) of detector panel for correction.
One of the most important features of TIGRE toolbox is that it contains many kinds of algorithms and extension algorithms: Analytic algorithms: FDK; Gradient descent algorithms: SART, SIRT; OS-SART; Krylov sub-space algorithms: CGLS; Statistical algorithms: MLEM; TV minimization: the projection onto convex sets (POCS) and the extension algorithms of POCS, SART-TV, Projection-controlled steepest descent (PCSD) and the extension of PCSD.
TIGRE toolbox provides filters: “ram-lak”, “shepp-logan”, “cosine”, “hann”. The toolbox implements the pixel-driven for back-projection. However, considering that in some iterative algorithms the forward projection may not match the back-projection calculation, the toolbox also implements the ray-driven back-projection calculation.
Two post-processing are available in TIGRE toolbox. One is denoising based on the TV of reconstructed image, the other is segmentation to the reconstructed volume (3D cutting). There are mainly four image quality metrics available in TIGRE for evaluation: Pearson correlation coefficient, RMSE, MSSIM, and UQI.
Since the scanners of modern CT are more complex, it is common that the rotation axis is arbitrary, which means the trajectory of the detectors and the source may be not on the same circle panel but on a “spherical” trajectory. The toolbox simulates the rotation axis of the devices by defining Euler angles, which are a set of independent angle parameters for uniquely determining the position of a fixed point rotating object.
Tomo series
TomoPy was developed by Argonne National Laboratory. TomoPhantom, ToMoBAR and Multi-Channel X-Ray CT package were developed and maintained by Daniil Kazantsev and his group in Manchester University. These four packages are feature connected and supported. In this part, they will be put together to introduce.
TomoPy is an open source package of Python. Most of the reconstruction algorithms in TomoPy were written in C compiling language for tomographic data processing and image reconstruction [77]. TomoPy is a framework for analyzing synchron tomographic data processing, dealing with image correction and reconstruction. The data analysis module is mainly composed of pre-processing module, reconstruction module and post-processing module [78]. Python front-end is used to connect these modules to the interface for user interaction. Any modules of the package can be connected to other language software separately.
TomoPhantom [79, 80] is a package that generates customizable 2D–4D phantoms (with a temporal capability). It provides a series of tools for testing reconstruction algorithms and simulating accurate projection. The core module of the software package is written by C-OpenMP and provides Python and MATLAB package interfaces, which are accessible by Cython and C-MEX respectively. The package relies on the ToMoBAR, ASTRA toolbox and TomoPy toolbox for reconstruction.
Tomographic model-based reconstruction (ToMoBAR) [81] is a library of direct and model-based regularized iterative reconstruction algorithms with a plug-and-play capability. The package uses TomoPhantom for simulations and ASTRA toolbox for projection operations. It also relies on the CCPi-Regularisation Toolkit for reconstruction. ToMoBAR can be used in MATLAB and Python.
The Multi-Channel X-ray CT package is used to provide numerical results for a new multi-channel reconstruction method for X-ray spectral CT proposed. The package relies on the ASTRA toolbox, compiled in C language and operate on the interface of MATLAB. It supports GPU acceleration to operate on the MATLAB interface.
TomoPhantom package provides a simple way to build 2D-4D analytical models, which are composed of piecewise constant, piecewise smooth and smooth analytic objects. TomoPhantom package contains the following scanned object models: Generate 2D models made of Gaussians, parabolas, ellipses, cones, rectangles; Generate 3D models and 4D (temporal) extensions of them; Calculate analytical Radon transform of 2D–4D models and their numerical projections.
TomoPy builds a framework for the analysis of synchrotron tomographic data. Synchronous accelerator data acquisition mode is to keep the position of detector and the source fixed, the object rotating in the middle. ToMoBAR uses ASTRA toolbox to define projection geometry and uses TomoPhantom to simulate for projection data. Multi-Channel X-ray CT also depends on ASTRA toolbox to define arbitrary projection geometry for spectral CT data.
TomoPhantom and ToMoBAR can simulate Gaussian noise and Poisson noise simulations. Besides, both of these packages provide the artifact simulations including some classical acquisition artifacts: Zingers, rings, shifts. The data generation script of Multi-Channel X-Ray CT software depends on the simulation of Photon Attenuation toolbox and Spektr software [82].
In the pre-processing part, TomoPy provides normalization, data padding, drift correction, zinger removal, FOV correction, strip removal and phase retrieval. TomoPy provides wavelet-Fourier filters to deal with stripes in projected sinusoidal images to remove annular artifacts from reconstructed images. TomoPy also provides Paganin-type and Bronnikov-type filters for single-distance phase restoration.
TomoPy provides many reconstruction methods for reconstruction module, including Gridrec, ART variant algorithms and MLEM variant algorithms. Taking into account the characteristics of synchronous accelerator scanners, TomoPy makes use of the characteristic of rotating center offset to generate artifacts, that is, the rotation center position corresponds to the image quality, which transforms the image reconstruction problem into the optimization problem of image quality, that is, the Centering reconstruction method in reconstruction module.
In ToMoBAR, simulated data reconstructed iteratively using FISTA or ADMM algorithms with multiple “plug-and-play” regularisers from CCPi-RegularisationToolkit. Researchers modified FISTA algorithm: convergence acceleration with ordered-subsets method, PWLS (penalized WLS), Huber, Group-Huber [83] and Student’ data fidelities [84, 85] to deal with noise and imaging artifacts. It also provides real data to demonstrate the reconstruction of the composite image. TomoPhantom uses ToMoBAR, ASTRA toolbox and TomoPy in the reconstruction module.
Multi-Channel X-ray CT package depends on CCPi-Regularisation Toolkit and Spot operator package (A linear-operator toolbox for MATLAB [86]) in the reconstruction module. It uses FISTA multi-channel reconstruction method [87]. The post-processing module in TomoPy contains segmentation and quantitative analysis.
Deep learning for CT image reconstruction
Over the past years, deep learning (DL) has been a hot topic and made a great difference in different research fields. In the perspective of imaging research, DL techniques are being actively developed for computed tomographic image reconstruction [88].
Development of DL
The development of modeling in image reconstruction was summarized to three stages by Zhang and Dong [89], starting from handcrafting modeling to DL modeling with a transition of handcraft plus data-driven modeling. The handcraft models may not be flexible and cannot fully deal with large data sets, however, these models are based on mathematical characterization with solid theoretical supports. Besides, DL models are in lack of theoretical foundations but much more effective and able to extract useful information from large data sets. It is concluded that combining the handcraft models with DL models is a major research trend in image reconstruction. However, there are still some challenges in the research. These challenges have been discussed in detail by Wang et al [90]. These challenges are also opportunities for researchers.
Techniques based on DL
Currently, many DL reconstruction techniques were developed for some specific problems. For instance, Shan et al [91] brought up a modularized neural network for low dose CT that performs an end-to-process mapping for denoising. Their study proves that the DL imaging approach has a better performance in denoising and structural fidelity. For low-dose CT imaging, Chen et al [92] combined the autoencoder, deconvolution network, and shortcut connections into the residual encoder-decoder convolutional neural network (RED-CNN), which achieved a good performance in noise suppression, structural preservation, and lesion detection. Chen et al [93] also brought up a deep convolutional neural network that was used to map low-dose CT images towards its corresponding normal-dose counterparts in a patch-by-patch fashion. The new noise reduction method demonstrated a great potential on artifact reduction and structure preservation and a better performance in computation time. Besides, Chen et al [94] constructed a learned experts’ assessment-based reconstruction network (LEARN) that was demonstrated the feasibility and merits in terms of artifact reduction, feature preservation, and computational speed, for sparse-data CT.
Wang [95] provides a comprehensive overview of DL-based tomographic imaging reconstruction techniques, including the concepts and components of representative deep neural networks, some open source toolkits, datasets for deep learning, as well as some hands-on neural network models for deep reconstruction. It can be also referred to [96]. In that paper, Shan et al reviewed deep learning methods for medical imaging that included CT image reconstruction, radiation therapy ranging from planning and verification to prediction, and the connection between them, based on the analysis of neural networks.
Open datasets
It is a challenging task to compare DL reconstruction approaches since these networks highly rely on the data that is used for training. Following are some freely-available CT datasets.
Low dose CT grand challenge
The low dose CT grand challenge [97] was to quantitatively evaluate the performance of denoising and iterative reconstruction techniques on common low-dose patient CT datasets. Projection and image data for a helical scan of the American College of Radiology (ACR) CT Accreditation [98] Phantom were provided. Although it is now complete, the phantom and patient training data are still being made available to interested researchers.
CQ500 dataset
To develop and validate a set of deep learning algorithms for detection of critical findings in head CT scans, the CQ500 dataset [99] of 491 scans with 193,317 slices are publicly available so that others can compare and build upon the results Chilamkurthy et al have achieved in the paper [100].
LIDC-IDRI
The lung image database consortium image collection (LIDC-IDRI) [101] consists of diagnostic and lung cancer screening thoracic CT scans with marked-up annotated lesions. It is a web-accessible international resource for development, training, and evaluation of computer-assisted diagnostic (CAD) methods for lung cancer detection and diagnosis.
The LoDoPaB-CT Dataset
The LoDoPaB-CT Dataset [102, 103] is a benchmark dataset for low-dose CT reconstruction methods. The dataset allows a fair comparison while training deep learning approaches. It contains over 40,000 scan slices from around 800 patients selected from the LIDC/IDRI Database.
ALERT Datasets
The purpose of ALERT Datasets [104] is to provide security-like data to academic researchers, to discover and evaluate the present state-of-art information, and to stimulate additional communication and research in the algorithm research community.
SPARE Challenge
The goal of the SPARE challenge [105] is to explore the possibility of high quality 4D CBCT while sparing patients the additional scan time and imaging dose that is currently used in the clinics. There are three main datasets in this challenge: Monte Carlo dataset. There are three different types of simulated projections: no scatter, scatter, and low dose. This allows the comparisons of both reconstruction algorithms and scatter correction. Clinical Varian dataset. This dataset contains oversampled clinical CBCT projections acquired from a Varian system using half-fan geometry. Clinical Elekta dataset. This dataset contains oversampled clinical CBCT projections acquired from an Elekta system using full-fan geometry.
Open X-ray Tomographic Datasets
The Open X-ray Tomography Datasets [106] include: Tomographic X-ray data of a lotus root slice filled with different chemical elements is an open-access dataset of tomographic X-ray data of a slice of a lotus root. The dataset consists of the X-ray sinogram of a single 2D slice of the lotus root with two different resolutions, and the corresponding measurement matrices modelling the linear operation of the X-ray transform. Tomographic X-ray data of a walnut consists of the X-ray sinogram of a single 2D slice of the walnut with three different resolutions and the corresponding measurement matrices modelling the linear operation of the X-ray transform. Each of these sinograms was obtained from a measured 120-projection fan-beam sinogram by down-sampling and taking logarithms. Tomographic X-ray data of a carved cheese consists of the X-ray sinogram of a single 2D slice of the cheese slice with three different resolutions, and the corresponding measurement matrices modelling the linear operation of the X-ray transform. Tomographic X-ray data of emoji phantom consists of the X-ray sinogram of a single 2D slice of 33 emoji faces (contains 15 different emoji faces) made by small squared ceramic stones, and the corresponding static and dynamic measurement matrices modelling the linear operation of the X-ray transform. Each of these sinograms was obtained from a measured 60-projection fan-beam sinogram by down-sampling and taking logarithms. Tomographic X-ray data of a 3D cross phantom. The dataset consists of the X-ray sinogram with 16- or 30-time frames (depending on the resolution) of 2D slices of the phantom, and the static and dynamic measurement matrices modelling the linear operation of the X-ray transform. Each of these sinograms was obtained from a measured 360-projection fan-beam sinogram by down-sampling and taking logarithms. Tomographic X-ray data of a dynamic agarose-gel phantom perfused with a liquid contrast agent consists of 17 consecutive X-ray measurements of a 2D slice of the gel phantom. This data is provided in two resolutions and stored in a special CT-data structure containing the metadata of the measurement setup. In addition, a denser projection sampling of the first-time step, and a very dense sampling of an addition 18th time step is provided as a complementary data.
Discussions and conclusion
In this review paper, fifteen open source reconstruction toolkits are presented from operating environment, CT system models, and algorithms. Based on the characteristics of each toolkit, considering the requirements of software development and research interaction, we will conclude from the following aspects: Some toolkits support GPU versions for calculating projection data. Based on hardware accelerator, reconstruction time can be greatly reduced, especially for iterative algorithms. This is essential for software development. The emergence of GPU acceleration technology creates more possibilities for iterative reconstruction techniques. Algorithms can be combined with current research frontier, including leveraging the noise and artifacts in real projection data. As artificial intelligence techniques have been actively developed for tomography reconstruction, especially DL techniques, it will be a trend to apply DL for image reconstruction. Many commercial toolkits have graphical interface. These toolkits provide great convenience for users to operate. Some toolkits provide open and flexible platforms that can be extended to current or future reconstruction algorithms. In addition, modularization connects features among different toolkits. This helps build a reconstruction system with better compatibility and stronger function.
With the development of CT technology, more and more CT image reconstruction algorithms and toolkits for different applications will be brought up, like multi-source and multi-detector structure image reconstruction algorithms, spectral CT image reconstruction algorithms. Considering low dose data, incomplete problem, and artifact correction, reconstruction methods based on prior information (regularization methods, DL techniques) need to be improved in toolkits. For example, under the support of NCI, GE global research, UMass Lowell and Duke University are developing an open-access simulation environment for X-ray based imaging for cancer [107, 108]. This simulation environment, the X-ray-based cancer imaging simulation toolkit (XCIST), will be physics-based, written in Python and C/C++, and will consist of four major subsets: digital patient phantoms with and without tumors, the simulator itself, image reconstruction algorithms, and a radiation dose estimation tool. To enable broad usage, XCIST will include a GUI similar to that of a scanner’s console, will be well documented and easy to use, and will be easily extendable by researchers who wish to model and evaluate new X-ray based imaging technologies.
In this review paper, we try to list CT image reconstruction toolkits as much as possible, but omissions cannot be avoided. We hope our work can be useful for researchers. Readers can choose whichever toolkits they need as a platform based on this review. In this way, researchers can conduct their research from a higher starting point.
Footnotes
Acknowledgments
Liu Shi was sponsored by National Key R&D Program of China No. 2017YFF0107201. Baodong Liu and Cunfeng Wei were sponsored by CAS Interdisciplinary Innovation Team, Project No. JCTD-2019-02. Long Wei was sponsored by the External Cooperation Program of BIC, Chinese Academy of Sciences No. 113111KYSB20150037. Li Zeng was sponsored by the National Natural Science Foundation of China No.61771003.