1. Introduction

The use of X-ray Computed Tomography (CT) is not limited to its most well-know application in medical diagnostics. Indeed, its ability to look at the interior of certain objects benefits many other application fields as well. Examples include biomedical research, that often relies on micro-CT scanning setups to investigate the effect of certain substances or medicines on small animals [1]; industrial applications, where tomographic reconstruction can be an invaluable tool for quality control of (semi-)automatic processes on a conveyor belt [2, 3]; chemistry, where it can be used to study and optimize reactions and their results (e.g., of polyurethane foams); and material science with various imaging modalities such as micro-CT, synchrotron radiation facilities, and electron microscopes [4]. While all of these applications share the same basic concept (multiple projection images are measured from different directions), the specific setup of the scanning system can differ from application to application.

Typically, reconstruction software comes bundled with the purchase of the scanner (e.g. NRecon for Bruker micro-CT systems, Phoenix datos|x for GE systems, and Inspect3D for FEI electron microscopes). A few commercial products are also available that are not tightly coupled to a specific scanning system, but can come at a premium rate (e.g., VG Studio Max, Octopus Reconstruction, Ikeda’s TomoShop). The main advantages of these commercial packages are their reliability due to their large user-base and their user-friendly graphical user interfaces. For the majority of CT users these tools will therefore be well sufficient. However, the implementations provided by these packages offer only a limited flexibility, restricting some use cases:

In Table 1, a short overview is given of different open source software packages for tomographic reconstruction that have recently been updated: pyHST2 [11], Ettention [12], and IRT [13]. While each of these packages provide excellent implementations of the most common tomographic routines, it can be noted that each package is focused on a specific target application. To the authors knowledge, there is currently no software package available that is more suited for the development of advanced algorithms that can be used in multiple application fields, using non-conventional geometries and various constraints.

Table 1. Overview of different open source toolboxes for tomographic reconstruction.

View Table | View all tables in this article

The ASTRA Toolbox (All Scales Tomographic Reconstruction Antwerp) introduced in this article is a software platform to address the need for a flexible, efficient, and easy to use development platform for tomographic algorithms. Since 2010, it has been developed by the Vision Lab, University of Antwerp, Belgium, and since 2014 jointly by Vision Lab and the Centrum Wiskunde & Informatica (CWI) in Amsterdam, The Netherlands. It is freely available as open source software under a GPLv3 license [14]. To achieve optimal efficiency, most of its code is written in C++, with its core computations offloaded to a GPU card using the CUDA language. Users with no access to GPU workstations can resort to CPU code, accelerated using OpenMP. The toolbox is accessible through a MATLAB and a Python interface, providing a fast and powerful platform for algorithm and application prototyping.

The ASTRA Toolbox is mainly intended for use in research (e.g., application prototyping and algorithm development), and less so for day to day usage by end-users. However, due to its open nature, the ASTRA Toolbox can be coupled with existing practical software packages. One example is its inclusion into TomoPy [15], which combines the user-friendliness and efficient pre-and postprocessing methods of TomoPy with the flexibility and optimized reconstruction methods of the ASTRA toolbox into a single framework. Also, for applications in the field of electron tomography, the ASTRA Toolbox is available through FEI’s Inspect3D software.

The main strength of the ASTRA Toolbox is its inherent flexibility, which comes in three forms:

This article provides an overview of the ASTRA Toolbox, its design, and some of its use-cases. Firstly, in Section (2), the key components of the toolbox are introduced. Then, Section (3) demonstrates how they can be used in a conventional cone beam setup. The subsequent two sections then delve deeper into some key advanced features of the toolbox that allow it to be used flexibly: Section (4) concerns the specification of non-conventional scanning setups, and Section (5) discusses how the ASTRA Toolbox can be easily integrated into external codebases. Finally, Section (6) concludes this work.

2. Framework concepts

There are seven different concepts and modules in the ASTRA Toolbox (Fig. 1). Volume and projection data objects (Section (2.1)) are linked to their corresponding geometry information (Section (2.2)), both of which are linked together in a projection algorithm building block (Section (2.3)), which is subsequently used together with the data in the reconstruction algorithm object (Section (2.4)).

 

Fig. 1 Schematic overview of all ASTRA Toolbox concepts.

