poly-DART: A discrete algebraic reconstruction technique for polychromatic X-ray CT

The discrete algebraic reconstruction technique (DART) is a tomographic method to reconstruct images from X-ray projections in which prior knowledge on the number of object materials is exploited. In monochromatic X-ray CT (e.g., synchrotron), DART has been shown to lead to high-quality reconstructions, even with a low number of projections or a limited scanning view. However, most X-ray sources are polychromatic, leading to beam hardening effects, which significantly degrade the performance of DART. In this work, we propose a new discrete tomography algorithm, poly-DART, that exploits sparsity in the attenuation values using DART and simultaneously accounts for the polychromatic nature of the X-ray source. The results show that poly-DART leads to a vastly improved segmentation on polychromatic data obtained from Monte Carlo simulations as well as on experimental data, compared to DART. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
X-ray computed tomography (CT) is a well known imaging method in which the interior of an object is reconstructed from a set of X-ray radiographs. High quality CT imaging generally requires hundreds or even thousands of radiographs acquired in a circular orbit with a large angular range. In many applications of X-ray CT, however, there is a need to lower the amount of radiographs to be acquired, either to reduce dose [1][2][3] or to shorten the acquisition time [4,5]. In time critical processes [6][7][8], for example, the object may change during the scan so reducing the amount of projections could lead to an improved image quality and/or a higher time resolution. Industrial CT for quality control, specifically in-line inspection, again has a need for reduced scanning times and in the case of larger parts, reconstruction from a limited angular range [1,9].
Unfortunately, lowering the number of projections to reduce dose or decrease acquisition time directly corresponds to a reduction of the sampling density in projection domain, which, if no precautions are taken, leads to undersampling artefacts in the reconstructed CT image. Indeed, when only few projections are available, the inverse problem of image reconstruction becomes severely underdetermined, causing streaking artefacts in the reconstructed image. By including prior knowledge about the sample, this underdetermination can be reduced to increase the quality of the reconstructed image.
One of the ways prior knowledge about the object to be imaged can be enforced, is Discrete Tomography, in which the number of possible materials in the sample is assumed to be known (and limited). For an overview of the mathematical background of discrete tomography, and several applications in medicine, we refer to [10]. The Discrete Algebraic Reconstruction Technique (DART) [11] is a practical discrete tomography algorithm which, if the sample consists of only a few materials, has been shown to generate reconstructions of much better quality than those from conventional reconstruction methods. The method has been applied in a variety of imaging domains such as X-ray CT [12][13][14][15], Electron tomography [16][17][18], Optical Diffraction Tomography [19], and Magnetic resonance imaging [20].
Many variations of the DART algorithm have been researched. SDART [21] and DART-ALBM [22], for example, focus on improving robustness against noise. rmwDART [23] specifically targets reduction of missing wedge artefacts. In PDART, the discreteness assumption on the whole sample is relaxed to local areas, making DART applicable to samples that are only partially discrete [18,24]. In PDM-DART, the grey values corresponding to the different materials are estimated in an automated way [25]. ADART [26] gradually reduces the number of unknowns in each iteration, and MDART [27] follows a multiresolution approach. EOD-DART [9], EDART [28] and TVR-DART [29] add additional prior knowledge about the sample to further improve the reconstruction quality.
All of the above DART based methods rely on the same linearized acquisition model as most classical techniques: it assumes that the log-corrected normalized projection is the sum of attenuation values along a ray. This is an adequate model for monochromatic but not polychromatic X-ray sources. As the X-ray beam travels through the sample, the low energy photons are attenuated more easily than those with high energy. This results in the absorption along a ray from a polychromatic source being a non-linear function of sample thickness. Reconstructing such a dataset with a linear reconstruction technique will lead to well known beam hardening (or cupping) artefacts [30]. At the same time, a single material no longer has a single defined attenuation value, as this value is energy dependent. Beam hardening artefacts and the inability to select a clear singular attenuation value for each material, lead to inaccurate reconstructions when DART is used in combination with polychromatic projections. One of the current methods for reducing beam hardening artefacts consists of placing metal filters in front of the source, to pre-harden the beam. However, the pre-hardened spectrum is still polychromatic. Moreover, the filtering decreases the number of photons available for imaging and hence lowers the SNR of the reconstructed image.
We propose a new discrete tomographic method for polychromatic X-ray CT, which incorporates a polychromatic forward model. Such an approach has previously successfully been introduced to extend a classical algebraic reconstruction method (SART [31]) to handle polychromatic X-ray data, referred to as pSART (polychromatic Simultaneous Algebraic Reconstruction Technique) [32,33]. In this paper, we propose a discrete reconstruction method that combines the best of both worlds: accounting for polychromaticity while exploiting the strength of DART to substantially reduce undersampling artefacts.
The paper starts with a short overview of DART. Next, we explain the poly-DART algorithm and the polychromatic model that is used. Finally, different discrete reconstructions from polychromatic datasets are compared to demonstrate the improvements of the proposed method. These datasets include Monte Carlo simulated data using the GATE framework [34], as well as an experimental dataset.

