Reduction of beam hardening artifacts on real C-arm CT data using polychromatic statistical image reconstruction

ZEMEDI-10802; No. of Pages 11

ARTICLE IN PRESS

ORIGINAL PAPER

Reduction of beam hardening artifacts on real C-arm CT data using polychromatic statistical image reconstruction Richard N.K. Bismark 1,∗ , Robert Frysch 2 , Shiras Abdurahman 2 , Oliver Beuing 3 , Manuel Blessing 4 , Georg Rose 2 1 STIMULATE, Otto-von-Guericke-University Magdeburg, Institute for Medical Engineering, Building 09, Universitaetsplatz 2, 39106 Magdeburg, Germany 2 Otto-von-Guericke-University Magdeburg, Institute for Medical Engineering, Germany 3 University Hospital Magdeburg, Institute for Neuroradiology, Germany 4 University of Heidelberg, Department of Radiation Oncology, University Medical Center Mannheim, Germany

Received 6 March 2019; accepted 7 October 2019

Abstract Purpose: This work aims at the compensation of beam hardening artifacts by the means of an extended three-dimensional polychromatic statistical reconstruction to be applied for flat panel cone-beam CT. Methods: We implemented this reconstruction technique as being introduced by Elbakri et al. (2002) [1] for a multi-GPU system, assuming the underlying object consists of several well-defined materials. Furthermore, we assume one voxel can only contain an overlap of at most two materials, depending on its density value. Given the X-ray spectrum, the procedure enables to reconstruct the energy-dependent attenuation values of the volume. Results: We evaluated the method by using flat-panel cone-beam CT measurements of structures containing small metal objects and clinical head scan data. In comparison with the water-corrected filtered backprojection, as well as a maximum likelihood reconstruction with a consistency-based beam hardening correction, our method features clearly reduced beam hardening artifacts and a more accurate shape of metal objects. Conclusions: Our multi-GPU implementation of the polychromatic reconstruction, which does not require any image pre-segmentation, clearly outperforms the standard reconstructions of objects, with respect to beam hardening even in the presence of metal objects inside the volume. However, remaining artifacts, caused mainly by the limited dynamic range of the detector, may have to be addressed in future work. Keywords: Beam hardening, C-Arm CT, Cone beam, Truncation, Polychromatic, Polyenergetic, Iterative Reconstruction

1 Introduction X-ray computed tomography imaging is a widely used technique, especially as a diagnostic method in daily clinical routine. Among others, the most common reconstruction method is the filtered backprojection based Feldkamp reconstruction (FDK). The FDK image suffers from well-known

beam hardening artifacts due to the nonlinear X-ray absorption property that is not taken into account during the reconstruction. Such artifacts may obscure subtle findings, like early infarct signs or small intracerebral bleedings, which are relevant for treatment decisions. In practice, X-ray sources produce a broad spectrum of photons of different energies which coins the term

∗ Corresponding author at: Richard N.K. Bismark, STIMULATE, Otto-von-Guericke-University Magdeburg, Institute for Medical Engineering, Building 09, Universitaetsplatz 2, 39106 Magdeburg, Germany. E-mail: [email protected] (R.N.K. Bismark). URL: http://www.forschungscampus-stimulate.de/ (R.N.K. Bismark).

Z Med Phys xxx (2019) xxx–xxx https://doi.org/10.1016/j.zemedi.2019.10.002 www.elsevier.com/locate/zemedi

ARTICLE IN PRESS 2

R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

‘polyenergetic’/‘polychromatic’. The higher the energy of Xray photons, the lower the corresponding attenuation value of matter. This causes an increase in the mean energy of a polyenergetic X-ray beam while propagating through matter. Thus, the X-ray radiation gets “harder” while propagating through an object. As a consequence, the effective attenuation coefficient of the medium decreases with the penetration depth, which is the main reason for the cupping artifact [2]. Additionally, bone can “overspill” in soft tissue areas [3]. Beam hardening effects become even more important at the edges between different material types and can lead to streak and shadow artifacts in the reconstructed image [4]. Another disadvantage of the FDK is its need for certain scanning geometries, such as circle trajectories. Iterative statistical methods can handle arbitrary scan trajectories and include prior knowledge, like object constraints and variance properties. Common methods to reduce beam hardening artifacts are based either on the pre-processing of the measured data, like a water correction [3,5] or a post-image processing of the reconstructed image [2,6,7]. Other approaches are based on ways to get different energy-dependent measurements of the same object. The most prominent approach, namely dual energy CT [8,9], requires two measurements with two different X-ray spectra. This enables the splitting of the attenuation coefficients into the main contributions, the Compton effect and the photoelectric effect, which is advantageous for quantitative CT and tissue characterization. The main drawback is the need for two X-ray spectra, thus the need for specialized hardware. Iterative reconstruction techniques were developed to solve the beam hardening issues and deliver promising results [10,11,33]. However, those methods rely on a proper segmentation and can suffer severely from wrong assumptions on their prior image and object constraints. A recent, very promising development of energy resolved photon-counting detectors has the potential to overcome the beam hardening problem [12]. However, the performance of these detectors, i.e. their energy resolution and counting rate, is still very poor and so this technique is not yet viable in real applications. Wu et al. developed a robust polychromatic statistical method, which can be of practical use for clinical applications. They state that they only need to approximate the polychromatic X-ray spectrum by three energy bins, which can be of importance for the computational speed of an algorithm [13]. Kachelrieß et al. developed an entirely empirical postprocessing cupping correction, which does not need any information about the X-ray spectrum or size and position of the reconstructed object [14]. Recent findings in the field of cone-beam consistency measures indicate that certain conditions in CT projection data can be used to compensate artifacts caused by motion as well as beam hardening [15,17,32].

