Local attenuation curve optimization framework for high quality perfusion maps in low-dose cerebral perfusion CT

Purpose: Cerebral perfusion x-ray computed tomography (PCT) is a powerful tool for noninvasive imaging of hemodynamic information throughout the brain. Conventional PCT requires the brain to be imaged multiple times during the perfusion process, and hence radiation dose is a major concern. The authors propose a PCT reconstruction algorithm that allows for lowering the dose while maintaining a high quality of the perfusion maps. It relies on an accurate estimation of the arterial input function (AIF), which in turn depends on the quality of the attenuation curves in the arterial region. Methods: The authors propose the local attenuation curve optimization (LACO) framework. It accurately models the attenuation curves inside the vessel and arterial regions and optimizes its shape directly based on the acquired x-ray projection data. Results: The LACO algorithm is extensively validated with simulation and real clinical experiments. Quantitative and qualitative results show that our proposed approach accurately estimates the vessel and arterial attenuation curves from only few x-ray projections. In contrast to conventional approaches, where the AIF is estimated based on the reconstructed images, our method computes an optimal AIF directly based on the projection data, resulting in far more accurate perfusion maps. Conclusions: The LACO algorithm allows estimating high quality perfusion maps in low dose scanning protocols. C 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1118/1.4967263]


INTRODUCTION
2][3][4] After an intravenous contrast bolus injection, the brain is scanned several times and from the series of brain volumes, a time concentration curve (TCC) describing the evolution over time of local contrast agent concentration can be extracted at the voxel-level.[8] Despite its clinical value, PCT suffers from a major drawback, which is the high radiation dose submitted to the patient caused by the requirement to scan the patient multiple times after the injection of the contrast bolus. 2,9everal approaches for reducing the radiation dose and/or increasing image quality have already been suggested.By lowering the current of the x-ray source, dose can be reduced in a straightforward way.Reduced beam intensity, however, directly lowers the SNR of the x-ray projections and hence also the SNR of the reconstructed images.1][12][13] This dose reduction technique has two additional advantages: First, it makes iterative reconstruction methods feasible in practice.Iterative reconstruction methods allow a more accurate modeling of the imaging physics, leading to improved reconstructions compared to filtered back projection (FBP).Moreover, they allow combining projection data from all time points to remove data redundancy.Despite these clear advantages, iterative reconstruction methods are infamous for their long computation times and high memory consumption.Reducing the number of projections is therefore highly beneficial both for lowering the computation time and memory consumption.Second, the negative effects of electronic noise increase with decreasing SNR in the projections.This effect has been carefully studied by Ma et al. 14 Opting for less but high SNR projections alleviates this problem.Reducing the number of projections, however, gives rise to limited data artifacts and a decreased SNR that may result in erroneous diagnosis.To avoid such artifacts, a number of tomographic image reconstruction methods were proposed that allow for reducing the number of x-ray projections per time frame (and thus also the radiation dose) without compromising image quality.
A first class of PCT reconstruction methods assumes that a high quality reconstruction is available a priori, either acquired by a nonenhanced prior CT scan or calculated from the time-series data.In this class of methods, the reconstruction is calculated by alternately enforcing data fidelity with the acquired projections and enforcing similarity with the prior image either based on nonlocal means 15,16 or on compressed sensing. 10,17,18n a second class of methods, perfusion specific model assumptions are enforced globally over the reconstructed time series of images.][21] Conventionally, the perfusion maps are calculated only after the series of brain volumes has been reconstructed and, optionally, postprocessed with a noise reduction filter.A well known PCT filter is the time-intensity profile similarity (TIPS) bilateral filter. 22To determine the perfusion maps, the arterial input function (AIF), i.e., the average TCC in the arterial region, is typically estimated from the reconstructed images, without any feedback mechanism to the original xray data.However, the estimation of attenuation curves in the vessels and arteries is often prone to reconstruction artifacts for low dose scanning protocols due to the large attenuation contrast difference with respect to neighboring tissue.This in turn results in substantial errors in the AIF estimation, which ultimately leads to inaccurate perfusion maps.By estimating the TCCs and AIF directly from the x-ray data, the vulnerability to reconstruction artifacts can significantly be reduced.
In this work, we propose the local attenuation curve optimization (LACO) framework.In this framework, the arterial TCCs are modeled by a linear combination of timeshifted gamma variate functions.Next, their coefficients are estimated directly from the measured projection data.As will be demonstrated, this results in a more accurate estimate of the AIF and ultimately leads to more accurate perfusion maps.
Our paper is organized as follows.In Sec. 2, a brief introduction to algebraic reconstruction methods and computed tomography of static and dynamic objects is given.In Sec. 3, the LACO framework is introduced.Simulation and experimental results are reported in Sec. 4. The results are discussed in Sec. 5 and the paper is summarized and concluded in Sec. 6.