Download Full Size | PDF

2.1. Projection and volume data

The data modules of the ASTRA Toolbox store and manage all datasets, both for volume data and for projection data. In order to perform actions on this data (e.g., a projection or a reconstruction) the toolbox first needs to have access to it, which is not the case if the data is stored only in the MATLAB or Python environment. Thus, memory first has to be allocated and data has to be copied from the interface layer to the actual toolbox. This operation then returns a unique identifier that can subsequently be used to refer to that particular set of data, much like how file input/output works in many programming languages.

2.2. Geometry

The experimental setup of the scanning system is specified by linking the data objects to a corresponding volume geometry or projection geometry. The volume geometry describes the pixel or voxel grid on which the reconstructed object is represented. A voxel is typically isotropic but can be specified to anisotropic as well. This is useful in cases where the effective resolution in a reconstruction is not the same in each dimension (e.g., tomosynthesis and laminography).

The projection geometry defines the trajectory of the X-ray source and detector relative to the volume geometry. This involves the number and size of the detectors pixels and the specification of the trajectory of the X-ray source and detector as they rotate around the object. Standard projection geometries that are available are parallel beam and fan beam for two-dimensional data, and parallel beam and cone beam for three-dimensional data.

Specifying a projection geometry in the ASTRA Toolbox can be done in two ways. In its default form, the source and detector follow a circular trajectory around the origin of the volume coordinate system and the detector plane lies perpendicular (and perfectly aligned) to the line through the source and the origin. Hence, for a circular cone-beam geometry, there are five sets of parameters to be specified: the number of detector pixels, the detector pixel size, the distance from the source to the centre of rotation, the distance from the centre of rotation to the centre of the detector and a list of projection directions (Fig. 2(a)). With this approach, the most widely used projection setups can be specified, but non-conventional, non-circular trajectories can not. Moreover, it is not possible to model common small detector misalignments, such as a detector tilt or detector shift (when the centre of the detector does not line up with line through the source and the centre of rotation), without resorting to interpolation in the projection domain, a process which will introduce inaccuracies in the reconstructed images and should thus at all times be avoided.

 

Fig. 2 Two approaches to the specification of a cone-beam projection geometry.

Download Full Size | PDF

Instead of specifying all these parameters in general with a fixed object source distance and detector size and instead of assuming that the X-ray source nicely rotates around the centre of rotation, the ASTRA Toolbox also allows users to specify each projection separately. This is called a vector-based geometry (Fig. 2(b)). As an example for a cone-beam projection model, each projection image is then defined by four 3D vectors in the volume coordinate system:

2.3. Projection operations

The most important building blocks of any reconstruction technique are the projection operations. Consider tomography in its algebraic form, the forward projection operator then simulates projection images from a 2D or 3D volume by computing a sparse matrix vector multiplication (SpMV):

The projection building blocks in the ASTRA Toolbox is where the values w_ij of W are calculated and used to compute projection and backprojection images. As these are the operations that typically make up most of the computation time and power of a tomographic algorithm, much effort was put into making these implementations as computationally efficient as possible. To achieve this, the ASTRA Toolbox uses the computational power provided by GPU cards through the NVIDIA CUDA language [16]. In case the user has no GPU cards available, OpenMP accelerated CPU implementations of the 2D versions of these building blocks are also provided.

2.4. Reconstruction algorithms

All tomographic reconstruction algorithms are built on top of the basic projection and back-projection building blocks. These building blocks are internally configured by the algorithm itself and thus need no intervention of the user. Reconstruction methods that are available in the ASTRA Toolbox include the analytical reconstruction methods FBP for parallel and fan beam and FDK for cone beam datasets, several algebraic reconstruction methods such as SIRT, SART, and the Krylov subspace method CGLS. All of these are available both for 2D and 3D datasets, and also allow for some optional customization such as reconstruction limited to pixels in a certain mask, or enforcing minimum/maximum constraints during iterations. Configuration of an algorithm can be done in the MATLAB or Python interface by linking it to the correct data identifiers and by setting some algorithm specific options. It is important to note that the forward-and backprojection building blocks themselves can also be directly accessed from within the external interfaces. This means that a user can easily develop new tomographic algorithms or prototypes in which a substantial portion of the computational burden is offloaded to a much faster GPU card.