In this paper, we present a multi-GPU implementation of a polyenergetic statistical reconstruction (PSR) approach for flat-panel cone-beam CT (FP-CBCT) with standard energyintegrating detectors. Such CT devices are commonly used in hospitals around the globe. Our implementation is based on the work of Elbakri and Fessler [1,18]. It utilizes a physical model of the forward projection to deal with beam hardening. Furthermore, the statistical character of the radiation is incorporated for noise reduction. It combines the properties of statistical weightings with the polyenergetic characteristic of the measured intensities. It is a versatile method and can be applied to CBCT. For the first time, we adapt this algorithm to real FP-CBCT data obtaining promising results for different objects, including phantoms as well as clinical head scan data.

2 Methods 2.1 Theoretical background The implemented reconstruction algorithm (PSR) is based on a penalized Poisson likelihood approach [1]. A cost function which is composed of a data consistency and a regularizer term is minimized by using quadratic surrogates with the optimization transfer principle [19–21]. The data consistency incorporates a polyenergetic X-ray projection model and assumes the measurements to be independently-distributed Poisson random variables. The attenuation μ of a voxel is represented by the product of its density ρ and the energy dependence m (E) of the corresponding material. The contributing materials have to be known to make an assignment from μ to ρ feasible. To speed up the monotonic minimization of the cost function, Elbakri and Fessler proposed to use ordered subsets and a “Precomputed Curvature”. Both proposals corrupt the monotonicity but drastically improve the convergence of the algorithm. Eventually, they obtain an iterative routine with an update step that includes a forward projection, which is then compared to the measurements. A backward projection of that comparison is then then performed to update the image. This step can incorporate a scatter estimation and a data regularization term. With [. . .]+ denoting the non-negativity operation, the update step of the density of each voxel j is written as ⎞⎤ ˆj + β ∂S M N ∂ρj ⎟⎥ ⎜ ⎢ ρj → ⎣ρj − Crelax ⎝ ⎠⎦ . ∂2 S dj + β ∂ρ2 ⎡

⎛

j

(1)

+

2

∂S Here, β ∂ρ and β ∂∂ρS2 are surrogate functions of a Huber j j

penalty regularization (smoothness constraints) regarding only each ρj voxel’s nearest neighbors, and β as a control parameter to adjust the penalty to the data consistency. M is a scalar which is correlated to the number of the used subset ξ of projections, in order to update the density ρj in a voxel j. We

ARTICLE IN PRESS R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

3

added a factor Crelax that operates as a relaxation parameter, in case we want to slow down or fasten up convergence speed. The default value is Crelax = 1. The “Precomputed Curvature” is defined as dj :=m2water (Eeff ) ·



(2)

aij γi Yi .

i

Eq. (2) basically represents an iteration-independent backprojection of the measured intensity Yi which can be calculated step (see “Precomputed Curvature” once as a pre-processing  in [1]). γ i = j aij denotes the path length through the object and mk (E) is the known energy dependence [22] of the kth material. Note that m2water (Eeff ) is a modification of the origi nal k m2k (Eeff ). This is viable due to the fact that our objects consist mostly of water. The effective energy of the X-ray beam with a measured energy spectrum I(E) is given by the centroid  Eeff =

E · I (E) dE  . I (E) dE

(3)

2.2 The approach and implementation

Nˆj contains a polychromatic forward projection Yi (Eq. (4)) of the current volume ρj which is compared with the measured intensities Yi and backprojected afterwards. The projection model for the ith ray is the energy-dependent Lambert–Beer law: 

 Yi =

dE Ii (E) exp −



 mk (E) sik (ρ) .

 Thereby, sik (ρ):= j aijk ρj is a classical forward projection of voxels consisting of material k. Ii (E) is a blank scan, containing information of the initial X-ray spectrum and the detector response. aijk = aij fjk are the geometry factors combined with the information about which material the voxel consists of. Elbakri et al. [1] used the following choice: 

1 0

voxel j belongs to the kth material else

(5)

Nˆj can be written as Nˆj =

 i∈ξ k

  Yi ∂Yi aijk 1 − . Yi ∂sik

We assume that the given set of K materials can be distinguished by their density ρk as ρ1 < ρ2 < · · · < ρK . To connect K materials as proposed, we need to define K − 1 cubic splines. With 0 < τ  ≤ 0.5 we define a smoothness parameter   τ k :=τ  ρk+1 − ρk

(7)

(4)

k

fjk =

