A Spatio-Temporal Deconvolution Method to Improve Perfusion CT Quantification

1182

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 5, MAY 2010

A Spatio-Temporal Deconvolution Method to Improve Perfusion CT Quantification Lili He*, Burkay Orten, Synho Do, W. Clem Karl, Avinish Kambadakone, Dushyant V. Sahani, and Homer Pien

Abstract—Perfusion imaging is a useful adjunct to anatomic imaging in numerous diagnostic and therapy-monitoring settings. One approach to perfusion imaging is to assume a convolution relationship between a local arterial input function and the tissue enhancement profile of the region of interest via a “residue function” and subsequently solve for this residue function. This ill-posed problem is generally solved using singular-value decomposition based approaches, and the hemodynamic parameters are solved for each voxel independently. In this paper, we present a formulation which incorporates both spatial and temporal correlations, and show through simulations that this new formulation yields higher accuracy and greater robustness with respect to image noise. We also show using rectal cancer tumor images that this new formulation results in better segregation of normal and cancerous voxels. Index Terms—Deconvolution, perfusion computed tomography (CT), singular value decomposition, spatial and temporal correlations.

I. INTRODUCTION ONVENTIONAL medical imaging, such as computed tomography (CT), magnetic resonance (MR), and X-ray radiography, provides anatomical structure information. Since many pathologies or injuries impact not only structure but also function, researchers have sought to expand the types of information that can be obtained using these modalities. One such approach is perfusion imaging—where hemodynamic parameters are obtained as the mixture of blood (endogenous or exogenous) and contrast agent is “perfused” throughout an anatomical

C

