Abstract

Deconvolution-based analysis of CT and MR brain perfusion data is widely used in clinical practice and it is still a topic of ongoing research activities. In this paper, we present a comprehensive derivation and explanation of the underlying physiological model for intravascular tracer systems. We also discuss practical details that are needed to properly implement algorithms for perfusion analysis. Our description of the practical computer implementation is focused on the most frequently employed algebraic deconvolution methods based on the singular value decomposition. In particular, we further discuss the need for regularization in order to obtain physiologically reasonable results. We include an overview of relevant preprocessing steps and provide numerous references to the literature. We cover both CT and MR brain perfusion imaging in this paper because they share many common aspects. The combination of both the theoretical as well as the practical aspects of perfusion analysis explicitly emphasizes the simplifications to the underlying physiological model that are necessary in order to apply it to measured data acquired with current CT and MR scanners.

1. Introduction

Tissue perfusion measurement from iodinated contrast agent enhancement on CT scans was first proposed by Axel in 1980 [1]; this was based on earlier developments by Meier and Zierler [2] for measuring blood flow and blood volume. At that time, the CT-based measurements were strictly limited to research because of the low speeds and narrow coverage of the existing CT scanners. However, the introduction of perfusion CT (PCT) helped expand the utility of CT significantly since it could now provide capillary level hemodynamic information. Within about a decade, perfusion imaging techniques were also adopted in MR [35].

With the advent of helical scanners and faster rotating gantries (0.33–0.5 s/rotation) in conjunction with multidetector geometries which provide larger coverage, PCT has now become part of the routine screening for many diseases.

Given the existing developments in perfusion imaging, the purpose of this paper is to focus on a detailed derivation of the theoretical model for deconvolution-based perfusion measurement. While the main equation of this model is well known, its derivation is spread over several publications.

We therefore first present a summary of the derivation, with the aim of fully explaining the parameters and the underlying assumptions that are made. Based on the main equation of the theoretical model, we also present a guideline for the algorithmic implementation of the deconvolution-based perfusion measurement. We discuss robust numerical deconvolution and discuss topics related to data pre-processing, providing references to the literature for each of the special topics. The overall aim of this paper is to provide an understanding of the underlying assumptions of the theoretical model and to show how the (simplified) model can be robustly implemented for clinical image analysis.

2. Clinical Applications of Perfusion Imaging

Perfusion imaging is most widely used in acute stroke and oncology [6]. When used in diagnosis of stroke, the purpose of perfusion imaging is to identify the extent of affected tissue and to delineate the ischemic tissue that can be reperfused. In oncology, perfusion imaging helps to identify angiogenetic tumors that alter the local tissue perfusion due to generation of neovasculature. Perfusion measurements are increasingly being used for assessment, staging, and monitoring posttherapy [6, 7].

Figure 1 shows common parameter maps based on a brain perfusion CT exam (Somatom Definition AS+, Siemens AG, Healthcare Sector, Forchheim, Germany) of a 69-year-old male stroke patient. The patient presented to the hospital with an acute high-grade hemiparesis on the right side. A CT angiography scan indicated an occlusion of the left middle cerebral artery. The time-to-peak (TTP) image shows a large lesion that illustrates the maximum affected tissue. In addition, the cerebral blood flow (CBF), cerebral blood volume (CBV), and mean transit time (MTT) images exhibit perfusion deficits in a smaller brain territory. In general, these perfusion CT maps are interpreted appropriately in order to guide the recanalization procedure of the occluded vessel.

Blood flow is critical to the functionality of any organ since it provides the essential nutrients and oxygen. In case of flow disruption, the body autoregulates the flow and pressure either by altering blood flow or volume or both. In the brain, there are some fairly well-defined thresholds for the cerebral blood flow in normal, reversibly damaged, and necrotic tissue. The normal value for the cerebral blood flow is between 50 and 60 mL/100 g/min for grey matter [8]. The average value decreases with age and is about 2 to 3 times lower in white matter compared to grey matter [9]. Any reduction in normal perfusion pressure results in vasodilation and hence an increase in blood volume and transit times. As the perfusion pressure falls lower, compensatory vasodilation is unable to offset the deficit. When this value falls below 20 mL/100 g/min (for grey matter), synaptic transmission ceases to function. When the flow is below 8–10 mL/100 g/min, the cell membrane pumps fail, causing irreversible damage to the cells [6].

2.1. Perfusion CT

In the acute stroke setting, conventional CT has been the imaging modality of choice for ruling out intracerebral hemorrhages (ICH). However, overall the sensitivity of CT for stroke detection is 60–65% [10, 11]. For ischemic stroke, which represents about 85% of all stroke cases, the inclusion of PCT along with CT angiography (CTA) can identify the subtle abnormalities in the cerebral tissue that can be missed on the noncontrast agent-enhanced scans. Most commonly, the perfusion scan consists of imaging one or two slices at the level of the basal ganglia. This allows inclusion of the branches of the carotid artery that are typically thrombosed. After approximately 7–10 s following an intravenous injection of iodinated contrast agent, continuous scanning is performed for about 50 s. Table toggling techniques are sometimes used to increase the coverage. More recent wide detector scanners allow whole brain coverage in each scan. The temporal scans are reconstructed and one of several approaches can be used to calculate the perfusion parameters.