Method
In this section, we first review the DART method. Next, since DART uses SIRT as an underlying reconstruction method, the pSIRT algorithm is introduced, followed by the proposed poly-DART algorithm.

Discrete reconstruction: DART
The DART algorithm is an efficient reconstruction method for discrete tomography, of which a complete discussion can be found in [11]. We briefly recall the different steps in the DART algorithm: 1. Initial reconstruction: An algebraic reconstruction method is performed to obtain an initial reconstruction V 0 .

Segmentation:
The reconstruction V 0 is segmented based on the a priori knowledge of the grey levels.
3. Masking: Sets of free and fixed pixels are selected. The free pixels encompass all boundary pixels, i.e. the pixels of which not all neighbouring pixels have the same grey value, and a small percentage of randomly chosen pixels. The fixed pixels are all pixels that are not free.

4.
Reconstruction: Perform a number of iterations of an algebraic reconstruction method (e.g., SIRT), on the set of free pixels only. The fixed pixels are kept at their segmented value. This results in an intermediate reconstruction V.

5.
Smoothing: Apply Gaussian smoothing to weaken the harsh fluctuations that can occur due to the masked reconstruction.
6. Iteration: Set V 0 = V and repeat steps 2-5 until a stopping criterion is met.
The DART algorithm has a number of parameters that will influence performance and are chosen heuristically. The most important ones are the percentage of random pixels chosen to be free in the masking step and the different amounts of iterations (initial, inner and outer).

Polychromatic SIRT (pSIRT)
The well-known Beer-Lambert law states that for a monochromatic beam the measured intensity due to absorption along a single ray L r , i.e. the projection data along the ray from the source to position r on the detector, can be modelled as follows: with I 0 the intensity of the source, µ the attenuation coefficient at every point in the object space, and x is a vector denoting a single point in the object space. This model is linearizeable by taking the natural logarithm of both sides of the equation. When a polychromatic source is assumed, the forward projection model is given by the integral over all possible energies the X-ray source emits, between some minimal and maximal energy, min and max respectively: where S( ) is the normalized source spectrum and D( ) is the detector response function. Define w( ) = I 0 D( )S( ), then discretization over E energy levels and M different materials yields the following projection value at detector pixel r [32]: with w(e) the effective spectrum, l r,m the distance that ray L r travels through material m, and µ m (e) the attenuation coefficient of material m at energy e. By log-correcting and normalizing the projectionP L r from Eq. (3) one can define a normalized polychromatic projection operator polyProj(·): with V the object that is projected.
Next, the necessary details of the update steps of pSIRT are explained. In matrix form, the kth update step of the SIRT algorithm is described as follows: with A the system matrix of the acquisition geometry, C and R diagonal matrices with the inverse sums of the columns and rows of A, respectively and s the sinogram. The multiplication AV (k) is the discretization of the monochromatic forward projection from Eq. (1), after taking the natural logarithm. By swapping out this monochromatic forward projection for the polychromatic projection in Eq. (4), the polychromaticity of the source and energy dependence of the attenuation coefficients is modelled, and one arrives at the update step of the pSIRT algorithm, which is the analogous extension of pSART as SIRT is to SART. The extra prior information needed for pSIRT involves the spectrum: E and w(·) and the different materials in the sample: µ m (·) and M.
Note that this prior knowledge is also a prerequisite of the DART algorithm.
For the forward projection of V (k) , the line-lengths l (k) r,m need to be calculated as follows. First, a reference energy ref is chosen, which provides reference attenuation values µ m = µ m ( ref ). In each step, the reconstruction represents the monochromatic attenuation map at energy level ref [32]. Next, for each pixel v in the current reconstruction V (k) , with value t v ∈ [µ i , µ i+1 ], the fractions in the following decomposition are calculated: This essentially models each pixel v as the mixture of the two materials whose attenuation values at the reference energy are closest to the current value of the pixel. These percentages are grouped, per material, in the masks M (k) m . The l (k) r,m are now found as the values in the product AM (k) m .

