Inferring CT perfusion parameters and their uncertainties using a Bayesian approach

— Objective: Perfusion CT imaging is commonly used for the rapid assessment of patients presenting with symptoms of acute stroke. Maps of perfusion parameters such as cerebral blood volume (CBV), cerebral blood flow (CBF), and mean transit time (MTT) derived from the perfusion scan data, provide crucial information for stroke diagnosis and treatment decisions. Most vendors implement singular value decomposition (SVD)-based methods on their scanners to calculate these parameters. However, some known issues are associated with conventional methods. Methods: In this work, we propose a Bayesian inference algorithm, which can derive both the perfusion parameters and their uncertainties. We apply the variational technique to the inference, which then becomes an expectation-maximization problem. The probability distribution (with Gaussian mean and variance) of each estimated parameter can be obtained, and the coefficient-of-variance is used to indicate the uncertainty. We perform evaluations in both simulations and patient studies. Results: In a simulation, we show that the proposed method has much less bias than conventional methods. Then, in separate simulations, we apply the proposed method to evaluate the impact of various scan conditions, i.e. with different frame intervals, truncated measurement, or patient movement, on the parameter estimate. In a patient study, we further test the case where patient movement is present; we use the proposed method to evaluate how well existing motion correction methods recover a motion-contaminated stroke scan. Conclusion: The proposed method can be used to evaluate confidence in parameter estimation and the scan protocol design. More clinically relevant studies are required to test the proposed method.

and recovery. To aid the diagnosis and selection of the appropriate treatment, the location and volume of the infarct core and penumbra are often calculated -infarcted core indicates the brain region that is not salvageable, while penumbra indicates the brain region that is at risk of progression to infarction. The penumbra usually surrounds the ischemic core and still may be salvageable. Classification of the infarct core and penumbra is often derived by thresholding certain brain hemodynamic parameters, e.g. cerebral blood volume (CBV), cerebral blood flow (CBF), mean transit time (MTT) and time to peak (TTP). For example, in [3], the ischemic core was defined as regions with relative MTT ≥ 145% of the normal tissue and absolute CBV <2.0 ml/100g.
Perfusion CT and dynamic contrast-enhanced MR are the primary imaging techniques to obtain the brain hemodynamic parameters. We are interested in perfusion CT in this study, which is performed after the administration of the intravenous iodine contrast and can monitor the first pass of a contrast bolus through the cerebral circulation. In a perfusion CT scan, mask image is first acquired at the beginning. The enhanced image at single time point can be obtained by subtracting the mask image from a later acquired image. Enhancement of tissues depends on iodine concentration which at one voxel g can be represented as: where h(t) is the flow-scaled residual functionthe fraction of contrast that remains in the voxel at a time t after its arrival. Ca (t) Inferring CT perfusion parameters and their uncertainties using a Bayesian approach Tao Sun * , IEEE Member, Zhanli Hu, IEEE Senior Member, Roger Fulton, IEEE Senior Member S Fig. 1 (a) The indicator dilution theory model describes the blood supply to tissue at one voxel in brain. After contrast injection, the bolus enters voxel via arterial inlet, is dilated into the capillary bed and leaves the voxel via the venous outlet. Due to the limited spatial resolution of the CT image, Ct(t) will contain numerous capillary beds as well as arterioles and venules. (b) h(t) is the residual function describing the hemodynamic functions in the physiological model. CBF, CBV and MTT can be derived from the h(t).
2 is the arterial input function (AIF) which indicates the contrast change in the arteries, and * indicates the convolution operation.
(1) describes the well-known indicator dilution theory ( [4], Fig.  1). Based on this theory, from a perfusion CT scan one can obtain the h(t), which can be used to derive the hemodynamic functional parameters related to the blood passage in the tissue, including CBF, CBV, MTT and TTP: MTT CBV / CBF = TTP arg max ( ) t gt = where ρ is the mean density in the capillary bed, which is set to 1.04 g/ml in this study [5]. A detailed explanation of how (2) is derived can be found in [6]. A most straightforward way to obtain the residual function h(t) is to perform deconvolution on (1) [7]. In practice, g(t) is sampled at discrete time points. Thus, the deconvolution problem is an inverse, algebraic problem, which can be solved by SVD. It was shown that such an inverse problem is ill-posed [6]. Thus, with SVD it is often necessary to truncate a range of singular values below a threshold to suppress the noise. In the clinical setting, perfusion analysis is usually done with SVD-based methods, the implementation of which varies across vendors.
On the other hand, SVD-based methods have some known problems. First, they make assumptions about the ideal underlying physiological model. For example, the residual function h(t) is assumed to be an impulse function, which is not valid for a real scan. Because of this, the estimated parameters can be biased. As shown in [8], [9], the SVD-based method can underestimate CBF and overestimate MTT. Several authors proposed advanced deconvolution methods, e.g. include regularization (e.g. Tikhonov) to address this problem [10]- [12]. The second issue is that these methods often cannot tolerate imperfect measurement well, e.g. relatively long interval scan intervals [13], [14] or truncated measurement in time [15], [16]. Thirdly, these methods are deterministic -they provide an estimate of the mean of the parameter of interest but not its probability distribution. Factors that affect the reliability of the measurement in a perfusion CT scan, e.g. protocol design and perturbance of the data, therefore cannot be accounted for in the estimation. For example, a too-long scan interval reduces the available samples for estimation and hence the accuracy; in other cases, the patient movement could produce inconsistent voxel attenuation measurements that would inevitably increase the uncertainty of the estimation; post-processing (e.g. spatial and temporal filtering) may further introduce uncertainty into the estimation. For quantification, it is desirable to account for such randomness with proper probability theory [17]. Therefore, a probabilistic inference method seems warranted to infer both the hemodynamic parameters and their associated uncertainties. And the method should provide less bias in estimation than conventional methods.
In this paper, we propose such a method and show its potential value in clinical application. We formulate the estimation problem as a Bayesian inference problem. A variational technique, mean-field approximation, was applied, which makes the problem more tractable as a classic expectation-maximization (EM) problem. The Gaussian distribution (with mean estimate and variance) of CBF and MTT can be obtained. Their coefficient-of-variance (CoV) is used to indicate the uncertainty in inference, which reflects the reliability of a given measurement. We will first demonstrate the superiority of the proposed method over conventional method in simulation studies where ground truth is available. We will also demonstrate the potential usage of the proposed method on evaluating different suboptimal scan conditions, including a patient scan contaminated with motion. The results are discussed at the end of this paper.