In animal studies, the product of CT cerebral blood volume (CBV) and flow (CBF) from CT measurements was found to have sensitivity of 90.6% and specificity of 93.3% (compared with histological measurements) for discerning ischemic and oligemic tissue [12]. One study that compared stroke diagnosis using CT perfusion plus angiography, against MRI, found good correlation and no significant prognosis differences [13]. Typical PCT scans add approximately 5 or less additional minutes to the scan time with around a 50 mL bolus of additional iodinated contrast agent. With regards to the X-ray dose in PCT, depending on the parameters, the effective dose is estimated to be between 1.2 mSv [6] and 3.4 mSv [14]. This is in the same range as the effective dose of a standard cerebral CT, which is reported to deliver about 2.5 mSv to the patient [14].

The gold standard for perfusion CT has been imaging with stable xenon as the contrast agent [7, 15]. This method involves inhalation of a mixture of stable xenon gas and oxygen followed by CT scanning. Because of the high atomic number of xenon, it serves as a radio-opaque contrast agent as it diffuses into the blood and neurons in a well-balanced manner. It has been proven to be accurate in quantifying perfusion by comparing the results with those obtained using radio-labeled microspheres.

2.2. Perfusion MR

Perfusion imaging in MR can be performed with or without contrast agent [16]. Noncontrast agent-enhanced perfusion imaging usually uses spin labeling of blood entering the imaging volume. This method is less commonly used because of the increased sensitivity to motion and related artifacts and low signal in case of slow flow. Gadolinium-based tracers such as Gd-DTPA are more commonly used for measuring perfusion derived from changes in the local susceptibility. Both spin echo (SE) and gradient echo (GRE) sequences have been applied successfully in perfusion MR. GRE sequences are most frequently used because they provide a better contrast-to-noise ratio for imaging of the contrast agent compared to SE sequences [1719]. However, GRE sequences have the disadvantage of disproportionately weighting the contribution of the contrast agent in relatively large vessels, whereas SE sequences provide a more accurate assessment of blood flow through vessels of all sizes [19]. After the 7–10 s interval that the gadolinium contrast agent takes to reach the brain following the intravenous injection, the signal in the cerebral tissue dips. The signal changes are most significant over about 15 s during which the change in T2* or equivalently the change in the associated relaxation rate R2* is monitored. Note that this also requires that the contrast agent is intravascular. Rapid imaging (interval less than 2 s) is required for accurate measurement of perfusion parameters. Typically echo-planar imaging (EPI) sequences are used for this purpose.

3. Theoretical Model

The aim of this section is to provide a compact outline of both some elementary as well as practically relevant theory of perfusion estimation based on previous work. In particular, we will introduce a theoretical physiological model of tissue perfusion for intravascular tracer systems and present the derivation of a deconvolution-based mathematical approach for the estimation of diagnostically important perfusion parameters. In addition, we will briefly describe alternative methods that do not require deconvolution.

3.1. Model of Microcirculation at the Tissue Level

For computing the tissue perfusion, we assume a physiological model of the blood supply to the tissue. Figure 2 shows this model that consists of a volume of interest covering the organ-specific parenchyma, the interstitial space, and the capillary bed. The volumes of the parenchyma and the interstitial space are denoted by , while the volume of the capillary bed is referred to as . The entire volume of interest shall be supplied with blood by a single arterial inlet and correspondingly drained by a single venous outlet. In general, it may have a different shape than the cuboid shown in Figure 2. A blood cell can take various paths through the capillary bed. The transit time it needs to pass through the capillary bed depends on the chosen path. We assume a stationary probability density distribution of transit times.

Once a contrast agent bolus has been injected, it enters the volume under consideration via the arterial inlet and is then diluted into the capillary bed. The local contrast agent concentrations and are measured directly adjacent to the capillary bed on the arterial and venous sides, respectively. Furthermore, the average contrast agent concentration within the volume of interest can also be measured. In perfusion CT, an iodinated contrast agent is used whereas, in perfusion MR, the measured signal difference is created by a paramagnetic contrast agent based on gadolinium (Gd) (see Section 2.2). The contrast agent concentration is defined as mass of iodinated contrast agent per volume (unit: g/mL) or amount of Gd-based contrast agent per volume (unit: mol/mL), respectively [20]. For the following analysis, we assume the contrast agent concentration to be measured as mass per volume, which can easily be related to amount per volume.

Figure 3 illustrates typical time-concentration curves , , and that may be measured in brain tissue, for example. For the sake of simplicity, the maximum contrast agent concentration has been normalized to 1. Note that the (average) enhancement within the volume of interest is commonly more than an order of magnitude below the enhancements of the feeding artery and the draining vein.

An additional important assumption is that the contrast agent remains in the intravascular space. For our case of cerebral perfusion, it should therefore not cross the blood-brain barrier (BBB). As a consequence, this means that all contrast agents entering from the arterial inlet will eventually leave the volume of interest at the venous outlet. A breakdown of the BBB may occur in tumor patients, in stroke patients, and in patients that suffer from inflammations or infections, for example. In these cases, the methods presented in this paper may lead to inaccurate perfusion estimates and particularly to an overestimation of the blood volume [21, 22]. Note that there exist other modelling approaches which do not assume that the contrast agent remains in the intravascular space. These models can be used for measuring tumor perfusion, for example [6, 20, 23].

Finally, we suppose that the contrast agent mixes perfectly with the blood and that the physical properties of the blood (its flow behavior, in particular) are not influenced by the contrast agent.

As we will see, only knowledge of the functions and is needed to compute the blood flow within the volume under consideration. In practice, the function —also known as the arterial input function (AIF)—is not measured directly at the respective volume of interest, but in a larger feeding artery in order to achieve a reasonable signal-to-noise ratio (SNR) (see Section 4.1).