Manuscript received January 11, 2010; revised February 09, 2010; accepted February 10, 2010. First published April 08, 2010; current version published April 30, 2010. This work was supported in part by grants from the Center for Integration of Medicine and Innovative Technology (CIMIT), U.S. Army Medical Research, and also in part by the National Science Foundation (NSF) Engineering Research Centers (ERC) grant to the Center for Subsurface Sensing and Imaging Systems (CenSSIS). Asterisk indicates corresponding author. *L. He is with the Laboratory for Medical Imaging and Computations, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Boston, MA 02114 USA (e-mail: [email protected]). B. Orten is with the Department of Electrical and Computer Engineering, Boston University, Boston, MA 02215 USA (e-mail: [email protected]). S. Do and H. Pien are with the Laboratory for Medical Imaging and Computations, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Boston, MA 02114 USA(e-mail: [email protected]; [email protected]). W. C. Karl is with the Department of Electrical and Computer Engineering and Biomedical Engineering, Boston University, Boston, MA 02215 USA (e-mail: [email protected]). A. Kambadakone and D. V. Sahani are with the Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Boston, MA 02114 USA (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMI.2010.2043536

region. Both MR and CT perfusion (MRp and CTp) have been found to be useful in a number of settings, including diagnosis, risk stratification, and therapeutic monitoring [1]–[5]. There are currently several different perfusion models in use, depending on the phase of contrast enhancement and the hemodynamic parameter to be estimated. In this paper, we focus on characterizing tissue perfusion by examining the initial passage of contrast (i.e., first-pass enhancement). Even with first-pass, several perfusion models have been developed; the two most widely studied classes of models are those that examine the rate of enhancement (maximum slope models), and those that seek to estimate a “residue” function convolution kernel that relates the input (arterial enhancement) and response (tissue enhancement) [1], [4], [6]. This second model forms the basis of the analysis of dynamic susceptibility contrast (DSC) MRI, is widely used in CTp analysis, and is the focus of our study. More specifically, the deconvolution method is based on the central volume principle and the assumption that the arterial input function (AIF) and tissue time-enhancement curves (TEC) are related through a convolution integral by the residue function, where the residue function describes the tracer remaining in the tissue at any instant of time [7]–[12]. As such, significant past research have focused on how best to deal with the inverse problem of estimating the residue function from the AIF and TEC. In this context, the most commonly used deconvolution method is truncated singular value decomposition (TSVD) and its extensions, such as circulant TSVD (cTSVD) [8], [9], [13]–[16]. Others, however, have noted the tendency of TSVD-based methods to introduce unwanted oscillations [17], [18]. This has led researchers to examine different regularization methods to stabilize the deconvolution, and have shown varying degrees of success in recovering the residue function [17], [19]–[21]. It is worth noting that these past studies have focused exclusively on imposing temporal regularization. Since perfusion images tend to be noisy, our aims are to develop a method to perform deconvolution-based first-pass hemodynamic parameter estimation that is more robust to noisy inputs, and to produce perfusion parameter maps with better signal-to-noise characteristics. To that end, we have developed a formulation which utilizes a Mumford–Shah-like functional to enforce both spatial and temporal correlations, and additionally terminate these correlations when a discontinuity between regions appears likely [22]–[24]. Because TSVD-based approaches estimate the residue function (and hence the perfusion parameter) for each voxel independently of its neighbors, our spatially- and temporally-correlated deconvolution approach mitigates the noise issues associated with traditional approaches. Furthermore, it is frequently the

0278-0062/$26.00 © 2010 IEEE Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

HE et al.: A SPATIO-TEMPORAL DECONVOLUTION METHOD TO IMPROVE PERFUSION CT QUANTIFICATION

case that an ad hoc moving average filter is used to postprocess the noisy TSVD-derived perfusion parameter maps; our formulation eliminates the need for such a step, potentially preserving the sharpness between adjacent regions. Although spatially-correlated image models have been used in several contexts (cf. [25]–[28]), to date we are not aware of any such work in the context of perfusion parameter estimation. In this paper, we utilize our spatio-temporal deconvolution method to estimate blood flow (BF). We utilize simulations both to select parameters in the algorithm, and to show that our approach leads to a more accurate estimate of the residue function. We further show that our approach generates perfusion parameter maps which are less noisy. Lastly, we use in vivo rectal tumor CTp data to demonstrate that our estimated BF values lead to better separation between cancerous tissue—which by its angiogenic nature tends to have greater blood flow—and normal tissue. This paper is organized as follows. In Section II, we provide the theoretical background for perfusion parameter estimation, and introduce our spatio-temporal deconvolution method. In Section III, we define our experiments and results from using both simulations and clinical datasets. Finally, we provide a discussion of our conclusions in Section IV.

1183

The concentration of contrast agent (or tissue time-enhancement , within the given VOI can be written as function), (6) Under this model, at time , a unit of contrast agent is injected as a bolus, and indicating that the entire mass of contrast agent is within the vascular network. After some time when all contrast have left the vascular network, . Note that in practice, contrast agents can recirculate back into the vascular network; therefore this model is valid only during the initial passage of contrast, not during the subsequent recirculation or equilibrium phases. and are measured with Let us assume that equally spaced time points , with time increment . The convolution can be discretized (7) where

.. .

II. THEORY

.. .

and

A. Perfusion Parameters Estimation Using the theoretical treatment provided by [13], we consider in a bolus of nondiffusible contrast agent injected at time a tissue volume of interest (VOI). The individual particles of the contrast agent follow different paths through the vascular structure of the tissue. The distribution of their transit time is , and defined by the transport function (1) is the concentration of tracer at the venous output, is the concentration at the arterial input (i.e., the AIF), and denotes convolution. The corresponding mean transit time (MTT) is defined as where

(2) and the blood volume (BV), represented by the amount of intravascular contrast agent in the tissue of interest, can be calculated by (3)

.. .

.. .

..

.

.. .

Note that if can be estimated, then can be computed , . As such, much of the remainder of since at in order to compute this paper deals with the estimation of . B. Circulant Truncated Singular Value Decomposition (cTSVD) is factored into In singular value decomposition, matrix and and a diagonal matrix , two orthogonal matrices , in descending order with singular values, along the diagonal (8) Equation (7) can then be rewritten as (9)

The central volume theorem [29] then states that (4) The residue function , which measures the mass of contrast media remaining in a given vascular network over time, is defined by (5)

Since smaller singular values are associated with higher frequency singular vectors, the reciprocal of these small singular values become large weighting coefficients of oscillatory singular vectors. TSVD regularizes the solution by truncating small singular values to zero according to a certain threshold and thus remove the corresponding terms from the solution. It has been noted that perfusion parameter estimations are sensitive to delays and dispersions between the AIF and the

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

1184

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 5, MAY 2010

Fig. 1. Performance of cTSVD with different thresholds  at PCNR

= 20.

tissue VOI. For example, in major vasculopathy cases, contrast agents may arrive earlier in the tissue than in the chosen AIF, leading errors in the estimation. Following previous studies [15], [16], [30], we use circulant (by means of a block-circulant matrix) instead of linear deconvolution to avoid version of such causality problems. The accuracy of the estimation of residue function by cTSVD is sensitive to the choice of the regularization parameter . Other researchers have variously set to 20% of the maximum ele[13], to 4% of ment in at signal-to-noise ratio [15], or adaptively the maximum element in at determined based on SNR [31], although questions about the conditions under which the adaptive technique may be useful has been raised [20]. In this paper, we utilize simulations to select a that yields accurate BF values across a range of noise levels. Specifically, we simulate two regions with different BF values—one at 70 ml/100g/min for abnormal tissue, the other at 30 ml/100g/min for normal tissue. Different levels of Gaussian noise are added to generate various levels of peak contrast-tonoise ratios (PCNRs, which are calculated by dividing the peak contrast difference between the two regions by the Gaussian noise standard deviation). A number of PCNRs are examined between the range of 20–70, with highly consistent results. Figs. 1 and 2 shows the performance of cTSVD with PCNRs of 20 and 60, where the dashed lines indicate the true BF values. The error bars in Figs. 1 and 2 are calculated from 100 random Monte Carlo realizations. According to our simulations, when is approximate 0.06 (6% of the maximum element in ), the most accurate BF estimates are consistently achieved. C. Proposed Spatio-Temporal Deconvolution Spatial correlations rest on the assumption that physiological changes to perfusion are regional effects, rather than singlevoxel effects. However, since different regions may be undergoing different perfusion changes, it is important to break that spatial continuity when the data indicates edges may be present. As such, we have chosen to implement a continuity constraint in the spirit of Mumford–Shah [22]. In our model, we estimate perfusion parameters by considering both spatial and temporal correlations.

Fig. 2. Performance of cTSVD with different thresholds  at PCNR

= 60.

1) Spatio-Temporal Regularization: To incorporate spatial continuity, we modify the terms in (7) such that there is both spatial and temporal dependency. Specifically, let denote contrast agent concentration in the tissue of the voxel and time point . Similarly, let at spatial location be the remaining contrast agent concentration and time point . The of the voxel at spatial location least-squares form of (7) is (10) However, noise in the data is usually magnified to such extent that the solution may have no practical value or meaning. In the spirit of Mumford–Shah [22], we incorporate a prior of not only temporal but also spatial correlations through the inclusion of two constraints to the original least-squares cost function. This results in the new cost function