II. METHODS
We first present the forward acquisition model of perfusion CT. Then we give the details of the proposed Bayesian inference algorithm and the quantitative analysis design.

A. Forward acquisition model
Starting from (1), the measurement process in a CT perfusion scan can be represented as: where Y(t) is the enhanced tissue time curve at voxel-of-interest, µ represent the hemodynamic parameters to be estimated, h(µ,t) is the flow-scaled residual function, ε(t) is the noise. The flow-scaled residual function h(µ,t) can be further written as a function of the CBF, tissue density ρ and residual function r(t): According to [18]- [20] r(t) can be modelled as a box shape function with exponential decay. Hence, h(µ,t) becomes: where T 0 = 0.632 × MTT . In this study, we consider µ contains two parameters -CBF and MTT. Other parameters, e.g. CBV, can be calculated as CBV=CBF×MTT, and TTP can be derived directly from the time-attenuation-curve (TAC) as the time when peak enhancement is reached as in (2). There is no bolus arrival delay modelled in this study. The convolution in (3) can be expanded by sampling the integral of the arterial enhancement curve: As explained in the introduction, Y(t) in practice is sampled at N discrete time points. Without knowledge of the measurement errors, ε is then modelled as zero-mean Gaussian distributed at each frame. Since the frames are temporally independent, we have the discrete version of (5) as: where the index of the time frame is i, and is proportional to the number of counts (duration of the ith frame as a surrogate). From (6) we can build the data mismatch term, and the inference problem becomes finding the estimate that maximizes the log-likelihood: In the case where the length of each frame is the same, one can set = 1. The above optimization problem is essentially equivalent to maximum likelihood estimation, which can be solved by nonlinear least squared fitting, where the way to perform the inference is similar to that of previous perfusion estimation methods [21].