TOMOGRAPHY MODEL
In this section, the tomography model used to describe PCT is explained.To that end, first algebraic reconstruction for stationary tomography is revised in Sec.2.A. Next, in Sec.2.B, the concepts of stationary CT are generalized to dynamic CT and the model assumptions for PCT are discussed.

2.A. Stationary tomography model
the pixel values of an (unknown) image of the scanned object.Furthermore, let p = (p i ) ∈ R M denote the log-corrected projection values for all angles θ, which we will refer to as the measured projection data.The projection data can be simulated as Wx, where W = (w i j ) ∈ R M ×N is a sparse matrix with weights w i j , where w i j represents the contribution of pixel value x j to the projection value p i (see Fig. 1).
Directly solving the system of linear equations Wx = p for an exact solution x is typically infeasible, since noise and discretization effects render the system of linear equations inconsistent.Therefore, algebraic methods typically minimize the projection distance ∥Wx − p∥ for some norm ∥ • ∥.The simultaneous iterative reconstruction technique (SIRT) algorithm is known to converge to a solution of argmin x ∥Wx − p∥ 2 R , where R = (r i j ) ∈ R M ×M is the diagonal matrix with inverse row sums of the projection matrix W with r ii = 1/  j w i j . 24Starting from an initial reconstruction x (0) = 0, SIRT iteratively updates the reconstruction as follows: F. 1. Illustration of the projection process.In this image, the contribution w i j of pixel j to the projection value with index i is represented as the ray-intersection length of projection line i with pixel j (picture taken from Ref. 23).
where C = (c i j ) ∈ R N ×N is defined as the diagonal matrix with the inverse column sums of W, that is c j j = 1/  i w i j .

2.B. Dynamic tomography model
In Sec.2.A, the conventional algebraic tomography model was described, which assumes the scanned object to remain unaltered throughout the entire data acquisition process.This assumption is no longer valid in PCT.Hence, algebraic reconstruction of stationary objects has to be generalized to that of dynamic objects.
A time-varying object is typically represented as a time series of images x r ∈ R N , where r ∈ {1,...,R} is the index referring to a particular time frame and R is the total number of time frames.During the acquisition, the gantry rotates multiple times around the object.Denote p r as the subvector of the measured projection data p which is acquired during a rotation of 180 • at time frame r and W r ∈ R M ×N as the corresponding forward projection matrix.
The reconstruction problem for the dynamic case can then be defined as the large system of linear equations where W represents the block diagonal matrix consisting of blocks W 1 ,W 2 ,...,W R and x ∈ R R N represents the vertical concatenation of x 1 , x 2 ,..., x R .In cerebral PCT, the time-varying object x is highly correlated over time.The time attenuation curve of vessels and arteries can typically be described by a linear combination of K functions: the constant function y K = 1 and K − 1 timeshifted gamma variate functions y 1 ,y 2 ,...,y K −1 . 20,25The latter are defined by their simplified form where κ and β are shape parameters and t k is a time shift.The shape parameters and the time shifts are chosen such as to place y 1 , y 2 ,..., y K −1 uniformly over the acquisition time interval. 20Hence we have x r ( j) =  K k=1 a k, j y k (t r ) with j ∈ V , a k, j ∈ R the coefficient corresponding to the kth basis function for the jth pixel and t r the time point corresponding to the rth time index.