Figure 1. The used parameter set (ρw = 1.0 g/cm3 , ρb = 1.92 g/cm3 , τ  = 0.25) for the reconstruction of the clinical head data Figure 7. The plot shows the resulting material fractions of water f w (red line) and bone fb (blue line) as a cubic polynom of the density ρ.

(6)

This type of polyenergetic reconstruction method needs a pre-segmentation of the reconstructed volume (see Eq. (5)). In a follow up, Elbakri et al. expanded the method with a “Displacement model” [18], which allows material mixtures in voxels via cubic splines.

and additionally the central point  1 k Rk := ρ + ρk+1 . 2 k :=Rk − τ k This enables us to define the threshold values ρmin k k k and ρmax :=R + τ . Finally we can write down the mixture of two materials as ⎧ k k , 1 ρmin < ρ < ρmax ⎪ ⎪ ⎨   k+1 k (8) f k (ρ) = f ρ; Rk , τ k ρmax < ρ < ρmin , ⎪ ⎪ ⎩ k+1 0 ρ > ρmin .

The fraction of the other material follows as f k+1 (ρ) = 1 − f k (ρ) in the corresponding   intervals. The function f ρ; Rk , τ k is used to smoothly connect the k+1 k edges of the fraction of the interval between ρmax and ρmin as shown in Figure 1 

k

f ρ; R , τ

k



3   ρ − Rk 3 ρ − Rk 1 :=  3 − + . k 4τ 2 4 τk 

(9)

ARTICLE IN PRESS 4

R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

We implemented a polychromatic ray-based forward projector and a voxel-driven footprint backprojector (TT-A1) [23,24] on multi-GPU hardware with OpenCL programming language. We do not observe any relevant mismatch between those forward and backprojection operators. The voxel-driven backprojector has been chosen for its advantage in computation performance compared to ray-based alternatives. This becomes even more important considering the options of parallelization. The spectrum is divided into 40 equally-spaced energy bins by linear interpolation of given sampling points with user-defined minimum energy value and bin width. Thus, the spectrum is represented by 40 values and two parameters. Moreover, we assume the object to consist of water, bone and a metal Depending on its density, the fraction  compound.  fjk = f k ρj of the materials in the jth voxel is approximated continuously by cubic splines (Eq. (8)) for material compositions (“Displacement Model” [18]). Thus, the material type is only identified by its density. Hence, it is permitted to change its type during the iterative process. The ray-based forward projection calculates the intersection points of a given ray with voxel borders. To distribute the computational load across multiple GPU devices, projections are divided into separate t-blocks (i.e. lines parallel to axial slices during a standard C-arm scan protocol) of the detector plane. For the backprojecting step, the volume is divided into y-blocks of the world coordinate system and split between the GPU devices [25]. We typically initialize the volume with zeros or with an FDK. The volume updates are iterated and start with a subset size of one angle per subset. After a user-defined number of iterations over all subsets, the angles per subset get doubled, thus halving the number of subsets. We repeat this procedure until one subset contains all angles [24]. The iteration with the biggest subset is repeated until the update step becomes small or a user-defined maximum number of iterations is achieved. To compare the results obtained by the PSR with a monochromatic approach, we realize a monochromatic statistical reconstruction (MSR) method by using the PSR algorithm, but treating all voxels as water with respect to their energy-dependent absorption. Additionally, the MSR routine is configured with a monoenergetic X-ray spectrum X(E) = δ (E − E0 ). In this setting, we use water corrected projections, which are also used for the FDK reconstruction. Thus, the MSR algorithm represents a statistical but monoenergetic iterative reconstruction technique that uses the identical forward and backward projection operators as the PSR. The consistency of the physical model of the implementation is illustrated, since the MSR results show similar beam hardening artifacts as the FDK. In contrast to Elbakri et al., Yi and its gradient are calculated on-the-fly from the forward projection sik , instead of

Figure 2. The X-ray spectrum [26] and the detector response function [29] normalized by its maximum value which are used for the PSR algorithm.

interpolating pre-calculated values. This means, the gradient is computed explicitly through derivation of (4): ∂Yi =− ∂sik







dE Ii (E) mk (E) exp −

 mk (E) sik (ρ)

.(10)

k

We estimated the X-ray spectrum X (E) by using the tool spekCalc [26], taking into account the 12◦ opening angle and 0.9 mm-Cu, 2.5 mm-Al pre-filtering of the incident X-rays (Figure 2). Furthermore, we used pre-processed projection data, which was normalized (provided as extinctions i = ln I0 /Yi ) and processed only with a kernel-based scatter correction routine [27]. Note that this means, in particular, no water correction has been applied to the data before reconstruction with the PSR. As mentioned above, standard water correction has been applied prior to reconstruction with MSR and FDK to provide reasonable comparison. The C-arm system is using an 1D anti-scatter grid and there was no bowtie filter used. We considered the possibility to leave out the provided kernel-based scatter correction, but can state that the resulting reconstructions were not enhanced by any means. We incorporated a detector response function D (E) into the forward projection (Eq. (4)) to model the measurement characteristics more accurately I (E) = X (E) · D (E) .

