tissue perfusion measurement from iodinated contrast agent enhancement on ct scans was first proposed by axel in 1980 ; this was based on earlier developments by meier and zierler   for measuring blood flow and blood volume . at that time , the ct - based measurements were strictly limited to research because of the low speeds and narrow coverage of the existing ct scanners . however , the introduction of perfusion ct ( pct ) helped expand the utility of ct significantly since it could now provide capillary level hemodynamic information . within about a decade ,
perfusion imaging techniques were also adopted in mr [ 35 ] . with the advent of helical scanners and faster rotating gantries ( 0.330.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 .
  perfusion imaging is most widely used in acute stroke and oncology . when used in diagnosis of stroke ,
the purpose of perfusion imaging is to identify the extent of affected tissue and to delineate the ischemic tissue that can be reperfused . in oncology , perfusion imaging helps to identify angiogenetic tumors that alter the local tissue perfusion due to generation of neovasculature .
perfusion measurements are increasingly being used for assessment , staging , and monitoring posttherapy [ 6 , 7 ] . 

figure 1   shows common parameter maps based on a brain perfusion ct exam ( somatom definition as+ , siemens ag , healthcare sector , forchheim , germany ) of a 69-year - old male stroke patient .
the patient presented to the hospital with an acute high - grade hemiparesis on the right side .
the time - to - peak ( ttp ) image shows a large lesion that illustrates the maximum affected tissue . in addition , the cerebral blood flow ( cbf ) , cerebral blood volume ( cbv ) , and mean transit time ( mtt ) images exhibit perfusion deficits in a smaller brain territory . in general , these perfusion
ct maps are interpreted appropriately in order to guide the recanalization procedure of the occluded vessel .
blood flow is critical to the functionality of any organ since it provides the essential nutrients and oxygen . in case of flow disruption ,
the body autoregulates the flow and pressure either by altering blood flow or volume or both . in the brain , there are some fairly well - defined thresholds for the cerebral blood flow in normal , reversibly damaged , and necrotic tissue .
the normal value for the cerebral blood flow is between 50 and 60  ml/100  g / min for grey matter .
the average value decreases with age and is about 2 to 3 times lower in white matter compared to grey matter .
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 810  ml/100  g / min , the cell membrane pumps fail , causing irreversible damage to the cells .   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 6065% [ 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 710  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 .
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 .
one study that compared stroke diagnosis using ct perfusion plus angiography , against mri , found good correlation and no significant prognosis differences .
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   and 3.4  msv .
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 .
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 .
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
gre sequences are most frequently used because they provide a better contrast - to - noise ratio for imaging of the contrast agent compared to se sequences [ 1719 ]
. however , gre sequences have the disadvantage of disproportionately weighting the contribution of the contrast agent in relatively large vessels , whereas se sequences provide a more accurate assessment of blood flow through vessels of all sizes . after the 710  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 .
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 . for computing the tissue perfusion , we assume a physiological model of the blood supply to the tissue .  
figure 2 shows this model that consists of a volume of interest 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
the transit time t it needs to pass through the capillary bed depends on the chosen path .
once a contrast agent bolus has been injected , it enters the volume vvoi 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 . for the following analysis
, we assume the contrast agent concentration to be measured as mass per volume , which can easily be related to amount per volume . 

figure 3 illustrates typical time - concentration curves cart(t ) , cvoi(t ) , and cven(t ) that may be measured in brain tissue , for example . for the sake of simplicity
note that the ( average ) enhancement within the volume of interest is commonly more than an order of magnitude below the enhancements of the feeding artery and the draining vein .
an additional important assumption is that the contrast agent remains in the intravascular space . for our case of cerebral perfusion
, it should therefore not cross the blood - brain barrier ( bbb ) . as a consequence ,
this means that all contrast agents entering from the arterial inlet will eventually leave the volume of interest at the venous outlet .
a breakdown of the bbb may occur in tumor patients , in stroke patients , and in patients that suffer from inflammations or infections , for example . in these cases , the methods presented in this paper may lead to inaccurate perfusion estimates and particularly to an overestimation of the blood volume [ 21 , 22 ] .
note that there exist other modelling approaches which do not assume that the contrast agent remains in the intravascular space .
these models can be used for measuring tumor perfusion , for example [ 6 , 20 , 23 ] .
finally , we suppose that the contrast agent mixes perfectly with the blood and that the physical properties of the blood ( its flow behavior , in particular ) are not influenced by the contrast agent .
as we will see , only knowledge of the functions cart(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 function h(t ) of the transit times , that is , furthermore , the residue ( or residual ) function r(t)compare represents an intermediate quantity of interest and is defined as 
 ( 2)r(t)={10th()d,for  t0,0,for  t<0 . 
 
the ( dimensionless ) residue function thus quantifies the relative amount of contrast agent that is still inside the volume vvoi 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 . using the parameters defined in table 1 ,
the accumulated masses of contrast agent that have entered and left the volume of interest during the time interval [ 0 , t ] , denoted as mc , voi , in(t ) and mc , voi , out(t ) , respectively , can be expressed as 
 ( 3)mc , voi , in(t)=f0tcart()d,mc , voi , out(t)=f0tcven()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 .
we can compute the mass mc , 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)=f0t(cart()cven())d. 
  the contrast agent concentration cven(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)=f0t(+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)=10th()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 = fvvoivoi . 
							
inserting this definition into ( 10 ) yields 
 ( 12)mc , voi(t)vvoi = cbfvoi+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)=cbfvoi+cart()r(t)d=cbfvoi(cartr)(t ) , 
							 where  denotes the convolution operator as usual , see also [ 21 , 26 ] .
a historical overview of the development of the indicator - dilution theory with numerous references to mathematical aspects can be found in .
note that the solution of ( 14 ) with respect to cbf and other clinically important perfusion parameters will be discussed in section 3.3 .
from a physiological point of view , it would be more meaningful to normalize cbf by the mass of the volume vvoi*. 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=fvvoivoi , 
							 would then lead to a corresponding alternative formulation of the indicator - dilution theory , 
 ( 16)cvoi(t)=cbfvoi(cartr)(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 section 3.1but 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   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 . 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)=cbfvoir(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)=(cartk)(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=1voimax  ( 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 function k(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=0r()dlim  (r()|0)=0r()d=1max  ( k())0k()d. 
  note that we have assumed that there is a constant t > 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 volume vvoi represents another diagnostically relevant perfusion parameter and is defined as 
 ( 23)cbv = vcapvoivvoi . 
							
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 .
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 = mttcbf=1voi0k()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 .
note that the definition of cbv that corresponds to the alternative definition of cbf in   is 
 ( 26)cbv=vcapvoivvoi. 
							 accordingly , this alternative definition relates the blood volume to the mass of the parenchyma ( and the interstitium ) only and explicitly omits the mass of the capillary bed itself .
furthermore , there are references in the literature that suggest measuring the blood volume in units of ml / ml . this alternative dimensionless quantity may therefore be considered as a measure of blood ( or vascular ) volume fraction . when relating the volume vcap 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   for both technical and clinical details . 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 .
firstly , there is an alternative method to compute the blood volume of the tissue volume under consideration .
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)=(voicbv)ccap(t ) . 
							
according to the principle of conservation of mass , it follows that 
 ( 28)mc , tot = f0cart()d=f0ccap()d=f0cven()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=1voi0cvoi()d0cart()d=1voi0cvoi()d0cven()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   for details and further references with a particular focus on mr perfusion measurements .
it is argued in   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)=f0tcart()d. 
							 recalling the cbf definition  compare ( 11)and that mc , voi(t ) = cvoi(t )  vvoi , we obtain 
 ( 32)cvoi(t)=voicbf0tcart()d. 
							 taking the derivative of ( 32 ) yields 
 ( 33)dcvoi(t)dt=voicbfcart(t ) , 
							 and since ( 33 ) must hold for all t , the blood flow is given by 
 ( 34)[dcvoi(t)dt]max  = voicbf[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 .
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 . according to ,
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 .
besides the aforementioned quantities cbv , cbf , and mtt , there are additional perfusion parameters such as the time - to - peak ( ttp ) of the time - concentration curve , the maximum contrast agent concentration cmax , 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 . 

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
a comparison of several perfusion parameters and their clinical impact on the treatment of stroke patients is given in .
in summary , table 2 covers the most common diagnostically relevant perfusion parameters and shows how they can be determined employing deconvolution - based and nondeconvolution - based methods . in principle , the central volume theorem  compare ( 24)may also be used to numerically estimate the mtt from the parameters cbf and cbv when the latter have been computed using nondeconvolution - based algorithms .
however , the authors are not aware of any reference that describes the application of this approach in clinical practice .
this section is devoted to the practical computer implementation of algorithms for perfusion image analysis .
first , we will discuss the necessary adaptations of the theoretical model from section 3 that are needed for its application to data from real ct and mr scanners .
afterwards , we will describe commonly used algebraic deconvolution methods and also give an overview of alternative approaches .
we will motivate the need for suitable regularization and discuss the influence of the regularization parameter on the resulting perfusion estimates . for the sake of completeness
, we will also address techniques for the pre - processing of the acquired perfusion data . in section 3.1
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 can not rely on these two aforementioned assumptions .
we will thus introduce two major adaptations of the physiological model which are necessary once it is to be applied to data from real scanners .
first , during a standard ct and mr perfusion exam , a volume of interest is scanned and the data is reconstructed on a grid of regularly spaced voxels . in the object domain , each voxel volume vvox ( 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 .
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 .
, 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
this approach leads to a traveling time of the contrast agent bolus from where the aif is measured to the location of the tissue volume where cvoi(t ) is measured .
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 ) .
first , the curve cvoi(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 .
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 function k(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 = arg  max  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 . in order to correct for bolus delay and dispersion several methods
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 . 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 .
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 .
we assume that the measured time - attenuation curves can be converted to time - concentration curves using a constant of proportionality of 1  g / ml / hu . details about the conversion , also discussing perfusion mr data , will be explained in section 4.5.4 . in practice ,
the time - concentration curves cart(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)=0cart()k(tj)dti=1ncart(ti)k(tji+1 ) , 
							 see .
we assume that the values of cart(t ) can be neglected for t > nt .
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)00cart(t2)cart(t1)0cart(tn)cart(tn1)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   , and k(tj ) and cvoi(tj ) represent the entries of the vectors k   and c   , respectively .
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 .
a standard approach to solve ( 37 ) for k is to use the singular value decomposition ( svd ) of a. for a matrix a   with r = rank  ( a ) linearly independent rows and columns , it is defined as 
 ( 39)a = uvt=i=1ruiivit , 
							 where u = [ u1 ,  , ur ] and v = [ v1 ,  , vr ] are unique orthogonal matrices composed of the left and right singular vectors ui and vi , respectively .
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=1ruitcivi , 
							 see again .
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 = arg  min  kn(||akc||22 ) . 
 
however , the least - square solution kls 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 .
as an example , figure 7 shows the solution kls of the example data from figure 6 .
it is obvious that this solution has nothing in common with the real physiological behavior of the flow - scaled residue function . in order to get a better understanding of why kls is not a meaningful solution and to motivate the regularization approach
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 ( uic)/i that occur in ( 40 ) .
these factors weight the right singular vectors vi of a. 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 uic 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   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 .
consequently , the absolute values of the ratios ( uic)/i  which represent the weighting factors of the right singular vectors vi  become very large . to obtain a numerically stable result ,
the filter should suppress the influences of small singular values i or , equivalently , the influences of high absolute values of the weighting factors ( uic)/i . the regularized solution k , where  is a regularization parameter , is given by 
 ( 42)k=i=1r(f,iuitci)vi . 
 
firstly , the filter factors f,i 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 are based on the tikhonov regularization approach and characterized by a smooth weighting function centered around  , 
 ( 44)f,i  ( tikh)=i2i2+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 , figure 9 shows a plot of the function f = /( +  ) which is  unlike ( 44)defined for a continuous range of . for determining f , we assumed 1 = 1 . it can be seen that , for increasing  ( i.e. , stronger regularization ) , the values of f decrease for all . interestingly , the solution k of ( 38 ) using the filter factors f,i is equivalent to minimizing the weighted sum of the squared euclidean residual norm ||ak  c||2 and the squared euclidean solution norm ||k||2 ; that is , 
 ( 46)k(tikh)=arg  min  kn(||akc||22+2||k||22 ) . 
 

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 ) , 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 .  
the algebraic deconvolution approach from section 4.2 is very commonly applied to analyze perfusion data . yet ,
deconvolution problems arise in many other applications , and numerous alternative algorithms to solve these problems have been developed . 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 , and it has also been evaluated to analyze perfusion data [ 45 , 5557 ] .
interestingly , the fourier transform - based deconvolution approach is mathematically equivalent to the svd - based deconvolution approach with a block - circulant matrix a , compare the appendix [ 47 , 5860 ] .
however , results obtained with svd - based and ft - based deconvolution can be different because the chosen regularization approaches for these two methods are usually not equivalent .
the regularization in the context of the fourier - based deconvolution approach can be implemented by means of a modified wiener filter , 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 .
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 .
however , if the underlaying residue function is different from the model the perfusion parameters may not be estimated correctly .
an iterative deconvolution algorithm based on maximum likelihood expectation maximization ( ml - em ) algorithm was proposed in .
the svd - based deconvolution approach , which is available in several software packages [ 6769 ] , is comparably simple to implement and can be considered as the current standard method in perfusion image analysis . 

figure 10(b ) demonstrated that ( the maximum of ) the solution k depends on the regularization parameter rel  .
consequently , the computed perfusion values  which can be derived from k according to table 2vary for different rel  . as an example , the cbf value will be underestimated systematically for large rel  .
this approach is often used in practice , and a typical value in brain perfusion ct is , for example , rel   = 0.2 .
however , there exist more sophisticated approaches as well to determine the values rel   independently for each voxel position .
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||2 of the solution versus the squared euclidean norm ||ak
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 .
furthermore , a parameter estimation method that uses a priori knowledge of the behavior of the residue function was proposed in .
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 .
this section gives an overview of pre - processing techniques that can be applied in order to enhance the quality of the estimated perfusion parameters .
pre - processing occurs prior to the deconvolution step which may be implemented as described in section 4.2 . a simple ,
yet mandatory , pre - processing step consists of the conversion to contrast agent concentration values , see section 4.5.4 . further pre - processing steps are used to enhance the image quality ( e.g. , noise reduction ) and to correct for artifacts ( e.g. , motion correction , partial volume correction ) and specific properties of the blood ( e.g. , correction of differences in hematocrit ) .
there are also pre - processing steps that can optimize the analysis of the perfusion value maps ( e.g. , segmentation of certain anatomic structures ) and the application workflow ( e.g. , automated aif estimation ) .
the order of the pre - processing steps presented in this section can act as a guideline for their practical implementation .
finally , this overview can not include all details regarding suitable pre - processing steps .
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 .
a 3d registration should be used because it can also correct motion that occurs perpendicular to the orientation of the reconstructed slices .
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 .
these approaches may be used for perfusion mr data as well since both types of data typically consist of t2 * -weighted   epi images . however , the dynamic signal changes are relatively higher in dsc - mr data when compared to fmri data 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 can not 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 .    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 .
both linear and nonlinear filtering in the temporal dimension  that is , between successive time frames  represent further methods for noise reduction . 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 .   
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 .
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 .    neither for the case of ct imaging nor for the case of mr imaging can the time - concentration curves cart(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 ,
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+b1)0),0=1b  i=1b  (ti ) . 
  in perfusion
mr , however , the contrast agent concentration value is not proportional to the received signal s(tj ) ( in one voxel ) .
instead , it can be determined using the following formula : 
 ( 48)c(tj)=kmrteln  ( s(tj+b1)s0),s0=1b  i=1b  s(ti ) , 
									 see . 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 ] .
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 .
note that if only one time frame is considered as the baseline ( i.e. , if b = 1 ) , then c(0 ) = 0 , and the matrix a defined by ( 37 ) and ( 38 ) will be rank deficient ,   compare section 4.2 .
hematocrit ( hct ) is a value that describes the proportion of the blood that consists of red blood cells .
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 .   
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 .
an interesting alternative approach is to estimate several local aifs , which would be better suited to the theoretical model that was introduced in section 3 [ 9496 ] .
for example , partial volume effects  compare section 4.5.7can be more severe when compared to choosing one global aif in a larger vessel .
perfusion analysis using local aifs is actually investigated in   and the authors state that it produced more useful cbf maps . besides the arterial input function cart(t )
knowledge about the venous outflow function could be used to automatically correct for partial volume effects which are described next .    due to limited spatial resolution in reconstructed perfusion ct and mr data ,
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 ,
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 .
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 .
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 .
 deconvolution - based analysis of ct and mr brain perfusion data is 
 widely used in clinical practice and it is still a topic of ongoing research activities . in this paper , we present a comprehensive derivation and explanation of the underlying physiological model for intravascular tracer systems . 
 we also discuss practical details that are needed to properly implement algorithms for perfusion analysis . 
 our description of the practical computer implementation is focused on the most frequently employed algebraic deconvolution methods based on the singular value decomposition . in particular , we further discuss the need for regularization in order to obtain physiologically reasonable results . 
 we include an overview of relevant preprocessing steps and provide numerous references to the literature . 
 we cover both ct and mr brain perfusion imaging in this paper because they share many common aspects . 
 the combination of both the theoretical as well as the practical aspects of perfusion analysis explicitly emphasizes the simplifications to the underlying physiological model that are necessary in order to apply it to measured data acquired with current ct and mr 
 scanners . 
1. Introduction
2. Clinical Applications of Perfusion Imaging
3. Theoretical Model
4. Practical Implementation
5. Summary