METHODS
In this section, the LACO method is described.The problem is formulated as a regularized least squares problem that enforces the model assumption of Sec.2.B in the artery and vessel regions of the reconstruction domain.The LACO framework can be integrated with any iterative reconstruction method, as is illustrated in the flowchart of Fig. 2.
It is assumed that a mask for the artery region, containing L arterial pixels, is known a priori.This mask can contain F. 2. The LACO method is applied at intermediate iterations of an iterative reconstruction algorithm of choice.In the flowchart, the intermediate optimization is performed every kth iteration.
vessel pixels as well, in order to speed up convergence and increase accuracy in those pixels as well.The mask can either by indicated manually or determined automatically. 17,19Since automatic artery/vessel detection is not the main topic of this paper, the artery and vessel region is indicated manually.
For each of the L arterial pixels and K basis functions, define ỹl,k [with (l,k) ∈ {1,...,L} × {1,...,K }] as the timevarying image that is zero everywhere, except in the lth arterial pixel, where the function values of the kth basis function are assigned.The basis functions were defined in Sec.2.B.Furthermore, define ỹ0 to be the time-varying image that is the same as the current reconstruction x in all pixels except in the arterial pixels, where the attenuation curves are replaced by zero.With these definitions, a time-varying parameterized reconstruction ỹ(a) can be defined as where a is a vector containing all coefficients a l,k , i.e., for all L arterial pixels and all K gamma variate basis functions.The intermediate optimization problem can now be written as where M(l) ⊂ {1,...,L} is the set of indices corresponding to the 8-connected neighborhood of the lth arterial pixel (considering only those pixels that belong to the artery region).The parameter µ controls the degree of similarity between coefficients corresponding to the same gamma variate basis function and to neighboring pixels.Let q 0 and q l,k be the simulated projections of ỹ0 and ỹl,k , respectively, i.e., q 0 = Wỹ 0 and q l,k = Wỹ l,k .This allows for rewriting the norm in Eq. ( 5) as The minimization problem in Eq. ( 5) can be solved by setting the partial derivatives (with respect to a l,k ) to zero.The partial derivative of the last term in Eq. ( 5) is given by Adding Eq. ( 7) to the partial derivative of Eq. ( 6) and setting the result to zero yields Since Eq. ( 8) holds for every (l,k) ∈ {1,...,L} × {1,...,K } and is linear with respect to the coefficients in a, this gives rise to the following matrix equation: where M ∈ R L K ×L K and b ∈ R LK are defined via the left-hand side and the right-hand side of Eq. ( 8), respectively.Note that M and all q l,k in Eq. ( 9) can be calculated in advance.This way, only the right hand side b of Eq. ( 9) needs to be updated if the LACO framework is applied at intermediate iterations of the reconstruction algorithm.To speed up the computations, Eq. ( 9) is solved for each connected component of the artery region containing a maximum of 30 pixels.The equations are solved by LU factorization with partial pivoting.

EXPERIMENTS AND RESULTS
The LACO framework was validated with both numerical and real data experiments based on different figures of merit that will be introduced in Sec.4.A.Next, in Sec.4.B, a simulation phantom is introduced and LACO is tested in various experimental settings.Finally, the LACO framework is validated on clinical perfusion CT data in Sec.4.C.

4.A. Figures of merit
Validation was performed based on both the reconstructed attenuation values and the derived CBF (ml/100 ml/min) and CBV (ml/100 ml) perfusion maps.The CBF and CBV maps were calculated with the deconvolution-based truncated singular value decomposition method with a fixed threshold value of 20% of the largest singular value. 5he relative root mean squared error (RRMSE) was utilized as a quality measure and is defined as where x denotes the calculated reconstruction (in which case the summation index i goes over all possible pixels at all possible time points) or perfusion map (in which case the summation index i goes over all possible pixels in the perfusion map) and x denotes the ground truth phantom or perfusion map.For some experiments, the RRMSE was inspected only in a region of interest (ROI), in which case the sum in Eq. ( 10) sums over all points in time and over all pixels in the specific ROI.