(11)

The used C-arm system features a CsI scintillator-based detector. Responses of NaI(Tl), CsI(Tl) and CsI(Na) detectors above 20 keV behave qualitatively similar [28]. Since the used energy bins lie in the interval from 20 to 110 keV, we used the response functions of Roberts et al. [29]. With the assumption equation (11) we are able to calculatethe intensities from the measured extinctions Yi = exp {− i } I (E) dE.

ARTICLE IN PRESS R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

5

Figure 3. Comparison of FDK (a), MSR (b) and PSR (c) reconstructions of the steel needle dataset.

2.3 Materials The first experiments were carried out by using a head of a pig cadaver, that contained a needle (steel) of 2.5 mm in diameter. The head was scanned on a C-arm system (Artis Zeego angiography system, Siemens Healthcare GmbH, Erlangen, Germany) using 616 × 480 detector resolution with a pixel size of 0.616 mm × 0.616 mm. The tube voltage was set to 102 kV. The angular resolution was 248 covering 200◦ of a circular short scan. Another measurement was done with a 5 euro cent coin, which mainly consists of steel (∼94%). The coin was placed inside a water basin. Additionally, we have performed scan simulations of a head phantom [30] including properties that cause artifacts, i.e. Xray scattering as well as limited dynamics of the detector. We motivate those simulations to study the remaining artifacts in the case of the 5 euro cent coin (see Figure 4). To simulate the saturation effects of the detector, we used a mapping of the form:  = a tanh a , with a = 4. To approximate the qualitative behavior of scattered radiation, we applied a Gaussian blur filter to the projections of the needle voxels. In a final evaluation step, we used real clinical data. Data were acquired from 7 patients by a C-arm system (Artis Q biplane angiography system, Siemens Healthcare GmbH, Erlangen, Germany) using a tube voltage of 108 kV, a 1240 × 960 pixel detector (resolution: 0.308 mm × 0.308 mm) and an angular resolution of 496 views for a 200◦ segment of a circular trajectory (short scan). We initialized the PSR with the uncorrected ML volume and performed only one relaxed (Crelax = 0.05 in Eq. (1)) iteration of the PSR over all views with single angle subsets. We found the PSR to be more sensitive to the truncated patient table. Usually, the truncation of the table does not corrupt the convergence of the other reconstruction techniques used in this work. Since the table is out of range of the reconstruction volume, the reconstruction methods tend to increase the attenuation at the bottom edge of the volume. The table is completely absent in around 15% of the projections. Due to a different mathematical structure, the ML is more robust in this setting [31] compared to the PSR, which can suffer from

truncation showing severe artifacts. Since the ML reconstruction does suffer a lot less from such truncation, we achieved proper beam hardening reduced reconstructions that are comparable to our PSR ‘early-stopping’ approach. We investigated this behavior and evaluated the methodology in [37]. We compared the results with iterative statistical but mono energetic maximum likelihood reconstructions (ML) [24] that were given the same data or bi-material beam hardening corrected projections (bmc) [32] of the same dataset. Note that ML reconstructions are used in the bmc-case instead of MSR, since the ML is more resistant to patient table truncation, therefore, providing a more reasonable comparison.

3 Results We reconstruct the head of the pig cadaver on a 256 × 256 × 199 grid with an isotropic voxel size of 1 × 1 ×1 mm3 . The density threshold values are set to ρw = 1.0 g/cm3 , ρb = 1.7 g/cm3 and ρm = 3.0 g/cm3 . The regularization parameter βpig is set to 10−1 (see Eq. (1)). We initialize the PSR with an iterative statistical (but monoenergetic) maximum likelihood reconstruction (ML [24]) and perform 20 iterations of ordered subsets. We compare the resulting volumes of different reconstruction methods at one specific axial slice containing the needle. The standard FDK image in Figure 3 has stronger shadow artifacts, caused by beam hardening, than the iterative MSR. Both display a similar representation of the needle. The PSR clearly outperforms the water corrected FDK and also the water corrected MSR, revealing low shadow and blurring artifacts as well as a sharper needle shape. The estimations of the needle diameter based on the PSR reconstruction reproduce the real numbers in good agreement (not shown here). The comparison of different detector response functions, i.e. constant over the energy, arbitrary linear function of the energy, and as defined by Roberts et al. [29] reveal similar image qualities (not shown here). This indicates a low sensitivity of the image quality to the detector response function, which is in accordance to what as Elbakri and Fessler observed for accuracy of the X-ray spectrum [1].

ARTICLE IN PRESS 6

R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

Figure 4. Comparison of the MSR (a) and PSR (b) reconstruction of the coin dataset (left). The blue line within image (a) indicates the data-range used for the contrast profile plots (right).