As a first diagnostically relevant perfusion parameter, the mean transit time (MTT) of the volume under consideration is defined as the first moment of the probability density function of the transit times, that is,

Furthermore, the residue (or residual) function —compare [24]—represents an intermediate quantity of interest and is defined as

The (dimensionless) residue function thus quantifies the relative amount of contrast agent that is still inside the volume of interest at time after an (idealized) delta-shaped contrast agent bolus has entered the volume at the arterial inlet at time ; that is, . Due to the various transit times within the capillary bed, the contrast agent will not leave the volume instantaneously, but gradually over time. In particular, this means that the residue function decreases continuously from to 0. Figure 4 shows typical examples of a distribution function of transit times as well as the corresponding residue function . In this example, the function is modeled by a gamma distribution [25].

3.2. Derivation of the Indicator-Dilution Theory

Using the parameters defined in Table 1, the accumulated masses of contrast agent that have entered and left the volume of interest during the time interval , denoted as and , respectively, can be expressed as The volume flow is assumed to be constant over time. The contrast agent concentrations and at the arterial inlet and the venous outlet, respectively, are time-dependent functions which we assume to be 0 for . These functions primarily depend on the parameters of the contrast agent injection and the patient's cardiac cycle.

We can compute the mass of a contrast agent within the volume of interest at time using the principle of conservation of mass as

The contrast agent concentration at the venous outlet can be computed from the contrast agent concentration at the arterial inlet by convolving it with the probability density function (). We therefore obtain

Note that throughout this paper, all integrals with infinite integration endpoints shall be interpreted as the limit of the integral when the respective endpoint approaches . Using (5), we can rewrite (4), by applying the delta function , as Changing the order of integration and rearranging this equation leads to By applying the substitution , recalling that, for , we have and considering that for , we obtain Equation (7) thus eventually reads We introduce the cerebral blood flow (CBF) as the blood volume flow normalized by the mass of the volume , Inserting this definition into (10) yields According to Table 1, we define the contrast agent concentration within the volume of interest as which finally leads to the following formulation of the indicator-dilution theory, where denotes the convolution operator as usual, see also [21, 26]. An alternative derivation of the same mathematical result is presented in [20]. A historical overview of the development of the indicator-dilution theory with numerous references to mathematical aspects can be found in [27]. Note that the solution of (14) with respect to CBF and other clinically important perfusion parameters will be discussed in Section 3.3.

From a physiological point of view, it would be more meaningful to normalize CBF by the mass of the volume . This volume contains the mass of the parenchyma (and the interstitium) only. In that case, CBF would be a local measure for the blood volume flow per mass of parenchyma (and interstitium) that actually requires blood supply for oxygen and nutrient delivery. In (11), however, the volume also contains the mass of the blood-filled capillary bed itself. Another aspect to consider is that the mean density of the volume, which influences the CBF value, actually depends on the (varying) mass of the contrast agent in the capillary bed. The alternative definition of CBF, would then lead to a corresponding alternative formulation of the indicator-dilution theory, From a practical perspective, however, it is more convenient to use the definition of CBF given by (11), see Section 4.1.

The derivation of the indicator-dilution theory in this section was focused on brain perfusion imaging. This theoretical model can be used in stroke patients if the BBB is intact—compare Section 3.1—but it is not suited for semi-permeable tumors, for example. With slight adaptations, this theoretical model can also be applied in other applications of perfusion imaging such as pulmonary perfusion imaging. See [28] for detailed discussions. A discussion of models in hepatic and renal perfusion imaging is given in [29, 30], respectively.

In the context of perfusion measurement, the term recirculation refers to the physiological phenomenon that, due to the patient's cardiac activity, the contrast agent passes through the volume under consideration multiple times. It can easily be shown, however, that there is no need to correct for recirculation when deconvolution methods are applied to determine perfusion parameters [31].

3.3. Computation of Perfusion Parameters Using Deconvolution

In (14), the variables and can be measured and have known values whereas the values of CBF, , and are unknown. In order to compute CBF as well as other diagnostically relevant tissue perfusion parameters, we first need to introduce an intermediate variable, the flow-scaled residue function , which is given in units of 1/s and can be determined directly from the measured data and . Using (17), (14) can be written as Hence, can be obtained from the measured data and using a deconvolution method. Since a fundamental property of the residue function is , we may then determine CBF as Using instead of has particular practical advantages that will be discussed in detail in Section 4.1.

The flow-scaled residue function can further be used to determine the MTT parameter of the tissue volume under consideration. From (2), it follows that, for , we have Equation (1) can thus be rewritten, and then using integration by parts and (17) and (19), we obtain

Note that we have assumed that there is a constant such that for . This assumption ensures that

The cerebral blood volume (CBV) corresponding to the tissue volume represents another diagnostically relevant perfusion parameter and is defined as It quantifies the blood volume normalized by the mass of and is typically measured in units of mL/100 g. The quantity CBV can be computed from the parameters CBF and MTT using the central volume theorem [22, 26], according to which holds for the perfused volume of interest. Interestingly, this theorem has been recognized for a long time and is already found in a historical publication from 1893 [32]. It states that the perfusion parameters CBV and CBF corresponding to the volume of interest are related by the respective temporal parameter MTT that quantifies the mean time that a blood cell needs to pass through its capillary bed. With (19) and (21), it follows from (24) that which demonstrates that the CBV parameter can be derived from the flow-scaled residue function as well.

