# Deconvolution-Based CT and MR Brain Perfusion Measurement: Theoretical Model Revisited and Practical Implementation Details

Authors: Andreas Fieselmann; Markus Kowarschik; Arundhuti Ganguly; Joachim Hornegger; Rebecca Fahrig
Journal: International Journal of Biomedical Imaging
 (2011)

---

## 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.


---

## Body

## 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 [3–5].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].Figure1  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.CT perfusion parameter maps of cerebral blood flow (CBF), cerebral blood volume (CBV), mean transit time (MTT), and time-to-peak (TTP). The ischemic stroke lesion is marked with arrows.
(a)
CBF map in mL/100 g/min(b)
CBV map in mL/100 g(c)
MTT map in s(d)
TTP map in sBlood 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 [17–19]. 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.


## 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 [17–19]. 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. Figure2 shows this model that consists of a volume of interest Vvoi covering the organ-specific parenchyma, the interstitial space, and  the capillary bed. The volumes of the parenchyma and the interstitial space are denoted by Vvoi*, while the volume of the capillary bed is referred to as Vcap. The entire volume of interest Vvoi=Vvoi*+Vcap 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 t it needs to pass through the capillary bed depends on the chosen path. We assume a stationary probability density distribution h(t) of transit times.Figure 2
Physiological model of the tissue perfusion. A blood cell can take several paths through the capillary bed. The variables are defined in Table1.Once a contrast agent bolus has been injected, it enters the volumeVvoi under consideration via the arterial inlet and is then diluted into the capillary bed. The local contrast agent concentrations cart(t) and cven(t) are measured directly adjacent to the capillary bed on the arterial and venous sides, respectively. Furthermore, the average contrast agent concentration cvoi(t) 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.Figure3 illustrates typical time-concentration curves cart(t), cvoi(t), and cven(t) 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.Examples of the time-concentration curvescart(t), cvoi(t), and cven(t) given in arbitrary units (a.u.). (b) Represents a zoomed view of (a) with a rescaled ordinate.
(a)
Typical time-concentration curves(b)
Zoomed view of (a)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 functionscart(t) and cvoi(t) is needed to compute the blood flow within the volume under consideration. In practice, the function cart(t)—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 functionh(t) of the transit times, that is,(1)MTT=∫0∞τh(τ)dτ.Furthermore, the residue (or residual) functionr(t)—compare [24]—represents an intermediate quantity of interest and is defined as(2)r(t)={1-∫0th(τ)dτ,fort≥0,0,fort<0.The (dimensionless) residue function thus quantifies the relative amount of contrast agent that is still inside the volumeVvoi of interest at time t after an (idealized) delta-shaped contrast agent bolus has entered the volume at the arterial inlet at time t=0; that is, cart(t)=δ(t). 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 r(0)=1 to 0. Figure 4 shows typical examples of a distribution function h(t) of transit times as well as the corresponding residue function r(t). In this example, the function h(t) is modeled by a gamma distribution [25].Examples of the distribution functionh(t) of transit times (the mean transit time is 4 s) and the corresponding residue function r(t).
(a)
Distributionh(t) of transit times(b)
Residue functionr(t)


### 3.2. Derivation of the Indicator-Dilution Theory

Using the parameters defined in Table1, the accumulated masses of contrast agent that have entered and left the volume of interest during the time interval [0,t], denoted as mc,voi,in(t) and mc,voi,out(t), respectively, can be expressed as(3)mc,voi,in(t)=F∫0tcart(τ)dτ,mc,voi,out(t)=F∫0tcven(τ)dτ.
The volume flow F is assumed to be constant over time. The contrast agent concentrations cart(t) and cven(t) at the arterial inlet and the venous outlet, respectively, are time-dependent functions which we assume to be 0 for t<0. These functions primarily depend on the parameters of the contrast agent injection and the patient's cardiac cycle.Table 1
Summary of parameters used to derive the indicator-dilution theory and to define clinically relevant tissue perfusion quantities.
VariableUnitDescriptionVvoimLTotal volume under considerationVcapmLVolume of the capillary bed within the volumeVvoiVvoi*mLVolumeVvoi without the volume of the capillary bed,Vvoi*=Vvoi-Vcapρvoig/mLMean density of the volumeVvoiρvoi*g/mLMean density of the volumeVvoi*mc,voi(t)gTotal mass of contrast agent in volumeVvoicart(t)g/mLLocal contrast agent concentration at the arterial inlet,cart(t)=dm/dV|t, measured at the arterial inletcven(t)g/mLLocal contrast agent concentration at the venous outlet,cven(t)=dm/dV|t, measured at the venous outletcvoi(t)g/mLAverage contrast agent concentration in the total volumeVvoi,cvoi(t)=mc,voi(t)/Vvoiccap(t)g/mLAverage contrast agent concentration in the capillary bed,ccap(t)=mc,voi(t)/Vcapcvoi*(t)g/mLAverage contrast agent concentration corresponding toVvoi*, cvoi*(t)=mc,voi(t)/Vvoi*FmL/sVolume flow at the arterial inlet and at the venous outleth(t)1/sProbability density function of the transit timesWe can compute the massmc,voi(t) of a contrast agent within the volume of interest at time t using the principle of conservation of mass as(4)mc,voi(t)=mc,voi,in(t)-mc,voi,out(t)=F∫0t(cart(τ)-cven(τ))dτ.The contrast agent concentrationcven(t) at the venous outlet can be computed from the contrast agent concentration cart(t) at the arterial inlet by convolving it with the probability density function h(t). We therefore obtain(5)cven(t)=∫-∞+∞cart(ξ)h(t-ξ)dξ.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 δ(t), as(6)mc,voi(t)=F∫0t(∫-∞+∞cart(ξ)δ(τ-ξ)dξ-∫-∞+∞cart(ξ)h(τ-ξ)dξ)dτ.
Changing the order of integration and rearranging this equation leads to(7)mc,voi(t)=F∫-∞+∞cart(ξ)(∫0t(δ(τ-ξ)-h(τ-ξ))dτ)dξ.
By applying the substitution τ′=τ-ξ, recalling that, for t≥0, we have(8)r(t)=1-∫0th(τ)dτ=∫0t(δ(τ)-h(τ))dτ,
and considering that h(t)=0 for t<0, we obtain(9)∫0t(δ(τ-ξ)-h(τ-ξ))dτ=∫-ξt-ξ(δ(τ′)-h(τ′))dτ′=r(t-ξ).
Equation (7) thus eventually reads(10)mc,voi(t)=F∫-∞+∞cart(ξ)r(t-ξ)dξ.
We introduce the cerebral blood flow (CBF) as the blood volume flow normalized by the mass of the volume Vvoi,(11)CBF=FVvoi⋅ρvoi.
Inserting this definition into (10) yields(12)mc,voi(t)Vvoi=CBF⋅ρvoi⋅∫-∞+∞cart(ξ)r(t-ξ)dξ.
According to Table 1, we define the contrast agent concentration cvoi(t) within the volume Vvoi of interest as(13)cvoi(t)=mc,voi(t)Vvoi,
which finally leads to the following formulation of the indicator-dilution theory,(14)cvoi(t)=CBF⋅ρvoi⋅∫-∞+∞cart(ξ)r(t-ξ)dξ=CBF⋅ρvoi⋅(cart*r)(t),
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 volumeVvoi*. This volume Vvoi* 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 Vvoi also contains the mass of the blood-filled capillary bed itself. Another aspect to consider is that the mean density ρvoi 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,(15)CBF*=FVvoi*⋅ρvoi*,
would then lead to a corresponding alternative formulation of the indicator-dilution theory,(16)cvoi*(t)=CBF*⋅ρvoi*⋅(cart*r)(t).
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 Section3.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 cart(t) and cvoi(t) can be measured and have known values whereas the values of CBF, r(t), and ρvoi 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 k(t),(17)k(t)=CBF⋅ρvoi⋅r(t),
which is given in units of 1/s and can be determined directly from the measured data cart(t) and cvoi(t). Using (17), (14) can be written as(18)cvoi(t)=(cart*k)(t).
Hence, k(t) can be obtained from the measured data cart(t) and cvoi(t) using a deconvolution method. Since a fundamental property of the residue function r(t) is r(0)=max(r(t))=1, we may then determine CBF as(19)CBF=1ρvoi⋅max⁡(k(t)).
Using max(k(t)) instead of k(0) has particular practical advantages that will be discussed in detail in Section 4.1.The flow-scaled residue functionk(t) can further be used to determine the MTT parameter of the tissue volume under consideration. From (2), it follows that, for t>0, we have(20)dr(t)dt=-h(t).
Equation (1) can thus be rewritten, and then using integration by parts and (17) and (19), we obtain(21)MTT=∫0∞τ(-dr(τ)dτ)dτ=∫0∞r(τ)dτ-lim⁡ξ→∞(τr(τ)|0ξ)=∫0∞r(τ)dτ=1max⁡(k(τ))⋅∫0∞k(τ)dτ.Note that we have assumed that there is a constantT>0 such that r(t)=0 for t>T. This assumption ensures that(22)lim⁡ξ→∞(τr(τ)|0ξ)=lim⁡ξ→∞(ξr(ξ))=0.The cerebral blood volume (CBV) corresponding to the tissue volumeVvoi represents another diagnostically relevant perfusion parameter and is defined as(23)CBV=Vcapρvoi⋅Vvoi.
It quantifies the blood volume normalized by the mass of Vvoi 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(24)CBF=CBVMTT
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 Vvoi 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 (25)CBV=MTT⋅CBF=1ρvoi⋅∫0∞k(τ)dτ,
which demonstrates that the CBV parameter can be derived from the flow-scaled residue function k(t) 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(26)CBV*=Vcapρvoi*⋅Vvoi*.
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 volumeVcap of the capillary bed to the entire volume Vvoi 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 cvoi(t) in the tissue volume can be related to the average contrast agent concentration ccap(t) in the capillary bed by(27)cvoi(t)=(ρvoi⋅CBV)⋅ccap(t).
According to the principle of conservation of mass, it follows that(28)mc,tot=F∫0∞cart(τ)dτ=F∫0∞ccap(τ)dτ=F∫0∞cven(τ)dτ,
where mc,tot is the total mass of contrast agent that has passed through the volume of interest. This results in an alternative expression for CBV,(29)CBV=1ρvoi⋅∫0∞cvoi(τ)dτ∫0∞cart(τ)dτ=1ρvoi⋅∫0∞cvoi(τ)dτ∫0∞cven(τ)dτ.
Hence, assuming a suitable correction for contrast agent recirculation [1, 34], CBV can be estimated from the integrals of either cvoi(t) and cart(t) or cvoi(t) and cven(t) 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(30)CBV=[cvoi(t)]max⁡[cven(t)]max⁡,
which avoids the computation of the integrals over time and only requires the maximum values of cvoi(t) and cven(t).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,(31)mc,voi(t)=mc,voi,in(t)=F⋅∫0tcart(τ)dτ.
Recalling the CBF definition—compare (11)—and that mc,voi(t)=cvoi(t)·Vvoi, we obtain(32)cvoi(t)=ρvoi⋅CBF⋅∫0tcart(τ)dτ.
Taking the derivative of (32) yields(33)dcvoi(t)dt=ρvoi⋅CBF⋅cart(t),
and since (33) must hold for all t, the blood flow is given by(34)[dcvoi(t)dt]max⁡=ρvoi⋅CBF⋅[cart(t)]max⁡,
which means that CBF can be estimated by dividing the maximum slope of the tissue time-concentration curve cvoi(t), shown as an example in Figure 5, by the maximum value of the contrast agent concentration cart(t) in the feeding artery.Figure 5
Perfusion parameters that are measured directly using the time-concentration curve. See Sections3.5 and 4.1 for explanations (BAT: bolus arrival time, TTP: time-to-peak, FM: first moment, AUC: area under the curve).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 concentrationcmax, 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.Figure5 illustrates the quantities cmax, TTP, and FM. The remaining parameter bolus arrival time (BAT) will be explained in Section 4.1. In practical measurements, the time point t=0 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, Table2 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.Table 2
Summary of perfusion parameters and how these parameters can be estimated using deconvolution-based and nondeconvolution-based methods.
Parameterw/Deconvolutionw/o DeconvolutionCBV(1/ρvoi)·∫0∞k(τ)dτ(1/ρvoi)·∫0∞cvoi(τ)dτ/∫0∞cart(τ)dτCBF(1/ρvoi)·max(k(t))(1/ρvoi)·[dcvoi(t)/dt]max/[cart(t)]maxMTT∫0∞k(τ)dτ/max(k(τ))see comment in Section3.5TTP—argmaxt(cvoi(t))FM—∫0∞cvoi(τ)τdτ/∫0∞cvoi(τ)dτ


## 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. Figure2 shows this model that consists of a volume of interest Vvoi covering the organ-specific parenchyma, the interstitial space, and  the capillary bed. The volumes of the parenchyma and the interstitial space are denoted by Vvoi*, while the volume of the capillary bed is referred to as Vcap. The entire volume of interest Vvoi=Vvoi*+Vcap 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 t it needs to pass through the capillary bed depends on the chosen path. We assume a stationary probability density distribution h(t) of transit times.Figure 2
Physiological model of the tissue perfusion. A blood cell can take several paths through the capillary bed. The variables are defined in Table1.Once a contrast agent bolus has been injected, it enters the volumeVvoi under consideration via the arterial inlet and is then diluted into the capillary bed. The local contrast agent concentrations cart(t) and cven(t) are measured directly adjacent to the capillary bed on the arterial and venous sides, respectively. Furthermore, the average contrast agent concentration cvoi(t) 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.Figure3 illustrates typical time-concentration curves cart(t), cvoi(t), and cven(t) 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.Examples of the time-concentration curvescart(t), cvoi(t), and cven(t) given in arbitrary units (a.u.). (b) Represents a zoomed view of (a) with a rescaled ordinate.
(a)
Typical time-concentration curves(b)
Zoomed view of (a)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 functionscart(t) and cvoi(t) is needed to compute the blood flow within the volume under consideration. In practice, the function cart(t)—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 functionh(t) of the transit times, that is,(1)MTT=∫0∞τh(τ)dτ.Furthermore, the residue (or residual) functionr(t)—compare [24]—represents an intermediate quantity of interest and is defined as(2)r(t)={1-∫0th(τ)dτ,fort≥0,0,fort<0.The (dimensionless) residue function thus quantifies the relative amount of contrast agent that is still inside the volumeVvoi of interest at time t after an (idealized) delta-shaped contrast agent bolus has entered the volume at the arterial inlet at time t=0; that is, cart(t)=δ(t). 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 r(0)=1 to 0. Figure 4 shows typical examples of a distribution function h(t) of transit times as well as the corresponding residue function r(t). In this example, the function h(t) is modeled by a gamma distribution [25].Examples of the distribution functionh(t) of transit times (the mean transit time is 4 s) and the corresponding residue function r(t).
(a)
Distributionh(t) of transit times(b)
Residue functionr(t)


## 3.2. Derivation of the Indicator-Dilution Theory

Using the parameters defined in Table1, the accumulated masses of contrast agent that have entered and left the volume of interest during the time interval [0,t], denoted as mc,voi,in(t) and mc,voi,out(t), respectively, can be expressed as(3)mc,voi,in(t)=F∫0tcart(τ)dτ,mc,voi,out(t)=F∫0tcven(τ)dτ.
The volume flow F is assumed to be constant over time. The contrast agent concentrations cart(t) and cven(t) at the arterial inlet and the venous outlet, respectively, are time-dependent functions which we assume to be 0 for t<0. These functions primarily depend on the parameters of the contrast agent injection and the patient's cardiac cycle.Table 1
Summary of parameters used to derive the indicator-dilution theory and to define clinically relevant tissue perfusion quantities.
VariableUnitDescriptionVvoimLTotal volume under considerationVcapmLVolume of the capillary bed within the volumeVvoiVvoi*mLVolumeVvoi without the volume of the capillary bed,Vvoi*=Vvoi-Vcapρvoig/mLMean density of the volumeVvoiρvoi*g/mLMean density of the volumeVvoi*mc,voi(t)gTotal mass of contrast agent in volumeVvoicart(t)g/mLLocal contrast agent concentration at the arterial inlet,cart(t)=dm/dV|t, measured at the arterial inletcven(t)g/mLLocal contrast agent concentration at the venous outlet,cven(t)=dm/dV|t, measured at the venous outletcvoi(t)g/mLAverage contrast agent concentration in the total volumeVvoi,cvoi(t)=mc,voi(t)/Vvoiccap(t)g/mLAverage contrast agent concentration in the capillary bed,ccap(t)=mc,voi(t)/Vcapcvoi*(t)g/mLAverage contrast agent concentration corresponding toVvoi*, cvoi*(t)=mc,voi(t)/Vvoi*FmL/sVolume flow at the arterial inlet and at the venous outleth(t)1/sProbability density function of the transit timesWe can compute the massmc,voi(t) of a contrast agent within the volume of interest at time t using the principle of conservation of mass as(4)mc,voi(t)=mc,voi,in(t)-mc,voi,out(t)=F∫0t(cart(τ)-cven(τ))dτ.The contrast agent concentrationcven(t) at the venous outlet can be computed from the contrast agent concentration cart(t) at the arterial inlet by convolving it with the probability density function h(t). We therefore obtain(5)cven(t)=∫-∞+∞cart(ξ)h(t-ξ)dξ.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 δ(t), as(6)mc,voi(t)=F∫0t(∫-∞+∞cart(ξ)δ(τ-ξ)dξ-∫-∞+∞cart(ξ)h(τ-ξ)dξ)dτ.
Changing the order of integration and rearranging this equation leads to(7)mc,voi(t)=F∫-∞+∞cart(ξ)(∫0t(δ(τ-ξ)-h(τ-ξ))dτ)dξ.
By applying the substitution τ′=τ-ξ, recalling that, for t≥0, we have(8)r(t)=1-∫0th(τ)dτ=∫0t(δ(τ)-h(τ))dτ,
and considering that h(t)=0 for t<0, we obtain(9)∫0t(δ(τ-ξ)-h(τ-ξ))dτ=∫-ξt-ξ(δ(τ′)-h(τ′))dτ′=r(t-ξ).
Equation (7) thus eventually reads(10)mc,voi(t)=F∫-∞+∞cart(ξ)r(t-ξ)dξ.
We introduce the cerebral blood flow (CBF) as the blood volume flow normalized by the mass of the volume Vvoi,(11)CBF=FVvoi⋅ρvoi.
Inserting this definition into (10) yields(12)mc,voi(t)Vvoi=CBF⋅ρvoi⋅∫-∞+∞cart(ξ)r(t-ξ)dξ.
According to Table 1, we define the contrast agent concentration cvoi(t) within the volume Vvoi of interest as(13)cvoi(t)=mc,voi(t)Vvoi,
which finally leads to the following formulation of the indicator-dilution theory,(14)cvoi(t)=CBF⋅ρvoi⋅∫-∞+∞cart(ξ)r(t-ξ)dξ=CBF⋅ρvoi⋅(cart*r)(t),
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 volumeVvoi*. This volume Vvoi* 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 Vvoi also contains the mass of the blood-filled capillary bed itself. Another aspect to consider is that the mean density ρvoi 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,(15)CBF*=FVvoi*⋅ρvoi*,
would then lead to a corresponding alternative formulation of the indicator-dilution theory,(16)cvoi*(t)=CBF*⋅ρvoi*⋅(cart*r)(t).
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 Section3.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 cart(t) and cvoi(t) can be measured and have known values whereas the values of CBF, r(t), and ρvoi 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 k(t),(17)k(t)=CBF⋅ρvoi⋅r(t),
which is given in units of 1/s and can be determined directly from the measured data cart(t) and cvoi(t). Using (17), (14) can be written as(18)cvoi(t)=(cart*k)(t).
Hence, k(t) can be obtained from the measured data cart(t) and cvoi(t) using a deconvolution method. Since a fundamental property of the residue function r(t) is r(0)=max(r(t))=1, we may then determine CBF as(19)CBF=1ρvoi⋅max⁡(k(t)).
Using max(k(t)) instead of k(0) has particular practical advantages that will be discussed in detail in Section 4.1.The flow-scaled residue functionk(t) can further be used to determine the MTT parameter of the tissue volume under consideration. From (2), it follows that, for t>0, we have(20)dr(t)dt=-h(t).
Equation (1) can thus be rewritten, and then using integration by parts and (17) and (19), we obtain(21)MTT=∫0∞τ(-dr(τ)dτ)dτ=∫0∞r(τ)dτ-lim⁡ξ→∞(τr(τ)|0ξ)=∫0∞r(τ)dτ=1max⁡(k(τ))⋅∫0∞k(τ)dτ.Note that we have assumed that there is a constantT>0 such that r(t)=0 for t>T. This assumption ensures that(22)lim⁡ξ→∞(τr(τ)|0ξ)=lim⁡ξ→∞(ξr(ξ))=0.The cerebral blood volume (CBV) corresponding to the tissue volumeVvoi represents another diagnostically relevant perfusion parameter and is defined as(23)CBV=Vcapρvoi⋅Vvoi.
It quantifies the blood volume normalized by the mass of Vvoi 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(24)CBF=CBVMTT
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 Vvoi 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 (25)CBV=MTT⋅CBF=1ρvoi⋅∫0∞k(τ)dτ,
which demonstrates that the CBV parameter can be derived from the flow-scaled residue function k(t) 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(26)CBV*=Vcapρvoi*⋅Vvoi*.
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 volumeVcap of the capillary bed to the entire volume Vvoi 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 cvoi(t) in the tissue volume can be related to the average contrast agent concentration ccap(t) in the capillary bed by(27)cvoi(t)=(ρvoi⋅CBV)⋅ccap(t).
According to the principle of conservation of mass, it follows that(28)mc,tot=F∫0∞cart(τ)dτ=F∫0∞ccap(τ)dτ=F∫0∞cven(τ)dτ,
where mc,tot is the total mass of contrast agent that has passed through the volume of interest. This results in an alternative expression for CBV,(29)CBV=1ρvoi⋅∫0∞cvoi(τ)dτ∫0∞cart(τ)dτ=1ρvoi⋅∫0∞cvoi(τ)dτ∫0∞cven(τ)dτ.
Hence, assuming a suitable correction for contrast agent recirculation [1, 34], CBV can be estimated from the integrals of either cvoi(t) and cart(t) or cvoi(t) and cven(t) 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(30)CBV=[cvoi(t)]max⁡[cven(t)]max⁡,
which avoids the computation of the integrals over time and only requires the maximum values of cvoi(t) and cven(t).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,(31)mc,voi(t)=mc,voi,in(t)=F⋅∫0tcart(τ)dτ.
Recalling the CBF definition—compare (11)—and that mc,voi(t)=cvoi(t)·Vvoi, we obtain(32)cvoi(t)=ρvoi⋅CBF⋅∫0tcart(τ)dτ.
Taking the derivative of (32) yields(33)dcvoi(t)dt=ρvoi⋅CBF⋅cart(t),
and since (33) must hold for all t, the blood flow is given by(34)[dcvoi(t)dt]max⁡=ρvoi⋅CBF⋅[cart(t)]max⁡,
which means that CBF can be estimated by dividing the maximum slope of the tissue time-concentration curve cvoi(t), shown as an example in Figure 5, by the maximum value of the contrast agent concentration cart(t) in the feeding artery.Figure 5
Perfusion parameters that are measured directly using the time-concentration curve. See Sections3.5 and 4.1 for explanations (BAT: bolus arrival time, TTP: time-to-peak, FM: first moment, AUC: area under the curve).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 concentrationcmax, 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.Figure5 illustrates the quantities cmax, TTP, and FM. The remaining parameter bolus arrival time (BAT) will be explained in Section 4.1. In practical measurements, the time point t=0 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, Table2 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.Table 2
Summary of perfusion parameters and how these parameters can be estimated using deconvolution-based and nondeconvolution-based methods.
Parameterw/Deconvolutionw/o DeconvolutionCBV(1/ρvoi)·∫0∞k(τ)dτ(1/ρvoi)·∫0∞cvoi(τ)dτ/∫0∞cart(τ)dτCBF(1/ρvoi)·max(k(t))(1/ρvoi)·[dcvoi(t)/dt]max/[cart(t)]maxMTT∫0∞k(τ)dτ/max(k(τ))see comment in Section3.5TTP—argmaxt(cvoi(t))FM—∫0∞cvoi(τ)τdτ/∫0∞cvoi(τ)dτ


## 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 Section3 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 Section3.1, we presented a model of microcirculation at the tissue level. We have assumed that we can measure the average contrast agent concentration cvoi(t) corresponding to a volume Vvoi under consideration which is supplied by one single capillary bed only. Furthermore, we have supposed that we can measure the contrast agent concentration cart(t) 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 volumeVvox (Vvox≫Vvoi) 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 Vvox is located completely within a larger artery or vein, there are of course no capillary beds located within Vvox.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 Vvox, we may adapt the model introduced in Section 3.1 as follows: the measured time-concentration curve cvoi(t) refers to the average perfusion from the arterioles through the capillary beds to the venules found in Vvox.The second adaptation of the model concerns the measurement ofcart(t). In reality, it is not possible to locally measure the concentration at the arterial inlet into the volume Vvox. 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 wherecvoi(t) 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 cvoi(t).The bolus delay has two implications. First, the curvecvoi(t) does not start to rise at the same time point as cart(t) 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 cvoi(t) 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 functionk(t) is equal to 0 from t=0 to t=BAT. In addition, due to the bolus dispersion, k(t) will not rise instantaneously to its maximum at t=BAT, but it will have a finite rise time. The time-to-maximum (TMAX) of the flow-scaled residue function, defined as(35)TMAX=argmax⁡t(k(t)),
has also been suggested as an additional perfusion parameter [41, 42]. Since the function k(t) can be 0 at t=0 (due to bolus delay), it is reasonable and recommended to estimate CBF as the maximum of k(t)—compare (19)—and not as the value of k(t) at time t=0.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 μart(tj) and μvoi(tj) shown in Figure 6 that were extracted from a real perfusion CT scan.Examples of measured time-attenuation curves in perfusion CT in (a) an arterial vessel and (b) in tissue. The time curves have been pre-processed by baseline subtraction and removal of the baseline time frames. The example data is measured atN=35 discrete time points.
(a)
Arterial time curveμart(tj)(b)
Tissue time curveμvoi(tj)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 Section4.5.4.In practice, the time-concentration curvescart(t) and cvoi(t) are sampled at discrete time points. We denote these time points as tj=(j-1)·Δt with j=1,…,N. A typical value of the sampling period Δt is 1 s, for example. We can discretize (18) as(36)cvoi(tj)=∫0∞cart(τ)k(tj-τ)dτ≈Δt∑i=1Ncart(ti)k(tj-i+1),
see [45]. We assume that the values of cart(t) can be neglected for t>NΔt. Since k(t)=0 for t<0, the end summation index could also be set to j instead of N. By rewriting this expression using matrix-vector notation, we obtain(37)Δt⋅(cart(t1)0⋯0cart(t2)cart(t1)⋯0⋮⋮⋱⋮cart(tN)cart(tN-1)⋯cart(t1))(k(t1)k(t2)⋮k(tN))=(cvoi(t1)cvoi(t2)⋮cvoi(tN)),
or shortly(38)Ak=c,
where Δt and cart(tj) are contained in the matrix A∈ℝN×N, and k(tj) and cvoi(tj) represent the entries of the vectors k∈ℝN and c∈ℝN, 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 A in order to reduce the influence of the bolus delay. See the appendix  for details.A standard approach to solve (37) for k is to use the singular value decomposition (SVD) of A. For a matrix A∈ℝN×N with r=rank(A) linearly independent rows and columns, it is defined as(39)A=U𝚺VT=∑i=1ruiσiviT,
where U=[u1,…,ur] and V=[v1,…,vr] are unique orthogonal matrices composed of the left and right singular vectors ui and vi, respectively [49]. The number of rows and columns in A that only contain zeros is determined by the number Nlz of leading zeros in the series cart(tj),j=1,…,N. Therefore, A has rank r≤N-Nlz. After the subtraction of the baseline, it may happen that the first entry cart(t1) is zero, see Section 4.5.4, and that A thus becomes rank-deficient. The diagonal matrix Σ=diag(σ1,…,σr) contains the singular values σi in nonincreasing order σ1≥σ2≥⋯≥σr>0. The least-squares solution kls of (38) using the SVD of A is given by(40)kls=∑i=1ruiTcσivi,
see again [49]. Note that the unique vector kls 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,(41)kls=argmin⁡k∈RN(‖Ak-c‖22).However, the least-square solutionkls does not represent a suitable solution of (38) if the matrix A is ill-conditioned. It can be shown that a matrix A 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 c (e.g., due to projection noise) can cause a large change in kls. The rate at which a change in c influences the solution kls is roughly proportional to the condition number of A, defined as σ1/σr [49].As an example, Figure7 shows the solution kls 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.Least-squares solution vectorkls of (38) using the example data from Figure 6. (kls)j denotes the jth entry of the vector kls. The plot shown in (a) illustrates the strong oscillations of kls. The plot given in (b) shows the amplitude |kls| of this function on a logarithmic scale.
(a)
kls on a linear scale(b)
|kls| on a logarithmic scaleIn order to get a better understanding of whykls 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 A and c. Figure 8 represents a plot of the absolute values of the expressions (uiTc)/σi that occur in (40). These factors weight the right singular vectors vi of A.Figure 8
SVD analysis of the matrixA constructed from the example data shown in Figure 6. The plot displays the absolute values of the weighting factors (uiTc)/σi and of their individual components |uiTc| and σi on a logarithmic scale.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 uiTc do not decay faster, on average, than the singular values σi 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 A 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 (uiTc)/σi—which represent the weighting factors of the right singular vectors vi—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σi or, equivalently, the influences of high absolute values of the weighting factors (uiTc)/σi. The regularized solution kλ, where λ is a regularization parameter, is given by(42)kλ=∑i=1r(fλ,iuiTcσi)vi.We will focus on two common definitions of the filter factorsfλ,i. Firstly, the filter factors fλ,i(tsvd) correspond to the truncated singular value decomposition (TSVD) approach and are defined with a sharp threshold at λ,(43)fλ,i(tsvd)={0,forσi<λ,1,forσi≥λ.
Secondly, the filter factors fλ,i(tikh) are based on the Tikhonov regularization approach and characterized by a smooth weighting function centered around λ,(44)fλ,i(tikh)=σi2σi2+λ2.The (absolute) regularization parameterλ is usually computed relative to the maximum singular value σ1, that is,(45)λ=λrel⁡σ1.
The relative regularization parameter λrel is supposed to lie in the interval between 0 and 1.In order to illustrate the Tikhonov filter factors, Figure9 shows a plot of the function fλ(tikh)=σ2/(σ2+λ2) which is—unlike (44)—defined for a continuous range of σ. For determining fλ(tikh), we assumed σ1=1. It can be seen that, for increasing λ (i.e., stronger regularization), the values of fλ(tikh) decrease for all σ.(a) Linear and (b) double logarithmic plot of the Tikhonov filter factorfλ(tikh) as a function of the singular value σ∈[10-5,1].
(a)
Linear plot offλ(tikh)(b)
Double logarithmic plot offλ(tikh)Interestingly, the solutionkλ(tikh) of (38) using the filter factors fλ,i(tikh) is equivalent to minimizing the weighted sum of the squared Euclidean residual norm ∥Ak-c∥22 and the squared Euclidean solution norm ∥k∥22; that is,(46)kλ(tikh)=argmin⁡k∈RN(‖Ak-c‖22+λ2‖k‖22).Figure10(a) shows the solution kλ(tikh) computed for two different regularization parameters. The solution for λrel=0.1 still shows some nonphysiological oscillations. However, the solution for λrel=0.3 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 max(kλ(tikh)), which is proportional to CBF (see Section 3.3 and Table 2), as a function of λrel. 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.Deconvolution with Tikhonov regularization: (a) Regularized solutionkλ(tikh) for two different regularization parameters λrel and (b) maximum of kλ(tikh) as a function of λrel. (kλ(tikh))j denotes the jth entry of the vector kλ(tikh).
(a)
Regularized solutionkλ(tikh)(b)
Dependency of max(kλ(tikh)) on λrel plotted on a logarithmic scale


### 4.3. Alternative Deconvolution Approaches

The algebraic deconvolution approach from Section4.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, 55–57]. Interestingly, the Fourier transform-based deconvolution approach is mathematically equivalent to the SVD-based deconvolution approach with a block-circulant matrix A, compare the appendix [47, 58–60]. 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 [67–69], 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

Figure10(b) demonstrated that (the maximum of) the solution kλ(tikh) depends on the regularization parameter λrel. Consequently, the computed perfusion values—which can be derived from kλ(tikh) according to Table 2—vary for different λrel. As an example, the CBF value will be underestimated systematically for large λrel.Therefore, an optimal choice ofλrel is crucial. A simple approach is to empirically determine a fixed value λrel. This approach is often used in practice, and a typical value in brain perfusion CT is, for example, λrel=0.2 [68]. However, there exist more sophisticated approaches as well to determine the values λrel 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 λrel) [31, 71, 72]. The L-curve is defined by a double logarithmic plot of the squared Euclidean norm ∥kλ∥22 of the solution versus the squared Euclidean norm ∥Akλ-c∥22 of the residual for a range of different λ values. The optimal regularization parameter λopt 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 λrel 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 Section4.2.A simple, yet mandatory, pre-processing step consists of the conversion to contrast agent concentration values, see Section4.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 curvescart(t) and cvoi(t) 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 cart(tj) and cvoi(tj) 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 kct = 1 g/mL/HU for the sake of simplicity. The baseline value μ0 can be computed as the mean of μ(tj) during the B acquired time frames before the contrast agent bolus arrives in the arterial input function. The conversion formula from an attenuation value μ(tj) (corresponding to a particular voxel) into the respective contrast agent concentration value c(tj) then reads as(47)c(tj)=kct(μ(tj+B-1)-μ0),μ0=1B∑i=1Bμ(ti).In perfusion MR, however, the contrast agent concentration value is not proportional to the received signals(tj) (in one voxel). Instead, it can be determined using the following formula:(48)c(tj)=-kmrTEln⁡(s(tj+B-1)s0),s0=1B∑i=1Bs(ti),
see [21]. Here, kmr is a constant of proportionality which—with a similar argument as for kct—can have a norm of 1 and TE is the echo time of the MR sequence. It must be noted, however, that the constant kmr 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., ifB=1), then c(0)=0, and the matrix A 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 cvoi(t) 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 [91–93].An interesting alternative approach is to estimate several local AIFs, which would be better suited to the theoretical model that was introduced in Section3 [94–96]. 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 functioncart(t) the venous outflow function cven(t) 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 [98–100]. 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].