B. Bayesian inference algorithm
Now we extend the above maximum likelihood estimation to a full Bayesian inference -which can provide the probability distribution rather than just the maximum likelihood estimate of a given parameter in (7). The derived distribution will inherently contain the measurement uncertainty of each parameter. According to Bayes' theorem, the posterior of the estimate µ over measurement Y is: where P(Y| µ) is the likelihood in (7), and P(µ) is prior. The fundamental difference from (7) is that µ no longer represents a maximum likelihood estimate but a full distribution. Estimating µ within the Bayesian framework becomes an inference problem, which involves maximizing the posterior probability P(µ|Y) of µ given the observed data Y. In practice, it is often difficult to compute the posterior P(µ|Y) in a closed form (or we say intractable). An approximation is therefore often used, for which the Markov Chain Monte Carlo (MCMC) and Variational Bayesian (VB) are two popular choices. MCMC is too computationally expensive for our problem given the volume size of the perfusion CT images. The approximate but efficient VB method is therefore used here. The basic idea is to find a simple analytical distribution q(µ|Y) to approximate the intractable posterior P(µ|Y), such that the Kullback-Leibler (KL) divergence of these two is minimized. Given the probabilistic model (8), the log evidence can be written as: where P(Y,µ) is the joint probability, E* denotes the expectation on q(µ|Y). ELBO is evidence lower bound while KL is the divergence between approximate q(µ|Y) and P(Y,µ). Because the KL-divergence is always positive, ELBO provides a lower bound on the log-likelihood evidence. So one can maximize ELBO by finding µ that approximates the posterior close to the true posterior (or minimizes KL): Now we have a second KL-divergence in the square bracket which is always positive, and the optimal solution can be found by equating the numerator and denominator.
Considering the numerator, we need first to ensure it is tractable. VB factorizes q(µ|Y) (in physics this is known as mean-field approximation), which collects the parameters into separate groups each with their own approximate posterior [22]. Note that this approximation does not imply that the parameters are uncorrelated for a given measurement. A simple choice of the distribution of the approximated true posterior is multivariate Gaussian. The mean-field theory assumes that the posterior distribution of each parameter is Gaussian. Supposing that we have a number, k, of parameters to be estimated, the approximate posterior has the analytical expression: Then for the denominator, we can insert the likelihood from (7). As for the prior, we choose the conjugate prior -the prior is said to be conjugate to the likelihood if and only if the factorized posterior has the same parametric form [23]. Therefore, the conjugate prior has multivariable Gaussian distribution with mean m0 and variance σ0. The numerator thus becomes: To ensure tractability and allow the method to generate to any nonlinear model, g(µ) is approximated by a first-order Taylor expansion about the mode (mean) of the posterior: where J is the Jacobian. With such linearization applied on g(µ), (11) becomes to: By equalizing (12) and (10), we have update equations for mean m and variance σ: The above derivation matches the description in [23], [24]: under certain conditions, a Bayesian inference problem reduces to an optimization problem, of which the optimal solution can be found by an algorithm similar to expectation maximization. In our case, because of the independent assumption of estimated parameters, each parameter in µ can be updated simultaneously, which iteratively maximizes the KL. Hence the parameter 2 k k k μ N(m ,σ ) that approximates q(µ|Y) to P(µ|Y) can be found, in which σk contains the uncertainty information in the k th parameter estimated from a given scan. Instead of directly using the variance σ 2 , we used the coefficient of variation, which has the advantage of being independent of the intensity scale. It can simply be represented as a ratio between the standard deviation and the mean of a given parameter: , which is unitless. One can conclude that the higher the CoV, the more uncertainty is associated with the parameter.

C. SVD methods
The deconvolution method was performed as illustrated in [25]. To suppress noise in the resulting residual function, the least significant eigenvectors were consecutively removed until h(t) had an oscillation index below a certain threshold. In this study, a threshold of 0.095 was used as suggested in [25]. CBF, MTT and CBV were calculated according to equations (2). An alternate of deconvolution methods with Tikhonov regularization was implemented as well. Tikhonov regularization was implemented using the Regularization Tools package [26]. The parameter selection and regularization setup was as suggested in [10].