(11) where introduces a penalty if the solution behaves undesirably according to the priori information, (i.e., we seek nonoscillating and smooth solutions). Various models have been proposed (cf [17]) and a first order 4-D derivative operator suffices for our purpose. denotes an edge field identifying discontinuities in the perfusion estimates, which is defined separately for the , , and components and its value ranges from 0 (no edge between voxels) to 1 (strong edge between voxels). We use the , and , where notation . The regularization terms and control the regularization strength, which must be optimized and are discussed in the following section. and to The minimizer of (11) requires both be 0, then we have (12) and

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

(13)

HE et al.: A SPATIO-TEMPORAL DECONVOLUTION METHOD TO IMPROVE PERFUSION CT QUANTIFICATION

The minimizer of (11) requires that (12) and (13) are simultaneously satisfied. The minimization is obtained iteratively by first fixing to solve for , and then fixing to solve for . The edge field is initialized to zero, and at each voxel is initialized to be a vector of zeros. Each iteration is solved by conjugate gradient method. 2) Method for Choosing Regularization Parameter: To examine the trade-off between data fidelity and regularization, we look at the L-curve characterizing versus [32]. Note in particular that we focus first on opsince (12) involves only . Subtimizing the selection of sequently, according to the (13), the ratio of serves as a threshold of whether the smoothing effect should turn on or not. We utilize the simulation of two regions with different BF values described previously. When the gradient across the boundary , between these two regions is much less than the ratio of , which means the smoothing effect should be turned on; on the other hand, when the gradient is much greater than , , the smoothing effect should be turned off. Using this approach, is calculated based on simulation data.

1185

^ k as a function of ^ k versus kC 0 C R Fig. 3. Example of the L-curve (kLR  ) for the synthetic data at PSNR = 10. Values of the regularization parameter  are listed along the curve and the optimum  = 4 is chosen as the corner of the curve.

the following, peak signal-to-noise ratio (PSNR) is calculated by dividing the peak value of tissue time-enhancement signal by the noise standard deviation . B. Simulation Results

