DIRA-3D—a model-based iterative algorithm for accurate dual-energy dual-source 3D helical CT

Quantitative dual-energy computed tomography may improve the accuracy of treatment planning in radiation therapy. Of special interest are algorithms that can estimate material composition of the imaged object. One example of such an algorithm is the 2D model-based iterative reconstruction algorithm DIRA. The aim of this work is to extend this algorithm to 3D so that it can be used with cone-beams and helical scanning. In the new algorithm, the parallel FBP method was replaced with the approximate 3D FBP-based PI-method. Its performance was tested using a mathematical phantom consisting of six ellipsoids. The algorithm substantially reduced the beam-hardening artefact and the artefacts caused by approximate reconstruction after six iterations. Compared to Alvarez-Macovski’s base material decomposition, DIRA-3D does not require geometrically consistent projections and hence can be used in dual-source CT scanners. Also, it can use several tissue-specific material bases at the same time to represent the imaged object.


Introduction
The model-based iterative reconstruction algorithm DIRA-3D was designed to address the need for better estimation of material parameters needed for radiation transport and dose calculations in brachytherapy with low energy photons and proton therapy. In the following paragraphs we explain why such a tool is needed and what solutions are already available.
In conventional external beam radiotherapy, high energy photons from linear accelerators are used. Compton scattering is the dominant interaction process in this energy range and thus CT images showing electron densities are sufficient. The situation is different when using low energy photons (<50 keV) for brachytherapy or with protons. With low energy photons, atomic composition is important since the photoelectric effect dominates even in soft tissues of low atomic number. With protons both the mean excitation energy (I-value) and the probability of nuclear reactions also depend on the atomic number. Today, in brachytherapy, dose calculations are approximate, since they are performed as if the body is a large water phantom without anatomical detail or size restrictions. All tissue absorbed doses are calculated as absorbed dose to water in water (Rivard et al 2004). The need for better characterization of the tissues has already been stated and is currently a bottle-neck for the use of model-based dose calculation methods for this modality, see e.g. (Beaulieu et al 2012, van Elmpt et al 2016. Computed tomography reconstructs linear attenuation coefficients (LACs). To suppress their strong dependence on photon energy, they are presented to the end-user in the form of CT numbers. The x-ray tube produces photon beams with a spectrum of energies and the scanned object normally contains multiple materials. The strong energy dependence of LACs causes beam-hardening artefacts.
In single-energy CT (SECT), the scan is performed using one x-ray tube voltage only. In dual-energy CT (DECT), two sets of data are obtained, each corresponding to a different tube voltage. The main advantage of DECT compared to SECT stems from the possibility to obtain more information about the material via material decomposition. Another advantage is a better suppression of the beam hardening artefact. Computer simulations have shown that DECT (compared to SECT) can improve the accuracy of radiation treatment planning (Bär et al 2017) and some clinical applications are already available (Wohlfahrt et al 2017). However, the need for image reconstruction algorithms that can accurately characterize material properties that are of interest in radiotherapy still exists; more information is in the following text.
Base material decomposition can be performed in the projection domain (PD) or the image domain (ID) (Heismann et al 2012). A widely used PD method is the one by Alvarez and Macovski (1976), which decomposes projections into two components, either two base materials or (approximate) photoelectric and Compton components. The decomposed projections are then reconstructed via filtered backprojection (FBP). This procedure inherently eliminates beamhardening artefacts. Disadvantages are that it requires geometrically consistent projections and the two base materials represent the whole scanned object; organspecific base materials can only be used during postprocessing.
The ID methods work with reconstructed images. They use two-material or three-material decomposition, see for example (Bazalova et al 2008, Liu et al 2009. The latter makes an extra assumption on the conservation of partial volumes of individual constituents in the mixture. The ID methods can use organ specific base materials. This concept was extended by Mendonca et al (2014) to multi-material decomposition (MMD). Their method first produces two monoenergetic images via the Alvarez Macovski method with the water-iodine base material doublet. Then a three-material decomposition is applied, where the pixel-specific material bases are selected according to the values of reconstructed LACs. Inspired by this method, Long and Fessler (2014) proposed a statistical image reconstruction method, a PL (penalized-likelihood) method, where the PL cost function contains a negative log-likelihood term and an edge-preserving regularization term for each basis material. Compared to the MMD method, the PL method better reduces noise, streak artefacts, and cross-talk. It is, however, very computationally demanding and thus it was implemented in 2D only.
The material decomposition from inconsistent rays (MDIR) (Maaß et al 2011) is an iterative image reconstruction method that produces 2 base material images. It uses FBP and the original measured projections inside the iterative loop. It can be used when geometrically consistent raw data projections are not available, for instance for a dual-source DECT scanner.
The original DIRA algorithm (Malusek et al 2017) is a model-based iterative reconstruction method which uses FBP inside its iterative loop. It works simultaneously on projections and images (a combined PD-ID method) and, as MDIR, it can work with inconsistent rays. Within its iterative loop, it uses twoand three-material decomposition to tissue specific, user defined material doublets and triplets for the characterization of body tissues. Soft tissue is decomposed for instance into lipid, protein and water as suggested by Woodard and White (1986). The algorithm effectively reduces beam hardening artefacts. The output from the original DIRA is a stack of images, each containing mass fractions of individual base materials, representing one slice. From this multidimensional image, virtual monoenergetic images at any photon energy can be computed if desired.
Modern clinical CT systems use helical scanning to speed up the data acquisition in 3D imaging, which is also used to delineate tumour and risk organs for treatment planning. The 3D weighted filtered backprojection (WFBP) (Stierstorfer et al 2004) is used on Siemens helical CT-scanners. This algorithm is approximate only and consequently produces small artefacts. To improve the reconstruction result and to suppress noise in helical CT, iterative algorithms based on backprojection are often used, see for instance (Beister et al 2012). There are many other examples on approximate helical 3D FBP algorithms, one being the PI-method (Turbell 2001), which similarly to the WFBP performs semi-parallel rebinning, horizontal ramp-filtering and cone-beam backprojection and thus produces similar artefacts. An exact helical 3D FBP algorithm was developed by Katsevich (2002). Both the PI and Katsevich's methods do not utilize projection data outside the Tam window (see figure 1(b)), which results in lower SNR compared to the WFBP method. This makes PI and Katsevich's methods less favourable in clinical environment.
Pure PD methods using geometrically consistent rays like the method by Alvarez and Macovski can simply be extended to 3D geometries. Nevertheless, the iterative algorithms MDIR and DIRA were for the sake of simplicity implemented for 2D geometries only. In the paper describing MDIR, a few simplified 3D-geometry cases are investigated, though. Since helical 3D FBP algorithms are approximate (apart from Katsevich algoritm), it is not obvious that an extension from 2D to 3D will work. The iterative algorithm should be able to correct artefacts caused by approximate reconstruction and beam-hardening artefacts at the same time. (For brevity, the former are called approximatereconstruction artefacts in this work.) Also, it should be able to handle the long object problem, see section 2.3.3.
Computer simulations demonstrated that DIRA had the potential to address these issues. To make DIRA clinically applicable, it has to be extended from 2D to 3D helical geometry. The aim of this work is to develop such an extension and test it in a proof of concept implementation. The extension should be able to handle geometrically inconsistent projections from dual-source CT scanners (e.g., Siemens SOMATOM definition Force). Since the PI-method was readily available to the authors, a solution based on this method was preferred.
2. Theory 2.1. Material decomposition Two-material decomposition (2MD) mathematically decomposes a mixture to two base materials. It determines mass fractions, w 1 and w 2 , of the two base materials and the mass density, ρ, of the mixture. Three-material decomposition (3MD) mathematically decomposes a mixture to three base materials. It determines mass fractions, w 1 , w 2 and w 3 , and the mass density of the mixture, which is calculated as , where ρ k and w k are the mass density and mass fraction, respectively, of the kth material. In both methods, the mass fractions are normalized so that å = w 1 k k . More information on the resulting systems of linear equations is in (Malusek et al 2017). When other than the base materials are present in the mixture, the mass fractions can become negative. In this case, w 1 , K, w 3 should be seen as linear coefficients that weight the mass attenuation coefficients of base materials so that the resulting calculated LACs at the energies E 1 and E 2 equal the reconstructed LACs of the mixture.

Forward projection generation
The logarithm of attenuation, here referred to as the polyenergetic projection P, is calculated as where I and I 0 are the detector responses with and without, respectively, the imaged object. The intensity I 0 , is calculated for an ideal energy integrating detector as where E is the photon energy and N(E) is the energy spectrum of photons emitted from the x-ray tube. The intensity I is calculated as at energy E and ò dl L is a line integral through the object. This calculation is time consuming since the line integrals must be calculated for all energies in the energy spectrum. When the material decomposition is used, the line integrals can be calculated through volume fractions of the different base materials and the intensity I is where μ k is the LAC of the kth base material and l k is computed as where ρ k is the tabulated density of the kth base material, ρ(x, y, z) is the calculated density in voxel (x, y, z) and w k (x, y, z) is the mass fraction of the kth material in voxel (x, y, z). The density ρ and the mass fractions w k are obtained from the two-material or three-material decomposition. Note that the demanding line integral computation is performed only once for each base material (equation (5)), not for every photon energy in the spectrum as in (3).
A monoenergetic projection P E i , where E i , i=1,2 is a specific energy, is calculated as The specific energy E i is the energy at which the LAC for water, μ w , equals the energy-fluence weighted mean LAC for water The line integral l k in equation (6) can in a voxelized geometry be calculated using Joseph's method (Joseph 1982) extended to three dimensions (Turbell 2001).
2.3. The PI-method, an FBP method for helical CT 2.3.1. Helical cone-beam geometry Dual-source DECT scanners typically use cone-beam projections and helical trajectory of the two sources, see figure 1(a). In this figure, κ is the cone-angle, γ is the fan-angle and P represents the pitch of the helix (the height of one complete helix turn, measured parallel to the axis of the helix). The maximum value of the cone-angle, κ max , is called the cone-beam angle and the maximum value of the fan-angle, γ max , is called the fan-beam angle. The s-axis starts at the point C, the cross-section between the central ray and the zaxis, and is aligned with the z-axis. The x-ray sources are visualized as red and blue spherical markers. The detectors are cylindrical, their axes are parallel with the z-axis and intersect corresponding x-ray source positions. The source-detector assemblies rotate around the z-axis in a helical trajectory.

The algorithm of the PI-method
In the PI-method (Turbell 2001), cone-beam projections are rebinned to semi-parallel projections, see figure 1(b). The mathematical relation between a cone-beam projection, p(β, γ, s), and the corresponding semi-parallel projection, p P (θ, t, s), is where the relations between the parameters (θ, t, s) and (β, γ, s) are In (9), R is the radius of the helix. The semi-parallel projection rays are parallel when viewed in the z-axis direction and fan-beam shaped when viewed in the tdirection, see figure 1(b). The PI-detector (associated with the Tam window (Tam et al 1998)) is located between two turns of the helix, see figure 1(b). Rebinned rays that fall outside the PI-detector are discarded. If the PI-detector is projected towards the (t, s)-plane, it forms a rectangular area. The s-axis has its origin at the centre of the PI-detector and is aligned with the z-axis. It has been shown that the PI-detector gives complete and non-redundant projection data, i.e. no redundant data are present during reconstruction (Turbell 2001). It has also been shown that the PIdetector restricts the illumination interval of a voxel (x, y, z) to exactly 180°. A basic outline for the simplest version of the PI-method is: (i) Obtain cone-beam projections from the CTscanner.
(iii) Discard projection data outside the PI-detector.
(iv) Pre-weight with k cos .
(v) Perform rampfiltering along rows in the t-direction of the virtual PI-detector.
(vi) Perform three-dimensional backprojection along the semi-parallel rays.
The PI-method is not an exact method and small artefacts can be seen in the reconstructed images. In general, the larger the cone-beam angle, the larger artefacts are observed (Sunnegårdh and Danielsson 2008).

The long object problem
The long object problem arises from the fact that only a limited range in the longitudinal (z-)direction of the patient is scanned. Since the projection angle interval required to reconstruct a voxel is at least 180 • , some voxels at the beginning and end of the z-range will lack projection data and become incompletely reconstructed, i.e. with errors.
In an iterative algorithm, projections are computed and compared with the measured projections. Some of the calculated projections are influenced by the incompletely reconstructed voxels. The technique by Magnusson et al (2006) handles the long object problem in single-energy CT volumes, but it has not been tested for dual-energy CT volumes yet.

DIRA-3D
DIRA-3D is an extension of the 2D algorithm DIRA presented in (Malusek et al 2017). The algorithm is illustrated in figure 2 and performs the following steps: 3. Perform automatic threshold segmentation on μ 1 and μ 2 .
4. Obtain base material reconstructed volumes by using the material decomposition methods (section 2.1).
6. Calculate monoenergetic forward projections P E i for parallel geometry, see (6).
Points 2-6 are repeated a predefined number of times. In mathematical terms, let denote the reconstructed images and the measured projections, respectively. Let the filtered backprojection operators be  PI and  ; II the former represents the PI-method working with semi-parallel projections, the latter represents ordinary parallel FBP. Let , which is the filtered backprojection result of the monoenergetic projections. The generation of monoenergetic projections followed by backprojection in DIRA serves as a regularization (Malusek et al 2017).

Methods
The new DIRA-3D algorithm was implemented according to the description in sections 2.4 and 3.2. It was tested in the geometry described in section 3.1.

Mathematical phantom and projection geometry
The mathematical phantom and projection geometry are shown in figure 3. The reconstructed area was a cylinder with diameter 353.3mm and height 132.5mm, which fitted into a box of size 353.3× 353.3×132.5mm 3 or 128×128×48 voxels. (Values in this section are rounded to 1 decimal digit.) The voxelsize was Δ x=Δ y=Δ z=2.8 mm. The phantom consisted of six ellipsoids with diameters 70.7mm, 70.7mm and 35.3mm. Two ellipsoids consisting of protein and water had their centers located at z=4.4 mm, which was within slice26. (The slices are numbered from 1 to 48). Four ellipsoids consisting of adipose tissue, lipid, muscle and compact bone had their centers located at z=−39.7 mm, which was within slice10. All the ellipsoids had their centers at the distance of 106.0mm from the rotation axis. The material composition of each ellipsoid is shown in table 1.
The x-ray source-to-rotation axis distance was 402.7mm. The number of projection angles was 800 (numbered from 0 to 799), the number of helix turns was 2, and the cone beam angle was ±κ max =± 3.14°. The pitch of the helix was P=88.3 mm or 32voxels. The corresponding spiral pitch factor defined as the ratio of the table feed per rotation to the total collimation width (NEMA PS3/ISO 12052) was 1 since the PI-detector is always located between two turns of the helix.
The detector size was N t ×N s =192×32 pixels. If the PI-detector is projected onto the (t, s)-plane, see figure 1(b), the detector pixel size becomes Δ t×Δ s=Δ x /1.5×Δ x /2. The two detectors (illustrated in figure 1(a)) had the same size and the phantom did not extend beyond their fields of view.
For parallel projection generation, the number of detector elements was 128 and the number of projection angles was 400 for the angular interval [0°,180°).
For the evaluation of DIRA-3D for the long object problem, the number of projections was reduced from 800 to 400, i.e. (0-399). Each cone-beam source made only one turn around the object volume, starting with projection number 0 and ending with projection number 399, see figure 3. As a consequence, the semiparallel projection at 400 was half filled only; the remaining data were missing. The procedure described in (Magnusson et al 2006) was used.

Implementation details
The rebinning step in the PI-method was omitted for simplicity. The 'measured' semi-parallel projections were simulated using the line integrals in (4) with , where m k are binary masks that indicate the position for the different ellipsoids. Energy spectra resembling Siemens SOMATOM definition Force 80 kV and Sn140kV spectra were used; their specific energies calculated from equation (7) were 50.0 and 88.5 keV, respectively. The algorithm is not very sensitive to the choice of these values. Quantum noise in the number of detected photons was not simulated.
At each iteration, reconstructed volumes μ 1 and μ 2 were threshold segmented into air, soft tissue and bone regions. The corresponding thresholds of 12 m −1 and 33 m −1 for the low-energy reconstructed LAC were chosen manually so that there were sufficient distances to the LACs of lipid (19.1 m −1 ) and protein (28.1 m −1 ), respectively.
Air was then decomposed into a (lipid, water) doublet, soft tissue was decomposed into a (lipid, protein, water) triplet and bone was decomposed into a (compact bone, bone marrow) doublet. As a consequence, the ellipsoid containing bone was decomposed into the (compact bone, bone marrow) doublet and the other ellipsoids were decomposed into the (lipid, protein, water) triplet.
The calculated semi-parallel polyenergetic projections were computed using (1), (2), (4) and the calculated parallel monoenergetic projections were computed using (6). The generation of parallel projections and following reconstruction was done slice by slice. The number of parallel detector elements was 128 and the number of projection angles in the interval [0°,179.55°] was 400.  Table 1. Elemental composition in mass fractions and corresponding density and tabulated LAC values μ(E 1 ) and μ(E 2 ) at E 1 =50 keV and E 2 =88.5 keV, respectively, for materials in the six ellipsoids constituting the phantom.

Ellipsoid
Material Mass fraction Density

Error calculation
The relative error of reconstructed μ 1 and μ 2 was estimated as d m m m m = - where μ t is the tabulated value (see table 1) and m is the average of the calculated LAC in a spherical region of interest (ROI). The ROI was defined as a sphere with a radius of one third of the evaluated ellipsoid radius (in the x, y plane) and positioned in the center of the evaluated ellipsoid.

Results
The ability of the iterative method to reconstruct accurate LAC values is demonstrated in figures 4, 5 and 6. Additional figures of reconstructed slices are shown in the supplementary material. Figure 4 shows how the reconstructed LACs for slice 10 (see figure 3) change with iteration number for the energy of E 1 =50 keV. The largest and the only clearly visible improvement was achieved for compact bone, where the LAC changed from approximately 59m −1 in iteration 1 to approximately 79.2m −1 in iteration 10, which is equal to the true value in table 1. Quantitative assessment of the improvement is presented in figure 5, which shows that relative errors (section 3.3) in LACs averaged over ROIs approached zero already after 5 iterations for the soft tissue ellipsoids and 10 iterations for the compact bone ellipsoid.
Beam hardening artefact suppression in the lipid and adipose tissue ellipsoids is illustrated in figure 6 for slice 10, iterations 1 and 10 and the photon energy of E 1 =50 keV. For adipose tissue, the LAC changed from from 19.3 to 20.2m −1 and for lipid the LAC changed from 18.3 to 19.1m −1 ; the resulting values agree with the true values in table 1. Figure 6 also  figure 3) shown for iterations 1, 2, 3 and 10, slice 10 and photon energy E 1 =50 keV.  shows the effect of the Hanning window included in the parallel FBP (and used in DIRA-3D by default). Compared to a situation without any window, the Hanning window reduced small overshoots and oscillations, but it slightly reduced the sharpness of the edges.
Suppression of the approximate-reconstruction artefacts with the increasing number of iterations is shown in figure 7 for slice14 (see figure 3). This position, close to the top of the ellipsoids, is known to be especially prone to these errors (Turbell 2001). The artefact is visible mainly in the air region where the LAC is close to 0. For this reason, the range of displayed LACs was restricted to the interval from −1m −1 to 1m −1 . Note that the shape of the artefacts differs between the low-energy and high-energy images due to the 90 • shift in the orientation of the x-ray sources. The convergence was very fast; visual differences between iterations 6 up to 25 (not shown) were small. Table 2 shows that there was a good agreement between mass fractions of lipid, protein and water calculated in the ellipsoids containing adipose tissue (region R 1 ) and muscle (region R 2 ) and the corresponding true values when the algorithm converged. In the ellipsoid containing compact bone, the calculated mass fractions of the compact bone and bone marrow were w 1 =0.43 and w 1 =0.57, respectively, for the first iteration and w 10 =1.0 and w 10 =0.0 for the tenth iteration. The relative error of the density decreased from 0.026 at iteration 1 to approximately 0 at iteration 10.
For the long object problem, figure 8 shows reconstructed LAC images similarly as in figure 7. The approximate-reconstruction artefact was decreased in a similar way. Areas approximated by red triangles have insufficient number of projections and will therefore never be perfectly reconstructed. The different positions of the red triangles in the top and bottom rows are caused by the 90 • offset of the second x-ray tube, see figure 1. Figure 9 shows that the relative errors for the long object experiment behave in a similar way as the relative errors in figure 5. Curves for the protein and water ellipsoids were not included; these two regions were poorly reconstructed owing to missing projection data.

Discussion
The presented algorithm represents a proof of concept. For this reason, a low resolution geometry where the computations can be done reasonably fast without requiring a specialized hardware was used. The phantom was designed to highlight artefacts caused by approximate reconstruction and the long object problem. A simple threshold segmentation was sufficient in this case. Some voxels were incorrectly segmented in the 0th iteration, but the accuracy of the segmentation increased with the number of iterations as the artefacts caused by approximate reconstruction and beam hardening artefacts disappeared, see also the discussion on the sensitivity to the selection of base materials below. For patient scanning, more advanced segmentation methods can be used, see for instance (Jeuthe 2017). In general, many voxels will be incorrectly segmented in the 0th iteration, but the accuracy of the segmentation increases with the number of iterations as the (mainly beam hardening) artefacts disappear. Alternatively, other segmentation methods may be used. We have experimented with hybrid segmentation methods MK2014 (Kardell et The first two terms in (12) and (13) corresponds to each other; the initial image is subtracted by an image reconstructed from forward poly-energetic Figure 8. The same as in figure 7 for the long object problem. The artefacts inside the red triangles are due to incomplete projection data. Figure 9. The same as in figure 5 for the long object problem. Curves for the water and protein ellipsoids are not plotted since they were affected by incomplete projection data.
projections of the ith image. The last term is different, however. For MDIR, the ith image is directly used in the update formula. In DIRA, an image reconstructed from forward mono-energetic projections of the ith image is used instead, i.e. a low-pass filtered version of the ith image. This procedure acts as a regularization that prohibit noise to grow by iteration. We have verified this in unpublished experiments with the original 2D DIRA. As mentioned in the introduction, MDIR was only implemented for simple 3D geometries. Especially, helical cone-beam geometry with approximate FBP, which are used in conventional medical CT-scanners, were not considered. The presented implementation of DIRA-3D focuses on the accuracy of LAC values, simplicity and speed. A better suppression of noise can be achieved by performing statistical optimization as in the ADMIRE algorithm (Ramirez-Giraldo et al 2018) or by replacing FBP with a statistical image reconstruction method inside the iterative loop of DIRA. The accuracy of LAC values close to object borders and the suppression of streak artefacts could be achieved by a more accurate modelling of the forward projections. For instance, subrays in the forward projections using the conebeam geometry followed with rebinning to semi-parallel projections could be used instead of the current simple approach with one ray per detector element.
DIRA-3D is conceptually similar to the original 2D version (Malusek et al 2017) except for the geometrical aspects. The main differences are: (i) the parallel FBP method was replaced with the approximate 3D FBPbased PI-method, (ii) the polyenergetic forward projections were modified to work in the helical 3D geometry, and (iii) the code for handling the long-object problem (Magnusson et al 2006) was added. It is reasonable to assume that both algorithms behave similarly in situations involving more realistic phantoms and statistical noise. Computer simulations performed with the 2D version showed that sufficiently good results were obtained after approximately 4-8 iterations while the noise did not increase with increasing number of iterations.
Of interest is whether the material decomposition to multiple doublets or triplets performed inside the iterative loop of DIRA brings any benefit over the material decomposition to just two bases in the iterative loop. Williamson et al (2006) showed that 2 material doublets (one for low-Z and the other for high-Z tissues) are sufficient to represent LAC curves of biological tissues with effective atomic numbers in the interval 2-20 (Calcium) and in the range 20-140keV. This was further confirmed for biological materials by Alvarez (2016). Nevertheless, Alvarez also found that a three or higher dimension basis set is needed if an externally administered high atomic number contrast agent is used (Alvarez 2016). Thus the benefit of the material decomposition inside the iterative loop may demonstrate itself in special cases only. The main application field of DIRA is radiation therapy with low energy photons. More work is needed to estimate the effects of base material selection on the accuracy of dose planning for instance when zinc is present in prostate calcifications Malusek et al (2018).
Similarly to some other algorithms like the one by Alvarez and Macovski (1976), DIRA-3D requires the knowledge of the x-ray spectra. In the present work, computed spectra were obtained from the manufacturer under a non-disclosure agreement. Alternatively, the spectra can be measured, see for instance (Matscheko and Alm Carlsson 1989, Duan et al 2011, Lin et al 2014.
The computation time was dominated by the calculation of forward projections and the reconstruction using the PI-method. When the resolution of the detector matches the resolution of the reconstructed volume, the computational complexity of the PI-method is  , where N is the number of voxels in the xor y-direction, N z is the number of voxels in the z-direction and N θ is the number of projection angles for one turn of the helix (Turbell 2001). Similarly, the complexities of monoenergetic and polyenergetic forward projections are  q ( ) N N N z 2 and  + q (( ) ) N K NN N z , respectively, where K is the number of energy bins of the x-ray spectrum. The relatively small increase in complexity of the polyenergetic projection compared to the monoenergetic one is caused by representation of the object materials with doublets and triplets.
The algorithm was implemented in MATLAB on a computer with 4 CPU cores, clock frequency of 3.8GHz and 16GB RAM. The computation time for one iteration was 145s. It was dominated by the polyenergetic forward projections (85 s) and the semi-parallel backprojections by the PI-method (45 s).

Conclusion
We have presented DIRA-3D, a model-based iterative reconstruction algorithm that estimates material composition of the imaged object from DECT projections obtained in helical geometries. Specifically, the algorithm determines mass fractions of components of user defined material doublets and triplets. The current proof of concept implementation was based on the PI-method, which was readily available to the authors. Nevertheless the WFBP should also work as both methods perform semi-parallel rebinning, 1D ramp-filtering and 3D backprojection; only the detector size and the weighting scheme differ.
The algorithm was evaluated using computer simulations with a simple phantom consisting of ellipsoids of different materials. In the studied geometry, the algorithm effectively removed artefacts caused by beam-hardening and approximate reconstruction and quickly converged. Also it was able to handle the long object problem. More work is needed to test the stability of the algorithm in the presence of quantum noise and scattered radiation.
The presented result indicate that the concepts used in DIRA can be extended to dual-sources, conebeams and helical scanning trajectories.