3. Conventional cone beam example

This section demonstrates the interface of the ASTRA Toolbox for a user wishing to perform a simple SIRT reconstruction of a set of conventional cone beam projection images. The code examples that are provided make use of the MATLAB interface to the ASTRA Toolbox. The Python interface uses the same design and interface principles and differs only in its syntax.

It is assumed that projection data is already stored as a 3D array in the MATLAB environment. As each CT machine utilizes different file formats to store projection images and geometry information, the ASTRA toolbox does not provide methods for reading in specific file formats, but instead relies on the user to be able to do this.

4. Non-conventional geometries

This section explores a few distinct applications of CT with a non-conventional projection geometry. The specification of each of these can be done in the ASTRA Toolbox using the vector based approach as described in Section (2.2) and Fig. 2(b). Specifically, with l the total number of projection images, an l × 12 matrix is constructed in which each row represents the geometry of a single projection. The first three columns represent the position of the X-ray source (S), the fourth to sixth column represent the position of the centre of the detector array (D), and columns seven to nine and ten to twelve represent both unit detector vectors (U and V). This matrix and the number of detectors are then supplied to the toolbox when building a cone beam projection geometry (Table 2).

Table 2. MATLAB code demonstrating how a conventional cone beam setup can be specified by a vector-based geometry object. The variable theta represents a l × 1 vector of projection directions; SOD and ODD represent the source-to-object and object-to-detector distance respectively; detWidth and detHeight represent the size of each detector pixel; and detCol and detRows represent the number of detectors in each projection image.

View Table | View all tables in this article

4.1. Rotary computed laminography

In conventional cone-beam CT, the object is imaged in a circular trajectory. However, for large flat objects, such as circuit boards and paintings, this is not possible as they might not physically fit longitudinally between the source and detector, or because the total X-ray attenuation in that direction might be too high. In those cases rotary computed laminography can offer a solution by tilting the object and rotation stage such that the rotation axis is no longer perpendicular to the plane defined by the point source and the central row of the detector (Fig. 3) [9]. This is equivalent to a setup with a fixed object, where the source and detector follow a circular trajectory parallel to the z = 0 plane. The radius of these circles depends on the angle α by which the object is tilted, and on the object-to-source or object-to-detector distance. The choice of α also represents a trade-off between a small field of view (FOV) (small α) and long intersections of the X-ray beam with the object and thus possibly noisy projection statistics (large α). For α = 90°, this setup is equivalent to a conventional cone-beam setup.

 

Fig. 3 In rotary computed laminography, the object is tilted such that its axis of rotation is not perpendicular to the plane defined by the X-ray source and the central row of the detector plane. This is equivalent to a setup with a fixed object where the source and detector follow a circular trajectory parallel to the z = 0 plane. (a) Side-view. Note how the rotation stage must have an opening for the X-ray beam to pass through. (b) Three-dimensional sketch of the projection geometry.

Download Full Size | PDF

Representing this type of projection geometry in the ASTRA Toolbox is relatively straightforward. Define Δt as the size of each detector pixel, Δs as the distance from the source to the centre of the coordinate system, and Δd as the distance from the centre of the detector array to the centre of the coordinate system. Let i ∈ {1, . . . , l} denote the projection index and let θ_i denote the viewing perspective corresponding to each projection. The full projection geometry is then:

In Fig. 4, an example is shown of a laminographic reconstruction of a large, yet flat, object. The CAD model depicted in Fig. 4(a) represents a 1000 × 1000mm circuit board with an array of holes with protruding edges drilled into it. In total, 270 800 × 800 projections (Fig. 4(b)) were simulated of this CAD model using the NVIDIA Optix raytracing framework. The laminographic tilt was set to α = 50°, the source-to-object distance to Δs = 750mm, and a object-to-detector distance to Δs = 1500mm. A realistic level of Poisson noise was applied to the projection images. A 300 iteration CGLS reconstruction was subsequently performed using the ASTRA Toolbox, of which two slices are shown in Figs. 4(c) and 4(d). It is clear that only a small portion of the object falls into the field of view.

 

Fig. 4 Simulation study of a laminographic reconstruction of a large flat object. The projection data (b) was simulated using the NVIDIA Optix raytracing toolbox.