4.B. Numerical simulations
In this section, the validation of the LACO framework with a series of numerical experiments is described.
The simulation phantom is a realistic digital brain phantom that was developed by Riordan et al. and extended by Manhart et al. 19,26,27 It is visualized in Fig. 3.
In contrast to classical digital CT phantoms, this phantom has reduced sparsity in the image domain, thereby not favoring reconstruction algorithms exploiting homogeneity (e.g., minimizing total variation).The simulation phantoms were defined on a 256 × 256 isotropic pixel grid, from which projections were generated with 384 detectors each.
For all experiments, the projections were sampled on angles defined by the golden ratio scanning protocol. 28Such angle selection protocol is more flexible than standard protocols in the sense that it allows the user to select an arbitrary number of projections per time frame after the data acquisition, while still approximately covering equiangular positions over the entire angular range.This allows the user to balance the temporal and spatial resolution a posteriori.
Projections were simulated using a fan beam geometry, 29 with a 12.23 • fan angle and a line detector with equidistant detector bins.To avoid the "inverse crime" of generating the data with the same model as the model for calculating the reconstruction, 30 each detector bin was sampled twice for simulating the projection data.Also, Poisson distributed noise was applied to the projection data assuming an incoming beam intensity of I 0 = 2 × 10 4 (photon count) per detector bin.
The LACO approach was implemented in two reconstruction algorithms: the SIRT algorithm and the recently proposed prior image constrained compressed sensing (PICCS) algorithm. 17,18The intermediate LACO optimization was applied every T = 20 iterations for the SIRT algorithm and every T = 2 iterations for the PICCS algorithm.The SIRT algorithm was implemented with a positivity constraint and 500 iterations for all subsequent experiments.The SIRT-LACO algorithm was implemented with a positivity constraint and 200 SIRT iterations.The PICCS algorithm was implemented as described in Refs.17 and 18.The chosen parameters associated with the PICCS algorithm are reported in Table I.
These parameters were optimized by selecting those parameters that yielded the lowest RRMSE with respect to the ground truth image.The prior image was calculated with the vessel-selective prior (see Nett et al. 17 ) with the same vessel mask that was used for the LACO framework.The total variation penalty was implemented as described in Chambolle. 31For comparison, the experiments in this section were also performed with filtered back projection (FBP), a standard noniterative reconstruction algorithm.
][34] Reconstructions were calculated on the same 256 × 256 pixel grid as the ground truth image in order to guarantee a nonsparse phantom and reconstruction.The phantom consists T I.The PICCS parameter values for the different experiments.The α parameter represents the relative weighting parameter in the objective function (Refs.17 and 18).The #iterations SART refers to the number of SART iterations within the PICCS algorithm to enforce data fidelity (Ref.17).The parameter numberSDIterations refers to the number of steps that are taken at each iteration in the steepest descent direction of the objective function (Ref.18).The size of the steepest descent direction (i.e., the negative gradient) is determined by the parameter steepest descent step size, which is multiplied with the negative gradient before adding it to the reconstruction.

Digital brain phantom
Real data α 0. In a first experiment, 50 projections were simulated per time frame.The numerical results for this experiment are summarized in Table II.The AIF estimate can be inspected in Fig. 4. The perfusion maps, obtained with and without the TIPS bilateral filter, 22 can be compared in Fig. 5.The performance of all reconstruction methods for a different number of projections per time frame and a varying amount of noise can be inspected in Figs. 6 and 7.
Furthermore, the perfusion map quality in function of the number of projections was studied for a given radiation dose per time frame.Instead of only varying the number of projections, the photon count was varied in such a way that the radiation dose is the same for each scan (radiation dose: 250 000 photons/time frame).No electronic noise was added.These results can be inspected in Fig. 8.