Polychromatic DART (poly-DART)
We propose a new algorithm based on the principles of DART and pSIRT, which aims at combining the benefits of both methods. In the first phase of the algorithm, an initial reconstruction V (0) is generated by performing a small number of pSIRT iterations. From this reconstruction and the attenuation values at the reference energy, the optimal grey values and thresholds are estimated by minimizing the projection difference (i.e., the difference between the polychromatic forward projection and the projection data) in a mean squares sense. This is a polychromatic extension of the projection distance minimization (PDM) algorithm of van Aarle et al [25]. It allows to correct for small errors in the assumed attenuation coefficients and/or in the estimated spectrum. Next, the initial reconstruction is segmented with the estimated segmentation parameters. Based on the segmented result, the pixels are classified into boundary or non-boundary pixels. All boundary pixels and a set of random pixels are marked 'free' and will be updated in the next reconstruction step. The other pixels remain fixed to their segmented value. Then, the line lengths through each material for the fixed pixels are precomputed to be used in the reconstruction, as these will not change until the next segmentation step. Masked pSIRT iterations are then performed, which only update the free pixels. Furthermore, these local iterations are relaxed by a parameter λ, equal to the percentage of free pixels, as the much lower amount of free pixels (compared to a conventional update over all pixels) would otherwise lead to high fluctuations in each step of the reconstruction. Hence the (internal) poly-DART update step becomes: Lastly, this reconstruction is segmented again and the previous steps are repeated until a stop criterion is met. A flowchart of poly-DART can be found in Fig. 1. The algorithm was implemented in Matlab, making use of the DART framework of the ASTRA toolbox [35].

Monte Carlo simulations
To validate the proposed algorithm, two simulation experiments were set up, with the phantoms shown in Figs. 2(a) and 2(b). Both phantoms consist of two materials, Plexiglass and aluminum, suspended in air. The attenuation values at different energy levels were used from the National Institute of Standards and Technology [36]. To avoid inverse crime, the phantom was analytically defined in the simulation framework, GATE [34]. For both phantoms, 300 fan-beam polychromatic projections were created using GATE. The simulated source was a 75 kVp X-ray source with a tungsten anode, of which the spectrum is shown in Fig. 2(c). Spekcalc was used to generate the source spectrum [37]. The projection angles were chosen in the interval [0, 2π] using golden angle sampling, which means that subsequent projections are 1+ √ 5 2 π radians apart from each other. The detector pixel size was 0.375 mm and the scanner set up had a magnification factor of 4. The reconstructed pixel size was 0.09375 mm. The detector was modelled as a photon counting detector with 400 pixels made of silicon, with a minimal activation energy of 5 keV. The detector was idealized in the sense that no cross-talk or blurring between pixels was present. As the projections were created using a Monte Carlo method, Poisson noise is present in the images. We simulated different noise levels by varying the number of photons from 4 * 10 5 to 2 * 10 7 photons for each projection, which is on average between 1000 and 50000 photons per pixel. For the experiments where varying noise levels were not taken into account, we chose the dataset with 10000 photons per detector pixel. To further simulate a real experiment, the tube spectrum (S) and detector response (D) were estimated with two step wedges. These wedges, one in Plexiglass and one in aluminium, were also defined in GATE. Multiple projections were taken orthogonal to the steps and then averaged to reduce the effects of noise in the measurements. We used a maximum likelihood expectation maximization (MLEM) algorithm to estimate the product S( )D( ) of the tube spectrum and the detector response from these measurements. The four characteristic peaks of a tungsten anode source spectrum were added manually to the initial guess for the MLEM iterations.
From the simulated polychromatic projection data, images were reconstructed with the following methods: pSIRT, SIRT, DART with manually optimized grey levels and poly-DART with automatically optimized grey levels. For the selection of the grey values for DART, the values with the lowest reconstruction error in the case of full angular sampling were chosen. Poly-DART was also compared to (segmented) SIRT and pSIRT. The segmentation for SIRT was performed with the same global threshold as DART. The segmentation for pSIRT was performed with the first set of estimated segmentation parameters that were obtained from the PDM step in the poly-DART algorithm. These automatically optimized grey values are only dependent on the initial pSIRT reconstruction and the projection data, and are more accurate than using the attenuation values at the reference energy.
All poly-DART reconstructions were performed with 20% random free pixels, 50 initial and 5 inner iterations of pSIRT. The same settings were used for DART. The percentage of free pixels was chosen as advised in [11]. For the polychromatic reconstruction methods, the spectrum shown in Fig. 2(c) was rebinned to 70 bins. A reference energy of 55 keV was chosen as optimal for pSIRT in terms of convergence speed an reliability [32]. The same reference energy was used for poly-DART. When testing the influence of the number of angles and the level of noise, 2000 (p)SIRT iterations were performed for all methods and the minimal achieved error over these iterations was then retained. That is, for both converging and diverging methods the best reconstruction (minimal error) is retained.