Download Full Size | PDF

4.2. Tomosynthesis

Tomosynthesis is an imaging technique mainly used in medical radiography when only a few limited-angle projection images can be taken and when the minimization on X-ray dose is of crucial importance. Examples include breast cancer screening, where tomosynthesis has a higher accuracy compared to conventional mammography [17]; and portable, in-bed, radiography for patients in intensive care units that can not be positioned in a CT or MRI scanner. Tomosynthesis can also be applied for non-medical purposes, such as baggage screening in airports [18].

In tomosynthesis, projections are acquired with an X-ray source which moves on a linear or circular trajectory above the object, and a stationary flat panel detector. An example is illustrated in Fig. 5(a). In this geometry, only a limited angle can be covered by the acquisitions resulting in a depth resolution that is lower than the in-plane resolution. Therefore, the reconstructed volume is typically formed by anisotropic voxels, with a voxel size which is several times higher in the in z-direction compared to the x- and y-direction.

 

Fig. 5 In tomosynthesis, a stationary detector is placed underneath the object. The X-ray source moves on a circular trajectory above the object, creating a set of projections from a limited angle.

Download Full Size | PDF

Let i ∈ {1, . . . , l} denote the projection index and let θ_i denote the acquisition angle corresponding to each projection. The detector plane is in a fixed position, centred underneath the object centre at distance Δd. The full projection geometry is then defined by:

Consider the example shown in Fig. 5(b). A total of 45 tomosynthesis projection images were acquired of a human hand. The projection angles θ_i were equally distributed over an angular range of [−20°, 20°]. The object centre was placed at a distance of Δd = 65mm above the detector panel and the object-to-detector distance was Δs = 1039mm. The detector counted 2208 × 2668 pixels of size 0.160μm. Images were acquired at 48kV and 0.32mAs. These projection images were reconstructed (Fig. 5(c)) using 100 SIRT iterations on a volume containing 552 × 667 × 30 anisotropic voxels of size (0.640, 0.640, 4.0)mm.

4.3. Conveyor belt tomography

In industrial applications, accurate quality control is often crucial in order to maximize profit and customer satisfaction. Non-destructive testing using computed tomography can be a valuable tool to achieve that [19]. In manufacturing, CT can be used to inspect the quality of welding seams [20], 3D printed objects [21], printed circuit boards, etc. In the food processing industry, using CT, foodstuffs that do not meet certain requirements (e.g., apples with early rotting on the inside) can be discarded before shipping [3], and the cutting of meat (e.g., pig carcasses) can be optimized when it is known where the meat, fat and bones are located [2].

For fast inspection of objects in a mass production setting (a production line), the use of a conventional CT is not convenient. Placing every object separately in a CT scanner is both expensive and very time consuming. In order not to interfere with the efficiency of the production line, inline tomography can be used. Here, an object moves on a conveyor belt where a fixed X-ray source is positioned on one side, and a fixed flat-panel detector is placed opposite to it. While the object translates on the conveyor belt with a speed v, it passes through the field of view of the source and detector, and multiple cone beam projection images can be recorded. To allow projections from the full angular range, also consider a rotation of the object around its central axis with rotation speed ω, perpendicular to the conveyor belt (Fig. 6). In the fixed coordinate system <x, y, z>, this geometry is defined by:

 

Fig. 6 Conveyor belt geometry of an object moving on a conveyor belt with a fixed source and flat-panel detector.

Download Full Size | PDF

In Fig. 7, an example is shown of an inline reconstruction of a toy. The inline projection data is derived from 493 projections of a conventional circular CT scan (Fig. 7(a)). After conversion, 100 524 × 3000 projections were obtained at a 10ms interval, with source-to-object distance Δs = 121mm, object-to-detector distance Δd = 40mm, detector pixel size Δu = 35.36μm, a conveyor belt velocity of v = 100μm/ms, an object rotation speed of ω = −0.36°/ms, and A₀ = −50mm. Fig. 7(b) shows a reconstruction of the central slice without the rotation of the object on the conveyor belt (i.e. ω = 0°/ms) and Fig. 7(c) shows the same reconstruction but with rotation ω = −0.36°/ms. It is clear that the extra rotation is required to obtain an accurate reconstruction.

 

