WO2005114229A2 - Perfusion weighted mri with local arterial input functions - Google Patents

Abstract

A method for perfusion weighted imaging in which a plurality of local arterial input functions is estimated at each of a plurality of voxels. At least in part of the basis of the local arterial input functions, a cerebral blood flow at a voxel associated with one of the local arterial input functions is estimated.

Description

PERFUSION WEIGHTED MRI WITH LOCAL

ARTERIAL INPUT FUNCTIONS

FIELD OF INVENTION

The invention relates to magnetic resonance imaging ("MRI"), and in particular, to perfusion weighted MRI. RELATED APPLICATIONS

This application claims the benefit of priority of U.S. Provisional Patent Application Serial No. 60/571,621, filed on May 14, 2004, the contents of which are incorporated herein by reference in their entirety.

BACKGROUND Perfusion- weighted magnetic resonance imaging is a common imaging technique used in the clinical treatment of patients with brain pathologies such as stroke or cancer. Perfusion-weiglited images are obtained by injecting a bolus of gadolinium chelate into a patient's bloodstream and imaging as it passes through the brain. The gadolinium acts as contrast dye due to its T2 and T2* effects, which cause a drop in transverse relaxation time. This signal drop can then be used to calculate the concentration of the dye in a given voxel of brain tissue over time. The resulting concentration-time curves are then used with standard tracer kinetic models to calculate perfusion metrics such as blood volume, blood flow, and mean transit time.

However, solving the tracer kinetic model equations to calculate blood flow requires the deconvolution of an arterial input function from the measured concentration-time curves. Since the arterial input function is not known explicitly, it must be estimated from measured data. In current practice, a trained specialist examines the measured data and selects a single estimate of the arterial input function. This single estimate is used for the entire brain.

The current practice of estimating an arterial input function relies on the assumptions that the contrast agent reaches all parts of the brain at nearly the same time and that the contrast agent does not disperse significantly on its path from the major arteries to the brain tissue. In many cases, these assumptions are incorrect. However, even if the assumptions were correct, the time required to manually select the arterial input function can make the current practice inconvenient or impractical. In an emergency, the extra time spent identifying a suitable arterial input function could spell the difference between saving and losing brain tissue.

SUMMARY

The invention is based on the recognition that a more accurate representation of cerebral blood flow can be obtained by replacing the global arterial input function with a spatially-variable local arterial input function.

In one aspect, the invention includes a method for perfusion weighted imaging in which a plurality of local arterial input functions is estimated at each of a plurality of voxels. At least in part of the basis of the local arterial input functions, a cerebral blood flow at a voxel associated with one of the local arterial input functions is estimated.

In one embodiment, estimating a plurality of local arterial input functions includes selecting a target voxel; defining a search neighborhood corresponding to the target voxel, the search neighborhood including a plurality of neighborhood voxels; and estimating a local arterial input function for the target voxel at least in part on the basis of measurements associated with each of the neighborhood voxels.

The search neighborhood can be defined as a search neighborhood centered on the target voxel. Alternatively, the search neighborhood can be defined as a search cube centered on the target voxel.

In another embodiment of the invention, estimating a local arterial input function includes evaluating a selected property of a concentration measurement associated with a neighborhood voxel; and determining that the selected property satisfies a criterion. The local arterial input function is then estimated at least in part on the basis of that concentration measurement.

Exemplary properties of the concentration measurement that can be evaluated include a peak amplitude, a full-width half maximum, a first moment, and a slope of a line extending between a baseline and the peak amplitude. In one embodiment of the invention, the inclusion of only those measurements for which the selected properties satisfy a criterion includes assigning, on the basis of values associated with the selected properties, an overall score to each of the neighborhood voxels; ranking the neighborhood voxels on the basis of their corresponding overall scores; and including only a selected number of neighborhood voxels, the selected neighborhood voxels being selected on the basis of their respective ranks.

Additional embodiments of the invention include those in which the included voxels are weighted on the basis of their respective distances to the target voxel; and a local arterial input function is estimated at least in part on the basis of a weighted combination of the included voxels.

In another aspect, this invention features an MRI ("Magnetic Resonance Imaging") system for perfusion weighted imaging. The system include an MRI scanner for collecting MRI data at each of plurality of voxels and a data processing system for controlling the scanner. The data processing system is configured to estimate a plurality of local arterial input functions at each voxel in the plurality of voxels; and, at least in part of the basis of the local arterial input functions, to estimate a cerebral blood flow at a voxel associated with one of the estimated local arterial input functions.

In some embodiments, the data processing system is configured to select a target voxel and to define a search neighborhood that includes neighborhood voxels and, that corresponds to the target voxel, and to estimate a local arterial input function for the target voxel at least in part on the basis of measurements associated with each of the neighborhood voxels. The search neighborhood can be centered on the target voxel. For example, the search neighborhood can be a search cube centered on the target voxel.

Other embodiments include those in which the data processing system is configured to evaluate a selected property of a concentration measurement associated with a neighborhood voxel, to determine that the selected property satisfies a criterion, and to estimate the local arterial input function at least in part or the concentration measurement. In yet another embodiments, the data processing system is configured to select the property of the concentration measurement from to be a peak amplitude, a full- width half maximum, a first moment, or a slope of a line extending between a baseline and the peak amplitude.

Additional embodiments include those in which the data processing system is configured to estimate the local arterial input function by assigning, on the basis of values associated with the selected properties, an overall score to each of the neighborhood voxels; ranking the neighborhood voxels on the basis of their corresponding overall scores; and including only a selected number of neighborhood voxels, the selected neighborhood voxels being selected on the basis of their respective ranks.