III. EXPERIMENTS AND RESULTS In this section, we first show our results from simulationbased experiments, followed by results from our clinical evaluations. A. Simulations 1) AIF Generation: The AIF is simulated using a gammavariate function [13], [33]; the analytic expression for is (14) where is bolus arrival time to any given region. Generally, , s are used to generate curves with much wider , since peaks for adults than for children the injected contrast agent dose is usually proportional to body weight and the amount of injected contrast agent is much larger , in adults than in children. Theses settings of and with can render an AIF with a shape and size that would typically be obtained using a standard injection scheme [34]. 2) Residue Function Generation: The transport function , where is a vector of model parameters, is simulated using the family of gamma distributions, as [18] (15) To ensure that the central volume theorem [(4)] is satisfied, we . The residue function can then be obset tained using (5). 3) Tissue Time-Enhancement Curve and Gaussian Noise , , , to generate Generation: We set , the AIF. Simulations are performed for in and BF values varied from 20 to increment. Correspondingly, MTT ranges to generate the from 3 to 12 s. We set residue function. Noise-free tissue time-enhancement curve can is be generated according to (7). Gaussian noise then added to the noise-free tissue time-enhancement signals. In

The L-curve used for optimizing the choice of is shown in Fig. 3. The vertical part of the curve corresponds to under-regularized estimations, where the solution is dominated by perturbation errors. On the other hand, the horizontal part of the curve corresponds to over-smoothed estimations, where the solution is dominated by residual fit errors. The optimum regularization is chosen as the corner between the horizontal parameter and vertical portions of the curve, which defines the transition between over- and under-regularization. It should be noted that is a relatively insensitive parameter, and that values near the “corner” of the L-curve yield very similar results. 1) Residue Function Recovery: The residue functions recovered by cTSVD and spatio-temporal deconvolution methods are and shown in Fig. 4. Under both circumstances of , the residue functions in Fig. 4(1c) and (2c) are in agreement with the true residual. On the other hand, the curves in Fig. 4(1b) and (2b) show relatively large baseline oscillations and less accurate peak values. 2) BF Estimation: From the recovered residue function, BF . We now examine the accuracy of our can be estimated at BF estimates at various noise levels. To do so, we generate a small region containing 40 40 voxels with the same perfusion characteristics and compute the mean and standard deviation over this region. The comparison is carried for different values of BFs and at both low and high noise levels. The simulations are repeated 25 times with random noise realizations. Shown in Fig. 5 are the BF values as estimated by cTSVD and spatio-temporal deconvolution, for two different noise levels. Spatio-temporal estimates can be visually observed to be more robust and accurate. To further summarize Fig. 5, we compute the mean squared error (MSE) as defined by (16) , are the true and estimated BF where and , values, respectively. The MSEs of cTSVD and spatio-temporal

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

1186

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 5, MAY 2010

Fig. 4. Recovery of residue functions by cTSVD and spatio-temporal methods. The parameters used for the simulation is BV = 4 ml=100 g . In the first row, PSNR = 20 and BF = 20 ml=100 g=min; in the second row, PSNR = 40 and BF = 80 ml=100 g=min. Column (a): True residue function. Column (b): Restored residue function by cTSVD. Column (c): Restored residue function by spatio-temporal deconvolution.

Fig. 5. Comparison of the accuracy in estimating BF by cTSVD and spatiotemporal deconvolution methods for different BF values. (a) PSNR = 10, MSE = 2:65 for spatio-temporal deconvolution and MSE = 134:52 for cTSVD; (b) PSNR = 40, MSE = 0:72 for spatio-temporal deconvolution and MSE = 13:1 for cTSVD.

deconvolution are 134.52 and 2.65, respectively, at , and 13.1 and 0.72 at . In addition, we are also interested in exploring how noise levels affect the performance in estimating BF of these two methods. For this purpose, we fix the BF and impose different levels of noise (PSNR varies from 5 to 60); the estimation results are shown in Fig. 6. It is shown that spatio-temporal deconvolution appears to provide more accurate estimates of BF than cTSVD across a broad range of noise levels. In addition, the accuracy of cTSVD in estimating degrades dramatically as a function of noise level, whereas the spatio-temporal deconvolution approach appears more stable. Furthermore, for the simulated uniform region, the ideal BF variability should be 0. Some variations are expected in the recovered BF maps and a better method should produce smaller variation. Examples of recovered BF maps are shown in Fig. 7. , cTSVD behaves When PSNR is extremely small poorly in recovering true object from background and defining

Fig. 6. Comparison of the accuracy in estimating BF via cTSVD and spatio-temporal deconvolution methods at different PSNRs. True BF = 30 ml=100 g=min.

edges, while in the BF map produced by spatio-temporal deconvolution method, the object can be easily distinguished from the background and the edge is well defined. When PSNR is , both methods can estimate BF values achigh curately, however, Fig. 7 (1a) contains more fluctuations than (1b). It is apparent that the estimated BF map produced by the spatio-temporal method is smoother than that by cTSVD. And further quantitative comparisons are shown in Fig. 8 (where BF is varied) and Fig. 9 (where PSNR is varied). In both figures, spatio-temporal deconvolution method produces lower BF variation than cTSVD does. 3) Contrast Preservation: Contrast is frequently considered a metric indicating how well different regions can be distinguished. For the problem at hand, one would expect that the contrast of the normal and abnormal tissue time-enhancement signals calculated using the recovered residue function should

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