## 4.1. Adaptations of the Model of Microcirculation

In Section3.1, we presented a model of microcirculation at the tissue level. We have assumed that we can measure the average contrast agent concentration cvoi(t) corresponding to a volume Vvoi under consideration which is supplied by one single capillary bed only. Furthermore, we have supposed that we can measure the contrast agent concentration cart(t) 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 volumeVvox (Vvox≫Vvoi) 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 Vvox is located completely within a larger artery or vein, there are of course no capillary beds located within Vvox.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 Vvox, we may adapt the model introduced in Section 3.1 as follows: the measured time-concentration curve cvoi(t) refers to the average perfusion from the arterioles through the capillary beds to the venules found in Vvox.The second adaptation of the model concerns the measurement ofcart(t). In reality, it is not possible to locally measure the concentration at the arterial inlet into the volume Vvox. 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 wherecvoi(t) 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 cvoi(t).The bolus delay has two implications. First, the curvecvoi(t) does not start to rise at the same time point as cart(t) 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 cvoi(t) 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 functionk(t) is equal to 0 from t=0 to t=BAT. In addition, due to the bolus dispersion, k(t) will not rise instantaneously to its maximum at t=BAT, but it will have a finite rise time. The time-to-maximum (TMAX) of the flow-scaled residue function, defined as(35)TMAX=argmax⁡t(k(t)),
has also been suggested as an additional perfusion parameter [41, 42]. Since the function k(t) can be 0 at t=0 (due to bolus delay), it is reasonable and recommended to estimate CBF as the maximum of k(t)—compare (19)—and not as the value of k(t) at time t=0.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 μart(tj) and μvoi(tj) shown in Figure 6 that were extracted from a real perfusion CT scan.Examples of measured time-attenuation curves in perfusion CT in (a) an arterial vessel and (b) in tissue. The time curves have been pre-processed by baseline subtraction and removal of the baseline time frames. The example data is measured atN=35 discrete time points.
(a)
Arterial time curveμart(tj)(b)
Tissue time curveμvoi(tj)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 Section4.5.4.In practice, the time-concentration curvescart(t) and cvoi(t) are sampled at discrete time points. We denote these time points as tj=(j-1)·Δt with j=1,…,N. A typical value of the sampling period Δt is 1 s, for example. We can discretize (18) as(36)cvoi(tj)=∫0∞cart(τ)k(tj-τ)dτ≈Δt∑i=1Ncart(ti)k(tj-i+1),
see [45]. We assume that the values of cart(t) can be neglected for t>NΔt. Since k(t)=0 for t<0, the end summation index could also be set to j instead of N. By rewriting this expression using matrix-vector notation, we obtain(37)Δt⋅(cart(t1)0⋯0cart(t2)cart(t1)⋯0⋮⋮⋱⋮cart(tN)cart(tN-1)⋯cart(t1))(k(t1)k(t2)⋮k(tN))=(cvoi(t1)cvoi(t2)⋮cvoi(tN)),
or shortly(38)Ak=c,
where Δt and cart(tj) are contained in the matrix A∈ℝN×N, and k(tj) and cvoi(tj) represent the entries of the vectors k∈ℝN and c∈ℝN, 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 A in order to reduce the influence of the bolus delay. See the appendix  for details.A standard approach to solve (37) for k is to use the singular value decomposition (SVD) of A. For a matrix A∈ℝN×N with r=rank(A) linearly independent rows and columns, it is defined as(39)A=U𝚺VT=∑i=1ruiσiviT,
where U=[u1,…,ur] and V=[v1,…,vr] are unique orthogonal matrices composed of the left and right singular vectors ui and vi, respectively [49]. The number of rows and columns in A that only contain zeros is determined by the number Nlz of leading zeros in the series cart(tj),j=1,…,N. Therefore, A has rank r≤N-Nlz. After the subtraction of the baseline, it may happen that the first entry cart(t1) is zero, see Section 4.5.4, and that A thus becomes rank-deficient. The diagonal matrix Σ=diag(σ1,…,σr) contains the singular values σi in nonincreasing order σ1≥σ2≥⋯≥σr>0. The least-squares solution kls of (38) using the SVD of A is given by(40)kls=∑i=1ruiTcσivi,
see again [49]. Note that the unique vector kls 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,(41)kls=argmin⁡k∈RN(‖Ak-c‖22).However, the least-square solutionkls does not represent a suitable solution of (38) if the matrix A is ill-conditioned. It can be shown that a matrix A 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 c (e.g., due to projection noise) can cause a large change in kls. The rate at which a change in c influences the solution kls is roughly proportional to the condition number of A, defined as σ1/σr [49].As an example, Figure7 shows the solution kls 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.Least-squares solution vectorkls of (38) using the example data from Figure 6. (kls)j denotes the jth entry of the vector kls. The plot shown in (a) illustrates the strong oscillations of kls. The plot given in (b) shows the amplitude |kls| of this function on a logarithmic scale.
(a)
kls on a linear scale(b)
|kls| on a logarithmic scaleIn order to get a better understanding of whykls 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 A and c. Figure 8 represents a plot of the absolute values of the expressions (uiTc)/σi that occur in (40). These factors weight the right singular vectors vi of A.Figure 8
SVD analysis of the matrixA constructed from the example data shown in Figure 6. The plot displays the absolute values of the weighting factors (uiTc)/σi and of their individual components |uiTc| and σi on a logarithmic scale.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 uiTc do not decay faster, on average, than the singular values σi 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 A 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 (uiTc)/σi—which represent the weighting factors of the right singular vectors vi—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σi or, equivalently, the influences of high absolute values of the weighting factors (uiTc)/σi. The regularized solution kλ, where λ is a regularization parameter, is given by(42)kλ=∑i=1r(fλ,iuiTcσi)vi.We will focus on two common definitions of the filter factorsfλ,i. Firstly, the filter factors fλ,i(tsvd) correspond to the truncated singular value decomposition (TSVD) approach and are defined with a sharp threshold at λ,(43)fλ,i(tsvd)={0,forσi<λ,1,forσi≥λ.
Secondly, the filter factors fλ,i(tikh) are based on the Tikhonov regularization approach and characterized by a smooth weighting function centered around λ,(44)fλ,i(tikh)=σi2σi2+λ2.The (absolute) regularization parameterλ is usually computed relative to the maximum singular value σ1, that is,(45)λ=λrel⁡σ1.
The relative regularization parameter λrel is supposed to lie in the interval between 0 and 1.In order to illustrate the Tikhonov filter factors, Figure9 shows a plot of the function fλ(tikh)=σ2/(σ2+λ2) which is—unlike (44)—defined for a continuous range of σ. For determining fλ(tikh), we assumed σ1=1. It can be seen that, for increasing λ (i.e., stronger regularization), the values of fλ(tikh) decrease for all σ.(a) Linear and (b) double logarithmic plot of the Tikhonov filter factorfλ(tikh) as a function of the singular value σ∈[10-5,1].
(a)
Linear plot offλ(tikh)(b)
Double logarithmic plot offλ(tikh)Interestingly, the solutionkλ(tikh) of (38) using the filter factors fλ,i(tikh) is equivalent to minimizing the weighted sum of the squared Euclidean residual norm ∥Ak-c∥22 and the squared Euclidean solution norm ∥k∥22; that is,(46)kλ(tikh)=argmin⁡k∈RN(‖Ak-c‖22+λ2‖k‖22).Figure10(a) shows the solution kλ(tikh) computed for two different regularization parameters. The solution for λrel=0.1 still shows some nonphysiological oscillations. However, the solution for λrel=0.3 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 max(kλ(tikh)), which is proportional to CBF (see Section 3.3 and Table 2), as a function of λrel. 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.Deconvolution with Tikhonov regularization: (a) Regularized solutionkλ(tikh) for two different regularization parameters λrel and (b) maximum of kλ(tikh) as a function of λrel. (kλ(tikh))j denotes the jth entry of the vector kλ(tikh).
(a)
Regularized solutionkλ(tikh)(b)
Dependency of max(kλ(tikh)) on λrel plotted on a logarithmic scale


## 4.3. Alternative Deconvolution Approaches

The algebraic deconvolution approach from Section4.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, 55–57]. Interestingly, the Fourier transform-based deconvolution approach is mathematically equivalent to the SVD-based deconvolution approach with a block-circulant matrix A, compare the appendix [47, 58–60]. 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 [67–69], 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