Fig. 7 Reconstruction of a toy in an inline scanning geometry.

Download Full Size | PDF

4.4. Adaptive zooming

In a cone beam geometry, the source-to-detector distance is sometimes referred to by the zooming factor, since the closer the source is to the object, the larger it appears on the detector. Ideally, it is chosen such that a projection of the object just fills the detector array. That way, the projection image has the highest achievable resolution without inducing object truncation and region-of-interest artifacts. Consider, however, an elongated object, such as in Fig. 8. In conventional setups, the source-to-object distance should, due to physical constraints, be chosen with respect to the orientation in which the longest axis of the object lines up with the source-to-detector line. This guarantees the best resolution in that direction (Fig. 8(a)), but leads to suboptimal detector usage in other projections (Fig. 8(b)).

 

Fig. 8 For elongated objects, a fixed source-to-detector distance leads to suboptimal use of the detector array in some viewing directions.

Download Full Size | PDF

With an adaptive zooming approach, the idea is to move the position of the object closer to or farther away from the source as the scan is being taken. The most optimal way of doing so is by precomputing the trajectory from the convex hull of the scanned object [10]. This approach allows to acquisition of more information about the object, thus increasing the reconstruction quality.

To mimic a tomographic system with variable source and detector position, an experiment was conducted on a desktop micro-CT system SkyScan-1172 (Bruker-MicroCT, Belgium). A piece of a pencil with a diameter of 7mm and a length of 15mm was used as an elongated object. Seven full-angle datasets were obtained, each containing 600 projections of 880 × 666, with the source-to-object distances ranging from 80.77 to 117.01mm. The source-to-detector distance was fixed at 216.392mm.

In Fig. 9(a), a 700 iteration SIRT reconstruction is shown from the largest source-to-object distance. Based on this reconstruction, a convex hull of the object was obtained. For each projection angle the closest possible source position was calculated according to the technique described in [10], and a projection was chosen from the dataset obtained from the smallest distance bigger than or equal to the distance from the calculated source position to the centre of rotation. The resulting trajectories are presented in Fig. 9(c), where is can be seen that the zooming factor is increased for the angles in which the object has a narrow shadow. The corresponding reconstruction is shown in Fig. 9(b), where the adaptive zooming approach results in a slightly improved contrast.

 

Fig. 9 Reconstructions of an elongated object (a piece of pencil) with a conventional and an adaptive zooming projection geometry. [10]

Download Full Size | PDF

4.5. Automatic projection alignment

Another advantage of considering the projection geometry in its vector form is its inherent flexibility when it comes to correcting for projection misalignment. Such misalignments may come in many forms. For example, the centre of rotation might not project exactly to the centre of the detector array and the detector plane may not be exactly perpendicular to the line connecting the source and the centre of rotation. Moreover, the rotation stage might have small inaccuracies, the focal spot size might not be constant, and the scanned object may move during the scan. All of these misalignments can be modelled simply by altering the vectors of the projection geometry. However, to do so, one first needs to know exactly which misalignments have occurred during the scan, and by how much. Recently, it has been shown that it is possible to automatically find these misalignment parameters by optimizing the vectors in the projection geometry such that the projection difference between the measured data and simulated data of the reconstructed images is minimal [22].

Here, we consider a simplified version of that work, focused on only a subset of misalignment factors: the detector offset δU and δV, and detector tilt ϕ (Fig. 10(a)). These misalignments generally result in double edges, such as in Fig. 11(b). When the offset and tilt factors are known, one can first translate and rotate the projection images, and subsequently continue with a reconstruction with a conventional geometry. However, instead of changing the data to fit the geometry, which requires interpolation in the projection domain, it is better to change the geometry to fit the data. Let R_i,ϕ denote the rotation matrix around the vector D_i − S_i with angle ϕ. The corrected projection geometry can then be expressed by:

 

Fig. 10 (a) A shifted and tilted detector array, seen from the point of view of the X-ray source. (b) Schematic overview of an iterative procedure to automatically estimate detector shift and tilt.

Download Full Size | PDF