Experimental dataset
Following the Monte Carlo simulations, a similar test was performed on an experimental dataset. The used phantom, called the Barbapapa phantom [38], consists of Plexiglass with inserted aluminum rods and is shown in Fig. 3. A phantom with the same materials as in the simulation study was used, to more easily compare the simulation results to the experimental results. A tube voltage of 130 kVp was employed. The phantom was scanned over 2400 equiangular projections in the interval [0, 2π).
To estimate the spectrum, the same technique as in the simulation case was used. For this experiment, a PVC step wedge with 11 steps, ranging from approximately 1 mm to 18 mm in thickness, was scanned. Next, the same MLEM spectrum estimation method as in the Monte Carlo case was performed, with 130 keV as maximum energy used as a boundary condition. From the central slice of the Barbapapa dataset, images were reconstructed with the same methods: pSIRT segmented, SIRT segmented, DART, and poly-DART. All poly-DART reconstructions were performed with 10% random free pixels, 50 initial and 10 inner iterations of pSIRT. The same settings were used for DART. A reference energy of 55 keV was chosen for pSIRT as well as for poly-DART. Since for real data no ground truth image is available, a pSIRT reconstruction from 2400 equiangular projection within [0, 2π] was segmented using Otsu thresholding [39] and the result was regarded as the ground truth reconstruction. When testing the influence of the variable number of angles, 2000 (p)SIRT iterations were performed for all methods and the minimal achieved error over these iterations was retained. The error measure is the rNMP calculated based on the golden standard segmented reconstruction. For both the simulated and real data, all parameters (e.g. percentage of free pixels, number of inner iterations and reference energy level) were either chosen empirically or in accordance with the literature. All reconstructions were performed using the ASTRA toolbox [35].

Results
To quantify the performance of poly-DART against DART, segmented SIRT and segmented pSIRT, we computed the relative number of misclassified pixels or rNMP [25]. The rNMP is equal to the amount of pixels in the reconstruction that have a different grey value compared to those of the corresponding ground truth image, divided by the total number of non-background pixels in the ground truth image.

Pie phantom
In Fig. 4(a), the rNMP of the different methods is shown as a function of the number of (p)SIRT iterations, for 300 projection angles between 0 and 2π. The rNMP for (p)SIRT was calculated at every iteration, whereas the rNMP for (poly-)DART was calculated at every DART iteration, i.e. every 5 (p)SIRT iterations. From Fig. 4(a), one can observe that both DART and poly-DART converge. However, poly-DART does so to a solution with a significantly lower error, whereas DART increases the error from its start point to its converged solution. The continuous methods, SIRT and pSIRT, have a best solution after a low amount of iterations and start to diverge from that point onwards. This is a known problem for continuous methods, as they tend to overfit to the noise in the projection data.
Next, the effect of varying amounts of projection angles on the rNMP was studied. These results are shown in Fig. 4(b). Again, poly-DART outperforms the other methods in terms of rNMP, reaching a lower rNMP for each number of projection angles. This suggests that poly-DART benefits from both the beam hardening correction property of pSIRT and the imposed discreteness. Note that because the minimal error for the methods over all iterations is compared, this is the best case scenario for SIRT and pSIRT, as they reach a minimum at some a priori unknown point. For down to about 150 projection angles, poly-DART and pSIRT have a comparable error. However, at lower amounts of projections, poly-DART outperforms pSIRT. A comparison of the different reconstruction methods for 40 projections can be found in Fig. 5, where each method was run for 250 (p)SIRT iterations. From this figure, it can be observed that the poly-DART algorithm shows the least reconstruction artefacts, with the other methods showing beam hardening artefacts and streaks due to the low amount of projection angles. Furthermore, the homogeneous interior of the object is well reconstructed with poly-DART, while pSIRT shows many misclassified (black) pixels. Finally, we tested the influence of noise on the reconstructions. Different levels of noise were simulated by varying the amount of photons emitted by the source. The amount of photons ranges from, on average, 1000 to 50000 photons per detector pixel. The minimal achieved rNMP over 2000 iterations is plotted in function of the number of photons. In Figs. 6(a) and 6(b) this rNMP is plotted for all methods for 10 and 100 projections, respectively. From these plots, it can be observed that poly-DART consistently outperforms the other techniques, especially for low amounts of photons.