A healthy human brain exhibits a CBV of about 4 mL/100 g for grey matter and a CBV of about 2 mL/100 g for white matter [8].

Note that the definition of CBV that corresponds to the alternative definition of CBF in [16] is Accordingly, this alternative definition relates the blood volume to the mass of the parenchyma (and the interstitium) only and explicitly omits the mass of the capillary bed itself.

Furthermore, there are references in the literature that suggest measuring the blood volume in units of mL/mL. This alternative dimensionless quantity may therefore be considered as a measure of blood (or vascular) volume fraction. When relating the volume of the capillary bed to the entire volume of interest, a typical average ratio of about 4% will result for the human brain. We refer to [33] for both technical and clinical details.

3.4. Overview of Nondeconvolution-Based Methods for Perfusion Imaging

For the sake of completeness, this section will briefly cover two alternative approaches for CBV and CBF estimation that are practical and relevant, and that do not involve deconvolution operations. Nondeconvolution-based methods for estimating perfusion parameters are also referred to as direct measurement-based approaches [26].

Firstly, there is an alternative method to compute the blood volume of the tissue volume under consideration [1]. This approach assumes that the average contrast agent concentration in the tissue volume can be related to the average contrast agent concentration in the capillary bed by According to the principle of conservation of mass, it follows that where is the total mass of contrast agent that has passed through the volume of interest. This results in an alternative expression for CBV, Hence, assuming a suitable correction for contrast agent recirculation [1, 34], CBV can be estimated from the integrals of either and or and over time. See [21] for details and further references with a particular focus on MR perfusion measurements.

It is argued in [34] that, particularly for the case of CT perfusion imaging of the brain, a physiologically reasonable approximation to (29) is given by which avoids the computation of the integrals over time and only requires the maximum values of and .

Secondly, there is a nondeconvolution-based approach to estimate the blood flow of the tissue volume under consideration; the maximum slope method [22, 34]. The derivation of this method is based on (4) and further assumes for simplicity's sake that there is no venous outflow from the tissue volume under consideration during the time of observation; that is, Recalling the CBF definition—compare (11)—and that , we obtain Taking the derivative of (32) yields and since (33) must hold for all , the blood flow is given by which means that CBF can be estimated by dividing the maximum slope of the tissue time-concentration curve , shown as an example in Figure 5, by the maximum value of the contrast agent concentration in the feeding artery.

An advantage of the maximum slope method is the shorter overall acquisition time. As a downside, however, it requires a faster contrast agent bolus injection rate in order to approximately fulfill the no-venous-outflow condition.

A more comprehensive discussion of the maximum slope method and a comparison with the deconvolution method is presented in [35]. According to [35], the clinical results based on these two approaches are generally of comparably high quality in CT imaging applications. However, in cases with insufficient data quality (e.g., in terms of noise, contrast agent concentration, bolus shape), deconvolution-based methods may lead to superior results. Moreover, violation of the aforementioned no-venous-outflow condition may yield incorrect perfusion estimates when the maximum slope method is employed. This can happen for penumbral regions of the brain which characterize the tissue at risk after an ischemic stroke.

3.5. Additional Perfusion Parameters

Besides the aforementioned quantities CBV, CBF, and MTT, there are additional perfusion parameters such as the time-to-peak (TTP) of the time-concentration curve, the maximum contrast agent concentration , as well as the first moment (FM) of the time-concentration curve, for example. The first moment can be computed by projecting the centroid of the area under the curve (AUC) of the time-concentration curve onto the time axis.

Figure 5 illustrates the quantities , TTP, and FM. The remaining parameter bolus arrival time (BAT) will be explained in Section 4.1. In practical measurements, the time point represents the start of the scanning. A detailed description and analysis of these additional quantities, however, is beyond the scope of this paper. A comparison of several perfusion parameters and their clinical impact on the treatment of stroke patients is given in [36].

In summary, Table 2 covers the most common diagnostically relevant perfusion parameters and shows how they can be determined employing deconvolution-based and nondeconvolution-based methods. In principle, the central volume theorem—compare (24)—may also be used to numerically estimate the MTT from the parameters CBF and CBV when the latter have been computed using nondeconvolution-based algorithms. However, the authors are not aware of any reference that describes the application of this approach in clinical practice.

4. Practical Implementation

This section is devoted to the practical computer implementation of algorithms for perfusion image analysis. First, we will discuss the necessary adaptations of the theoretical model from Section 3 that are needed for its application to data from real CT and MR scanners. Afterwards, we will describe commonly used algebraic deconvolution methods and also give an overview of alternative approaches. We will motivate the need for suitable regularization and discuss the influence of the regularization parameter on the resulting perfusion estimates. For the sake of completeness, we will also address techniques for the pre-processing of the acquired perfusion data.

4.1. Adaptations of the Model of Microcirculation

In Section 3.1, we presented a model of microcirculation at the tissue level. We have assumed that we can measure the average contrast agent concentration corresponding to a volume under consideration which is supplied by one single capillary bed only. Furthermore, we have supposed that we can measure the contrast agent concentration locally at the arterial inlet into the capillary bed. However, real CT and MR scanners are characterized by limited spatial (and contrast) resolution and, in reality, one cannot rely on these two aforementioned assumptions. We will thus introduce two major adaptations of the physiological model which are necessary once it is to be applied to data from real scanners.

First, during a standard CT and MR perfusion exam, a volume of interest is scanned and the data is reconstructed on a grid of regularly spaced voxels. In the object domain, each voxel volume () contains numerous capillary beds as well as arterioles and venules that supply and drain these capillary beds, respectively. For the particular case when the volume is located completely within a larger artery or vein, there are of course no capillary beds located within .