HE et al.: A SPATIO-TEMPORAL DECONVOLUTION METHOD TO IMPROVE PERFUSION CT QUANTIFICATION

1187

Fig. 7. Examples of recovered BF maps via (a) cTSVD and (b) spatio-temporal deconvolution methods. The leftmost image is true BF map with BF = 60 ml=100 g=min. The PSNRs in first and second rows are 80 and 7, respectively.

Fig. 10. Simulated CTp data. First row and second row are 1-D time series (residue function) of a voxel which is spatially located in the regions I and II, respectively. The parameters used to generate residue function: a = 1, b = 3, c = 1:5 and t = 0 ; Voxels in region I have same perfusion characteristic, BF = 70 ml=100 g=min, and BV = 4 ml=100 g. Similarly, voxels in region II have same perfusion characteristic, BF = 30 ml=100 g=min and BV = 4 ml=100 g. Fig. 8. Comparison of cTSVD and spatio-temporal deconvolution methods in reducing variations over homogeneous region for different BF values. The ideal variation is 0. (a) PSNR = 10. (b) PSNR = 40.

Fig. 11. Comparison of cTSVD and spatio-temporal deconvolution methods in preserving contrast over different levels of noise.

Fig. 9. Comparison of cTSVD and spatio-temporal deconvolution methods in reducing variations over homogeneous region at different PSNRs. The ideal variation is 0. True BF = 30 ml=100 g=min.

comparable to that of the noise-free residue function. To compare the performance of cTSVD and spatio-temporal deconvolution methods in preserving contrast, we generate synthetic

CTp data spatially containing two 40 40 uniform regions with , different perfusion characteristic ( ), temporal resolution of 1 s, and total duration of 60 s, as shown in Fig. 10. A plot of true PCNRs against estimated PCNRs is shown in Fig. 11, from which PCNRs estimated by spatio-temporal deconvolution method can be observed to be consistently higher

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

1188

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 5, MAY 2010

scan parameters were used for the cine acquisition: 1 s gantry rotation time, 100–120 kVp, 200–240 mA, 10 s scanning delay from initiation of contrast media injection and the total duration of scanning was 45 s (four sections per gantry rotation of 5 mm thickness). Given the slice thickness, we processed our data in instead of four, since significantly three dimensions different anatomy can be observed in consecutive -slices. D. Clinical Experiments Results

Fig. 12. Comparison of (a) cTSVD and (b) spatio-temporal deconvolution methods in preserving edges between two adjacent regions over different levels of noise. PCNR = 2, 5, 10, and 20 at the first, second, third, and fourth row, respectively.

and more accurate than those by cTSVD, when . , cTSVD produces marginally higher PCNRs At than spatio-temporal deconvolution. However, from Fig. 12, it clearly shows that spatio-temporal deconvolution performs better than cTSVD in preserving edges bestween two adjacent regions consistently over different levels of noise, especially . when PCNRs are low C. Clinical Studies In this section we describe the results from comparing our approach with cTSVD on four subjects—two with rectal tumors and two normal rectums. The presence and precise location of tumors were identified by board-certified radiologists. 1) Data Acquisition: Perfusion CT was performed with a four slice MDCT scanner (LightSpeed QX/I, GE Medical Systems, Milwaukee, WI). 500 ml of saline was instilled rectally for adequate distension of the rectum and proper delineation of the tumor. An unenhanced scan of the pelvis (5 mm thickness) was performed first to localize the tumor. A 2 cm region of the tumor at its maximum dimension was selected for perfusion imaging. Dynamic scanning was then performed at a static table position in the region of the tumor following intravenous injection of 125 mL of iopamidol (Isovue 300, Bracco) at an injection rate of 7 cc/s through a 18 guage IV cannula. The following