4.C. Clinical perfusion data
To test the performance of the LACO method on clinical data, a cerebral perfusion CT dataset was acquired with a Discovery CT750 HD (GE Healthcare) scanner.The scan was performed with a source voltage of 80 kVp and xray tube current of 500 mA.Following a bolus injection of 50 ml, with injection rate of 4 ml/s, 912 projections were acquired per 180 • rotation of source and detector.Based on this large amount of projection data, high quality images were reconstructed on a 512×512 pixel grid at 24 time frames with a 2.8094 s interframe temporal distance.These high quality reconstructions are utilized as ground truth images in the subsequent experiment.
Based on the ground truth images, 150 fan beam projections per time frame were simulated.Poisson noise was applied to the projections, assuming a beam intensity of I 0 = 8 × 10 5 (photon count) at the source.Reconstructions were calculated on a 256 × 256 pixel grid assuming a linear projection kernel.The artery region was indicated manually and is displayed in Fig. 9.A vessel mask, also displayed in Fig. 9, was defined as all pixels that are close to the arterial pixels and have a peak attenuation value above a certain threshold.
The numerical results are summarized in Table III and the estimated AIF is displayed in Fig. 10.
The CBF perfusion maps can be compared in Fig. 11.
F. 8.The RRMSE in the dynamic ROI is plotted for the CBF map as a function of the number of projections per time frame for the digital brain phantom.In this experiment, projections were simulated with a radiation dose of 250 000 photons/time frame.

5.A. Brain phantom
From Table II and Fig. 5, it is clear that the LACO framework substantially improves image quality for all reconstruction methods.The RRMSE values in Table II indicate that the LACO framework improves the quality of the TCCs in the vessel region and tissue region (third and last row, respectively) and hence also the quality of the perfusion maps (row 1 and 2).The improved quality of the perfusion map can be attributed to two contributions: First, the more accurate AIF estimate, which is displayed in Fig. 4. Second,  the improvement of the TCCs in the tissue region.Note that the LACO framework does not optimize the tissue TCCs directly.The improvement of these TCCs is the result of more accurate artery TCCs at intermediate iterations, these will guide the whole reconstruction domain to a more accurate reconstruction.Another important observation is that the AIF of the SIRT reconstruction is more severely underestimated in comparison to the FBP reconstruction (see Fig. 4).Note that this effect cannot be reduced by increasing the number of SIRT iterations.It is inherent to SIRT in the sense that the frequency response of this algorithm on undersampled data is very low for the mid and high frequency range. 35This causes underestimation of small high intensity regions, but also reduces the amount of streak artifacts compared to FBP.The positive effect of the LACO framework on the perfusion maps is also visually noticeable in Fig. 5. Again, the same observation can be made for different number of projections per time frame or for different noise levels, of which the results are plotted in Figs. 6 and 7.
Figure 5 shows both the results with and without the TIPS filter.Although the TIPS filter improves the results significantly (especially for FBP and SIRT), the PICCS reconstruction algorithm and LACO framework improve the results even more.If the TIPS filter is combined with the SIRT-LACO, PICCS, and PICCS-LACO methods an increased SNR, though with a slight decrease of spatial resolution, can be observed.
From Fig. 8, it is clear that conventional algorithms such as FBP and SIRT benefit from more projections (with F. 10.The AIF extracted from the different reconstructions in the real data experiment. a lower photon count).Between 10 and 50 projections, the RRMSE of SIRT-LACO and PICCS improves for an increasing number of projections.For 50 projections or more, the RRMSE remains stable.These results indicate that the number of projections can be lowered to 50 projections (while increasing the photon count) without significant quality loss.Note that electronic noise was not considered during this experiment.Adding electronic noise would have a bigger negative impact on the results with low photon count per projection (and thus many projections).These observations show that careful undersampling can be beneficial for the reconstruction quality of low dose CT with the additional benefit of shorter computation times and lower memory consumption.The same conclusions can be drawn from the results of PICCS-LACO.Here the number of projections can even be lowered to only ten projections without compromising the reconstruction quality.

5.B. Clinical perfusion data
As can be seen from the numerical results in Table III, the LACO framework greatly improves image quality.This can also be seen on the CBF maps in Fig. 11, where the SIRT-LACO reconstruction can discriminate between the two different structures that are indicated by the gray arrow in Fig. 11, whereas the SIRT reconstruction cannot.The AIF's estimate, displayed in Fig. 10, shows again a great improvement of the SIRT-LACO reconstruction in comparison to the SIRT reconstruction.As the AIF's estimate derived from the PICCS reconstruction is already quite accurate, the PICCS-LACO reconstruction produces approximately the same result.