In some embodiments, the data processing system is configured to weight the included voxels on the basis of their respective distances to the target voxel; and to assign a local arterial input function at least in part on the basis of a weighted combination of the included voxels.

Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. Although methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present invention, suitable methods and materials are described below. All publications, patent applications, patents, and other references mentioned herein are incorporated by reference in their entirety. In case of conflict, the present specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and not intended to be limiting.

Other features and advantages of the invention will be apparent from the following detailed description, and from the claims.

BRIEF DISCRIPTION OF THE FIGURES

FIG. 1 is a typical arterial input function. FIG. 2 is a flowchart of a method for estimating local arterial input functions.

DETAILED DISCRIPTION

To monitor cerebral blood flow in a subject, one injects a contrast dye (e.g. gadolinium chelate) into the subject. This contrast dye is selected for its conspicuousness in a magnetic resonance image of the brain.

The contrast dye enters the brain through a number of arteries. The cerebral blood flow eventually distributes the dye throughout the brain. By using magnetic resonance imagining to observe spatial and temporal evolution of dye concentration throughout the brain, one can infer the cerebral blood flow that drives the distribution of the dye.

At any voxel in the brain, the concentration of dye entering that voxel, shown in FIG. 1, takes the form of a time-varying function having a rapid rise to an initial narrow peak D, followed by an extended dilution period. The extended dilution period is followed by a second, somewhat broader peak F, which corresponds to the re-circulation of blood containing residual contrast dye back into the voxel. The curve shown in FIG. 1 is referred to as an "arterial input function."

The concentration of dye in a particular voxel depends on the temporal convolution of the arterial input function at that voxel and a residency function, also evaluated at that voxel. As noted above, the arterial input function describes the rate at which the contrast dye enters the voxel. The residency function describes how long the contrast dye remains in the voxel before dissipating outward. A typical residency function has an initial value of unity that diminishes with time. This is consistent with the notion that a contrast dye in a voxel will gradually diffuse outward in response to a concentration gradient.

When multiplied by the cerebral flow rate at a particular voxel, the convolution of the arterial input function at that voxel and the residency function yields the measured concentration of contrast dye at that voxel. Thus, to determine the unknown cerebral flow rate from the measured dye concentration, one must deconvolve the arterial input function from the residency function. This requires knowledge of the arterial input function at that voxel, i.e., the "local" arterial input function.

Referring to FIG. 2, a method for estimating a local arterial input function can be divided into three stages: a preparatory stage 10, a searching stage 12, and a deconvolution stage 14.

The preparatory stage 10 begins by converting the signal intensity curves into concentration-time curves (steps 16, 18). In doing so, a linear relationship between concentration and the change in transverse relaxation is assumed. A baseline intensity for each voxel is obtained by first finding an average peak time for the whole brain by averaging the time-to-peak for every voxel in the brain. Following that, a baseline intensity is calculated on a voxel-by- voxel basis to be the average of the intensities from time zero to a time that is twelve seconds before the average peak time. The concentration-time curves are then calculated. For each voxel, the cerebral blood volume for that voxel is calculated as the area under the concentration curve for that voxel (step 20).

The next step is to exclude from further consideration those voxels that are clearly unfit for arterial input function selection (step 22). This is done by first excluding voxels having a cerebral blood volume that is either too low (i.e., less than 5% of the maximum cerebral blood volume) or too high (i.e., greater than 60% of the maximum cerebral blood volume). Following this, the first moment of the concentration curve for each voxel is calculated (step 24) and a 30 second arrival window centered on the average peak time is defined. The remaining voxels are inspected to exclude voxels corrupted by noise, vessel pulsation, patient motion, susceptibility artifacts or other sources of artificial signal drops. In particular, a voxel is excluded if the measured concentration of dye entering that voxel has one or more of the following properties:

• A first moment that is more than 7.5 seconds before the average peak time. (Because having a majority of the concentration before the expected arrival of the bolus of dye is indicative of noise.) • A maximum value outside of the arrival window. (Because having maximum values well before or after the expected bolus arrival time is indicative of artificial signal drops.)

• Two points at 60% of the maximum value separated by more than 26 seconds. (Because a 26 second interval has been found to be long enough for bolus passage.)

• A "negative concentration" with an absolute value more than half of its maximum. (Since noise is usually distributed evenly around zero, this will exclude noisy voxels.)

• An occurrence, before the arrival window, of a negative concentration change that is greater than 40% of the maximum concentration change. (A large negative concentration change is expected only after the bolus arrives.)

• An occurrence, after the arrival window, of a positive concentration change that is greater than 60% of the maximum concentration charge. (A large positive concentration change is expected only when the bolus arrives.)

• A mean concentration after the arrival window that is less than two standard deviations above the mean concentration before the arrival window. (Valid arterial input functions have a re-circulation artifact that will increase the post-bolus mean concentration.) For each of the remaining voxels (hereafter referred as "target" voxels), one then defines a search cube of neighboring voxels centered on a target voxel (step 26). The measured concentration function at each of the neighboring voxels is then evaluated on the basis of how likely it is that that concentration function corresponds to an arterial input function.

In particular, a concentration function for a neighboring voxel is assigned a score on the basis of certain properties. Referring back to FIG. 1 , these properties include the amplitude of the first peak D, the first moment of the concentration function, the width of the peak E, and the slope C of a line connecting the peak to the beginning of the rise in concentration B. The four properties noted above are then combined into a single scalar score associated with that neighboring voxel (step 28). The three highest-scoring neighboring voxels are then identified.