Figure10(b) demonstrated that (the maximum of) the solution kλ(tikh) depends on the regularization parameter λrel. Consequently, the computed perfusion values—which can be derived from kλ(tikh) according to Table 2—vary for different λrel. As an example, the CBF value will be underestimated systematically for large λrel.Therefore, an optimal choice ofλrel is crucial. A simple approach is to empirically determine a fixed value λrel. This approach is often used in practice, and a typical value in brain perfusion CT is, for example, λrel=0.2 [68]. However, there exist more sophisticated approaches as well to determine the values λrel 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 λrel) [31, 71, 72]. The L-curve is defined by a double logarithmic plot of the squared Euclidean norm ∥kλ∥22 of the solution versus the squared Euclidean norm ∥Akλ-c∥22 of the residual for a range of different λ values. The optimal regularization parameter λopt 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 λrel 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 Section4.2.A simple, yet mandatory, pre-processing step consists of the conversion to contrast agent concentration values, see Section4.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 curvescart(t) and cvoi(t) 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 cart(tj) and cvoi(tj) 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 kct = 1 g/mL/HU for the sake of simplicity. The baseline value μ0 can be computed as the mean of μ(tj) during the B acquired time frames before the contrast agent bolus arrives in the arterial input function. The conversion formula from an attenuation value μ(tj) (corresponding to a particular voxel) into the respective contrast agent concentration value c(tj) then reads as(47)c(tj)=kct(μ(tj+B-1)-μ0),μ0=1B∑i=1Bμ(ti).In perfusion MR, however, the contrast agent concentration value is not proportional to the received signals(tj) (in one voxel). Instead, it can be determined using the following formula:(48)c(tj)=-kmrTEln⁡(s(tj+B-1)s0),s0=1B∑i=1Bs(ti),
see [21]. Here, kmr is a constant of proportionality which—with a similar argument as for kct—can have a norm of 1 and TE is the echo time of the MR sequence. It must be noted, however, that the constant kmr 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., ifB=1), then c(0)=0, and the matrix A 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 cvoi(t) 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 [91–93].An interesting alternative approach is to estimate several local AIFs, which would be better suited to the theoretical model that was introduced in Section3 [94–96]. 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 functioncart(t) the venous outflow function cven(t) 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 [98–100]. 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].