For cTSVD, a threshold of 6% of the maximum singular value is used, in accordance with our simulation studies (Section II-B). Similarly, the regularization parameters for the spatio-temporal deconvolution method are chosen using the L-curve method discussed in Section II-C-2. Processing is performed on a computer with Intel Core 2 Quad CPU, [email protected] GHz and 2.39 GHz, 3 GB of RAM. It takes approximately 700 s to process a clinical data set of by spatio-temporal deconvolution method with 5 iterations, and approximately 50 s by the cTSVD method. 1) BF Estimation: BF maps were calculated for each data set using cTSVD and spatio-temporal deconvolution. We first compare two methods by visually observing the estimated BF maps. As shown in Figs. 13 and 14, variations in estimated blood flow maps are reduced greatly by our method. The tumor region are more evidently defined, while the normal region are more smooth. 2) Cancerous Voxels Clustering: By aggregating all voxels (within the ROI) from the two normal patient data sets into a “normal” cluster, and two tumor patient data sets into a “abnormal” cluster, we have two clusters of voxels—one consamples of BF values from normal rectal voxels, the taining other containing samples of BF values from cancerous rectal voxels. In our case, and . To quantify the separability between normal and tumor BF values, we define the distance between these two clusters as (17)

, are the means, and , are the standard deviwhere ations of BF in the normal and tumor clusters, respectively. We hypothesized that our spatio-temporal deconvolution algorithm to produce larger distance as defined in (17), that is, to more definitively differentiate between normal and tumor patients. Figs. 15 and 16 show scatter plots of normal versus abnormal cluster. The coordinate value of each point is “the number of pixels”—increasing pixel number moves from the top-left to the bottom-right of the region of interest as delineated by a radiologist. It is apparent that the two clusters are more separatable in data processed via spatio-temporal deconvolution than cTSVD . 3) Diagnostic Test: Let us define sensitivity as the proportion of samples with abnormal BF values which test positive, and specificity as the proportion of samples with normal BF values which test negative, at a specific threshold. Fig. 17 is the receiver operating characteristic (ROC) curve drawn based on 1247 normal samples and 2269 abnormal samples in which we examine a spectrum of thresholds. The plot shows the tradeoff between true positive rate (sensitivity) and false positive rate

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

HE et al.: A SPATIO-TEMPORAL DECONVOLUTION METHOD TO IMPROVE PERFUSION CT QUANTIFICATION

1189

Fig. 13. BF maps and zoomed-in regions of two rectal tumor patients estimated by (a) cTSVD and (b) spatio-temporal deconvolution. Tumor regions are delineated in red.

Fig. 14. BF maps and zoomed-in regions of two patients with normal rectums estimated by (a) cTSVD and (b) spatio-temporal deconvolution.

Fig. 15. Two clusters of normal vs. abnormal generated by cTSVD method. The distance d between two clusters is 31.37.

Fig. 16. Two clusters of normal versus abnormal generated by our spatial-temporal deconvolution method. The distance d between two clusters is 81.23.

(1-specificity). The closer the curve is to the upper left corner, the more accurate the test. Fig. 17 shows that spatio-temporal

deconvolution appears to be considerably more accurate than cTSVD, leading to more efficient diagnostics.

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

1190

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 5, MAY 2010

especially high noise levels. The proposed spatio-temporal deconvolution method improves the characterization of residue function by both producing a more accurate representation (which results in improved estimation of BF) and reducing the baseline oscillations. Spatio-temporal deconvolution also reduces unexpected variations within a local uniform region of interests; preserves the contrast of tissue-enhancement signals between different types (i.e., between normal and tumor) tissues, and appears to facilitate clinical diagnosis.

Fig. 17. ROC curves generated by spatial-temporal deconvolution and TSVD. Area under curve (AUC) of spatio-temporal deconvolution is 0.9800 and the AUC of cTSVD is 0.7987.

IV. DISCUSSION AND CONCLUSION In this paper, a spatio-temporal deconvolution method for improving characterization of residue function and for quantifying perfusion parameters was introduced. The proposed method was compared to the cTSVD, and numerical simulations with known AIF and residue function were conducted to validate the proposed methodology. For our simulation studies, a gamma-variate probability density function was used, since the corresponding residue functions presented a wide range of shapes, from exponential to box-car (with increasing ), representative of normal to ischemic tissues. Our experiments were conducted using different shapes of residue functions and results were highly consistent. The L-curve method was chosen for determining the regularin our work as noted in Section II-C-2. It ization parameter would be interesting to compare L-curve with other approaches, such as discrepancy principle, generalized cross-validation, and other statistical approaches [32]. It is also noted that the first-order derivative operator was used in the presented work, although other higher order derivative operators (e.g., , ) may also be effective [17]. It has been found that the presence of bolus delay and dispersion between the AIF and the tissue of interest can introduce considerable errors in the quantification of BF [34]. The problem can generally be minimized by either correcting the dispersion or using a local AIF [35]. The proposed spatio-temporal deconvolution method made an effort to avoid such problem via the usage of a block-circulant version of , as we did with the cTSVD method as well. In summary, by exploiting contextual information from neighboring voxels both spatially and temporally, a spatio-temporal regularizer has been introduced. It has been demonstrated by both simulated and clinical data that the proposed method appears to be more robust than cTSVD with respect to noise,