Each of these four properties is normalized to fall between 0 and 1. This is done by subtracting out the minimum value of each property in the search cube followed by dividing by the maximum value of each property in the search cube. The score "5" for each of the neighboring voxels in the search cube is given by

S = FM+FWHM-PV-AS where FM is the first moment, FWHM (which stands for "full width half max") is the width of the first peak at the point at which it has reached its maximum value, PV is the peak value, and AS is the slope of the line extending between the peak value and the baseline of the peak. A low value of E is indicative of minimal delay. A low value of FWHM, together with high values of PV and AS is indicative of minimal dispersion. The score S falls between -2 and 2, with the best voxels having the lowest scores.

Referring back to FIG. 2, the three lowest-scoring voxels are then weighted on the basis of their location relative to the target voxel (step 30). In particular, the closer a neighboring voxel is to the target voxel, the higher will be the weight associated with that voxel. This is done by interpolating the three best voxels with a 3-D 27mm FWHM Gaussian kernel. This result is stored as the arterial input function for the target voxel centered in the search cube. The foregoing procedure is then repeated for additional target voxels (step 32) until all voxels in the brain have been processed as described above.

To prevent a bias on the edges of the brain, the arterial input functions are extended out by 14mm in all directions at the brain periphery, thereby avoiding edge effects during the smoothing process (step 34). This is accomplished by adding extra slices above and below the brain, as well as by using morphological operators to replicate the values on the edge of the brain out an extra 14mm. The final deconvolution stage begins by smoothing the local arterial input functions (for continuity) with the same 3-D Gaussian kernel used for interpolation (step 36). Finally, the local arterial input functions are de- convolved from the concentration curves (step 38). The deconvolution is carried out by singular value decomposition on a voxel-by- voxel basis.

Exemplary code for implementing the foregoing method using MATLAB (TM) is given in Appendix A.

It is to be understood that while the invention has been described in conjunction with the detailed description thereof, the foregoing description is intended to illustrate and not limit the scope of the invention, which is defined by the scope of the appended claims. Other aspects, advantages, and modifications are within the scope of the following claims.

Having described the invention, and a preferred embodiment thereof, what is claimed as new, and secured by letters patent is:

Appendix A:

Local AIF Algorithm Code

Main Function function [AIF, cbv, CBF,MTT] =

LocalAIFvl_3 (inDir,outDir, Prefix,Nslice, Postfix, Type, TR,TE, thick,gap, FOV x,FOVy, smooth, exclude, post_ex,kk) ;

% LocalAIF_vl_3 Computes a local AIF, CBV, CBF, and MTT given a PWI data set

% [AIF, cbv, CBF,MTT] = LocalAIF_vl_3 (inDir, outDir, Prefix,Nslice, Postfix,

%

Type, TR, TE,Thick, Gap, FOVx, FOVy, Smooth, Exclude, Post_ex,kk)

% Where:

% "inDir" is the directory where the input data is located.

% "outDir" is the directory where the output data will be stored.

% "Prefix" is the prefix before the dot, without slice numbers.

% "Nslice" is the number of slices.

% "Postfix" is after and including the dot.

% "Type" specifies whether the data is in one big chunk or in slices.

% For instance, if you had inOO.bshort ... inlO.bshort as input,

% you would have 'in' as prefix, '11' as nslice, '.bshort¹ as postfix, % and 0' as type. % However, if you just have input.bshort in 12 sequential slices

% and time interleaved, you would have 'input' as the prefix, ^λ 12' as

% nslice, '.bshort' as postfix, and *1' as type. % "TR" is the repition time (time between samples)

% "TE" is the echo time.

% "Thick" is the slice thickness, in mm.

% "Gap" is the spacing between the slices, in mm.

% "FOVx" is the Field of View in the x dimension (in cm) . % "FOVy" is the Field of View in the y dimension (in cm) .

% "Smooth = 0" means do not smooth input signal in time.

% "Exclude" is the number of initial time points to exclude.

% "Post_ex" is the number of end points to exclude.

% "AIF" is the locally defined AIF, (not written) . % "CBV" is the Cerebral Blood Volume (written to OutDir) .

% "CBF" is the Cerebral Blood Flow (written to OutDir) .

% "MTT" is the Mean Transit Time (written to OutDir) .

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Start

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Read in the Dataset in one of two ways . %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cd(inDir); if (Type == 0) i = 0;

[temp, nr, nc, nsamp, e] = read_mgh_tb (strca (Prefix, ' 00 ' , Postfix) ) ; Sig = zeros (nr, nr, Nslice, nsamp-exclude-post_ex) ; Sig ( : , : , i+l, : ) = temp ( : , : , 1+exclude :nsamp-post_ex) ; for i = l :Nslice-l , if (i <10) num = strcat ( ' 0 ' ,num2str(i) ) ; else num = num2str(i) ; end [temp , nr, nc , nsamp , e] = read_mgh_tb ( strcat (Prefix, num,

Postf ix) ) ;

Sig ( : , : , i+1 , : ) = temp ( : , : , 1+exclude : nsamp-post_ex) ; end else [Sig , nr, nc , total , e] = read_mgh_tb ( strcat ( Prefix, Postfix) ) ;

[Sig , nsamp] = temporal ( Sig , Nslice ) ; Sig = Sig ( : , : , : , 1+exclude : nsamp-post_ex) ; end clear temp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Extract out the brain and calculate certain statistics for later use %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [nr,nc,nsl,nsamp] = size(Sig);

%% Noise Threshold in STDs kk=8;