A. Simulation study
We simulated 3D (2D+t) perfusion CT scans using a modification of the brain phantom described in [27] (downloadable at https://www5.cs.fau.de/research/data/digitalbrain-perfusion-phantom/index.html). We defined two annotated tissue classes: healthy tissue, and tissue with reduced CBF and increased MTT (Fig. 2a). Perfusion parameters, as defined in Table I were assigned to the annotated tissues. From the perfusion phantom, we could generate dynamic projection data based on the forward model defined in (4) and (5). A total of 50 1-s frames were simulated, for which TACs are shown in Fig. 2b. Data were generated by forward projecting each frame, with an acquisition protocol (Table II) based on a simulated scanner model. Poisson noise was added to the projections by setting the number of detected counts in a blank scan to 1e5. Noise-free and noisy projections were reconstructed to generate two sets of acquired dynamic images, from which perfusion parameters were inferred. In all simulations, reconstruction was performed with a helical Feldkamp-Davis-Kress (FDK) algorithm [28], [29]. A total number of four simulations were performed to evaluate the proposed method: (1) In the first simulation, we compared the proposed method with the conventional methods with ideal simulated measurements. The noise-free dynamic images were used to compute the perfusion maps with both methods, and the bias was assessed by comparing image profiles across the lesion in perfusion maps generated from the SVD-based methods, the proposed method, and the truth. Note that only this simulation uses the noise-free images, and the rest used noisy simulated images. (2) In the second simulation, we evaluated the scan protocols with different intervals using the proposed method. It was shown previously that scan interval variation could have an impact on perfusion imaging. Altering the acquisition interval may introduce a bias in the perfusion parameters [13], [14]. We varied the scan intervals by undersampling the simulated noisy dynamic frames. For a normal setup, the scan interval was zero; for long-interval setups, the scan intervals were either 1 sec or 3 sec, reducing the number of frames to 25 and 12, respectively. The hypothesis was that Bayesian method can tolerate better with large interval than the SVD-based method does. Moreover, the measurement would be less reliable when fewer data were available, i.e. that the proposed Bayesian method would identify larger uncertainties in parameters. We also calculated the predicted enhanced TAC for an ROI in the lesion for Bayesian methods by applying the estimated perfusion parameters to forward model in (6). Then this enhanced TAC was compared with the true one. (3) In the third simulation, we evaluated the effect of truncated measurement on inference using the proposed method.
Truncation is a common issue in perfusion CT when patients have poor cardiac output or large vessel occlusion, in which case the acquisition is too short of capturing the complete attenuation change [30]. It has been shown that incomplete acquisition could have a negative impact on perfusion estimation. Decreasing scan duration falsely reduced perfusion estimates, hence may alter patient management [15], [16]. To simulate such an effect, we made the last 20 or 30 s of measurement unavailable for the parameter inference. The hypothesis was that the proposed Bayesian method can better handle missing data than the SVD-based method does. Moreover, the measurement would be less reliable because of the missing data, i.e. that the proposed Bayesian method would identify larger uncertainties in parameters. We also calculated the predicted enhanced TAC for an ROI in the lesion for Bayesian methods and compared them to the true one. (4) In the fourth simulation, we evaluated the effect of the motion on parameter inference using the proposed method. Involuntary patient motion can create inconsistencies in temporal measurements and hence have a negative impact on perfusion map estimation [31], [32]. We examined the effect of motion on the uncertainty of parameter estimates. For this investigation, we generated simulated motion-corrupted projections of the phantom by applying the segment of motion in Fig. 2c during forward projection, and reconstructed the motion-corrupted images. Perfusion maps derived from both the motion-free and motion-corrupted dynamic images were compared. The hypothesis was that the measurement would be less reliable in the motion case because of inconsistent data, i.e. that the proposed Bayesian method would identify larger uncertainties in parameters estimated from motion-corrupted data.

B. Patient study
The anonymized raw data of a patient who had previously undergone head CT scan in the Department of Radiology at Westmead Hospital, Sydney, Australia were collected with the approval of the Human Research Ethics Committee of the Western Sydney Local Health District. The scan was acquired on a Siemens Force scanner (Definition Healthineers, Forcheim, Germany). Before the scan, 45 mL of iodine contrast agent was injected intravenously at a rate of 6 mL/s, followed by a 40 mL saline flush. After that, a total of 33 frames (26×1.5s, 7×3s) were acquired in axial (shuttle) mode. The first acquired frame was labelled as the baseline frame (mask). The reconstruction was initially performed with the vendor's filtered back-projection algorithm with a smoothing kernel of H40f. The pixel dimensions of the reconstructed image were 0.4×0.4×5 mm 3 . DICOM volumes were exported without further post-processing.
The exported dynamic reconstructed frames presented both inter-and intra-frame motion artifacts. These images were labelled as NMC. The patient motion would be expected to bring the uncertainty of parameter estimates, while the reduction in motion artifacts and blurring following application of a motion correction algorithm may recover the reliability of parameter estimates. As a test of the proposed method, we aim to apply the proposed method to evaluate the effectiveness of motion correction methods, in terms of the perfusion parameters.  Two motion correction methods were implemented in this study. The first only compensated for inter-frame motion, by directly operated on the exported images [33], [34]. All dynamic frames were rigidly registered to the first (baseline) frame to remove the inter-frame motion. The corrected images were labelled as MC1. The second performed an additional step before inter-frame motion compensation, which compensated for intra-frame motion using the method described in [35]. This iterative data-driven method jointly estimated the six degree-of-freedom rigid motion of the head at each projection and the motion-corrected reconstruction from the raw projection data. At each iteration, head pose estimates were updated analytically at multiple projection angles by adjusting head pose to reduce differences between forward and measured projections. The motion-corrected image was updated by taking the estimated motion into account during reconstruction with a fully 3D Maximum Likelihood Expectation Maximization (MLEM) or FDK algorithm. Motion correction was applied to each frame independently. Then an inter-frame motion correction was followed as we did for the first correction method. The corrected images were labelled as MC2.
After motion correction, perfusion analysis was performed using the proposed Bayesian inference algorithms for NMC, MC1 and MC2 images. Unlike simulation study, AIF was unknown, and it was derived directly from the images as follows. We first applied a threshold to the summation of all frames and drew a 2D region-of-interest over the internal carotid region. We then corrected the AIF for the partial volume effect due to the small diameter of the target arteries -a venous output function (VOF) was derived similarly in a straight sinus region, and the AIF was scaled to have the same area under the curve as the VOF to account for the underestimation [36]. After inference, the hematocrit correction was performed on CBV and CBF (multiplying with a constant 0.733 [5]) to account for the different hematocrit values between arteries and capillaries.

C. Quantification
All generated perfusion parameter images are compared with the truth in Fig. 2a by calculating means-square-error (MSE) difference and the bias. Assuming IM0 is the truth image and IM1 is the image to compare with, we have MSE between these two images: where N is the total number of voxels, j is the index of the voxel. The bias between two images can be represented as: where calculation was only performed in a region-of-interest (ROI) drew inside the lesion. Average CBF and MTT intensity distribution were also calculated for that ROI, which has Gaussian shape. The mean of the distribution was the mean at all voxel in ROI; the standard deviation of the distribution was the mean of the square root of the variance at all voxels in ROI.
The wider the distribution is, the greater the uncertainty of the estimated parameters in ROI is.

A. Simulation study
In the first simulation, we compared the proposed method with the deconvolution methods in terms of the accuracy of the derived perfusion maps (Fig. 3a). MSE was measured for the comparison. One ROI was labelled in lesion regions in the estimated MTT images. Bias analysis was performed for this ROI. Compared to the SVD-based method, the proposed Bayesian method resulted in lower bias for CBF and MTT (hence CBV). The profile plots corresponded well to the visual appearance and quantification (Fig. 3b).
In the second simulation, we assessed the impact of varying the scan interval on estimation. Quantification measurement of MSE and bias were assessed. The SVD-based method cannot tolerate these improper protocols, but the proposed method can still produce relatively unbiased estimates. With larger scan intervals, we observed that the mean CBF and MTT maps were noisier than with normal sampling (Fig. 4). The CoV maps and the CBF distribution plots indicated differences in CBF mean and variance with different scan intervals. We also calculated the predicted enhanced TAC sampled in a lesion ROI for all methods in Fig.  5a.
In the third simulation, our observations in Fig. 6 were similar to those of Fig. 4. The SVD-based method cannot tolerate these imperfect scans, but the proposed method can still produce relatively unbiased estimates. With severely-truncated measurement data, the perfusion estimate was biased and its    The average CBF intensity distribution of a sampled white matter region. MTT has similar results which were not shown. (c) A segment of realistic rigid motion [46] that was applied. Note that only in-plane motion was simulated (one rotation, two translations) in this 2D+t simulation.  variance higher, as indicated by the CoV images, compared to the case with full measurement data. We also calculated the predicted enhanced TAC sampled in a lesion ROI for all methods in Fig. 5b.
In the fourth simulation with simulated motion, we observed that the mean CBF and MTT map had noticeable visual artifacts (Fig. 7a). The CBF_CoV images indicated that motion indeed made the scan less reliable. The high CoV values corresponded to regions with artifacts in the mean estimate maps, and CoV values were generally lower in the motion-free image. This suggests that CoV mapping is useful for identifying parts of the brain where parameter estimates are problematic due to motion-induced uncertainty.

B. Patient study
MC1 and MC2 were generated as in Methods section. The goal is to apply the proposed method to evaluate the effectiveness of two motion correction methods, in terms of recovering the perfusion parameters. For both CBF and MTT, mean estimates and CoV were derived with MC1, MC2 and noMC. Fig. 8a indicates that CBF from NMC has severe artifacts due to motion presented in dynamic frames, both by observing the mean and CoV images. MC1 successfully recovered the quality of the dynamic scan and the inferred parameters were more reliable than with ones noMC, although there are still some residual artifacts. MC2 further recovered the images, and the reliability of the estimation was the highest among these three. The results correspond well to the hypothesis of progressive recovery of the correction methods (from noMC to MC1 and MC2).

V. DISCUSSION
CT perfusion imaging results can be affected by many factors in practice, including the selection of scanner acquisition protocols, the choice of post-processing technique, data truncation, and data perturbance caused, for example, by patient movement. All may affect the reliability of hemodynamic parameter estimation, which is crucial for the diagnosis and treatment plan of stroke patients. Conventional inference methods cannot tolerate this imperfection well, and cannot assess the reliability of the data although such information could be useful diagnostically. We proposed a Bayesian method which can calculate CT perfusion maps and estimate the uncertainty at each voxel. Compared with the conventional SVD-based methods, this approach has the advantages of (1) estimating parameters with less bias when comparing to conventional methods, and (2) providing probability of an estimate. It can be used to evaluate a scan protocol design or to compare different data processing or artifacts reduction methods. We applied the proposed method to three simulated abnormal scan scenarios, i.e. with increased frame intervals, truncated measurement, and patient movement. The results indicated that the uncertainty measure provided by the method allows scan protocols of varying reliability to be differentiated. In a patient study, we evaluated the proposed method on a stroke scan in which the patient exhibited severe motion. Two motion correction techniques were compared in terms of how much they improved the scan reliability, as assessed by the proposed Bayesian method.
The proposed method successfully provided mean and uncertainty estimates for CBF and MTT. However, it is also possible to infer other parameters such as bolus arrival time delay and permeability, if the associated acquisition model is assumed. For example, a delay exists between the injection time and the contrast agent arrival time in the tissue of interest. To include this delay, the convolution model in (1) could be modified in future work, and the inference proceeds with minor change. Permeability is an important indicator of the degree of blood-brain barrier damage. The Patlak model is the common choice to infer this parameter and could be incorporated with the proposed Bayesian method. However, a longer scan time (~200 s) is often needed to provide enough information for the Patlak model, for which the scan protocol is not available at this moment in our center.
Another possible way to infer the distribution of the perfusion parameters is the frequentist method, which often requires repeatable measurements (which is not feasible on patients) or simulated measurements using techniques such as bootstrapping. We chose a Bayesian method in this study instead, since we believe it is more appropriate for our application where the measured data are fixed, and parameters are random rather than the opposite. Also, the Bayesian method has the advantage that prior knowledge can easily be included in the inference process, while it is not obvious how to do this for bootstrapping [37]. The ability to include spatial-temporal regularization would further enable the proposed method to overcome the high radiation dose in a perfusion CT scan.
In this work, we inferred the parameter distributions from dynamic reconstructed images. In principle, it is possible to infer the desired parameters directly from the measured raw projections without reconstruction. Applications of directly estimating the kinetic parameters exist in CT [38], PET [39], [40] and MR [41], [42]. However, these approaches can only provide the maximum likelihood estimate of parameters, not the full distribution. Another concern is that it is difficult to separate the coupling between the attenuation and kinetic distribution since a perfect forward model can often not be guaranteed. For example, if there is patient motion during a scan, simultaneous motion correction and perfusion map inference can be problematic because errors from any side will likely propagate to another side.
A limitation of our method is, in our measurement model, the observation noise is assumed to be an i.i.d. Gaussian variable. Within each frame, the noise level was scaled according to the duration of that frame. This assumption leads to great convenience for computation. If the CT scan dose is very low, the actual noise model can be complicated [43], [44] and the one assumed here may yield poor inference results. Other inference techniques may be required to consider the correct noise model. For example, techniques approximating the Poisson to Gaussian distribution may be required before applying Bayesian inference. Another limitation is we did not evaluate the estimated parameter maps for determining the infarct core and penumbra region, which could be more relevant to a clinical decision. A more detailed investigation of the clinical impact of the method is needed.

VI. CONCLUSION
We propose and apply a Bayesian approach which infers the hemodynamic parameters in perfusion CT imaging. The method can yield mean parameter estimates and further provide uncertainty estimates of such parameters under imperfect scan scenarios. The proposed method can be used to assess the confidence in parameter estimation under imperfect scan conditions and used to evaluate the protocol design. Future clinically relevant studies are required to evaluate the proposed method.