CONCLUSION
In general, the radiation dose in PCT experiments can be decreased by lowering the number of acquired projections per rotation of the gantry.Standard approaches that reconstruct the object independently at different time frames, however, result in images containing artifacts that are introduced by the lack of projection data.
In this paper, the LACO framework was presented.It is an iterative approach tailored specifically to the reconstruction of cerebral PCT images.The algorithm exploits prior knowledge about the attenuation curves by modeling it as a linear combination of time-shifted gamma variate functions and optimizing its corresponding coefficients directly based on the projection data.This results in a significantly increased accuracy of the vessel's attenuation curves, of the derived AIF, the tissue attenuation curve, and consequently also in increased quality of the estimated perfusion maps.
The LACO framework was validated with simulation experiments and clinical data (based on a high quality reconstruction of a Discovery CT750 HD scan).These experiments illustrate that, in comparison to standard approaches, the LACO framework significantly improves image quality for the same radiation dose and has similar image quality for a lower radiation dose.

F. 3 .
The simulation phantom is adopted from Manhart et al. (Refs.19 and 26).Some exemplary TCCs are displayed in (b), corresponding to the indicated pixels in (a).However, all TCCs vary from pixel to pixel.This way, this phantom does not favor sparsity exploiting algorithms.It was created by starting from a CBF and CBV map (with manually indicated reduced and severely reduced tissue regions) and calculating the TCCs based on these maps.By adding random perturbations on the CBF and the extracted MTT (=CBV/CBF) map, the final TCCs have reduced sparsity characteristics.In (c), the dynamic ROI is indicated in blue and the artery/vessel ROI is indicated in red.(See color online version.)Medical Physics, Vol.43, No. 12, December 2016

F. 4 .
The AIF extracted from the different reconstructions of the brain phantom.Physics, Vol.43, No. 12, December 2016 F. 5. Comparison of the extracted CBF maps (ml/100 ml/min) for the different reconstructions of the brain phantom.The top row contains the CBF maps and the bottom row the difference with respect to the ground truth image.In this experiment, 50 projections were simulated per time frame assuming an incoming beam intensity of I 0 = 2 × 10 4 (photon count) per detector bin.F. 6. RRMSE values as a function of the number of projections per time frame for the digital brain phantom.The RRMSE in the dynamic ROI is plotted for CBF map (a) and CBV map (b), the vessel/artery ROI for the TTCs (c), and the dynamic ROI without the vessel/artery ROI (the tissue ROI) for the TTCs (d).These ROIs are indicated in Fig. 3(c).In this experiment, projections were simulated assuming an incoming beam intensity of I 0 = 2 × 10 4 (photon count) per detector bin.F. 7. RRMSE values as a function of the incoming beam intensity I 0 for the digital brain phantom.The RRMSE in the dynamic ROI is plotted for CBF map (a) and CBV map (b), the vessel/artery ROI for the TTCs (c), and the dynamic ROI without the vessel/artery ROI (the tissue ROI) for the TTCs (d).These ROIs are indicated in Fig. 3(c).In this experiment, 50 projections were simulated per time frame.

F. 9 .
The masks that were utilized in the clinical perfusion data experiment.Red indicates the artery ROI, yellow the artery/vessel mask, and light blue indicates the dynamic ROI.(See color online version.)Medical Physics, Vol.43, No. 12, December 2016 T III.RRMSE in the dynamic region for the clinical perfusion data experiment for perfusion maps and attenuation values (rows) of different reconstruction algorithms (columns

F. 11 .
Comparison of the CBF maps for the different reconstructions in the real data experiment.The top row contains the full CBF maps and the bottom row a zoomed in portion of the CBF maps.Medical Physics, Vol.43, No. 12, December 2016 T II.RRMSE values for perfusion maps evaluated in the dynamic region of the brain phantom (row 1 and 2), RRMSE of the TTCs evaluated in the artery/vein region and tissue region (third and last row) and for different reconstruction algorithms (columns).In this experiment, 50 projections were simulated per time frame assuming an incoming beam intensity of I 0 = 2 × 10 4 (photon count) per detector bin.The numbers in bold indicate the lowest achieved RRMSE value for the given perfusion map and ROI.
).The numbers in bold indicate the lowest achieved RRMSE value for the given perfusion map and ROI.