computer programs
Foamlike phantoms for comparing tomography algorithms
^{a}LIACS, Leiden University, Leiden, The Netherlands, and ^{b}Computational Imaging Group, CWI, Amsterdam, The Netherlands
^{*}Correspondence email: d.m.pelt@liacs.leidenuniv.nl
Tomographic algorithms are often compared by evaluating them on certain benchmark datasets. For fair comparison, these datasets should ideally (i) be challenging to reconstruct, (ii) be representative of typical tomographic experiments, (iii) be flexible to allow for different acquisition modes, and (iv) include enough samples to allow for comparison of datadriven algorithms. Current approaches often satisfy only some of these requirements, but not all. For example, realworld datasets are typically challenging and representative of a category of experimental examples, but are restricted to the acquisition mode that was used in the experiment and are often limited in the number of samples. Mathematical phantoms are often flexible and can sometimes produce enough samples for datadriven approaches, but can be relatively easy to reconstruct and are often not representative of typical scanned objects. In this paper, we present a family of foamlike mathematical phantoms that aims to satisfy all four requirements simultaneously. The phantoms consist of foamlike structures with more than 100000 features, making them challenging to reconstruct and representative of common tomography samples. Because the phantoms are computergenerated, varying acquisition modes and experimental conditions can be simulated. An effectively unlimited number of random variations of the phantoms can be generated, making them suitable for datadriven approaches. We give a formal mathematical definition of the foamlike phantoms, and explain how they can be generated and used in virtual tomographic experiments in a computationally efficient way. In addition, several 4D extensions of the 3D phantoms are given, enabling comparisons of algorithms for dynamic tomography. Finally, example phantoms and tomographic datasets are given, showing that the phantoms can be effectively used to make fair and informative comparisons between tomography algorithms.
Keywords: tomography; phantom; simulation; opensource; experiment design.
1. Introduction
In tomographic imaging, an image of the interior of a scanned object is obtained by combining measurements of some form of penetrating wave passing through the object. Tomographic imaging is routinely used in a wide variety of application fields, including medical imaging (Goo & Goo, 2017), materials science (Salvo et al., 2003), biomedical research (Metscher, 2009), and industrial applications (De Chiffre et al., 2014). To extract relevant information from the acquired data, the measurements are often processed by several mathematical algorithms in a processing pipeline. Common processing steps include tomographic reconstruction (Kak et al., 2002; Marone & Stampanoni, 2012; Ravishankar et al., 2020), artifact removal (Barrett & Keat, 2004; Münch et al., 2009; Miqueles et al., 2014), and image segmentation (Iassonov et al., 2009; Foster et al., 2014; Perciano et al., 2017). Because of the importance of tomography in practice, a wide variety of algorithms have been developed for these processing steps, and tomographic algorithm development remains an active research field. In addition to classical image processing algorithms, the use of datadriven machine learning algorithms has become popular in tomography in recent years (Jin et al., 2017; Yang et al., 2017; Pelt et al., 2018; Adler & Öktem, 2018; Liu et al., 2020; Yang et al., 2020).
To properly assess the available algorithms, it is essential to compare them with each other in a fair, reproducible, and representative way. Such comparisons are important for algorithm developers to understand how newly developed algorithms compare with existing approaches. Proper comparisons are also important for end users of tomographic imaging to learn which algorithms to use for certain experimental conditions, and to know what results to expect from each available algorithm. A common approach to compare tomography algorithms is to take a set of tomographic datasets, apply several algorithms to the data, and compare results. To make this approach as informative as possible, the chosen datasets should ideally satisfy the following requirements:
(i) The datasets should be challenging: it should not be trivial to obtain accurate results for them.
(ii) The datasets should be representative of typical objects, experimental conditions, and data sizes that are used in practice.
(iii) The datasets should be flexible with respect to object complexity and experimental properties, making it possible to explore the capabilities and limitations of each algorithm for different acquisition modes, experimental conditions, and object complexities.
(iv) The datasets should include enough samples to allow for the comparison of datadriven algorithms that require a large number of similar samples for training and testing.
The datasets that are used to compare algorithms in the current literature typically satisfy some of the requirements above, but not all. For example, realworld datasets from public databases (Hämäläinen et al., 2015; McCollough et al., 2017; Jørgensen et al., 2017; Singh et al., 2018; De Carlo et al., 2018; Der Sarkissian et al., 2019) are often used for comparison. While these datasets are both challenging and representative since they are obtained in actual tomographic experiments, they are typically not flexible as it is impossible to change the acquisition mode and experimental conditions that were used in the experiment. In addition, while some datasets are specifically designed for datadriven applications (McCollough et al., 2017; Der Sarkissian et al., 2019), other realworld datasets are often not suitable for datadriven approaches, since the number of scanned samples is often limited.
A common alternative to comparing results for realworld datasets is to use computergenerated phantom images for which virtual tomographic datasets are computed. One advantage of using such mathematical phantoms is that the true object is readily available, allowing one to compute accuracy metrics with respect to an objective ground truth. Another advantage is that this approach is flexible: since the tomographic experiment is performed virtually, different acquisition modes and experimental conditions can be easily simulated.
Popular examples of phantoms used in tomography include the Shepp–Logan head phantom (Shepp & Logan, 1974), the FORBILD head phantom (Yu et al., 2012), and the MCAT phantom (Segars & Tsui, 2009). In addition to predefined phantom images, several software packages have been recently introduced that allow users to design their own custom mathematical phantoms and generate simulated tomography datasets for them (Ching & Gürsoy, 2017; Faragó et al., 2017; Kazantsev et al., 2018). The main disadvantage of popular mathematical phantoms is that they typically consist of a small number of simple geometric shapes (i.e. less than 100). As a result, the phantoms are often not representative of realworld objects, which typically contain a much larger number of more complicated features. Several oftenused phantoms, e.g. the Shepp–Logan head phantom, consist of large uniform regions and can therefore be relatively easy to reconstruct accurately using certain algorithms, making it difficult to compare algorithms using these phantoms. Finally, predefined phantoms usually consist of only a single sample and manually defined phantoms require considerable time to design, making it difficult to effectively use them for datadriven approaches that require multiple samples for training.
Because of the aforementioned disadvantages of both realworld datasets and mathematical phantoms, a hybrid approach is used in practice as well (Adler & Öktem, 2018; Leuschner et al., 2019; Hendriksen et al., 2020). In this approach, reconstructed images of realworld tomographic datasets are treated as phantom images for which virtual tomographic datasets are computed. In this way, different acquisition modes and experimental conditions can be simulated for realistic phantom images. However, the approach has several disadvantages when comparing algorithms with each other. First, the reconstructed images have to be represented on a discretized voxel grid and often include various imaging artifacts, resulting in inaccurate representations of the actual scanned objects. Second, the approach can lead to the `inverse crime', i.e. when the same image formation model is used for both data generation and reconstruction, which can lead to incorrect and misleading comparison results (GuerquinKern et al., 2012). Finally, since imaging artifacts such as noise and data sampling artifacts are present in the phantom images, artifactfree objective ground truth images with which to compute accuracy metrics are not readily available.
To summarize, new datasets are needed that satisfy all requirements given above for improved comparisons between tomography algorithms. In this paper, we present a family of mathematically defined phantoms that aim to satisfy all requirements. The phantoms consist of threedimensional foamlike structures and can include more than 100000 features. Since foamlike objects are often investigated using tomography (Babin et al., 2006; Roux et al., 2008; Hangai et al., 2012; Raufaste et al., 2015; Evans et al., 2019), the proposed phantoms are representative of a popular class of objects. Furthermore, foamlike objects are typically challenging to accurately reconstruct and analyze due to the fact that they exhibit both largescale and finescale features (Brun et al., 2010), making them well suited for comparing tomography algorithms. Tomographic datasets can be computed for the proposed phantoms for a wide variety of experimental conditions and acquisition modes and with data sizes that are common in realworld experiments, making the approach both flexible and representative. Finally, an effectively unlimited number of random variations of samples can be generated, enabling comparisons of datadriven algorithms that require multiple samples for training.
The proposed family of simulated foam phantoms has already been used for comparing algorithms in several papers from various research groups (Pelt et al., 2018; Hendriksen et al., 2019, 2020; Liu et al., 2020; Etmann et al., 2020; Marchant et al., 2020; Renders et al., 2020; Ganguly et al., 2021). In this paper, a formal definition of the phantoms is given, and mathematical and computational details about both the phantom generation and tomographic experiment simulation are discussed. This paper is structured as follows. In Section 2.1 a mathematical description of the foam phantoms is given, and in Section 2.2 we introduce an algorithm that can compute such phantoms efficiently. We explain how, given a generated foam phantom, tomographic projections can be computed in Section 2.3. Several 4D (i.e. dynamic) variations of the proposed phantoms are introduced in Section 2.4. In Section 3, several experiments are performed to investigate the influence of various parameters on the generated phantoms, the computed projection data, and the final tomographic reconstruction. Furthermore, we discuss the required computation time for generating phantoms and computing projection data. In Section 4, we give a few concluding remarks.
2. Method
In this section, we give the mathematical definition of the proposed family of phantoms, and describe how such phantoms can be efficiently generated. In addition, we explain how projection images can be computed for both parallelbeam and conebeam geometries. As explained above, the main inspiration for the design of these phantoms is the continued popularity of investigating a wide variety of realworld foam samples using tomography. In Fig. 1, two examples are given of such samples, in addition to an example of a foam phantom from the proposed family of phantoms, showing the similarities in features between the proposed foam phantoms and realworld foam samples.
2.1. Mathematical description
In short, each phantom from the foam phantom family consists of a singlematerial cylinder with a large number (e.g. >100000) of nonoverlapping spheres of a different material (or multiple different materials) inside. A more detailed explanation follows. Each phantom is defined in continuous 3D space . In all phantoms, a cylinder is placed in the origin, parallel to the zaxis (the rotation axis). This cylinder has an infinite height and a radius of 1, with all other distances defined relative to this unit radius. Inside the main cylinder, N nonoverlapping spheres are placed, which will be called voids in the rest of this paper. Each void i has a radius and a position with p_{i} = (x_{i}, y_{i}, z_{i}). The area outside the main cylinder consists of a background material that does not absorb any radiation, i.e. all positions (x, y, z) with (x^{2} + y^{2})^{1/2} > r_{C} where r_{C} is the radius of the main cylinder (defined to be r_{C} = 1). The foam itself, i.e. all positions that are within the main cylinder but not inside a void, consists of a single material with an of 1, with all other attenuation coefficients defined relative to this. Each separate void in the phantom can be filled with a different material, each with its own userdefined . In the default settings, all voids are filled with the background material. To summarize, each void i can be completely characterized by a vector with five elements: its position x_{i}, y_{i}, and z_{i}, its radius r_{i}, and its c_{i}. Similarly, each foam phantom is completely characterized by the set of s_{i} vectors of all voids S = {s_{1}, s_{2},…, s_{N}}. The definition of a foam phantom is shown graphically in Fig. 2.
The vertical size of a phantom is controlled by ensuring that the z_{i} position of each void i satisfies , with a maximum position . In addition, no part of any void is allowed to exist outside the main cylinder, i.e. #!(x_{i}^{2}+y_{i}^{2})^{1/2} + #!r_{i} 1 for all voids. Also, no part of any void is allowed to overlap with any other void, i.e. d(p_{i}, p_{j}) ≥ r_{i} + r_{j} for all voids i and j, where d(a, b) is the Euclidian distance between points a and b. Finally, the size of the voids is controlled by choosing a maximum radius r_{max} and ensuring that r_{i} ≤ r_{max} for all voids.
2.2. Phantom generation
Foam phantoms are generated by starting with the main cylinder and repeatedly adding voids until N voids are placed. When placing the ith void at position p_{i}, the radius of the new void is limited by three considerations: (1) the distance to the outside of the main cylinder, , (2) the distance to the closest edge of any existing void, i.e. , and (3) the maximum allowed radius r_{max}. The final radius r_{i} is chosen to be as large as possible, meaning that
Since placed voids are made as large as possible, the size of newly placed voids naturally becomes smaller during the generation of the phantom: at the end of phantom generation, not much room is left for any new voids, resulting in smaller sizes. Another consequence is that each void i either touches the outside of the main cylinder, i.e. #!(x_{i}^{2}+y_{i}^{2})^{1/2}+r_{i} = 1, or touches at least one other void, i.e. d(p_{i}, p_{j}) = r_{i} + r_{j} for some j. As a result, the radius of none of the voids can be increased without either making the void overlap another void or having part of the void outside the main cylinder.^{1}
We found empirically that realistic looking phantoms (see Fig. 1) are obtained when new voids are placed in positions that allow for the largest possible void size out of all possible positions. Finding such optimal positions given a partial set of voids is not trivial, since the number of possible positions is, in theory, infinite. We note that it might be possible to deterministically find an optimal position in a computationally efficient way by using a plane sweep algorithm approach (Nievergelt & Preparata, 1982). However, we propose to use a much simpler approach: a set of N_{p} randomly picked trial points is available at all times, in which each trial point is a valid position where a void could be placed (i.e. inside the main cylinder but not inside an existing void). Then, out of all trial points, a void is placed at the point that results in the largest void. There are several advantages to this approach: (1) it is relatively simple to implement, (2) it is random in nature, enabling the generation of an infinite number of different phantoms, and (3) it is computationally efficient, allowing generation of phantoms with many voids within reasonable time.
To summarize, the algorithm to generate foam phantoms works as follows:
(1) Create a list of N_{p} trial points, randomly placed within the main cylinder and satisfying the maximum height z_{max}.
(2) Pick the trial point that results in the largest void (if multiple exist, randomly select one) and remove it from the list.
(3) Add a void at the picked trial point with the largest possible radius [equation (1)].
(4) Remove trial points that are inside the newly placed void from the list.
(5) Add new trial points to the list until the list has N_{p} points again, each randomly placed at a valid position (inside the main cylinder, outside any existing void, and satisfying the maximum height z_{max}).
(6) Repeat steps (2) to (5) until N voids are placed.
Note that each foam phantom can be recreated deterministically given the following values: the number of voids N, the number of trial points N_{p}, the maximum void size r_{max}, the maximum height z_{max}, and the random seed used for the random number generator.
There are several implementation tricks that can improve the computational performance of generating phantoms in practice. For example, the maximum possible radius [equation (1)] of each trial point can be precomputed and stored in a sorted data structure such as a skip list (Pugh, 1990) to enable fast access to the trial point with the largest possible radius. After placing a new void, the maximum radii have to be updated and reinserted in the sorted data structure, which can be efficiently done during step (4) above. For more details about these implementation tricks, we refer to the computer code that is available under an opensource license (Pelt, 2020).
2.3. Computing projections
The foam phantoms presented in this paper were developed for use in tomography research. As such, it is important that tomographic projections of these phantoms can be computed accurately and efficiently. Here, we assume that the projections are formed by the Radon transform: a measurement is computed by taking a line integral of the attenuation coefficients of the sample over the virtual Xray i. The orientation and direction of the virtual ray depends on the tomographic acquisition geometry that is simulated. Measurements collected by 2D pixels with a certain area, which often represent realworld experiments better than individual rays, can be approximated by supersampling, i.e. averaging the measurements of multiple rays within a single pixel.
In many tomographic experiments, projections are formed by rotating the sample in front of a 2D detector (or, equivalently, rotating the detector around the sample) and acquiring separate 2D projection images at different angles. In these cases, the projection data are naturally described by a set of 2D projection images, each taken at a specific angle . Depending on the experimental setup, incoming rays of a single projection image are often assumed to be either parallel to each other (parallelbeam geometries) or to originate from a point source (conebeam geometries).
In many existing comparisons between algorithms in which tomographic experiments are simulated, projections are formed by first discretizing the object on a discrete voxel grid and computing line integrals for the discrete object afterwards. As mentioned above, this approach can lead to the `inverse crime', which can produce incorrect and misleading comparison results (GuerquinKern et al., 2012). For the proposed foam phantoms, we prevent the inverse crime by computing projections analytically in the continuous domain, using the exact intersection between a ray and the phantom. Specifically, simulated Xray projections of the proposed foam phantom are computed raybyray. The line integral of the sample over a certain ray is computed by first computing the intersection of the main cylinder with that ray and subsequently subtracting the intersections with all voids, taking into account their attenuation factors. If we denote the intersection of ray i with the main cylinder by L_{i} and the intersection of ray i with void j by l(i, j), we can describe the projection P_{i} of ray i mathematically by
where c_{j} is the attenuation factor of void j, as described above. Note that for each void we have to both subtract the intersection that was counted in L_{i} and add the attenuation of the void itself, resulting in a factor of (1 − c_{j}).
In parallelbeam geometries, the intersection L_{i} of the main cylinder with ray i can be computed by
where dz_{i} is the shortest distance between any point along ray i and the zaxis (the rotation axis). The computation of this intersection is shown graphically in Fig. 3. For conebeam geometries, it is also possible to analytically compute the intersection between a ray and the main cylinder, although it is more complicated than equation (3). For the sake of brevity, we refer to the computer code (Pelt, 2020) for more details about this computation.
The computation of intersections between rays and voids is similar to that of the main cylinder,
where dv(i, j) is the shortest distance between the center of void j and any point along ray i, and r_{j} is the radius of void j, as described above. The derivation of equation (4) is based on Fig. 3, extended to three dimensions. For more details about the analytic expression of the Radon transform of a sphere, we refer to Toft (1996). Note that shortest distances dv(i, j) can be computed efficiently by projecting the center p_{j} of void j along the direction of ray i. For parallelbeam geometries, all rays of a single projection image are parallel to each other, which enables precomputing the projections of all void centers for each projection image, significantly reducing the required computation time.
Detector noise can be simulated by applying Poisson noise to each measurement. First, a virtual photon count I_{i} for measurement i is computed using the Beer–Lambert law: #!I_{i} = , where I^{0} is the number of incoming photons and γ is a factor that controls the amount of radiation absorbed by the phantom. Afterwards, a noisy photon count is computed by sampling from a with I_{i} as the expected value. The noisy photon count is transformed back to a noisy measurement = . In realworld tomographic experiments, other artifacts are often present in the measured data in addition to Poisson noise, for example due to source characteristics [e.g. beam hardening (Barrett & Keat, 2004)], additional photon interactions [e.g. free space propagation (Moosmann et al., 2011)], optical effects (Ekman et al., 2018), and detector defects (Miqueles et al., 2014). Simulating such additional artifacts is not yet supported in the current version of the computer code. However, we note that it might be possible to include such artifacts in the future, either during the computation of projections within the code or as a postprocessing step afterwards, possibly taking advantage of existing software packages that support them (Allison et al., 2016; Faragó et al., 2017).
2.4. 4D extensions
In recent years, improvements in radiation sources and detector equipment have increased interest in dynamic tomography of timeevolving samples (dos Santos Rolo et al., 2014; Maire et al., 2016; GarcíaMoreno et al., 2018). In these applications, samples are fourdimensional (4D) in nature, consisting of three spatial dimensions and one time dimension. To enable quantitative comparisons between algorithms for dynamic tomography (Kazantsev et al., 2015; Mohan et al., 2015; Van Nieuwenhove et al., 2017; Nikitin et al., 2019), 4D phantoms are needed. Similar to 3D phantoms, these phantoms should be challenging, representative, flexible, and suitable for datadriven applications. Here, we introduce such 4D phantoms by adapting the 3D foam phantoms described above, adding timeevolving aspects in different ways. Currently, the computer code includes three types of 4D extensions, which are described below. Additional 4D extensions are planned for future inclusion.
The first 4D extension is a moving phantom, in which the voids move along the zaxis. All voids move with the same velocity, but the velocity changes randomly during the experiment. The second extension is an expanding phantom, in which the size of all voids increases during the experiment. The third extension is an infiltration phantom, in which the voids are slowly filled by a material with a different than the initial void material. Specifically, all voids at a certain chosen height are filled at the start of the experiment. Then, each unfilled void with an edge close to a filled void is filled after a randomly chosen interval. This process is repeated until all voids are filled. Example phantoms for the three 4D extensions are shown in Fig. 6.
For all 4D extensions, several parameters can be chosen to adjust the time evolutions of the generated phantoms. After generating each 4D phantom, a dynamic tomography experiment can be simulated by virtually rotating the phantom during its time evolution, and computing projections as described in Section 2.3. Changes within the sample that happen during the acquisition of a single projection can be modeled by supersampling in time.
3. Experiments
3.1. Implementation details
Computer code to generate the proposed foam phantoms and simulate tomographic experiments is available as the opensource foam_ct_phantom software package (Pelt, 2020). The code is available for the Windows, MacOS, and Linux operating systems, and can be installed using the Conda package management system:
The code is implemented in the Python 3 (Van Rossum & Drake, 2009) programming language. Parts of the code with a high computational cost are implemented in the C programming language (Kernighan & Ritchie, 1988), using OpenMP (Dagum & Menon, 1998) for parallelization. Projections for conebeam geometries can also be computed using NVidia Graphic Processor Units (NVidia, Santa Clara, CA, USA), which significantly reduces the required computation time. The GPU code is implemented using the Numba package (Lam et al., 2015).Generated phantoms and projection datasets are stored in HDF5 file containers (Folk et al., 2011), using a simple custom data format that includes metadata about how the phantom or dataset was generated. A skip list (Pugh, 1990) is used to enable fast access to the trial point with the largest possible radius during generation of phantoms, and random numbers are generated using the Mersenne Twister algorithm (Matsumoto & Nishimura, 1998). The experiments in this paper were performed on a workstation with an AMD Ryzen 9 3900X CPU (AMD, Santa Clara, CA, USA), running the Fedora 32 Linux operating system. Experiments involving GPU computations were performed using a server with four NVidia GeForce RTX 2080 Ti GPUs (NVidia, Santa Clara, CA, USA), running the Fedora 30 Linux operating system.
3.2. Code examples
Below, we give a few code examples to show how the computer code can be used in practice to generate new phantoms, simulate tomographic experiments, and reconstruct projection data. First, new foam phantoms can be generated using the following Python code:
Here, the five parameters that determine the phantom shape (see Section 2.2) are given by #!nspheres, #!ntrial _ points, #!random _ seed, #!rmax, and #!zmax.
Once a phantom has been generated, parallelbeam projection data can be computed by the following Python code:
Here, the #!supersampling parameter controls the number of rays that are simulated within each pixel. Specifically, #!supersampling^{2} (i.e. #!supersamplingsquared) rays are simulated within each pixel, evenly distributed over the area of the pixel in a #!supersampling × #!supersampling grid. The measured projection of a pixel is then the average value of the measurements of all rays within that pixel. For conebeam projection data, only the geometry specification has to be changed to
Here, #!sod and #!odd denote the source–object distance and object–detector distance, respectively.
The computer code also includes utility functions to assist in reconstructing the generated projection data using existing tomography toolboxes such as the ASTRA toolbox (Van Aarle et al., 2016) and TomoPy (Gürsoy et al., 2014). For the ASTRA toolbox, functions are included to convert defined geometries to equivalent ASTRA geometries:
More code examples are included in the source code of the foam_ct_phantom package.
3.3. Phantom examples
In this section, we present several examples of generated phantoms, and investigate the effect of the various generation parameters on the phantom characteristics. As explained in Section 2.2, each phantom is defined by the number of voids N, the number of trial points N_{p}, the maximum void size r_{max}, and the maximum height z_{max}. In the following, the values used for generating phantoms are N = 150000, N_{p} = 10^{6}, r_{max} = 0.2, and z_{max} = 1.5 and all voids are filled with the background material, unless stated otherwise. Note that the code supports filling voids with other materials as well, making it possible to simulate objects with various characteristics, e.g. with lowcontrast features. In Fig. 4, generated phantoms are shown for various numbers of included voids N. Since the other parameters are identical for all shown phantoms, the figure also shows how a phantom is generated by increasing the number of included voids. As expected, phantoms with a small number of voids mostly include relatively large voids, while the void size decreases with increasing numbers of included voids. In addition, the figures show both the largescale and finescale features that are present phantoms with relatively large numbers of voids.
In Fig. 5, generated phantoms are shown for three maximum void sizes r_{max} and N = 150000. The results show that the phantom features depend significantly on the choice of r_{max}: for a relatively large maximum void size (r_{max} = 0.8), there are a few large voids present in the phantom and a large number of relatively small voids, as the large voids restrict the available space for the remaining voids. For a relatively small maximum void size (r_{max} = 0.05), most voids in the phantom have a similar size. The phantom with an intermediate maximum void size (r_{max} = 0.2) exhibits both characteristics to a lesser degree. These results show that the proposed phantom family can be used to simulate a wide variety of foam structures. Examples of the phantoms generated by the 4D extensions described in Section 2.4 are shown in Fig. 6.
3.4. Projection data and reconstruction examples
In this section, we present several examples of generated projection data and compare reconstruction results using several popular tomographic reconstruction algorithms. In all cases, we use the foam phantom generated with N = 150000, N_{p} = 10^{6}, r_{max} = 0.2, and z_{max} = 1.5. Parallelbeam projections are computed for a detector with 2560 × 2160 pixels and 16 rays per pixel (i.e. 4 × 4 supersampling), mimicking a PCO.edge 5.5 sCMOS detector (PCO, Kelheim, Germany) that is commonly used at synchrotron tomography beamlines (Mittone et al., 2017). The width and height of a detector pixel was set to 3/2560, resulting in a detector width of 3, with the sample (which has a fixed radius of 1) projecting on twothirds of the detector width. Projections were computed for four imaging scenarios: `highdose', with a large number of noisefree projections; `noise', with a large number of projections with a significant amount of Poisson noise applied; `few projections', with a relatively low number of noisefree projections; and `limited range', with a large number of noisefree projections acquired over less than 180°. Specific details about scenarios are given in Table 1, and example projection data are shown in Fig. 7.
We compare results for several popular tomographic reconstruction algorithms: the filtered backprojection method (FBP) (Kak et al., 2002), SIRT (Kak et al., 2002), CGLS (Scales, 1987), SART (Andersen & Kak, 1984), and SIRT and SART with additional nonnegativity constraints on the pixel values (Elfving et al., 2012). All reconstruction images were computed using the optimized GPU implementations of the ASTRA toolbox (Van Aarle et al., 2016). We compare the reconstructed images using three popular image quality metrics, the root mean square error (RMSE), peak signaltonoise ratio (PSNR), and the multiscale structural similarity index (MSSSIM) (Wang et al., 2003). We also compare the images using two segmentation metrics, in which the images are segmented using thresholding, and Dice scores (Bertels et al., 2019) are computed for voxels inside large voids (with radii r_{i} ≥ 0.1) and small voids (with radii r_{i} < 0.05). All metrics are computed with respect to a ground truth image that consists of a discretization of the foam phantom with the same number of voxels as the reconstructions and using 64 (4 × 4 × 4) sampling points per voxel. For the iterative algorithms, the number of iterations that minimizes the RMSE is used, which is determined using the Nelder–Mead method (Nelder & Mead, 1965) for CGLS and a simple grid search for all other algorithms.
In Table 2, the quality metrics are given for the central slice of the phantom and the four projection data scenarios given above. The results show that in most cases FBP produces images with the highest RMSE and lowest MSSSIM values, the three unconstrained iterative methods (SIRT, CGLS, and SART) produce images with lower RMSE and higher MSSSIM than FBP, and the iterative methods with nonnegativity constraints produce images with the lowest RMSE and highest MSSSIM. However, the segmentationbased metrics show more nuanced results. For example, in the `limited range' scenario, the Dice score for large voids of the FBP reconstruction is close to the Dice scores of the iterative algorithms, even though the RMSE is significantly higher and MSSSIM significantly lower. This shows that, if the specific application of tomography would require only the analysis of large voids, the FBP algorithm would be sufficient, even though its image metrics are significantly worse than other methods. Similar results are shown in Fig. 8, in which a few selected reconstructed images are shown. Such results partly explain the continued popularity of FBPlike methods in practice (Pan et al., 2009).