The density threshold values of the 5 euro cent coin reconstruction are set to ρw = 1.0 g/cm3 and ρm = 3.0 g/cm3 . The regularization parameter βcoin is set to zero to emphazize the effectiveness of the polychromatic nature of the measurements and the resulting reconstructions. We initialize the PSR with the FDK and perform 20 iterations of ordered subsets. We repeat this procedure with the MSR in the same manner. The following figures visualize slices of reconstructed volumes. The density values are given in g/cm3 or density Hounsfield units (HU = 1000 · (ρ − ρw ) /ρw ) respectively. Since the used detector is intended to measure extinction values that are common in head scans, the dynamic range of the detector is not suitable for high attenuation of metal objects like the coin. The thickness-to-diameter ratio of a 5 mm euro cent coin is 21.25 1.67 mm ≈ 12.7. The extinctions found in a frontal and a side view projection should have a similar ratio (as the path length through the coin increases correspondingly). Figure 5 shows the measured data and reveals that the actual ratio of the corresponding measured extinc4.8 tions is approximately 4.2 ≈ 1.14 and by that, far away from the expected value. This property, among others (e.g. approximated I0 , insufficient compensation of scatter), prevents the estimation of quantitatively correct Hounsfield units. However, even under this circumstances, the PSR remains a robust reconstruction routine. As shown in Figure 4c, the remaining blurring yields an upper bound for the resulting thickness of the coin by 5 voxels (=5 ˆ mm) and a lower bound by 2 voxels (=2 ˆ mm). The estimated thickness and the diameter of the PSR coin representation are in good agreement with the true values of 1.67 and 21.25 mm. The results of the simulations (Figure 6) reveal that scattering causes smaller attenuation values around the metal needle. On the contrary, the limited detector dynamics lead to brighter regions around the needle as well as shadow artifacts at both of its ends. This indicates that the appearance of the resulting artifacts depends strongly on the shape of the metal object. The results of the clinical head scan data are presented in Figures 7 and 8. In the first case we also show the resulting FDK that was given the in-build water correction of the C-arm and performed a ring artifact correction as a post-processing

step similar to the standard method described in [27]. Such ring artifacts are common in C-arm imaging. They are caused by defect detector pixels or insufficient homogeneity of detector dynamics. The bi-material beam hardening corrected ML still suffers from remaining beam hardening artifacts caused by the dense bone of the skull, whereas our PSR method significantly compensates for these artifacts, i.e. the “spill-over” effects of the bones are drastically reduced. The reconstruction volumes have a size of 464 × 464 × 176 voxels with 0.5 × 0.5 × 1.0 mm3 resolution. The density threshold values were set to ρw = 1.0 g/cm3 and ρb = 1.92 g/cm3 . The regularization parameter βpatient was set to 1.0. We observed that βpatient > βpig works out better (probably due to the different spatial resolution), while testing different values by trial-and-error. For all obtained results, the smoothness parameter τ  (Eq. (7)) was set to 0.5. 3.1 Clinical relevance As beam hardening causes the bone to spill over, physicians cannot access important information in those regions. Figure 9 shows the impact of this behavior on the clinical data and the benefit of this preliminary promising result. The green arrow marks a small bleeding. The white arrows point out sulci of the brain reconstruction. The brain in the marked area is swollen, which is indicated by the missing sulci (blue arrows). Such an absence is an important indicator for the radiologist. Thus, the beam hardening induced cupping on the left decreases the diagnostic value of the reconstruction and the PSR on the right is increasing the relevant information for the radiologist.

4 Discussion We have implemented and evaluated an extended polychromatic statistical image reconstruction (PSR) for flat panel cone-beam CT. For the first time, we used the PSR on clinical data measured with a common C-arm device. The results are promising, as we have shown that the PSR can drastically reduce typical beam hardening induced cupping on several

ARTICLE IN PRESS R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

7

Figure 5. The measured frontal and side view projections of the coin show an inconsistency in the extinction values from those two directions. The relation of the extinction value is expected to be proportional to its geometric size (except for beam hardening effects), which is above 4.8 ≈ 1.14. a factor of the order 10. However, the measured ratio lies around 4.2

Figure 6. Density maps [g/cm3 ] showing one slice of the head phantom (a) containing a metal object, (b) PSR result of projections including a scatter simulation and (c) PSR result of projections including a limited detector dynamics range.

Figure 7. Comparison of the reconstructions using the ML method with bi-material beam hardening corrected (bmc) projections of the head (mid row), PSR (bottom row) and water corrected (wc) FDK of the same measurement of a patient with intracranial hemorrhage (visible in the fourth row).

ARTICLE IN PRESS 8

R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

Figure 8. Comparison of uncorrected ML, bi-material beam hardening corrected ML and PSR on different clinical head scan data sets. The blue lines in the first column indicate the data-range used for the contrast profile plots in the right column. They emphasize the compensation of the beam hardening induced cupping. Note that there was no ring artifact correction used.

ARTICLE IN PRESS R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

9

Figure 9. An enlarged view of the last slices shown in Figure 7. With PSR (right) and bi-material corrected ML (left), circumscribed mild brain edema is depicted, but barely visible. The blue arrows show the absence of sulci, which should be visible similar to the white arrow indications.