To demonstrate this method, consider Fig. 11. A plastic box containing candy and a highly attenuating QRM resolution phantom was scanned in the TomoLab microCT at the SYRMEP beamline at Elettra, Trieste, Italy. In total 1444 projections of 890 × 1336 were acquired, of which only 288 are considered here. The source-to-object distance was 6000mm, the object-to-detector distance was 1000mm, and the detector size was 0.375 μm. A single projection image is shown in Fig. 11(a). The uncorrected reconstruction shown in Fig. 11(b) clearly shows severe misalignment artifacts. After 4 iterations of the automatic alignment method presented above, these artifacts have disappeared (Fig. 11(c)). Figs. 11(e) and 11(f) show the corresponding projection difference.

 

Fig. 11 Automatic detector offset and tilt estimation of a dataset containing a box of candy with a highly attenuating resolution phantom.

Download Full Size | PDF

5. Algorithm prototyping

Suppose a user has a custom MATLAB library with an implementation of some special advanced reconstruction algorithm. One example might be some form of total variation minimization (TVmin), for example using Chambolle-Pock [23], which currently does not have a native implementation inside the ASTRA framework. Custom libraries are typically built for general use, and their forward- and backprojection operations are likely performed as sparse matrix vector multiplications (SpMV). In order to use these libraries, the user is thus responsible for the construction of the projection matrix. While the ASTRA Toolbox provides the functionality to turn any 2D geometrical setup into a sparse projection matrix (see Table 3), this can only be advised for very small problems. For realistic data sizes, the projection matrix is simply too large to fit into the main memory space. Moreover, as the performance in tomographic reconstructions is latency-bound rather than computation-bound [24], it is best to always recompute the projection matrix on the fly, i.e. during each forward- or backprojection. The ASTRA Toolbox has very efficient implementations for these projection operators, exploiting the power of the GPU hardware. Ideally therefore, the projection operations of the ASTRA Toolbox are used inside the advanced logic of the custom reconstruction library.

Table 3. MATLAB code demonstrating how to extract a projection matrix for a certain projection and volume geometry from the ASTRA toolbox, and how to use this matrix inside an external function performing a SIRT reconstruction. Due to memory constraints, this approach is feasable only for small problems.

View Table | View all tables in this article

The naive approach to this would be to open up the external library and replace all the SpMV occurrences with function calls to the appropriate elements of the ASTRA Toolbox (Table 4). While possible, this is not advisable as the user might not be familiar with the internal code of the custom library and maintaining its code can be quite difficult.

Table 4. MATLAB code demonstrating a naive approach to include ASTRA projector operators in a SIRT function.

View Table | View all tables in this article

A better approach is to use a recent MATLAB library called the Spot-Toolbox [25]. This allows one to overwrite operations like matrix-multiplications with custom logic. The ASTRA Toolbox for MATLAB comes combined with such a spot-operator called opTomo [26]. This operator can be configured with the projection and volume geometry of the setup, after which MATLAB will regard it in many ways as a normal matrix. It can thus be passed to the external library as a normal function parameter. However, each time MATLAB detects a multiplication of this matrix and a vector or a multiplication of the transposed matrix and a vector, respectively the ASTRA forward projection or backprojection code is called (Table 5).

Table 5. MATLAB code demonstrating how the opTomo object is configured and simply passed on to the same SIRT-function as in Table 3.

View Table | View all tables in this article

As an example, consider the Simultaneous Iterative Reconstruction Technique (SIRT). This algebraic reconstruction method regards the reconstruction problem as the solving of the system

In Fig. 3, a MATLAB function is shown that performs a SIRT reconstruction given a projection matrix (W), projection data (p), an initial estimation (v⁽⁰⁾), and a number of iterations. In practice, this script would not be useful as SIRT also has a native implementation inside the ASTRA Toolbox, which is much faster. However, this does provide the flexibility to quickly prototype new ideas. To use the SIRT function, a projection matrix is extracted from the ASTRA Toolbox and passed to the function as a parameter. This method requires a lot a memory, and does not exploit the efficient implementation of the ASTRA Toolbox. With a naive linking of the ASTRA Toolbox and the script in Table 4, the projection operations run on the GPU using the ASTRA Toolbox implementation, but the code turns from 6 elegant lines to 10 lines that are a lot less readable. Also note that the function parameter list is slightly different, which might break other pieces of code. With the usage of the Spot tools, in Table 5, the best of both worlds are combined.