The measured signal (X-ray attenuation or MR relaxation rate) in a voxel is thus a combination of the signals from both the capillary beds as well as the arterial and venous vessels [37]. The perfusion parameters that are computed from the voxel's time-concentration curve are therefore not true parameters of the capillary perfusion. If no larger artery or vein is located inside the volume , we may adapt the model introduced in Section 3.1 as follows: the measured time-concentration curve refers to the average perfusion from the arterioles through the capillary beds to the venules found in .

The second adaptation of the model concerns the measurement of . In reality, it is not possible to locally measure the concentration at the arterial inlet into the volume . Instead, it is common practice that a global arterial input function (AIF) is chosen in a large arterial vessel. In brain perfusion imaging, for example, the anterior cerebral artery is often selected [38].

This approach leads to a traveling time of the contrast agent bolus from where the AIF is measured to the location of the tissue volume where is measured. We will refer to this traveling time as the bolus delay. Another physical effect that needs to be taken into consideration is bolus dispersion [39]. It appears as a widening of the shape of the bolus that is caused during the flow from the remote AIF location to the measurement site of .

The bolus delay has two implications. First, the curve does not start to rise at the same time point as starts to rise. The difference between these two time points can be defined as the bolus arrival time (BAT), which may be considered as an additional perfusion parameter [40]. Alternatively, the BAT can be defined as the time interval between the start of the scanning and the time when begins to rise, see Figure 5. The results obtained with this alternative definition differ from the results obtained with the first definition by a constant value only.

Second, the flow-scaled residue function is equal to 0 from to . In addition, due to the bolus dispersion, will not rise instantaneously to its maximum at , but it will have a finite rise time. The time-to-maximum (TMAX) of the flow-scaled residue function, defined as has also been suggested as an additional perfusion parameter [41, 42]. Since the function can be 0 at (due to bolus delay), it is reasonable and recommended to estimate CBF as the maximum of —compare (19)—and not as the value of at time .

Bolus delay and dispersion may lead to an underestimation of CBF [39]. In order to correct for bolus delay and dispersion several methods have been proposed [43, 44]. The use of local arterial input functions could also reduce the effect of bolus dispersion, see Section 4.5.6. On the other hand, new perfusion parameters (BAT, TMAX) are motivated by these two effects and can be defined accordingly. They represent perfusion characteristics related to the flow of the contrast agent bolus from the selected feeding artery to the respective tissue site, see again Figure 5.

4.2. Deconvolution Using Algebraic Methods

In this section, we will discuss the robust numerical solution of the main equation of the indicator-dilution theory—(18)—by means of algebraic deconvolution methods. An overview of further deconvolution methods will then be given in Section 4.3. We will introduce the discretization of (18) and show that its solution without regularization leads to nonphysiological results. We will explain and motivate suitable regularization approaches by a singular value decomposition-based analysis. To illustrate the mathematical concepts, we will provide examples using the time-attenuation curves and shown in Figure 6 that were extracted from a real perfusion CT scan.

We assume that the measured time-attenuation curves can be converted to time-concentration curves using a constant of proportionality of 1 g/mL/HU. Details about the conversion, also discussing perfusion MR data, will be explained in Section 4.5.4.

In practice, the time-concentration curves and are sampled at discrete time points. We denote these time points as with . A typical value of the sampling period is 1 s, for example. We can discretize (18) as see [45]. We assume that the values of can be neglected for . Since for , the end summation index could also be set to instead of . By rewriting this expression using matrix-vector notation, we obtain or shortly where and are contained in the matrix , and and represent the entries of the vectors and , respectively. Different ways to discretize (18) are investigated in [46]. For example, it was suggested in [47, 48] to use a discretization method with a block-circulant matrix in order to reduce the influence of the bolus delay. See the appendix for details.

A standard approach to solve (37) for is to use the singular value decomposition (SVD) of . For a matrix with linearly independent rows and columns, it is defined as where and are unique orthogonal matrices composed of the left and right singular vectors and , respectively [49]. The number of rows and columns in that only contain zeros is determined by the number of leading zeros in the series . Therefore, has rank . After the subtraction of the baseline, it may happen that the first entry is zero, see Section 4.5.4, and that thus becomes rank-deficient. The diagonal matrix contains the singular values in nonincreasing order . The least-squares solution of (38) using the SVD of is given by see again [49]. Note that the unique vector is referred to as least-squares solution since determining it from (40) is equivalent to minimizing the squared Euclidean residual norm of the linear systems given by (37) and (38); that is,

However, the least-square solution does not represent a suitable solution of (38) if the matrix is ill-conditioned. It can be shown that a matrix with a structure as shown in (37) or (A.3), also known as a Toeplitz matrix, is in fact ill conditioned [50, 51]. In that case, a small change in (e.g., due to projection noise) can cause a large change in . The rate at which a change in influences the solution is roughly proportional to the condition number of , defined as [49].

As an example, Figure 7 shows the solution of the example data from Figure 6. The solution is strongly oscillating and even has a rising amplitude. It is obvious that this solution has nothing in common with the real physiological behavior of the flow-scaled residue function.

In order to get a better understanding of why is not a meaningful solution and to motivate the regularization approach, we will investigate the individual terms of (40). We use the data shown in Figure 6 to obtain and . Figure 8 represents a plot of the absolute values of the expressions that occur in (40). These factors weight the right singular vectors of .