head scans (Figures 7 and 8) as well metal artifacts (Figures 3 and 4). A comparison between FDK (with simple water correction) and PSR shows that the dark regions caused by beam hardening and metal artifacts are strongly reduced by the PSR and the structure of the surrounding volume became clearer. However, in the presence of such metal objects, shadow-like corruptions of the image are still present, probably caused by photon starvation due to the high photon attenuation (Figure 4). The areas perpendicular to those dark regions are brighter which corresponds to higher density values. This qualitative behavior can be observed in the simulations as well (Figure 6c) and is caused by the limited detector dynamics. As the main result, we have shown that PSR clearly outperforms the other methods, which do not take into account the energy-dependent absorption of matter directly. The cupping as well as streak artifacts caused by the beam hardening of metal objects are significantly reduced. This can be observed in the contrast profile plots shown in Figure 8. Furthermore, the edges of the structures are less blurred, thus the shapes of the metal objects are largely preserved. Looking on the performance of the PSR with respect to the detector response function, we observed a low sensitivity on this property. We have shown in Figure 5 that the limited detector dynamic of C-arm devices are a limiting factor in metal artifact reduction. Similar behavior is assumed for CT devices, where the detectors have higher bit depth but still suffer from photon starvation. If we initialize the PSR routine with a monochromatic spectrum and treat all voxels as water (MSR), we observe beam hardening artifacts similar to the FDK (Figure 3). This observation supports the underlying physical model of the PSR. The comparison of MSR/ML and PSR is more meaningful since those reconstruction methods are statistical and iterative algorithms.

We defined one density value per material as threshold for the displacement model (Eq. (7) to (9)). Those thresholds need to be sufficiently different from each other to produce stable results. Thus, this method might not perform as well in the case of contrast agent (CA) in blood vessels, since there is a very broad density range of blood-CA mixtures, i.e. their density ranges may overlap. A similar problem would occur for a complex multi-metal object containing metal compounds of similar densities. Such conditions could be overcome by a proper segmentation. For example, in perfusion imaging, one can find the voxels that contain CA by subtracting the mask from the volume. For this case, the “Solution model” could be applied, as proposed in the original source [18] for the case of a mineral solution. Available prior knowledge, for example, from CAD models of metallic implants could lead to a nearly perfect segmentation. Such a segmentation could also open up the opportunity to approximate a good scatter estimation of the metallic device or even fill spots of the sinogram that are corrupted by photon starvation. In principle, the PSR is capable of including scatter estimations as additive term to the forward projections. Further, advanced pre-processing steps of the projections should be a compatible modification as well (e.g. primary X-ray modulation that was already shown to work on C-arm devices [36]). This will be subject of future work. Xu et al. [33] proposed a method that does not depend on prior knowledge about the highly absorbing materials and instead just assumes unknown polynomial absorption properties. They overcome the segmenting problem by using 2D–3D registration of the implant scan in air or by a given CAD model. This can limit the feasibility in some cases in daily practice (e.g. unknown implants). On the other hand, their method is supposed to work in VOI imaging, which can be of high value in other cases of the clinical routine. A similar idea of an

ARTICLE IN PRESS 10

R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

approximated outer volume for the PSR needs to be implemented and evaluated first to compare the performance in a VOI imaging setting. Since the PSR is already sensitive to the truncation of the patient table outside of the reconstruction volume, we initialized the reconstruction with a non-beam hardening corrected ML reconstruction that does suffer less from such conditions [31]. In the case of the clinical head scan datasets we performed only one relaxed iteration of PSR that includes all views once. Alternatively, one could increase the size of the reconstruction volume to get a better approximation of the patient table. Then, we can initialize the volume with zeroes and perform full ordered subsets iterations of the PSR. We have shown that such an approach can produce similar results in [37], but needs further algorithmic development to achieve the same spatial resolution. The independence of a good initialization is the motivation of such an approach and will be subject of future work. Humphries et al. [34,35] developed a method similar to the PSR. They introduced the polyenergetic forward model into an iterative weighted least squares reconstruction that also features TV-superiorization. Their method shows promising results on simulated 2D data even when confronted with a limited angle setting. A direct comparison of the PSR with this method is of interest and will be subject of future work as well. To improve the image quality further, more information about the measuring process could be taken into account. For instance, in most FP-CBCT systems the X-ray intensity of the tube is adapted during the scanning procedure, in order to save dose as well as to optimally utilize the dynamic range of the detector. Hence, the approximation that I(E) in (Eq. (11)) is constant with respect to the recording angle does not hold true, which corrupts statistical assumptions as well as the polyenergetic specifications and might be considered in the future. Since the PSR produces good approximations of edges of highly absorbing materials, one could consider the PSR as a practical tool for segmentation, which is a bottleneck of a lot of methods in computed tomography.

5 Conclusion We present a cone-beam flat panel CT implementation of the PSR routine of Elbakri and Fessler [18]. For the first time, this routine is used on clinical cone-beam CT data. The results on measured data for an angiographic C-arm device show a clear image quality improvement. Beam hardening artifacts were reduced while the surrounding tissue was less blurred, allowing a better recognition of a hemorrhagic stroke lesions. When metal objects are present, artifacts remain due to the limited dynamic range of the detector as well as due to scattered radiation, which were not yet considered in this work. For the latter, given a good scatter approximation, the PSR is capable of taking this information into account.