Derenzo phantom
The rNMP is not a suitable error measure for the Derenzo phantom, as every misclassified pixel contributes equally to the error. Of more interest in this phantom is to what extent the different cylinders can be discerned. Therefore, we visually show reconstructions for all methods as a function of different numbers of projections to illustrate the level of visibility. These reconstructions, all after 450 (p)SIRT iterations, are shown in Fig. 7. From Fig. 7, similar reconstruction results of the Derenzo phantom as for the Pie phantom can be observed. Even when a large number of projections is available, two artefacts arise as before, beam hardening and overfitting to noise. Both SIRT and DART suffer from the beam hardening artefacts. For 50 projections and less, these beam hardening artefacts, coupled with undersampling, make it increasingly more difficult to see the cylinders at any resolution. Neither the poly-DART nor pSIRT reconstructions show significant beam hardening artefacts. Even with a large number of projections, pSIRT shows artefacts due to overfitting to noise, though the cylindrical features can still be discerned at every resolution. As the number of projections decreases, both the size and the placement of the cylinders become harder to distinguish from the pSIRT reconstruction. The poly-DART reconstructions are more accurate for each number of projections, most notably when less than 50 projections are present. Both position and size of the cylinders more closely match those of the phantom image in Fig. 2(b).

Experimental dataset
Based on the experimental data, two experiments were performed. First, the effect of varying amounts of projection angles on the rNMP was studied. As can be observed from Fig. 8(a), poly-DART outperforms the other methods in terms of misclassified pixels, reaching a lower rNMP at any number of projection angles. The plot looks similar to the one in Fig. 4(b), indicating  Secondly, we studied the effect of a missing wedge in the projection data, i.e. a dense sampling of the object, but over a limited range [α, π − α] ∪ [π + α, 2π − α], α ∈ R. For the full sampling, 400 equiangular projections were taken from the dataset. From the previous experiment, it is clear that this is a sufficiently dense sampling for all methods to generate high quality reconstructions. Next, a number of projections were deleted symmetrically around both 0 and π. For the created missing wedge datasets in this experiment, α ∈ [0, π 4 ] was used. As before, the minimal error over 2000 iterations was plotted in Fig. 8(b). The DART algorithm outperforms SIRT from about 15 • of missing wedge and outperforms pSIRT after around 25 • . The proposed poly-DART algorithm reaches the lowest reconstruction error of all the tested methods at any of missing wedge angle. This shows that the imposed prior knowledge is able to compensate for the missing data.

Conclusion
Many objects consist of a limited number of materials. Using discrete tomography, this prior knowledge can be exploited in the reconstruction of images from X-ray projection data to reduce undersampling artefacts. Current discrete tomography methods, however, do not account for polychromaticity of X-ray sources, leading to various reconstruction artefacts and limiting their applications. In this paper, poly-DART was proposed, a discrete tomography method that exploits sparseness in attenuation values, while taking a polychromatic projection model into account. Reconstruction experiments on both simulated and experimental data from polychromatic sources revealed that poly-DART results in substantially improved image reconstruction quality compared to DART or segmented versions of SIRT or pSIRT for polychromatic X-ray data. This allows DART, which has been successfully used in a monochromatic setting in different applications, to be extended to data acquired with polychromatic lab sources.