REFERENCES [1] K. Miles and M. Griffiths, “Perfusion CT: A worthwhile enhancement?,” Br. J. Radiol., vol. 76, pp. 220–231, 2003. [2] K. Miles, “Perfusion CT for the assessment of tumour vascularity: Which protocol?,” Br. J. Radiol., vol. 76, pp. S36–S42, 2003. [3] M. Konig, “Brain perfusion ct in acute stroke: Current status,” Eur. J. Radiol., vol. 45, pp. S11–S22, 2003. [4] E. Hoeffner, I. Case, R. Jain, S. Gujar, G. Shah, J. Deveikis, R. Carlos, B. Thompson, M. Harrigan, and S. Mukherji, “Cerebral perfusion CT: Technique and clinical applications,” Radiology, vol. 231, pp. 632–644, 2004. [5] D. Torigian, S. Huang, M. Houseni, and A. Alavi, “Functional imaging of cancer with emphasis on molecular techniques,” CA: Cancer J. Clinicians, vol. 57, pp. 206–224, 2007. [6] M. Harrigan, J. Leonardo, K. Gibbons, L. Guterman, and L. Hopkins, “CT perfusion cerebral blood flow imaging in neurological critical care,” Neurocritical Care, vol. 2, pp. 352–366, 2005. [7] K. Miles, C. Charnsangavej, F. Lee, E. Fishman, K. Horton, and T. Lee, “Application of CT in the investigation of angiogenesis in oncology,” Acad. Radiol., vol. 7, pp. 840–850, 2000. [8] A. Cenic, D. Nabavi, R. Craen, A. Gelb, and T. Lee, “Dynamic CT measurement of cerebral blood flow: A validation study,” Am. J. Neuroradiol., vol. 20, pp. 63–73, 1999. [9] A. Cenic, D. Nabavi, R. Craen, A. Gelb, and T. Lee, “A CT method to measure hemodynamics in brain tumors: Validation and application of cerebral blood flow maps,” Am. J. Neuroradiol., vol. 21, pp. 462–470, 2000. [10] M. Wintermark, J.-P. Thiran, P. Maeder, P. Schnyder, and R. Meuli, “Simultaneous measurement of regional cerebral blood flow by perfusion CT and stable xenon CT: A validation study,” Am. J. Neuroradiol., vol. 22, pp. 905–914, 2001. [11] J. Eastwood, M. Lev, T. Azhari, T. Lee, D. Barboriak, D. Delong, C. Fitzek, M. Herzau, M. Wintermark, R. Meuli, D. Brazier, and J. Provenzale, “CT perfusion scanning with deconvolution analysis: Pilot study in patients with acute middle cerebral artery stroke,” Neuroradiology, vol. 222, pp. 227–236, 2002. [12] L. Ostergaard, “Cerebral perfusion imaging by bolus tracking,” Topics Magn. Reson. Imag., vol. 15, pp. 3–9, 2004. [13] L. Ostergaard, R. Weisskoff, D. Chesler, C. Gyldensted, and B. Rosen, “High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis,” Magn. Reson. Med., vol. 36, pp. 715–725, 1996. [14] L. Otergaard, A. Sorensen, K. Kwong, R. Weisskoff, C. Gyldensted, and B. Rosen, “High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. part II: Experimental comparison and preliminary results,” Magn. Reson. Med., vol. 36, pp. 726–736, 1996. [15] O. Wu, L. Ostergaard, R. Weisskoff, T. Benner, B. Rosen, and A. Sorensen, “Tracer arrival timing-insensitive technique for estimating flow in MR perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix,” Magn. Reson. Med., vol. 50, pp. 164–174, 2003. [16] H.-J. Wittsack, A. Wohlschlager, E. Ritzl, R. Kleiser, M. Cohnen, R. Seitz, and U. Modder, “CT-perfusion imaging of the human brain: Advanced deconvolution analysis using circulant singular value decomposition,” Computerized Med. Imag. Graphics, vol. 32, pp. 67–77, 2008. [17] F. Calamante, D. Gadian, and A. Connelly, “Quantification of bolustracking MRI: Improved characterization of the tissue residue function using tikhonov regularization,” Magn. Reson. Med., vol. 50, pp. 1237–1247, 2003. [18] K. Mouridsen, K. Friston, N. Hjort, L. Gyldensted, L. Ostergaard, and S. Kiebel, “Bayesian estimation of cerebral perfusion using a physiological model of microvasculature,” NeuroImage, vol. 33, pp. 570–579, 2006.

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.