% We gain information on the noise by looking at the corner of the set mMean = mean(Sig (nr-floor (nr/2) ,nc-5,nsl-l, :)) ; mSTD = std(Sig (nr-floor (nr/2) ,nc-5,5, : ) ) ;

%Find the AIFs Using an approximately 28x28x28 mm^Λ3 search, ni = floor ( (nr/FOVx)*2.8/2) *2 + 1; nj = floor ( (nc/FOVy) *2.8/2) *2 + 1; nk = floor (12/ (thick+gap) )*2+l; ext = max(ni,nj);

% Initializes a simple mask to sort out the noise. pre_mask = zeros (nr,nc,nsl) ; pre_maskp = pre_mask; extra_mask = zeros (nr,nc,nsl) ; warning off; for i = l:nsl,

% Look at the first image to determine brain or not. temp = Sig( : , : , i, 1) ;

% This is a simple threshold based on user input.

Q = find(temp <= mMean + kk*mSTD) ;

% Now, morphological operators are used to extract brain temp2 = ones(nr,nc); temp2 (Q) = 0; temp2 = erode (temp2, ones (3,3)) temp2 = erode(temp2, ones(3,3)) temp2 = erode (temp2, ones (3, 3)) temp2 = erode (temρ2, ones (3, 3)) temp2 = erode (temp2, ones (3, 3)) temp2 = erode (temp2, ones (3,3)) tetnp2 = dilate (temp2, ones (3,3)) temp2 = dilate (temp2, ones (3, 3)) temp2 = dilate (temp2, ones (3, 3)) temp2 = dilate (temp2, ones (3, 3)) temp2 = dilate (temp2, ones (3, 3)) temp3 = dilate(temp2,ones(ext,ext) ) ; pre_mask ( : , : , i) = temp2 ; shell (:,:,i) = temp2 - erode (temp2,ones (3,3) ) ,- extra_mask ( : , : , i) = temp3-temp2 ; end

% This is then an indication of what is brain and what isn't. Relevant = find(pre_mask) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Smooth the input data in space %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% scon = 3*22/128; dimpri = ma (nr, nc) /ma (FOV , FOVy) *scon; sdim = 2*floor (dimpri/2) + 1; for i=l:nsamp, Sig(:,:,:,i) = Gsmoot (Sig (:,:,:, i) , sdim, .8*floor (dimpri/2) ) ; end

% Find where largest signal drop is, but only in brain. [Y, I] = min(Sig, [] ,4); Relevant_I = I (Relevant);

%This calculates average arrival time, (in samples) arrive_time = mean (Relevant_I ( : ) ) ; display (arrive_time) ; clear Relevant_I Relevant;

%% Calcualte baseline to be twelve seconds before the arrival time Base = mean (Sig ( : , : , : , 1 :round(arrive_time-12/TR) ) ,4) ;

% Calculate the Concentration and subject it to the masking.

% Replicate Baseline

Baseline = padarray (Base, [0 0 0 nsamp-1], 'replicate', 'pre');

% Make sure not to divide by zero Q = find(IMG_in==0) ;

Q2 = find(Baseline==0) ;

Q3 = union (Q,Q2) ;

IMG_in(Q3)=l;

Baseline (Q3) = 1; % Calculate concentration

Cone = -log(Sig. /Baseline) ; clear Sig, Baseline;

%% Discard extraneous non-brain voxels mask = padarray (pre_mask, [0 0 0 nsamp-1], 'replicate', 'pre'); Cone = Conc./TE .* mask; clear Base;

% Calculate CBV. cbv = trapz (Cone, 4) ; Q = find (cbv <0) ; cbv(Q) = 0; cbv = cbv . * mask (:,:,:, 1) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Smooth the input data in time for low SNR data sets %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if not (smooth == 0) temp_Conc = Cone; Conc2 = zeros (size (temp_Conc) ) ; Conc3 = zeros (size (temp_Conc) ) ; Concl = 2*temp_Conc/3;

Conc2 (:,:,: ,2 -.nsamp) = temp_Conc (:,:,:, 1 :nsamp-l) /6; Conc3 (:,:,:, 1 :nsamp-l) = temp_Conc (:,:,:, :nsamp) /6 ; temp_Conc = Concl + Conc2 + Conc3; Cone = temp_Conc ; clear Concl temp_Conc Conc2 Conc3; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q = find (Cone < 0) ; Conc(Q) = 0; • jui* (XSU'DU'JU) sauo = εuoτ-ta-iτJD

.τeaχo

^•rjuauiow sjτa sτ uoτja^τjco -isατ^ % Q

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

£9

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