It is known from numerical analysis that the discrete Picard condition represents a means to analyze discrete ill-conditioned problems [50, 51]. This condition is violated, if the expressions do not decay faster, on average, than the singular values until a threshold value is reached where the singular values level off. The reader is referred to [51] for a more detailed explanation of the discrete Picard condition and its relation to the Picard condition from which it is derived. A usual reason for the violation of the discrete Picard condition is noise in the measured data that the matrix is based on. We can see that the discrete Picard condition is actually violated in the example shown in Figure 8 [52]. Consequently, the absolute values of the ratios —which represent the weighting factors of the right singular vectors —become very large.

To obtain a numerically stable result, a filter is used for regularization. The filter should suppress the influences of small singular values or, equivalently, the influences of high absolute values of the weighting factors . The regularized solution , where is a regularization parameter, is given by

We will focus on two common definitions of the filter factors . Firstly, the filter factors correspond to the truncated singular value decomposition (TSVD) approach and are defined with a sharp threshold at , Secondly, the filter factors are based on the Tikhonov regularization approach and characterized by a smooth weighting function centered around ,

The (absolute) regularization parameter is usually computed relative to the maximum singular value , that is, The relative regularization parameter is supposed to lie in the interval between 0 and 1.

In order to illustrate the Tikhonov filter factors, Figure 9 shows a plot of the function which is—unlike (44)—defined for a continuous range of . For determining , we assumed . It can be seen that, for increasing (i.e., stronger regularization), the values of decrease for all .

Interestingly, the solution of (38) using the filter factors is equivalent to minimizing the weighted sum of the squared Euclidean residual norm and the squared Euclidean solution norm ; that is,

Figure 10(a) shows the solution computed for two different regularization parameters. The solution for still shows some nonphysiological oscillations. However, the solution for can in fact be interpreted as a flow-scaled residue function in the presence of bolus delay and dispersion, compare Section 4.1. Figure 10(b) illustrates a plot of , which is proportional to CBF (see Section 3.3 and Table 2), as a function of . Apparently, CBF depends on the choice of regularization parameter. Choosing an optimal regularization parameter that will lead to physiologically reasonable estimates will be discussed in Section 4.4.

4.3. Alternative Deconvolution Approaches

The algebraic deconvolution approach from Section 4.2 is very commonly applied to analyze perfusion data. Yet, deconvolution problems arise in many other applications, and numerous alternative algorithms to solve these problems have been developed [53]. In this section, we provide a brief overview of alternative deconvolution approaches that have also been applied to perfusion data.

The Fourier transform represents a standard method to solve deconvolution problems [54], and it has also been evaluated to analyze perfusion data [45, 5557]. Interestingly, the Fourier transform-based deconvolution approach is mathematically equivalent to the SVD-based deconvolution approach with a block-circulant matrix , compare the appendix [47, 5860]. However, results obtained with SVD-based and FT-based deconvolution can be different because the chosen regularization approaches for these two methods are usually not equivalent. The regularization in the context of the Fourier-based deconvolution approach can be implemented by means of a modified Wiener filter [55], for example. The reader is referred to [60, 61] for a detailed analysis of the equivalence of SVD-based and Fourier-based regularization approaches.

In contrast to the model-independent deconvolution approaches also model-dependent approaches exist. Model-dependent approaches assume a certain shape of the residue function. For example, in [45, 62] a decaying exponential function was used which makes the deconvolution more stable since it reduces the degrees of freedom of the residue function [45]. However, if the underlaying residue function is different from the model the perfusion parameters may not be estimated correctly.

Deconvolution using orthogonal polynomials was investigated in [63]. An iterative deconvolution algorithm based on maximum likelihood expectation maximization (ML-EM) algorithm was proposed in [64]. An approach using Gaussian processes was evaluated in [65]. The deconvolution algorithm in [66] uses a nonlinear stochastic regularization method.

A comprehensive comparison of all available deconvolution methods has not been carried out yet. The SVD-based deconvolution approach, which is available in several software packages [6769], is comparably simple to implement and can be considered as the current standard method in perfusion image analysis.

4.4. Determination of the Regularization Parameter

Figure 10(b) demonstrated that (the maximum of) the solution depends on the regularization parameter . Consequently, the computed perfusion values—which can be derived from according to Table 2—vary for different . As an example, the CBF value will be underestimated systematically for large .

Therefore, an optimal choice of is crucial. A simple approach is to empirically determine a fixed value . This approach is often used in practice, and a typical value in brain perfusion CT is, for example, [68]. However, there exist more sophisticated approaches as well to determine the values independently for each voxel position [70]. Since the required amount of regularization depends roughly on the signal-to-noise ratio (SNR), these approaches can be more flexible when the SNR is spatially variant.

In [45, 47, 48], an oscillation index (OI) was defined to determine the intensity of oscillations of the flow-scaled residue function. The regularization can then be varied until the OI value falls below a certain threshold.

The L-curve criterion represents a model-independent method to determine (and ) [31, 71, 72]. The L-curve is defined by a double logarithmic plot of the squared Euclidean norm of the solution versus the squared Euclidean norm of the residual for a range of different values. The optimal regularization parameter can be found at the location of the characteristic corner point of the L-curve.

Another method to determine an appropriate regularization parameter is generalized cross-validation as described in [50, 73]. An implementation of the L-curve method and the generalized cross-validation can be found in [52].

Furthermore, a parameter estimation method that uses a priori knowledge of the behavior of the residue function was proposed in [74].