HE et al.: A SPATIO-TEMPORAL DECONVOLUTION METHOD TO IMPROVE PERFUSION CT QUANTIFICATION

[19] N. Pack, E. DiBella, T. Rust, D. Kadrmas, C. McGann, R. Butterfield, P. Christian, and J. Hoffman, “Estimating myocardial perfusion from dynamic contrast-enhanced CMR with a model-independent deconvolution method,” J. Cardiovas. Magn. Reson., vol. 10, p. 52, 2008. [20] I. Andersen, A. Szymkowiak, C. Rasmussen, L. Hanson, J. Marstrand, H. Larsson, and L. Hansen, “Perfusion quantification using Gaussian process deconvolution,” Magn. Reson. Med., vol. 48, pp. 351–361, 2002. [21] K. Wong, C. Tam, M. Ng, S. Wong, and G. Young, “Improved residue function and reduced flow dependence in MR perfusion using least-absolute-deviation regularization,” Magn. Reson. Med., vol. 61, pp. 418–428, 2009. [22] D. Mumford and J. Shah, “Boundary detection by minimizing functionals,” in IEEE Conf. Comput. Vis. Pattern Recognit., 1985, pp. 22–26. [23] B. Orten, W. Karl, D. Sahani, and H. Pien, “A new deconvolution approach to perfusion imaging exploiting spatial correlation,” Proc. SPIE Med. Imag., 2008. [24] L. He, B. Orten, S. Do, W. Karl, A. Kambadakone, D. Sahani, and H. Pien, “Spatio-temporal deconvolution of perfusion CT data in rectal tumor patients,” in IEEE Int. Symp. Biomed. Imag., 2009. [25] X. Descombes, F. Kruggel, and D. V. Cramon, “Spatio-temporal fMRI analysis using Markov random fields,” IEEE Tran. Med. Imag., vol. 17, no. 6, pp. 1028–1039, Dec. 1998. [26] K. Katanoda, Y. Matsuda, and M. Sugishita, “A spatio-temporal regression model for the analysis of functional MRI data,” NeuroImage, vol. 17, pp. 1415–1428, 2002. [27] L. He and I. Greenshields, “An MRF spatial fuzzy clustering method for fMRI SPMs,” Biomed. Signal Process. Control, vol. 3, pp. 327–333, 2008.

1191

[28] V. Schmid, B. Whitcher, A. Padhani, N. Taylor, and G. Yang, “Bayesian methods for pharmacokinetic models in dynamic contrast-enhanced magnetic resonance imaging,” IEEE Trans. Med. Imag., vol. 25, no. 12, pp. 1627–1636, Dec. 2006. [29] P. Meier and K. Zierler, “On the theory of the indicator-dilution method for measurement of blood flow and volume,” J. Appl. Physiol., vol. 6, pp. 731–744, 1954. [30] O. Wu, L. Ostergaard, and A. Sorensen, “Technical aspects of perfusion-weighted imaging,” Neuroimag. Clin. N. Am., vol. 15, pp. 623–637, 2005. [31] H.-L. Liu, Y. Pu, Y. Liu, L. Nickerson, T. Andrews, P. Fox, and J.-H. Gao, “Cerebral blood flow measurement by dynamic contrast MRI using singular value decomposition with an adaptive threshold,” Magn. Reson. Med., vol. 42, pp. 167–172, 1999. [32] W. Karl, “Regularization in image restoration and reconstruction,” in Handbook of Image and Video Processing, A. Bovik, Ed. New York: Elsevier, 2005, pp. 183–202. [33] M. Rausch, K. Scheffler, M. Rudin, and E. Radu, “Analysis of input functions from different arterial branches with gamma variate functions and cluster analysis for quantitative blood volume measurements,” Magn. Reson. Imag., vol. 18, pp. 1235–1243, 2000. [34] F. Calamante, D. Gadian, and A. Connelly, “Delay and dispersion effects in dynamic susceptibility contrast MRI: Simulations using singular value decomposition,” Magn. Reson. Med., vol. 44, pp. 466–473, 2000. [35] F. Calamante, “Bolus dispersion issues related to the quantification of perfusion MRI data,” J. Magn. Reson. Imag., vol. 22, pp. 718–722, 2005.

Authorized licensed use limited to: BOSTON UNIVERSITY. Downloaded on July 12,2010 at 14:58:44 UTC from IEEE Xplore. Restrictions apply.