We show that the PSR significantly enhances the image quality, in both low and high contrast scenarios, making the clinical diagnosis more reliable.

Compliance with ethical standards The authors have no relevant conflicts of interest to disclose. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. For this type of study formal consent is not required.

Acknowledgments The work of this paper is partly funded by the German Ministry of Education and Research (BMBF) within the Forschungscampus STIMULATE under grant number ‘ 13GW0095A’. We want to thank Thomas Hoffmann, Sebastian Gugel and Axel Boese for the smooth execution of the measurement tasks. Furthermore we thank Steffen Serowy and Martin Skalej for providing the clinical data sets and Tim Pfeiffer for some helpful discussions.

Appendix A Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.zemedi.2019.10.002.

References [1] Elbakri I, Fessler J. Statistical image reconstruction for polyenergetic Xray computed tomography. IEEE Trans Med Imaging 2002;21(2):89–99, http://dx.doi.org/10.1109/42.993128. [2] Brooks RA, Chiro GD. Beam hardening in X-ray recontomography. Phys Med Biol 1976;21(3):390 structive http://stacks.iop.org/0031-9155/21/i=3/a=004. [3] Joseph PM, Spital RD. A method for correcting bone induced artifacts in computed tomography scanners. J Comput Assist Tomogr 1978;2:100–8, http://dx.doi.org/10.1097/00004728-197801000-00017. [4] Slaney ACKM. Principles of computerized tomographic imaging. New York: IEEE Press; 1988. [5] Nalcioglu O, Lou RY. Post-reconstruction method for beam hardening in computerised tomography. Phys Med Biol 1979;24(2):330 http://stacks.iop.org/0031-9155/24/i=2/a=009. [6] Hsieh J, Molthen RC, Dawson CA, Johnson RH. An iterative approach to the beam hardening correction in cone beam CT. Med Phys 2000;27(1):23–9, http://dx.doi.org/10.1118/1.598853. [7] Meagher JM, Mote CD, Skinner HB. CT image correction for beam hardening using simulated projection data. IEEE Trans Nucl Sci 1990;37(4):1520–4, http://dx.doi.org/10.1109/23.55865. [8] Alvarez RE, Macovski A. Energy-selective reconstructions in Xray computerised tomography. Phys Med Biol 1976;21(5):733 http://stacks.iop.org/0031-9155/21/i=5/a=002.

ARTICLE IN PRESS R.N.K. Bismark et al. / Z Med Phys xxx (2019) xxx–xxx

[9] Fessler JA, Elbakri IA, Sukovic P, Clinthorne NH. Maximum-likelihood dual-energy tomographic image reconstruction, vol. 4684; 2002. p. 38–49, http://dx.doi.org/10.1117/12.467189. [10] Yan CH, Whalen RT, Beaupre GS, Yen SY, Napel S. Reconstruction algorithm for polychromatic CT imaging: application to beam hardening correction. IEEE Trans Med Imaging 2000;19(1):1–11, http://dx.doi.org/10.1109/42.832955. [11] Lemmens C, Faul D, Nuyts J. Suppression of metal artifacts in CT using a reconstruction procedure that combines MAP and projection completion. IEEE Trans Med Imaging 2009;28(2):250–60, http://dx.doi.org/10.1109/TMI.2008.929103. [12] Gonzales B, Lalush D. Full-spectrum CT reconstruction using a weighted least squares algorithm with an energyaxis penalty. IEEE Trans Med Imaging 2011;30(2):173–83, http://dx.doi.org/10.1109/TMI.2010.2048120. [13] Wu M, Yang Q, Maier A, Fahrig R. A practical statistical polychromatic image reconstruction for computed tomography using spectrum binning. In: SPIE, editor. Proc SPIE medical imaging 2014. 2014. p. 9033–126 http://www5.informatik.uni-erlangen.de/Forschung/Publikationen/2014/ Wu14-APS.pdf. [14] Kachelrieß M, Sourbelle K, Kalender WA. Empirical cupping correction: a first-order raw data precorrection for conebeam computed tomography. Med Phys 2006;33(5):1269–74, http://dx.doi.org/10.1118/1.2188076. [15] Aichert A, Berger M, Wang J, Maass N, Doerfler A, Hornegger J, et al. Epipolar consistency in transmission imaging. IEEE Trans Med Imaging 2015;34(11):2205–19, http://dx.doi.org/10.1109/TMI.2015.2426417. [17] Frysch R, Rose G. Rigid motion compensation in interventional C-arm CT using consistency measure on projection data. In: Proceedings of the 18th international conference on medical image computing and computer-assisted intervention – MICCAI 2015 – vol. 9349. New York, NY, USA: Springer-Verlag New York, Inc.; 2015. p. 298–306, http://dx.doi.org/10.1007/978-3-319-24553-9 37. [18] Elbakri IA, Fessler JA. Segmentation-free statistical image reconstruction for polyenergetic X-ray computed tomography with experimental validation. Phys Med Biol 2003;48(15):2453 http://stacks.iop.org/0031-9155/48/i=15/a=314. [19] Erdogan H, Fessler JA. Monotonic algorithms for transmission tomography. IEEE Trans Med Imaging 1999;18(9):801–14, http://dx.doi.org/10.1109/42.802758. [20] Pierro ARD. On the relation between the ISRA and the EM algorithm for positron emission tomography. IEEE Trans Med Imaging 1993;12(2):328–33, http://dx.doi.org/10.1109/42.232263. [21] Pierro ARD. A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography. IEEE Trans Med Imaging 1995;14(1):132–7, http://dx.doi.org/10.1109/42.370409. [22] Hubbell JH, Seltzer SM. Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients from 1 kev to 20 mev for elements z = 1 to 92 and 48 additional substances of dosimetric interest (1989, 1990, 1996). [23] Long Y, Fessler JA, Balter JM. 3D forward and back-projection for X-ray CT using separable footprints. IEEE Trans Med Imaging 2010;29(11):1839–50. [24] Frysch R, Pfeiffer T, Bannasch S, Serowy S, Gugel S, Skalej M, et al. C-arm perfusion imaging with a fast penalized maximum-likelihood approach, vol. 9033; 2014, http://dx.doi.org/10.1117/12.2043450. p. 90332M-8.