Kudo et al. [68] reported that two manufacturers applied a fixed threshold value in their perfusion analysis software. Unfortunately, the clinical use of methods with adaptive threshold values is rarely described in the current literature.

4.5. Perfusion Data Preprocessing

This section gives an overview of pre-processing techniques that can be applied in order to enhance the quality of the estimated perfusion parameters. Pre-processing occurs prior to the deconvolution step which may be implemented as described in Section 4.2.

A simple, yet mandatory, pre-processing step consists of the conversion to contrast agent concentration values, see Section 4.5.4. Further pre-processing steps are used to enhance the image quality (e.g., noise reduction) and to correct for artifacts (e.g., motion correction, partial volume correction) and specific properties of the blood (e.g., correction of differences in hematocrit). There are also pre-processing steps that can optimize the analysis of the perfusion value maps (e.g., segmentation of certain anatomic structures) and the application workflow (e.g., automated AIF estimation).

The order of the pre-processing steps presented in this section can act as a guideline for their practical implementation. However, a different ordering can of course be reasonable as well. Finally, this overview cannot include all details regarding suitable pre-processing steps. The reader is referred to the available literature for in-depth discussions.

4.5.1. Motion Correction

Patient motion (e.g., due to head movement or breathing) can result in a sudden change of the attenuation values at the fixed (stationary) voxel positions. Since this change in the attenuation value is caused by motion and not by contrast agent flow, the computed perfusion values can be severely biased. A practical approach for motion correction is to register all time frames of the reconstructed data set onto the first time frame [75]. A 3D registration should be used because it can also correct motion that occurs perpendicular to the orientation of the reconstructed slices. For a brain perfusion scan, a rigid registration may be sufficient. Conversely, in abdominal perfusion imaging, a non-rigid registration may be better suited to compensate for the deformations due to breathing.

As an alternative to registration, use of groupwise motion correction based on an optimization of a global cost function has been suggested [76]. There are also several approaches for motion correction in fMRI data [77]. These approaches may be used for perfusion MR data as well since both types of data typically consist of T2*-weighted EPI images [78]. However, the dynamic signal changes are relatively higher in DSC-MR data when compared to fMRI data [78]—in particular during the contrast agent bolus passage—which must be taken into account when adapting fMRI-based motion correction algorithms to DSC-MR data.

A related issue is streak artifact in reconstructed perfusion CT images that are caused by patient motion that occurs while the projection data corresponding to a single time frame is acquired. In perfusion MR images, ghosting artifacts can arise if the patient moves during the data acquisition. These kinds of artifact cannot be corrected by inter-frame motion correction. Instead, dedicated reconstruction algorithms would be required. As a practical alternative, time frames that exhibit severe reconstruction artifacts may simply be removed from the data set (i.e., from the series of successive time frames), which corresponds to the elimination of invalid sampling points of the voxel-specific time-concentration curves.

4.5.2. Noise Reduction

In the course of a perfusion exam, the measured signal in tissue that is caused by the contrast agent flow can be very low. For the case of perfusion CT, for example, tissue enhancements of less than 10 HU are measured. Hence, noise in the reconstructed images can be of a similar order of magnitude as the signal in tissue itself. Consequently, noise reduction should be taken into consideration in order to improve the accuracy of the estimated perfusion parameters.

Noise reduction can be implemented as a spatial smoothing of the data. Using a basic approach, each time frame can be filtered independently of the other time frames, and linear isotropic filters (e.g., based on a Gaussian filter kernel) may be applied. Alternatively, anisotropic filters that preserve edges and avoid blurring of large vessels can also be employed [79].

Both linear and nonlinear filtering in the temporal dimension—that is, between successive time frames—represent further methods for noise reduction [80]. It should be noted, however, that the regularization during the deconvolution step is equivalent to linear filtering in the temporal domain.

Recently, sophisticated 4D filtering techniques have been proposed that perform filtering in both the spatial and the temporal dimension and that are optimized for perfusion data [81, 82]. Fitting of the time-concentration curves to a model function such as a gamma-variate function is also a means for noise reduction [75].

4.5.3. Segmentation

A segmentation of certain anatomic structures in the reconstructed data set can optimize the perfusion image analysis [69, 83]. For example, the time-concentration curves could then be analyzed only in regions of interest where blood flow is actually expected [84]. Other regions such as air, bone, cerebrospinal fluid (CSF), and calcifications can be neglected. A segmentation and the subsequent removal of vessels is useful in order to optimize the quantitative analysis of perfusion parameters in tissue. Such a vessel segmentation can be performed prior to the deconvolution step, but it can also be implemented as a postprocessing step as described in [85].

4.5.4. Conversion to Contrast Agent Concentration

Neither for the case of CT imaging nor for the case of MR imaging can the time-concentration curves and be measured directly. Instead, the measurement is a superposition of the signal from the tissue itself and the contrast agent. Since the deconvolution approach presented in Section 4.2 expects that the functions and only refer to the signal caused by the contrast agent, the tissue signal must be subtracted. Furthermore, the measured signal must be converted to a contrast agent concentration value.

In perfusion CT, it is assumed that the (underlying) contrast agent concentration value is proportional to the (measured) X-ray attenuation value [86, 87]. Since deconvolution is a linear operation, the constant of proportionality does not influence the computed flow-scaled residue function. It can also be seen that the additional perfusion parameters from Section 3.5 are independent of this constant. Therefore, this constant is usually set to = 1 g/mL/HU for the sake of simplicity. The baseline value can be computed as the mean of during the acquired time frames before the contrast agent bolus arrives in the arterial input function. The conversion formula from an attenuation value (corresponding to a particular voxel) into the respective contrast agent concentration value then reads as