## 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 curvescart(t) and cvoi(t) 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 cart(tj) and cvoi(tj) 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 kct = 1 g/mL/HU for the sake of simplicity. The baseline value μ0 can be computed as the mean of μ(tj) during the B acquired time frames before the contrast agent bolus arrives in the arterial input function. The conversion formula from an attenuation value μ(tj) (corresponding to a particular voxel) into the respective contrast agent concentration value c(tj) then reads as(47)c(tj)=kct(μ(tj+B-1)-μ0),μ0=1B∑i=1Bμ(ti).In perfusion MR, however, the contrast agent concentration value is not proportional to the received signals(tj) (in one voxel). Instead, it can be determined using the following formula:(48)c(tj)=-kmrTEln⁡(s(tj+B-1)s0),s0=1B∑i=1Bs(ti),
see [21]. Here, kmr is a constant of proportionality which—with a similar argument as for kct—can have a norm of 1 and TE is the echo time of the MR sequence. It must be noted, however, that the constant kmr 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., ifB=1), then c(0)=0, and the matrix A 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 cvoi(t) 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 [91–93].An interesting alternative approach is to estimate several local AIFs, which would be better suited to the theoretical model that was introduced in Section3 [94–96]. 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 functioncart(t) the venous outflow function cven(t) 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 [98–100]. 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.



---

*Source: 467563-2011-08-28.xml*