11

[25] Pfeiffer T, Frysch R, Gugel S, Rose G. ML reconstruction of conebeam projections acquired by a flat-panel rotational X-ray device. In: Proceedings of SPIE medical imaging conference, vol. 8668. 2013. p. 86682V, http://dx.doi.org/10.1117/12.2043450. [26] Poludniowski G, Landry G, DeBlois F, Evans PM, Verhaegen F. Spekcalc: a program to calculate photon spectra from tungsten anode X-ray tubes. Phys Med Biol 2009;54(19):N433 http://stacks.iop.org/0031-9155/54/i=19/a=N01. [27] Zellerhoff M, Scholz B, Ruehrnschopf E-P, Brunner T. Low contrast 3D reconstruction from C-arm data; 2005, http://dx.doi.org/10.1117/12.593433. [28] Aitken D, Beron B, Yenicay G, Zulliger H. The fluorescent response of NaI(Tl), CsI(Tl), CsI(Na) and CaF2(Eu) to X-rays and low energy gamma rays. IEEE Trans Nucl Sci 1967;14(1):468–77, http://dx.doi.org/10.1109/TNS.1967.4324457. [29] Roberts D, Hansen V, Niven A, Thompson M, Seco J, Evans P. A low Z linac and flat panel imager: comparison with the conventional imaging approach. Phys Med Biol 2008;53(22):6305–19 http://epubs.surrey.ac.uk/738905/. [30] Aichert A, Manhart MT, Navalpakkam BK, Grimm R, Hutter J, Maier A, et al. A realistic digital phantom for perfusion C-arm CT based on MRI data. In: 2013 IEEE nuclear science symposium and medical imaging conference (2013 NSS/MIC). 2013. p. 1–2, http://dx.doi.org/10.1109/NSSMIC.2013.6829168. [31] Frysch R, Bismark R, Maier A, Rose G. Ray-density weighted algebraic reconstruction for volume-of-interest CT. In: The 14th international meeting on fully three-dimensional image reconstruction in radiology and nuclear medicine. 2017. [32] Abdurahman S, Frysch R, Bismark R, Melnik S, Beuing O, Rose G. Beam hardening correction using cone beam consistency conditions. IEEE Trans Med Imaging 2018, http://dx.doi.org/10.1109/TMI.2018.2840343. [33] Xu S, Uneri A, Jay Khanna A, Siewerdsen JH, Stayman JW. Polyenergetic known-component CT reconstruction with unknown material compositions and unknown X-ray spectra. Phys Med Biol 2017;62(8):3352–74, http://dx.doi.org/10.1088/1361-6560/aa6285/pdf. [34] Humphries T, Winn J, Faridani A. Superiorized algorithm for reconstruction of CT images from sparse-view and limited-angle polyenergetic data. Phys Med Biol 2017;62(16):6762–83, doi:10.1088%2F13616560%2Faa7c2d. [35] Humphries T, Gibali A. Superiorized polyenergetic reconstruction algorithm for reduction of metal artifacts in CT images. In: 2017 IEEE nuclear science symposium and medical imaging conference (NSS/MIC). 2017. p. 1–6. [36] Bier B, Berger M, Maier A, Kachelrieß M, Ritschl L, Müller K, et al. Scatter correction using a primary modulator on a clinical angiography C-arm CT system. Med Phys 2017;44(9):e125–37, http://dx.doi.org/10.1002/mp.12094. [37] Richard NK, Bismark, Oliver Beuing, Georg Rose, International Society for Optics and Photonics. Truncation artifacts caused by the patient table in polyenergetic statistical reconstruction on real C-arm CT data. In: Samuel Matej, Scott D, Metzler, editors. 15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 11072. SPIE; 2019. p. 425–9, http://dx.doi.org/10.1117/12.2534435.

Available online at www.sciencedirect.com

ScienceDirect