In perfusion MR, however, the contrast agent concentration value is not proportional to the received signal (in one voxel). Instead, it can be determined using the following formula: see [21]. Here, is a constant of proportionality which—with a similar argument as for —can have a norm of 1 and TE is the echo time of the MR sequence. It must be noted, however, that the constant can be different for blood and tissue due to differences in T2* relaxivities [37, 88]. This complicates absolute quantification of cerebral perfusion as discussed in [89]. Furthermore, studies have shown that fully oxygenated blood, for example, demonstrates a nonlinear relationship between the measured difference in T2* relaxation rate and contrast agent concentration [90].

Note that if only one time frame is considered as the baseline (i.e., if ), then , and the matrix defined by (37) and (38) will be rank deficient, compare Section 4.2.

4.5.5. Correction of Hematocrit Differences

Hematocrit (Hct) is a value that describes the proportion of the blood that consists of red blood cells. Hct is higher in arteries than in capillaries. Consequently, the proportion of the plasma in the blood, given by the difference (1-Hct), has a higher value in capillaries than in arteries. Since the contrast agent is distributed in the plasma only, the amount of plasma has a direct influence on the measured Hounsfield value or MR relaxation rate.

If the Hct difference is not corrected, it may bias the absolute quantification of the contrast agent concentration. A constant dimensionless correction factor , derived from the known Hct values in arteries and capillaries (often set to = 0.73) has been proposed [22, 85]. The measured time-concentration curve is then multiplied with to avoid the bias due to different Hct values.

4.5.6. Automated AIF Estimation

The total time for the perfusion image analysis can be shortened and the analysis can be made user independent by an automated estimation of the arterial input function. Several methods have been proposed that detect one global AIF [9193].

An interesting alternative approach is to estimate several local AIFs, which would be better suited to the theoretical model that was introduced in Section 3 [9496]. Since the local arteries are often small, this approach can have several disadvantages [89]. For example, partial volume effects—compare Section 4.5.7—can be more severe when compared to choosing one global AIF in a larger vessel. Perfusion analysis using local AIFs is actually investigated in [97] and the authors state that it produced more useful CBF maps.

Besides the arterial input function the venous outflow function could also be detected automatically. Knowledge about the venous outflow function could be used to automatically correct for partial volume effects which are described next.

4.5.7. Correction of Partial Volume Effects in the AIF

Due to limited spatial resolution in reconstructed perfusion CT and MR data, the AIF can suffer from partial volume effects [26]. This effect can lead to an underestimation of the AIF and consequently to incorrect perfusion values. To correct for partial volume effects in the AIF, several methods have been proposed [98100]. Commonly, the peak concentration value within a larger venous vessel or the area under the curve of a large venous vessel is used to rescale the AIF [31].

5. Summary

We have presented an overview of algorithms for the estimation of the most prominent perfusion parameters from CT or MR measurements that play an essential role in the assessment of flow altering diseases such as stroke, for example. In particular, we have emphasized the class of deconvolution-based methods that result from the application of the indicator-dilution theory, which is also derived in detail. Alternative approaches that do not use a deconvolution method are addressed briefly as well. The robust numerical solution of the resulting system of linear equations represents the second major topic of this paper. We have included a detailed discussion regarding the application of the singular value decomposition method as well as the practically relevant introduction of a suitable regularization technique in order to avoid physiologically unrealistic behavior of the estimated solution. Since this paper is intended to provide an introduction both to the underlying theory and to implementation-relevant aspects, we have provided a survey of preprocessing techniques that should be considered when designing a clinically useful tool for CT or MR perfusion analysis.

The novel contribution of this paper is to present the fundamental model, the mathematical deconvolution with regularization, and the practical pre-processing steps in one place. For a thorough understanding of perfusion image analysis, knowledge of all of these aspects is important and we have elaborated several links between these topics.

Appendix

The matrix in (38) can be replaced by a block-circulant matrix to reduce the influence of the bolus delay, compare Section 4.1, and thus to become independent of time shifts in the tissue time-concentration curve. Several studies actually exhibited an improvement of the accuracy of the perfusion estimates when using this alternative discretization method compared to the approach given by (36) [47, 48, 101]. On the other hand, in a receiver operating characteristics analysis—concerning infarct prediction in acute stroke patients—both discretization methods led to almost equal results [36].

The elements of —with denoting the row index () and denoting the column index () as usual—are defined as see (37). In order to assemble the block-circulant matrix , the size of the time series must be increased from to () using zero padding. We denote the new zero-padded time series as . The size of must be changed accordingly in order to retain consistency in (38).

The elements of the block-circulant matrix can then be defined as As an example, for , the matrix has the following structure:

The horizontal and vertical lines drawn in (A.3) subdivide the matrix into four quadrants. As can be seen, the matrix is a submatrix of , and it appears in the upper left and lower right quadrant.

Acknowledgments

The authors gratefully acknowledge funding of the Erlangen Graduate School in Advanced Optical Technologies (SAOT) by the German Research Foundation (DFG) in the framework of the German excellence initiative. Financial support was also provided through NIH 1K99EB007676 and the Lucas Foundation. Furthermore, the authors wish to give thanks to T. Struffert, MD, Department of Neuroradiology, Friedrich-Alexander University of Erlangen-Nuremberg, for providing the clinical data.