• (Hl'

^•ai(τχa

pue χaxoΛ tpea 30 .SSΘU^IV_I 3t{4 -IOJ jpaq %

.'J_UI = (_ζO) 4tiauroui ςς •'(Hi/S'Δ - auιτ_f~ aΛτ ze > :}uauιouι)puτj = εo

(^• aΛτ-tJE snχoq atπ aero jag "3*T) % • oχ oo_t axe s-iuauιouι ΘSOUVΛ s-tuτod p-tsosfα %

.'(0 == 3uauιouι)puτj = 0 ■auoχB sχaxθΛ _upooB„ BAEBX pue 'JUI Λrjτu juy. 03 sχaxθΛ uP^ ,, ^S %

.'JUI = (θ)-iuauιouι £f

^■aβuBJt aumχoΛ aq^ uτ τj ?oπ op _n_tπ s^uatuoui asoτπ - τuτjuτ 03 ^as %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%% sχaxoΛ peq pjB sτp 03 si au 30 saτ;tas qBnottn 03 %%

ζ£

.'douoo jteaχo

•' (ϊ-'douoo) z -eα-i = rπiauioui

-luauiow -isjτa stR

%

.' (XSU'OU'-tU) SOI9Z = 3UΘIUOU1

.' ζβταo Λqo βτ:co Λqo siunα αeaχo oε .' sumu* ^■ douoo = douoo

.'i[SBUi αeaχo

^■S.3IV jo uoτ^oaχas aip αo sχoΛ «oχ puε u;^βτt{ aτj-t auτ ap am ajaH %

-' (

Λ90. )Λeχdsτp

.' (ouoo)ΛSD = Bτj:o-Λqo g t^?99ΪO/£OOZSfl/I3d 6ZZHΪ/S00Z OΛV criterion3 = ones (nr,nc, nsl) *Inf; criterion4 = ones (nr,nc, nsl) *Inf;

% Third Criterion is highest peak. [Y,I] = max(Conc, [] ,4) ; criterion3 = Y; criterion3 (Q) = Inf ;

% Second Criterion is Average slope, run = zeros (nr,nc, nsl) ,- rise = zeros (nr,nc, sl) ; for i=l:nr*nc*nsl, if not (criterionl (i) == Inf) if (Hi) < round(18/TR)) before = I (i) ; else before = round(18/TR) ; end

% This computes the average Slope.

1 = ( (I (i) -before) *nr*nc*nsl:nr*nc*nsl : (I (i) -1) *nr*nc*nsl) look = squeeze (Cone (i+1) ) ;

Qp = finddook < .05); if not (isempty(Qp) )

[Y2,I2] = max(Qp) ; run(i) = before-Y2; rise(i) = Y(i) -Cone (i+nr*nc*nsl* (I (i) -run(i) ) ) ; % Fourth Criterion is the FWHM.

1 = (0 :nr*nc*nsl : (nsamp-1) *nr*nc*nsl) ; criterion4 (i) = fwhm(squeeze (Cone (i+1) )) ; end end end warning off ; criterion2 = rise. /run; warning on; criterion2 (Q) = Inf; criterion4 (Q) = Inf; clear rise run; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%

%% Create the matrix that will accomplish the interpolation and filtering.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% sigma = sqrt(.5 * ( (ni-1) /2) ^Λ2/log (2) ) ;

G2 = fspecial ( 'gaussian' , [ni,nj] , sigma) ;

% This also gives a half max at 28mm. alph = (nk/2-l/2)^Λ2/log(2) ;

Gl = 1:1 :nk; Gl = Gl- (nk/2+1/2) Gl = exp(-Gl.^A2/alph) ;

G2p = padarray(G2, [0 0 nk-1] , 'replicate', 'pre'); % Replicates G2 in 3rd Dim.

G3 = zeros (ni,nj ,nk) ; for i = l:ni, for j = l:nj ,

G3(i,j,:) = Gl' .* squeeze (G2 (i, j, :)) _; end end

G3 = G3*1000;

G4 = padarray(G3, [0 0 0 nsamp-1], 'replicate¹, 'pre'); ^•'(I)frD *- (^: 'i{I+^'Cl+C'TI+T)^ouoO = x^χ •-Ij -mBτa ue ssnBO pus saτjtas auιτ_ι s τ 533 %%

OL

^■ axoΛ S9 punoj au; jo a_iBuτpj:ooo ^o xa aq^ a_i χn %%

3sχa

•τ = Λ-tp -'i = εi -τ =

^•βuτq_tauιos puτj BM aαns

jτ 09

.' JUI = ( 1 ) aqno qojtBas •'((:) aqno- u jteas ) uτui = [I'A] • (Hi'ϊ_'^τα 'ε^τα 'ε-iτj 'i-iτα 'aqno^'- uoo) εs^c s = aqno- qoαeas ςς

^• apnχoxa uaqq puB 'xapuτ :}aβ 'a_ιoos umuiτuτui mjiΛ XΘXOΛ puτ,£[ % pua (: 'sttj+jfi-fn-jj' Cq+C: Cq-C' q+τ:τq-τ) ouoo = aqno- uoo oς (3[t[+3[:5[t[-5[' q+C: pq-C'τq+τ:τq-τ)

(5(u;+i(:i(t[-i(' Cq+C: q-C 'τti+τiτn-τ) εuoταa-iτα = ε^ϊ-^ (i[q+5[:i(t(-i[' Cq+C: Cq-C 'τq+τi q-τ) suoτja_ι j =

( q+5[: j[u;-i[' Ct+C: Cq-C'τq+τ:τu;-τ) xuoταa^τα = 1-iτ.to ζ

•0 = >1PPB •saoτχs axppτw % asxa

(: 'is r-fq - ' Cq+ C'τq+τ:τq-τ) uoo = aqno- uoo (χstι:3[q-3t τq+τ:τq-τ) f.uoταa^τ.t = f xo 017 (χsu:.fq-3[ Pn-f τq+τ : τq-τ) εuoταa^τ = ^ αo

(ISU:3pi-- ' t+C τq+τ:τq-τ) zuoταa-iταD = ε3 _to (χsu:i(q-5[ τq+τcτq-τ) lucc-ta^f-t = τ^τxo .'0 = yppv aoτχs -ISBT % xsu < 5[q+5[ j"fasχa

(: '5[iι+3t:χ'Ctl+ C'τq o- uoo sε '3[t[+3[:i' fq+ : Cq-C τq+τ

ϊ'-tτao 'i[t[+i[:χ ' Cq+p: τq+τ:τq-τ) εuoταa-tτα = ε^T-³⁰ 'i[q+ :χ' Cq+C: fq-f τq+τ : τq-τ) guoτja-iτj = S-iτj 'j[u;+3[:T' Cq+C: fH-f τq+τ :τq-τ) xuo αa-iτj = T. _Ϊ-30

.'t+3[-3[q = ifPPB oε aoτχs

% τ > tq- JT q JBas iχnj B op .[outiBO qoτqΛ SΘUXBΛ auiacrpca aq .toj ^sa-i asaqj, % •uτ αq uτ 30U sχaxoΛ JO a-)BχnoχB ,uocι % X == (-(' C ' )i{SBUι^— a:πd jτ sz

'χsu:χ = ^■% JtOJ ' .' (ε/Cu) j:ooχj-ou: (ζ/Cu) χτao = C αo

' (2/τu) j:ooχ -_ιu: (s/ u) χτao = τ αoj dooi q jBas afiJBq; %% oz

•' (ιsu'ou';tu) soαaz = dswT. barj •' (Hi,'ϊ_'Uoταa-(τΛ ' uoτja-i α 'suo j:a-i -tD'ιuo -ιa_ιτα ) ^aαoos = aqno- axoos%

^• sajtoos O-iuτ ταa-iταo

%

.' (duiBsu'χsu'ou'α

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 01

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

•ε = Λτp

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t99ΪO/£OOZSΛ/13d 6ZZMΪ/S00Z OΛV freq_map (i+Ii, j+Ij ,k+Ik) = freq_map(i+Ii, j+Ij , k+Ik) +1;

% Find Next lowest score value, get its index and then exclude. [Y,I2] = min(search_cube( .- ) ) ; search_cube (12) = Inf; if (Y == Inf) % make sure we find something. x2 = 0; x3 = 0; 13 = 1; div = 1; else

%% Calculate the coordinate of the found voxel . Ip = 12-1;

Ik = floor (Ip/ (ni*nj) ) -hk+addk; I2p = Ip- (Ik+hk-addk)*ni*nj;

Ij = floor (I2p/ni) -hj ; Ii = I2p- (Ij+hj ) *ni-hi;

%% Get its time series and Gaussian weight it. x2 = Conc(i+Ii, j+Ij,k+Ik, : ) .* G4(I2); freq_map (i+Ii, j+Ij ,k+Ik) = freq_ma (i+Ii, j+Ij ,k+Ik) +1;

% Find last lowest score value, and its index. [Y,I3] = min(search_cube ( : ) ) ; if (Y == Inf) x3 = 0; div = 2; else

%% Calculate the coordinates Of voxel. Ip = 13-1; Ik = floor (Ip/ (ni*nj) ) -hk+addk;

I2p = Ip- (Ik+hk-addk) *ni*nj ; Ij = floor (I2p/ni) -hj ; Ii = I2p- (Ij+hj ) *ni- hi;

%% Get its time series and Gaussian weight it. x3 = Conc(i+li, j+Ij,k+Ik, :) .* G4(I3); freq_map(i+Ii, j+Ij ,k+Ik) = ... freq_map (i+Ii, j+Ij ,k+Ik)+l; end end

%% Sum up time series, normalize and store AIF_int(i, j,k, :) = (xl+x2+x3) / (G4 (I) +G4 (12) +G4 (13) ) ; end end end end end clear G4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%

%% Extend the AIF out to avoid tailing off at the edges. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AIF_ext = zeros (nr,nc,nsl,nsamp) ; extp = (ext+l)/2;

%% Shell defines the outside shell of the brain, we extend this in time shelled = padarray (shel1, [0 0 0 nsamp-1], 'replicate', 'pre'); %% Get this outside shell from AIF above, shelled = shelled .* AIF_int; dime = 2*extp+l; for i = l+extp:nr-extp, for j = l+extp:nc-extp, for k = 1 :nsl, if extra_mask (i, j ,k) ==1 %% Find all elements in shell within a certain box around the voxel

Q = find (shelled (i-extp: i+extp, ... j-extp: j+extp,k, floor (arrive_time) ) ) ; num = length (Q) ; for h = l:num,

%% Find where these elements exist Qp = Q(h)-1; Ij = floor(Qp/dime) -extp; Ii = Qp- (Ij+extp) *dime- extp;

%% Take their average as the resulting edge AIF AIF_ext (i, j ,k, : ) = AIF_ext (i, j ,k, : ) + ... shelled(i+Ii, j+Ij ,k, : ) /num; end end end end end

%% Add in actual AIF to get final AIF. AIF_ext = AIF_ext + AIF_int; clear AIF_int; %% Extend AIF to an extra slice above and an extra slice below to avoid a bias on the slices.

AIF_int3 = zeros (nr,nc,nsl+2*hk,nsamp) ;

AIF_int3 ( : , : , 1+hk:nsl+hk, : ) = AIF_ext ;

AIF_int3 (: , :,l:hk, :) = ... padarray (AlF_ext(:, :,1, :) , [0 0 hk-1 0] ,' replicate ', 'pre ') ;

AIF_int3 (: , : ,nsl+hk+l:nsl+2*hk, :) = ... padarray (AIF_ext (: , : ,nsl, :) , [0 0 hk-1 0] ,' replicate ', 'pre' ) ; clear AIF_ext; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Spatial filtering of the AIF. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AIF_intf = zeros (nr,nc,nsl+2*hk,nsamp) ; for i = 1:nsamp,

AIF_intf (:,:,: ,i) = imfilter (AIF_int3 (:,:,:, i) , G3); end

AIF = AIF_intf ( : , : , l+hk:nsl+hk, : ) ; clear AIF_int3 AIF_intf ;

AIF = AIF .* padarray(pre_mask, [0 0 0 nsamp-1], 'replicate', 'pre'); display ( ' Filtering Done ' ) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Calculate the CBF. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %CBF_Ai = CBF_Deconv(Conc, AIF,pre_mask) ; Old Function call warning off;

CBF_Ai = zeros (nr,nc, nsl) ;

%% Get indices for diagonal elements.

W = eye (nsamp) ;

Q = find(W) ;

R = zeros (1, nsamp) ; for i=l:nr for j=l:nc, for k=l:nsl, if pre_mask(i, j , k) ==1 %% Don't calculate for non-brain temp_AIF = squeeze (AIF (i, j ,k, :)) ;

%% this linearizes the data.

AIF2 = zeros (nsamp, 1) ;

AIF3 = zeros (nsamp, 1) ;

AIF1 = 2*temp_AIF/3;

AIF2 (2 :nsamp) = temp_AIF (1 :nsamp-l) /6 ;

AIF3 (l:nsamp-l) = temp_AIF (2 :nsamp) /6 ; temp_AIF = AIF1 ₊ AIF2+AIF3;

% R is used to create the block-matrix for the SVD R(l) = temp_AIF(l) ; [U,S,V] = svd(toeplitz(temp_AIF,R) ) ,- % Q is defined earlier. W(Q) = l./S(Q);

%Find those diagonal entries that are less than some threshold

% and then set them to zero. Hi = max(S(Q) ) ; Ex = find(S(Q) < Hi*.2); W(Q(Ex) ) = 0;

%% Multiply everything out to get final deconvolution CBF_Ai(i, j ,k) = max(V*W * (transpose (U) * ... squeeze (Cone (i, j ,k, : ) ) ) ) ; end end end end warning on; ϊϊScScScϊK^MJiStSrMSc ns M^nscS^StSr&MScStS M

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Final Post-processing: calculate MTT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Remove Bad Values that can appear. Q = find(isnan(CBF_Ai) ) ; CBF_Ai (Q) = 0; Q = find(CBF_Ai == Inf); CBF_Ai (Q) = 0; Q = find (isnan (cbv) ) ; cbv(Q) = 0 ;

% Smooth CBF and calculate MTT. CBF_Ais = Gsmooth (CBF_Ai ,3,1.1); warning off; MTT = cbv./CBF_Ais; warning on; % Remove bad values and scale MTT. Q = find (isnan (MTT) ) ; MTT(Q) = 0 ; Q = find(MTT == Inf); MTT(Q) = 0; mMTT = max (MTT ( : ) ) ; MTT = 50 * MTT/mMTT;

%% Normalize CBF by its median. Q = find (CBF_Ai > 0) ; mF = median (CBF_Ai (Q) ) ; CBF = CBF_Ai . /mF;

%% Write Final Output cd (outDir) ; write_mgh_mcc ( 'UCBV.bfloat ' , cbv) ; write_mgh_mcc ( ' CBF.bfloat ' , CBF) ; write_mgh_mcc ( 'MTT.bfloat ' ,MTT) ;

display ( ' Finished ' ) ; Support Functions

Score3 function score_cube = ... score3 (conc_cube, eriterionl, criterion2 , criterion3 , criterion4,TR) % Score3 computes the score for a given set of 4 criteria.

% Specifically, it just normalizes then to be between 0 and 1.

%

% score_cube = score3 (conc_cube, eriterionl, criterion2, criterion3,criterion4,TR) ; % Where:

% "score_cube" is the score for every voxel.

% "eriterionl" is the first moment

% "criterion2" is the average slope

% "criterion3" is the highest peak % "criterion4" is the FWHM

% "TR" is the repetition time. warning off MATLAB : divideByZero; [ni,nj ,nk, nsamp] = size (conc_cube) ; nQ = find (eriterionl < Inf); score_cube = ones (ni,nj ,nk) *Inf ; if not (isempty (nQ) )

% Now we normalize all of these, mcl = min(eriterionl (:)) ; eriterionl = criterionl-mcl; mc2 = min(criterion2 ( : ) ) ; criterion2 = criterion2-mc2; mc3 = min(criterion3 ( : ) ) ; criterion3 = criterion3-mc3 ; mc4 = min(criterion ( : ) ) ; criterion4 = criterion4-mc4;

Mcl = max(eriterionl (nQ) ) ; if (Mcl == 0) eriterionl (nQ) = eriterionl (nQ) +1; else eriterionl = criterionl/Mcl; end

Mc2 = max(criterion2 (nQ) ) ; if (Mc2 == 0) criterion2 (nQ) = criterion2 (nQ) +1; else criterion2 = criterion2/Mc2 ; end Mc3 = max(criterion3 (nQ) ) ; if (Mc3 == 0) criterion3 (nQ) = criterion3 (nQ) +1; else criterion3 = criterion3/Mc3 ; end nQ = find(criterion4 < Inf) ; if isempty (nQ) score_cube (nQ) = eriterionl (nQ) - criterion2 (nQ) criterion3 (nQ) ; else

Mc4 = max(criterion4 (nQ) ) ; if (Mc4 == 0) criterion4 (nQ) = criterion4 (nQ) +1; else criterion4 = criterion4/Mc4 ; end score_cube = ones (ni,nj ,nk) *Inf ; score_cube (nQ) = eriterionl (nQ) - criterion2 (nQ) - ... criterion3 (nQ) + criterion4 (nQ) ;

Q = find(criterion4 > (15*1.5/TR) score__cube(Q) = Inf; end end warning on MATLAB : divideByZero;

Gsmooth function IMG = Gsmooth (IMG_in,N, sigma) ;

%Gsmooth Takes a 4-D MRI Dataset and smooths it with an NxN Gaussian kernel .

% [IMG] = Gsmooth (IMG_IN,N, SIGMA)

% Where "IMG_IN" is the 4-D Dataset, "N" is the dimension of the % Gaussian kernel, and "SIGMA" its std dev.

G = fspecial ( 'gaussian' , [N N] , sigma) ; [nr,nc,nsl,nsm] = size (IMG_in) ,- for i = l:nsl, for j=l:nsm,

IMG(:,:,i,j) = filter2 (G, IMG_in ( -. , : , i, j ) ) ; end end

FWHM function fwhmax = fwhm(data) ;

% FWHM computes the FWHM of a time series.

% fwhmax = fwhm(data) ;

% Where:

% "data" is the time series % "fwhmax" is the FWHM value mdata = max (data) ;

I = find (data > mdata/2) ; = length (I) ; if (L == 0) fwhmax = Inf ; elseif (I( ) == length(data) | 1(1) == 1) fwhmax = Inf ; else width = I(L) -1(1) ; bO = data (I (1) ) ; aO = data(Kl)-l) ; xO = (mdata/2-aO)*(l/(bO-aO) ) ; bl = data (I (L) ) ; al = data (I (L)+l) ; xl = (mdata/2-bl)*(l/(al-bl) ) fwhmax = width + xl-xO+1; end

Claims

1. A method for perfusion weighted imaging, the method comprising:

estimating a plurality of local arterial input functions at each of a plurality of voxels; and

at least in part of the basis of the local arterial input functions, estimating a cerebral blood flow at a voxel associated with one of the estimated local arterial input functions.

2. The method of claim 1, wherein estimating a plurality of local arterial input functions comprises:

selecting a target voxel;

defining a search neighborhood corresponding to the target voxel, the search neighborhood including a plurality of neighborhood voxels; and

estimating a local arterial input function for the target voxel at least in part on the basis of measurements associated with each of the neighborhood voxels.

3. The method of claim 2, wherein defining a search neighborhood comprises defining a search neighborhood centered on the target voxel.

4. The method of claim 2, wherein defining a search neighborhood comprises defining a search cube centered on the target voxel.

5. The method of claim 2, wherein estimating a local arterial input function comprises:

evaluating a selected property of a concentration measurement associated with a neighborhood voxel;

determining that the selected property satisfies a criterion; and estimating the local arterial input function at least in part on the basis of the concentration measurement.

6. The method of claim 5, further comprising selecting the property of the concentration measurement from the group consisting of: a peak amplitude, a full-width half maximum, a first moment, and a slope of a line extending between a baseline and the peak amplitude.

7. The method of claim 5, wherein estimating the local arterial input function comprises:

assigning, on the basis of values associated with the selected properties, an overall score to each of the neighborhood voxels;

ranking the neighborhood voxels on the basis of their corresponding overall scores; and

including only a selected number of neighborhood voxels, the selected neighborhood voxels being selected on the basis of their respective ranks.

8. The method of claim 7, further comprising:.

weighting the included voxels on the basis of their respective distances to the target voxel; and

estimating a local arterial input function at least in part on the basis of a weighted combination of the included voxels.

9. A computer-readable medium having encoded thereon instructions that, when executed, cause a computer to practice the method recited in claim 1.

10. An MRI system for perfusion weighted imaging, the system comprising:

an MRI scanner for collecting MRI data at each of plurality of voxels; and a data processing system for controlling the scanner, the data processing system being configured to

estimate a plurality of local arterial input functions at each voxel in the plurality of voxels; and

at least in part of the basis of the local arterial input functions, to estimate a cerebral blood flow at a voxel associated with one of the estimated local arterial input functions.

11. The system of claim 10, wherein the data processing system is configured:

to select a target voxel;

to define a search neighborhood corresponding to the target voxel, the search neighborhood including a plurality of neighborhood voxels; and

to estimate a local arterial input function for the target voxel at least in part on the basis of measurements associated with each of the neighborhood voxels.

12. The system of claim 11, wherein the data processing system is configured to define a search neighborhood centered on the target voxel.

13. The system of claim 11, wherein the data processing system is configured to define a search cube centered on the target voxel.

14. The system of claim 11, wherein the data processing system is configured:

to evaluate a selected property of a concentration measurement associated with a neighborhood voxel;

to determine that the selected property satisfies a criterion; and

to estimate the local arterial input function at least in part on the basis of the concentration measurement.

15. The system of claim 14, wherein the data processing system is configured to select the property of the concentration measurement from the group consisting of a peak amplitude, a full- width half maximum, a first moment, and a slope of a line extending between a baseline and the peak amplitude.

16. The system of claim 14, wherein the data processing system is configured to estimate the local arterial input function by:

assigning, on the basis of values associated with the selected properties, an overall score to each of the neighborhood voxels;

ranking the neighborhood voxels on the basis of their corresponding overall scores; and

including only a selected number of neighborhood voxels, the selected neighborhood voxels being selected on the basis of their respective ranks.

17. The system of claim 16, wherein the data processing system is configured:

to weight the included voxels on the basis of their respective distances to the target voxel; and

to assign a local arterial input function at least in part on the basis of a weighted combination of the included voxels.

18. The system of claim 10, further comprising a computer-readable medium having encoded thereon instructions that, when executed, cause the data processing system to execute the method recited in claim 1 .