In terms of performance, Fig. 12(a) shows the timings of these different SIRT implementations for 100 iterations on a 128 × 128 volume with 180 projections, as performed on a system containing an Intel Core i7-6700HQ CPU running at 2.60GHz, an NVIDIA GeForce GTX 960M GPU, and 32GB of system memory. Of the three implementations presented in this section, the original script (Fig. 3) is the slowest as it does not make use of the efficiently implemented ASTRA Toolbox projection operations. Both the naive (Table 4) and the opTomo (Table 5) code are substantially faster, with the opTomo code being only very slightly the weaker of the two. Also included in Fig. 12 are two native ASTRA implementations of SIRT. The one that fully runs on the CPU, turns out to be substantially slower than any of these MATLAB implementations. This is due to MATLAB also using GPU computations for performing SpMV’s. The ASTRA GPU based implementation is slightly more than a factor two faster than the fastest MATLAB code. This is mainly due to the constant copying back and forth of the data that is going on in these MATLAB versions. In Fig. 12(b), the same experiment is repeated with a 100 iteration CGLS reconstruction algorithm.

 

Fig. 12 Timings for 100 iterations of various implementations of the (a) SIRT and (b) CGLS reconstruction algorithm on a 100 × 100 volume. The first three bars correspond to Table 3, 4, and 5, respectively. The final two bars refer to native ASTRA implementations of the SIRT and CGLS algorithm.

Download Full Size | PDF

In Fig. 13, reconstruction times of 100 SIRT iterations are plotted as a function of the volume size (both 2D and 3D). In each, a total of 180 parallel-beam projection angles were considered. This experiment was performed on two different systems: System A, as described in the previous paragraph; and System B, a 12-core Intel Xeon E5-2630 computation server running at 2.30GHz and fitted with 64GB of memory and a powerful Tesla K20 GPU card.

 

Fig. 13 Reconstruction times of 100 SIRT iterations as a function of the volume size (both 2D and 3D), on two distinct systems.

Download Full Size | PDF

6. Discussion and conclusions

In this article, the ASTRA Toolbox has been introduced and several use-cases were presented showing off its main features. These include its free, open source nature [14], its ability to describe virtually any 2D and 3D projection setup using vector geometries, and its flexibility regarding the easiness with which it can be included in other frameworks. Combined, they make the ASTRA Toolbox very suitable for fast prototyping of new applications, new geometry setups, and new reconstruction methods. Furthermore, due to its efficiently implemented building blocks, these prototypes can straightforwardly be upscaled to realistic data sizes and usages.

However, it should be noted that the ASTRA Toolbox is by no means the first and only tool to consider when dealing with tomographic reconstruction. Many users might prefer commercial software packages as they provide an easy-to-use graphical user interface, which the ASTRA Toolbox does not. Also, because of its flexibility regarding application fields, scanning devices and protocols, the ASTRA Toolbox can not provide out-of-the box support for various file formats. In order to use the toolbox, the user thus requires knowledge of these file formats, and the skill to parse and process them in the MATLAB or Python layer. Moreover, the algorithms bundled in the ASTRA Toolbox are limited to reconstruction methods, and do not include typical preprocessing (e.g., flat field correction, phase retrieval) and postprocessing (e.g., segmentation, morphological operations, mesh generation) algorithms. For practical use, the toolbox must thus be linked with other tools such as TomoPy [27]. These disadvantages limit the target audience of the ASTRA Toolbox mainly to researchers and users with a expertise in computer science.

While this may seem a huge issue, these research communities are in fact relatively large, and more of them are appearing in the literature. The list of applications that can benefit from the ASTRA Toolbox is extensive, and includes (i) material sciences, through electron tomography [4], diffraction contrast tomography with synchrotron radiation [28], microCT systems, laminography systems [9], etc.; (ii) biomedical sciences, through microCT [1], tomosynthesis setups, etc.; (iii) industrial purposes, through conveyor belt systems [3], microCT systems, etc.; and (iv) security purposes [18].

In summary, the ASTRA toolbox is an excellent platform to bridge the large gap between application scientists and researchers in the field of numerical mathematics and linear solvers [29]. That way, new and advanced numerical solvers can be tested on realistic data, benefiting both communities.