In Fig. 9, the PSNR of FBP, SIRT, and SIRT with a nonnegativity constraint is given as a function of both the number of projection angles and the number of voids in the phantom. In each case, data were generated with a low amount of Poisson noise (I_{0} = 10^{7}) and a foam material that corresponds to an average absorption of 10% of virtual photons for a phantom with 150000 voids. The results show how the behavior of each reconstruction algorithm depends on the complexity of the scanned sample and the amount of acquired data. For example, the results show that the accuracy of FBP does not depend significantly on the complexity of the phantoms, while the accuracy of the SIRT algorithms is significantly improved for lowcomplexity samples compared with highcomplexity samples. We note here that the proposed family of foam phantoms is especially well suited for performing such detailed comparisons, as the phantoms exhibit features at multiple scales, the level of complexity is tunable, and complete groundtruth information about the void positions and sizes is known. For another example of a detailed taskbased analysis that uses the foam phantoms, we refer to Marchant et al. (2020).
3.5. Computation time
In this section, we present results for measurements of the required computation time for generating a foam phantom and computing projection data. A theoretical analysis of the required computation time for phantom generation is technically complicated due to the random nature of placing trial points. However, we hypothesize that, for large numbers of voids, the most timeconsuming part is step (5) of the algorithm described in Section 2.2, and that the required computation time scales with , where N is the number of voids and N_{p} the number of trial points used during generation. The various terms come from the fact that the required time for inserting an item in a skip list scales with , the number of new trial points that have to be placed scales with N_{p}, the number of voids that have to be checked for overlap for each new trial point scales with N, the number of required random trials until a valid position is found scales with N (since the available space decreases when more voids are placed), and step (5) has to be evaluated N times. The required computation time for simulating a projection scales with N_{r}N, where N_{r} is the number of simulated rays, which depends on the size of the detector and the amount of supersampling used.
In Fig. 10, the computation time required to generate a phantom with r_{max} = 0.2 is given as a function of the number of included voids. The results show that phantoms with a large number of voids, e.g. 100000 voids, can be generated within a few minutes. The results also show that using multiple CPU cores can significantly reduce the required computation time, especially for large numbers of voids. It is interesting to note that three different phases can be identified in Fig. 10. We hypothesize that different parts of the algorithm are dominant within each phase. During generation of the first 100 voids, we expect that most time is spent inserting newly placed trial points in the linked list data structure, which is not parallelized in the current implementation. Between 100 and around 10^{5} voids, we expect that most time is spent updating the maximum possible radius of each trial point, which is highly parallelizable. Finally, for more than 10^{5} voids, we expect that most time is spent finding valid positions while randomly placing new trial points, which is not parallelized in the current implementation as well. It may be possible to use such observations to reduce computation time for generating phantoms in the future.
In Fig. 11, the computation time for generating projections is given as a function of the number of rows and columns in each projection, for a phantom with 150000 voids and r_{max} = 0.2. The results show that it is possible to compute a parallelbeam projection with common highresolution numbers of pixels, e.g. 1024 × 1024 pixels, in less that a tenth of a second using a modest multicore CPU system. This computational efficiency makes it possible to generate full tomographic datasets within a few minutes. As explained in Section 2.3, computing conebeam projections is more computationally demanding than computing parallelbeam projections. This is indeed visible in Fig. 11, which shows that even when using multiple CPU cores, generating a conebeam projection can take considerable time. However, multiple GPUs can be used to significantly speed up these computations, reducing the required computation time to a few seconds per projection for common detector sizes.
4. Conclusion
In this paper, we introduced a family of foamlike phantoms for comparing the performance of tomography algorithms. The generated phantoms are challenging to reconstruct, representative of typical tomography experiments, and flexible, as projections can be calculated for various acquisition modes. In addition, an unlimited number of varying foamlike phantoms can be generated, enabling comparisons of datadriven algorithms. The phantoms consist of a main cylinder with a large number of randomly placed voids (e.g. more than 100000), resulting in foamlike structures with both largescale and finescale features. We also introduced several 4D extensions to the static 3D phantoms, resulting in timeevolving phantoms for comparing algorithms for dynamic tomography.
Computationally efficient ways of both generating a phantom and simulating projection data for a given phantom were discussed, and a software package that implements these algorithms, the foam_ct_phantom package (Pelt, 2020), was introduced. Experimental results show that it is possible to generate phantoms on a modest workstation within a few minutes, and that projection data can be simulated for common highresolution detector sizes within a few minutes as well. Comparisons between common reconstruction algorithms for several experimental settings show that it is possible to perform detailed analyses of algorithm performances using the proposed phantom family. These results show that the phantoms can be effectively used to make fair and informative comparisons between tomography algorithms.
Footnotes
^{1}In rare cases, voids with the maximum size r_{max} can exist that touch neither the outside of the main cylinder nor another void. However, the radius of these voids cannot be increased as well, since their radius would become larger than the maximum allowed radius.
Funding information
Funding for this research was provided by: Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant No. 016.Veni.192.235 to Daniël M. Pelt; grant No. 639.073.506 to Kees Joost Batenburg).
References
Aarle, W. van, Palenstijn, W. J., Cant, J., Janssens, E., Bleichrodt, F., Dabravolski, A., De Beenhouwer, J., Joost Batenburg, K. & Sijbers, J. (2016). Opt. Express, 24, 25129–25147. Web of Science PubMed Google Scholar
Adler, J. & Öktem, O. (2018). IEEE Trans. Med. Imaging, 37, 1322–1332. Web of Science CrossRef PubMed Google Scholar
Allison, J., Amako, K., Apostolakis, J., Arce, P., Asai, M., Aso, T., Bagli, E., Bagulya, A., Banerjee, S., Barrand, G., Beck, B., Bogdanov, A., Brandt, D., Brown, J., Burkhardt, H., Canal, P., CanoOtt, D., Chauvie, S., Cho, K., Cirrone, G., Cooperman, G., CortésGiraldo, M., Cosmo, G., Cuttone, G., Depaola, G., Desorgher, L., Dong, X., Dotti, A., Elvira, V., Folger, G., Francis, Z., Galoyan, A., Garnier, L., Gayer, M., Genser, K., Grichine, V., Guatelli, S., Guèye, P., Gumplinger, P., Howard, A., Hřivnáčová, I., Hwang, S., Incerti, S., Ivanchenko, A., Ivanchenko, V., Jones, F., Jun, S., Kaitaniemi, P., Karakatsanis, N., Karamitros, M., Kelsey, M., Kimura, A., Koi, T., Kurashige, H., Lechner, A., Lee, S., Longo, F., Maire, M., Mancusi, D., Mantero, A., Mendoza, E., Morgan, B., Murakami, K., Nikitina, T., Pandola, L., Paprocki, P., Perl, J., Petrović, I., Pia, M., Pokorski, W., Quesada, J., Raine, M., Reis, M., Ribon, A., Ristić Fira, A., Romano, F., Russo, G., Santin, G., Sasaki, T., Sawkey, D., Shin, J., Strakovsky, I., Taborda, A., Tanaka, S., Tomé, B., Toshito, T., Tran, H., Truscott, P., Urban, L., Uzhinsky, V., Verbeke, J., Verderi, M., Wendt, B., Wenzel, H., Wright, D., Wright, D., Yamashita, T., Yarba, J. & Yoshida, H. (2016). Nucl. Instrum. Methods Phys. Res. A, 835, 186–225. Web of Science CrossRef CAS Google Scholar
Andersen, A. H. & Kak, A. C. (1984). Ultrason. Imaging, 6, 81–94. CrossRef CAS PubMed Google Scholar
Babin, P., Della Valle, G., Chiron, H., Cloetens, P., Hoszowska, J., Pernot, P., Réguerre, A.L., Salvo, L. & Dendievel, R. (2006). J. Cereal Sci. 43, 393–397. Web of Science CrossRef Google Scholar
Barrett, J. F. & Keat, N. (2004). RadioGraphics, 24, 1679–1691. Web of Science CrossRef PubMed Google Scholar
Bertels, J., Eelbode, T., Berman, M., Vandermeulen, D., Maes, F., Bisschops, R. & Blaschko, M. B. (2019). Proceedings of the International Conference on Medical Image Computing and ComputerAssisted Intervention, pp. 92–100. Springer. Google Scholar
Brun, F., Mancini, L., Kasae, P., Favretto, S., Dreossi, D. & Tromba, G. (2010). Nucl. Instrum. Methods Phys. Res. A, 615, 326–332. Web of Science CrossRef CAS Google Scholar
Ching, D. J. & Gürsoy, D. (2017). J. Synchrotron Rad. 24, 537–544. Web of Science CrossRef IUCr Journals Google Scholar
Dagum, L. & Menon, R. (1998). IEEE Comput. Sci. Eng. 5, 46–55. Web of Science CrossRef Google Scholar
De Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, K. J., Ludwig, W., Mancini, L., Marone, F., Mokso, R., Pelt, D. M., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004. Web of Science CrossRef Google Scholar
De Chiffre, L., Carmignato, S., Kruth, J.P., Schmitt, R. & Weckenmann, A. (2014). CIRP Annals, 63, 655–677. Web of Science CrossRef Google Scholar
Der Sarkissian, H., Lucka, F., van Eijnatten, M., Colacicco, G., Coban, S. B. & Batenburg, K. J. (2019). Sci. Data, 6, 215. Web of Science CrossRef PubMed Google Scholar
Ekman, A., Weinhardt, V., Chen, J.H., McDermott, G., Le Gros, M. A. & Larabell, C. (2018). J. Struct. Biol. 204, 9–18. Web of Science CrossRef PubMed Google Scholar
Elfving, T., Hansen, P. C. & Nikazad, T. (2012). SIAM J. Sci. Comput. 34, A2000–A2017. Google Scholar
Etmann, C., Ke, R. & Schönlieb, C.B. (2020). Proceedings of the 2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP), 21–24 September 2020, Espoo, Finland, pp. 1–6. IEEE. Google Scholar
Evans, L. M. (2019). Xray tomography image data of a graphite foam block (KFoam) and tortuosity analysis, https://doi.org/10.5281/zenodo.3532935. Google Scholar
Evans, L. M., Margetts, L., Lee, P., Butler, C. & Surrey, E. (2019). Carbon, 143, 542–558. Web of Science CrossRef CAS Google Scholar
Faragó, T., Mikulík, P., Ershov, A., Vogelgesang, M., Hänschke, D. & Baumbach, T. (2017). J. Synchrotron Rad. 24, 1283–1295. Web of Science CrossRef IUCr Journals Google Scholar
Folk, M., Heber, G., Koziol, Q., Pourmal, E. & Robinson, D. (2011). Proceedings of the EDBT/ICDT 2011 Workshop on Array Databases (AD'11), 21–25 March 2011, Uppsala, Sweden, pp. 36–47. Google Scholar
Foster, B., Bagci, U., Mansoor, A., Xu, Z. & Mollura, D. J. (2014). Comput. Biol. Med. 50, 76–96. Web of Science CrossRef PubMed Google Scholar
Ganguly, P. S., Pelt, D. M., Gürsoy, D., de Carlo, F. & Batenburg, K. J. (2021). arXiv:2103.08288. Google Scholar
GarcíaMoreno, F., Kamm, P. H., Neu, T. R. & Banhart, J. (2018). J. Synchrotron Rad. 25, 1505–1508. Web of Science CrossRef IUCr Journals Google Scholar
Goo, H. W. & Goo, J. M. (2017). Kor. J. Radiol. 18, 555–569. Web of Science CrossRef Google Scholar
GuerquinKern, M., Lejeune, L., Pruessmann, K. P. & Unser, M. (2012). IEEE Trans. Med. Imaging, 31, 626–636. Web of Science CAS PubMed Google Scholar
Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188–1193. Web of Science CrossRef IUCr Journals Google Scholar
Hämäläinen, K., Harhanen, L., Kallonen, A., Kujanpää, A., Niemi, E. & Siltanen, S. (2015). arXiv:1502.04064. Google Scholar
Hangai, Y., Takahashi, K., Yamaguchi, R., Utsunomiya, T., Kitahara, S., Kuwazuru, O. & Yoshikawa, N. (2012). Mater. Sci. Eng. A, 556, 678–684. Web of Science CrossRef CAS Google Scholar
Hendriksen, A. A., Pelt, D. M. & Batenburg, K. J. (2020). IEEE Trans. Comput. Imaging, 6, 1320–1335. Web of Science CrossRef Google Scholar
Hendriksen, A. A., Pelt, D. M., Palenstijn, W. J., Coban, S. B. & Batenburg, K. J. (2019). Appl. Sci. 9, 2445. Web of Science CrossRef Google Scholar
Iassonov, P., Gebrenegus, T. & Tuller, M. (2009). Water Resour. Res. 45, w09415. Web of Science CrossRef Google Scholar
Jin, K. H., McCann, M. T., Froustey, E. & Unser, M. (2017). IEEE Trans.Image Processing, 26, 4509–4522. Web of Science CrossRef Google Scholar
Jørgensen, J. S., Coban, S. B., Lionheart, W. R., McDonald, S. A. & Withers, P. J. (2017). Meas. Sci. Technol. 28, 124005. Google Scholar
Kak, A. C., Slaney, M. & Wang, G. (2002). Med. Phys. 29, 107–107. CrossRef Google Scholar
Kazantsev, D., Pickalov, V., Nagella, S., Pasca, E. & Withers, P. J. (2018). SoftwareX, 7, 150–155. Web of Science CrossRef Google Scholar
Kazantsev, D., Thompson, W. M., Lionheart, W. R. B., Van Eyndhoven, G., Kaestner, A. P., Dobson, K. J., Withers, P. J. & Lee, P. D. (2015). Inverse Probl. 9, 447–467. Web of Science CrossRef Google Scholar
Kernighan, B. W. & Ritchie, D. M. (1988). The C Programming Language, 2nd ed. Prentice Hall. Google Scholar
Lam, S. K., Pitrou, A. & Seibert, S. (2015). Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC (LLVM '15), 15 November 2015, Austin, TX, USA. New York: Association for Computing Machinery. Google Scholar
Leuschner, J., Schmidt, M., Baguer, D. O. & Maaß, P. (2019). arXiv:1910.01113. Google Scholar
Liu, Z., Bicer, T., Kettimuthu, R., Gursoy, D., De Carlo, F. & Foster, I. (2020). J. Opt. Soc. Am. A, 37, 422–434. Web of Science CrossRef Google Scholar
Maire, E., Le Bourlot, C., Adrien, J., Mortensen, A. & Mokso, R. (2016). Int. J. Fract, 200, 3–12. Web of Science CrossRef CAS Google Scholar
Marchant, D., Munk, R., Brenne, E. O. & Vinter, B. (2020). Proceedings of the 2020 IEEE/ACM 2nd Annual Workshop on ExtremeScale ExperimentintheLoop Computing (XLOOP), 12 November 2020, Atlanta, GA, USA, pp. 23–28. Google Scholar
Marone, F. & Stampanoni, M. (2012). J. Synchrotron Rad. 19, 1029–1037. Web of Science CrossRef CAS IUCr Journals Google Scholar
Matsumoto, M. & Nishimura, T. (1998). ACM Trans. Model. Comput. Simul. 8, 3–30. CrossRef Google Scholar
McCollough, C. H., Bartley, A. C., Carter, R. E., Chen, B., Drees, T. A., Edwards, P., Holmes, D. R. III, Huang, A. E., Khan, F., Leng, S., McMillan, K. L., Michalak, G. J., Nunez, K. M., Yu, L. & Fletcher, J. G. (2017). Med. Phys. 44, e339–e352. Web of Science CrossRef CAS PubMed Google Scholar
Metscher, B. D. (2009). BMC Physiol. 9, 11. Google Scholar
Miqueles, E. X., Rinkel, J., O'Dowd, F. & Bermúdez, J. S. V. (2014). J. Synchrotron Rad. 21, 1333–1346. Web of Science CrossRef CAS IUCr Journals Google Scholar
Mittone, A., Manakov, I., Broche, L., Jarnias, C., Coan, P. & Bravin, A. (2017). J. Synchrotron Rad. 24, 1226–1236. Web of Science CrossRef IUCr Journals Google Scholar
Mohan, K. A., Venkatakrishnan, S., Gibbs, J. W., Gulsoy, E. B., Xiao, X., De Graef, M., Voorhees, P. W. & Bouman, C. A. (2015). IEEE Trans. Comput. Imaging, 1, 96–111. Google Scholar
Moosmann, J., Hofmann, R. & Baumbach, T. (2011). Opt. Express, 19, 12066–12073. Web of Science CrossRef PubMed Google Scholar
Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591. Web of Science PubMed Google Scholar
Nelder, J. A. & Mead, R. (1965). Comput. J. 7, 308–313. CrossRef Web of Science Google Scholar
Nievergelt, J. & Preparata, F. P. (1982). Commun. ACM, 25, 739–747. CrossRef Web of Science Google Scholar
Nikitin, V. V., Carlsson, M., Andersson, F. & Mokso, R. (2019). IEEE Trans. Comput. Imaging, 5, 409–419. Web of Science CrossRef Google Scholar
Pan, X., Sidky, E. Y. & Vannier, M. (2009). Inverse Probl. 25, 123009. Web of Science CrossRef Google Scholar
Pelt, D. M. (2020). dmpelt/foam_ct_phantom, https://doi.org/10.5281/zenodo.3726909. Google Scholar
Pelt, D. M., Batenburg, K. J. & Sethian, J. A. (2018). J. Imaging, 4, 128. Web of Science CrossRef Google Scholar
Perciano, T., Ushizima, D., Krishnan, H., Parkinson, D., Larson, N., Pelt, D. M., Bethel, W., Zok, F. & Sethian, J. (2017). J. Synchrotron Rad. 24, 1065–1077. Web of Science CrossRef IUCr Journals Google Scholar
Pugh, W. (1990). Commun. ACM, 33, 668–676. CrossRef Web of Science Google Scholar
Raufaste, C., Dollet, B., Mader, K., Santucci, S. & Mokso, R. (2015). EPL (Europhysics Lett), 111, 38004. Google Scholar
Ravishankar, S., Ye, J. C. & Fessler, J. A. (2020). Proc. IEEE, 108, 86–109. Web of Science CrossRef Google Scholar
Renders, J., Sijbers, J. & De Beenhouwer, J. (2020). Proceedings of the 6th International Conference on Image Formation in Xray Computed Tomography, 3–7 August 2020, Regensburg, Germany, pp. 154–157. Google Scholar
Roux, S., Hild, F., Viot, P. & Bernard, D. (2008). Composites Part A, 39, 1253–1265. Web of Science CrossRef Google Scholar
Salvo, L., Cloetens, P., Maire, E., Zabler, S., Blandin, J. J., Buffière, J.Y., Ludwig, W., Boller, E., Bellet, D. & Josserond, C. (2003). Nucl. Instrum. Methods Phys. Res. B, 200, 273–286. Web of Science CrossRef CAS Google Scholar
Santos Rolo, T. dos, Ershov, A., van de Kamp, T. & Baumbach, T. (2014). Proc. Natl Acad. Sci. USA, 111, 3921–3926. Web of Science PubMed Google Scholar
Scales, J. A. (1987). Geophysics, 52, 179–185. CrossRef Web of Science Google Scholar
Segars, W. P. & Tsui, B. M. (2009). Proc. IEEE, 97, 1954–1968. Web of Science CrossRef Google Scholar
Shepp, L. A. & Logan, B. F. (1974). IEEE Trans. Nucl. Sci. 21, 21–43. CrossRef Web of Science Google Scholar
Singh, K., Menke, H., Andrew, M., Rau, C., Bijeljic, B. & Blunt, M. J. (2018). Sci. Data, 5, 180265. Web of Science CrossRef PubMed Google Scholar
Toft, P. A. (1996). The Radon TransformTheory and Implementation. PhD thesis, Technical University of Denmark, Denmark. Google Scholar
Van Nieuwenhove, V., De Beenhouwer, J., Vlassenbroeck, J., Brennan, M. & Sijbers, J. (2017). Opt. Express, 25, 19236–19250. Web of Science CrossRef PubMed Google Scholar
Van Rossum, G. & Drake, F. L. (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace. Google Scholar
Wang, Z., Simoncelli, E. P. & Bovik, A. C. (2003). Proceedings of the ThirtySeventh Asilomar Conference on Signals, Systems and Computers, 2003, 9–12 November 2003, Pacific Grove, CA, USA, Vol. 2, pp. 1398–1402. IEEE. Google Scholar
Yang, X., De Carlo, F., Phatak, C. & Gürsoy, D. (2017). J. Synchrotron Rad. 24, 469–475. Web of Science CrossRef IUCr Journals Google Scholar
Yang, X., Kahnt, M., Brückner, D., Schropp, A., Fam, Y., Becher, J., Grunwaldt, J.D., Sheppard, T. L. & Schroer, C. G. (2020). J. Synchrotron Rad. 27, 486–493. Web of Science CrossRef IUCr Journals Google Scholar
Yu, Z., Noo, F., Dennerlein, F., Wunderlich, A., Lauritsch, G. & Hornegger, J. (2012). Phys. Med. Biol. 57, N237–N252. Web of Science CrossRef PubMed Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.