Fast Estimation of Diffusion Tensors under Rician noise by the EM algorithm

0.61 MB·English

by  Jia Liu

Fast Estimation of Diffusion Tensors under Rician noise by
the EM algorithm
JIA LIU ∗DARIO GASBARRA †AND JUHA RAILAVO ‡.
December2014
Abstract
Diffusiontensorimaging(DTI)iswidelyusedtocharacterize,invivo,thewhitematterofthecentral
5 nervesystem(CNS).Thisbiologicaltissuecontainsmuchanatomic,structuralandorientationalinforma-
1
tionoffibersinhumanbrain.Spectraldatafromthedisplacementdistributionofwatermoleculeslocated
0
2 in the brain tissue are collected by a magnetic resonance scanner and acquired in the Fourier domain.
AftertheFourierinversion,thenoisedistributionisGaussianinbothrealandimaginarypartsand,asa
n
a consequence,therecordedmagnitudedataarecorruptedbyRiciannoise.
J
Statistical estimation of diffusion leads a non-linear regression problem. In this paper, we present a
8
fast computational method for Maximum Likelihood estimation (MLE) of diffusivities under the Rician

] noise model, based on the Expectation Maximization (EM) algorithm. By using data augmentation, we
O
areabletotransformanon-linearregressionproblemintotheGeneralizedLinearModeling(GLM)frame-
C
work,reducingdramaticallythecomputationalcost.TheFisher-scoringmethodisusedforachievingfast
.
t convergenceofthetensorparameter. Thenewmethodisimplementedandappliedusingbothsynthetic
a
t andrealdatainawiderangeofb-amplitudesupto14000s/mm2. Higheraccuracyandprecisionofthe
s
[ Ricianestimatesareachievedcomparedwithotherlog-normalbasedmethods.Inaddition,weextendthe

1 MLframeworktothemaximumaposterior(MAP)estimationinDTIundertheaforementionedscheme
v byspecifyingthepriors. Wewilldescribehowclosenumericallyaretheestimatorsofmodelparameters
8
obtainedthroughMLandMAPestimation.
9
8
1
Keywords dataaugmentation,Fisherscoring,maximumlikelihoodestimator,maximumaposterioresti-
0
. mator,RicianLikelihood,reducedcomputation
1
0
5
1 1 Introduction
:
v
i
X Diffusion tensor imaging (DTI) is a powerful tool to detect, in vivo, the white matter anatomy and struc-
r turesofthebrain. TherawMR-dataarecollectedbyamagneticresonancescannerandconsistofspectral
a
measurement from the displacement distribution of water molecules constrained into cellular structures.
Diffusionanisotropycharacterizesthenervousfibers.
∗Correspondingauthor,DepartmentofMathematicsandStatistics,UniversityofJyväskylä,P.O.Box(MaD)FI-40014Finlande-
mail:[email protected]
†Department of Mathematics and Statistics, University of Helsinki P.O. Box 68 FI-00014 Finland e-mail:
[email protected]
‡HUSe-mail:[email protected]
1
After the Fourier inversion, the MR-signals are corrupted by a complex Gaussian noise, and conse-
quently,therecordedmeasurementmagnitudes,referredasdiffusionweightedmagneticresonanceimag-
ing (DW-MRI) data, will follow the Rician distribution. The noise distribution, however, will still stay
Gaussian in both real and imaginary components. The simplest method for diffusion tensor estimation
(DTE) is based on the linearized log-normal regression model, where the residual variance is assumed to
beeitherconstant(LeastSquares)ordependingonthesignalamplitude(WeightedLeastSquares). These
Gaussiannoisemodelsfailtofitthehighfrequencydata,whichcarryinformationaboutthehigherorder
diffusioncharacteristics. IntheexistingliteratureRajan, J.etal.(2011);Veraart, J.etal.(2011);Andersson
J.L.R.(2008)ontheML-estimationofdiffusiontensorsundertheRiciannoise,themaximizationalgorithm
involves repeated computation of modified Bessel functions. By using data augmentation we are able to
replacetheRicianlikelihoodbyaPoissonlikelihoodwhichisstandardintheframeworkofGLM.
SuchsimplificationreducesdramaticallythecomputationalburdenoftheFisher-scoringmaximization
algorithm. Thisappliesalsoathighb-amplitudes, whereinthelowsignalregimemeasurementsbelowa
thresholdarecustomarilycodedaszeros. InthestandardLSorWLSapproaches,zero-measurementsare
problematicsincetheycannotbefittedbyalog-normaldistribution,andsimplydiscardingtheminduces
selection bias. The appropriately modeled noise level provides capability of data correction in further
insights,e.g. removingartefactsfromtherawdata.
Thispaper isstructuredas follows. Section2 describesdataaugmentation andspecifiesthe statistical
model for DTE. In Section 3 we discuss the implementation of the EM and the Fisher-scoring algorithms
intheDTIcontext. Inaddition, wealsospecifypriorsfortheparametersanddiscussthecomputationof
theMaximumaPosterioriEstimator(MAPE)underthesamescheme. Section4illustratestheresultsfrom
bothsyntheticandrealdata.InSection5weconcludewithanoverviewofthemethodsandtheundergoing
developments.
2 GLM for MRI observations
2.1 RiciannoiseinMRI
In magnetic resonance imaging (MRI), we usually need to take the noise in the raw MR-acquisitions into
account. Thecomplexvaluednoise(cid:15)iscomposedoftwoi.i.d.Gaussianrandomvariableswithzeromean
andvarianceσ2,onefortherealandtheotheronefortheimaginarycomponent. AftertheFourierinver-
sion, the signal intensity S ≥ 0 is corrupted by the the complex Gaussian noise, and Y = |S +(cid:15)| will be
observed.
Consequently, the observed MR-signal magnitudes follow a Rician distribution resulting in the likeli-
hoodfunction
y (cid:18) y2+S2(cid:19) (cid:18)yS(cid:19)
pS,σ2(y)= σ2 exp − 2σ2 I0 σ2 , (2.1)
whereI istheα-ordermodifiedBesselfunctionoffirstkind. Forα=0ithasalsothefollowingrepresen-
α
tationintermsofGaussianhypergeometricseriesGradshteyn,I.S.,Ryzhik,I.M. (2007):
(cid:88)∞ τ2n
I (2τ)= F (1,τ2)= . (2.2)
0 0 1 (n!)2
n=0
2
Lett=S2/(2σ2),thenEq. (2.1)gives
y (cid:18) y2 (cid:19) (cid:18)y√ (cid:19)
Pt,σ2(Y ∈dy)= σ2 exp −t− 2σ2 I0 σ 2t dy (2.3)
√
withτ =yS/(2σ2)= 2ty/(2σ).
2.2 Dataaugmentation
We follow the strategy presented in Gasbarra D. et al. (2014) introducing an augmented data N from a
Poisson distribution with mean t > 0. The likelihood of the observed data can be transformed from the
RicianlikelihoodEq. (2.3)toajointaugmenteddensity
Pt,σ2(N =n,Y2 ∈dy2)=Pt,σ2(N =n,X ∈dx) (2.4)
(tx)n (cid:18) x (cid:19)
=P (N =n)P (X ∈dx|N =n)= exp −t− dx,
t σ2 (n!)2(2σ2)n+1 2σ2
whereX isfromtheconditionaldistributionGamma(N +1,1/(2σ2))givenN. Eq. (2.4)providesatrans-
formationfromanon-linearregressionproblemtotheGLMframework
(cid:18) (cid:19)
zξ−a(ξ)
f (z)=c(z,φ)exp (2.5)
ξ,φ φ
withzcorrespondingtotheresponseingeneral,seeMcCullagh,P.,Nelder,J.A.(1989)formoredetails.
3 Method
3.1 DW-MRIandparametrization
InDW-MRI,thesignalismodeledasthefirstequality
(cid:0) (cid:1) (cid:0) (cid:1)
S(q)=S exp −bd(g) =S exp Zθ ,
0 0
where the control vector q ∈ R3 is determined by the sequence of gradient pulses, b = |q|2, and g =
q/|q| ∈ S2 isavectorofunitlength. TheMR-signaldecaysexponentiallywithrespecttotheb-amplitude.
Dependingonthegradientdirectiongthedecayismodeledbythereflectionsymmetricdiffusivityfunction
d:S2 →R+.
Greateffortshavebeendevotedtomodelingthediffusivity,andingeneralwecanhaveparametrization
as the second equality. In the simplest model the diffusivity is expressed by a symmetric and positive
definiterank-2tensorD ∈R3×3,giving
logS(q)=logS −bg(cid:62)Dg=logS +Zθ ,
0 0
whereinthelefthandsidethediffusiontensorisparametrizedas
θ =(θ ,...,θ )(cid:62) :=(cid:0)D ,D ,D ,D ,D ,D (cid:1)(cid:62)
1 6 xx yy zz xy xz yz
withadesignmatrix
Z =Z(q)=−b(cid:0)g2,g2,g2,2g g ,2g g ,2g g (cid:1).
x y z x y x z y z
3
Inhighangularresolutionmodels(HARDI)(seee.g. BarmpoutisA.etal.(2009)),thediffusivityismodeled
withatotallysymmetriccartesiantensorDofordern∈N,as
3 3 3
(cid:88) (cid:88) (cid:88)
d(g):= ··· D g g ···g .
(cid:96)1,(cid:96)2,...,(cid:96)n (cid:96)1 (cid:96)2 (cid:96)2n
(cid:96)1=1(cid:96)2=1 (cid:96)2n=1
3.2 EMinMLE
In the optimization of the likelihood, we employ the EM (Expectation - Maximization) algorithm, which
is one among the iterative methods in the MLE or in the Maximum a Posterior Estimation (MAPE). The
EMalgorithmproceedsintwostepsandshortensthecomputationalcomplexitybyusingaugmenteddata.
In terms of our case, in the E-step we calculate the expectation of the log likelihood w.r.t the conditional
distributionofN givenbytheobservationsandotherparameterswithfixedvalues. IntheM-step,wefind
theMLparameterofS2andσ2bymaximizingtheaugmentedloglikelihoodquantities.Thecomputational
0
detailsarelistedinAppendixB.
TheloglikelihoodfromEq. (2.4)isexpressedas
Q:=log(cid:0)pt,σ2(N =n,Y)(cid:1)=c(Y,N)+Nlog(t)−(N +1)log(σ2)−t− 2Yσ22, (3.6)
wherec(Y,N)=Nlog(Y2)−2log(N!)−(N+1)log(2)doesnotdependon(t,σ2)whichwillbeomittedin
theM-step. FromSection3.1,wehavet=S2exp(2Zθ)/2σ2.
0
IntheEM-iteration,giventhecurrentparameterestimates(θ(k),S2(k),σ2(k)),weupdatetheconditional
0
expectationoftheaugmenteddataby
(cid:104)N(cid:105)(k) :=Et(k),σ2(k)(cid:0)N(cid:12)(cid:12)Y(cid:1)= τ(kI)0(cid:0)I21τ(cid:0)2(kτ)((cid:1)k)(cid:1) with τ(k) = YS0(k)2eσx2p(k(Z) θ(k))
IntheM-stepweupdateσ2andS2bytherecursions
0
(cid:18) m (cid:19)(cid:30)(cid:18) m (cid:19)
(σ(k+1))2 = (cid:88)(cid:0)(S(k))2exp(cid:0)2Z θ(k)(cid:1)+Y2(cid:1) 2m+4(cid:88)(cid:104)N (cid:105)(k) (3.7)
0 i i i
i=1 i=1
and
(cid:18) m (cid:19)(cid:30)(cid:18) m (cid:19)
(S(k+1))2 =2(σ(k))2 (cid:88)(cid:104)N (cid:105)(k) (cid:88)(cid:0)exp(cid:0)2Z θ(k)(cid:1) , (3.8)
0 i i
i=1 i=1
wheremisthenumberofacquisitionsateachvoxel.
Forthetensorparameterθ,weemployastabilizedFisherscoringmethod: giventhestabilizingparam-
eterα∈[0,1],weiteratetherecursion
(cid:18) (cid:19)−1
θ →θ+ (1−α)J(θ)+αS(θ)(cid:62)S(θ) S(θ), (3.9)
untilconvergencetoafixedpointLangeK.(2013). InEq. (3.9)thescoreS(θ)isgivenby
m m
S(θ)=2(cid:88)Z (cid:104)N (cid:105)(k)−(cid:0)S(k)/σ(k)(cid:1)2(cid:88)exp(2Z θ)Z(cid:62) ,
i i 0 i i
i=1 i=1
andthecorrespondingFisherinformationis
m
J(θ)=2(cid:0)S(k)/σ(k)(cid:1)2(cid:88)exp(2Z θ)Z(cid:62)Z .
0 i i i
i=1
4
TheinitialsoftheEMalgorithmcanbeobtainedthroughtheleastsquare(LS)fromatruncateddataset
with the diffusion weighting ranging from 0 ∼ 1000s/mm2 in order to fit the Gaussian model (see Jones
D.K.,BasserP.J.(2004),Barber,P.A.etal.(1998)). Topursuehigherqualityoftheinitials,wecouldfurther
apply the weighted least square (WLS) described in Zhu H. et al. (2007). In the Appendix C we compare
thedifferencesbetweenourEMalgorithmandthedirectoptimizationoftheRicianlikelihoodinEq. (2.1),
whichiscommonlyusedtocomputetheMLEinDTI.ItshouldbenotedthattheEMalgorithmisneeded
because of the latent augmented variables; it does not decrease the marginal likelihood of the data, see
AppendixAfortheproof.
3.3 EMinMAPE
IntheBayesianframework,theMaximumaPosteriorEstimation(MAPE)aimstoobtainthepointestimates
bymaximizingtheposteriordensity.ThedifferencebetweenMLEandMAPEinthisscenarioisintheprior
probabilityπ(ξ).Giventhedatay,thenormalizingconstantintheposteriordensityπ(ξ|y)doesnotdepend
ontheparameterξ. WefindtheMAPEbymaximizingthejointdensityπ(ξ)p (y),andthisisachievedby
ξ
iteratingtheEM-recursionwiththepenalizationlogπ(ξ)
(cid:26) (cid:27)
ξ(k+1) =argmax Eξ(k)(cid:0)logpξ(z,y)(cid:12)(cid:12)y(cid:1)+logπ(ξ) (3.10)
ξ∈Ξ
untilconvergencetoafixedpoint.Thelog-priorpenalizationtermhasaregularizingeffect,whichvanishes
asymptoticallyasthesamplesizegrowsAnderssonJ.L.R.(2008).
InDTE,wecanassignconjugatepriorsinlightofSection3.2forσ2 andS2. Sincewehavelittleknowl-
0
edge of the tensor parameter θ, we may choose non-informative priors which are either scale- or shift-
invariantJaynesE.T.(2002). AsimpleBayesianhierarchicalmodelisobtainedafterthefollowingchoices:
• σ2hasscaleinvariantimproperpriorwithdensityπ(σ2)∝1/σ2,
• S2 ∼Gamma(c ,c ),wherec ,c areverysmall.
0 1 2 1 2
• θ ∈RdhastheisotropiccenteredGaussianpriorN(0,Ω−1),whereΩisad×dprecisionmatrix.
ThepenalizedEM-updatesforMAPEaregivenby
(cid:18) m (cid:19)(cid:30)(cid:18) m (cid:19)
(σ(k+1))2 = 1(cid:88)(cid:0)(S(k))2exp(cid:0)2Z θ(k)(cid:1)+Y2(cid:1) (cid:88)(2(cid:104)N (cid:105)(k)+1)+1 (3.11)
2 0 i i i
i=1 i=1
and
(cid:18) m (cid:19)(cid:30)(cid:18) m (cid:19)
(S(k+1))2 = (cid:88)(cid:104)N (cid:105)(k)+c 1 (cid:88)(cid:0)exp(cid:0)2Z θ(k)(cid:1)+c . (3.12)
0 i 1 2(σ(k))2 i 2
i=1 i=1
Additionally,thisgivesthemodifiedscoreandFisherscoring
S˜(θ)=S(θ)−Ωθ, andJ˜=J(θ)+Ω, respectively.
UnderourBayesianmodelwithweakpriorstheMAPestimationEq. (3.11)andEq. (3.12)aresimilar
m
(cid:80)
as the ML updates Eq. (3.7) and Eq. (3.8). Indeed, usually (cid:104)N (cid:105) (cid:29) 1, and we can omit the difference
i
i=1
betweenEq.(3.7)andEq.(3.11).Thenwhenc andc aresmallenough,thedifferencebetweenthelikelihood
1 2
andposteriormodeofS , expressedinEq.(3.8)andEq. (3.12)respectively, canalsobeignored. Theonly
0
5
difference when updating θ, is that we have considered the correction between the elements of a tensor
representedbythepriordistribution,theinversecovariancematrix,Ω. Suchcorrectionmaybeignorable.
Remark: Sometimes the MLE can be treated as a special case of the MAPE where the precision of the
parametersdependonthechosenprior. Iftheeffectsofthepriorsareweakenoughtobeignored,thenthe
posteriordistributionisasymptoticallyapproximatedbythelikelihood. Theconsequenceisthatnumeri-
callytheMAPtendtotheMLestimatesnumerically. Suchremarkisnotunusual(seeSparacino, G.etal.
(2000))butnearlyhasneverappearedintheDTIliterature.
4 Results
4.1 SyntheticData
Syntheticdatasetsweresimulatedbychoosingapositivetensorof2nd-orderandof4th-orderwithfixed
S and the noise variance σ2. The simulated data sets in the experiments arise from models with param-
0
etervaluesresemblingtherealscenario. Everydatasetcontains1440measurementswhichweresampled
from 32 distinct gradients and 15 distinct increasing b values (knots) up to 14000s/mm2, each repeated
three times. The ground truth (GT) of high (H-) and low (L-) Rician noise, σ2, are 93,0405 and 12,8821,
respectively. Tocomparetheperformance,wefirstplottheMLestimatedRiciansignalsandtheGTshown
in Fig. 1. Fig.2 shows the empirical signal to noise ratio (SNR) of the distinct b values from the first 480
measurementsasanillustration. Forcomparisonofthemethods,wesimulate100datasetsfromhighnoise
caseunderthe4th-ordertensorandcomparethefirst480measurementsofthesamplemeansofSNRand
the GT under different methods in Fig.3, where "∗" denotes that only the low frequencies (b values less
than 1000s/mm2) are considered in the estimation. This figure reveals that our MLE and the WLS under
atruncateddatasetfittheGTwell. However,whencomparingthemeansquareerrors(MSE)ofthenoise
variance,theWLS*hasahugebias,54.777,comparedwiththeMSEoftheMLestimates,10.358. Thefitted
Rician signals are depicted in Fig. 4, where the estimated signals are retrieved from a low b and a high
b-value casesestimated fromthe first480measurements. Itreveals thatin thelowb-value caseboth the
WLS* and the MLE perform well. But our ML estimates show advantages in the high b-value case. The
reasonisthatfortheWLS*weusedatainformationwhichdonotconsiderthehighfrequencies,butonly
fit back to fetch reliable estimates of signals. We further check the MSE on tensor coefficients from the
15 distinct b values. The results are described in Fig.5, where we can conclude that the MLE is the best
onecomparedwithothermethods. TheaveragecomputationaltimeoftheaforementionedMLEmethod
underthe4th-ordertensormodelis0.4868seconds, whichisextremelyshorterthantheminutesrunning
timepervoexlfromthecurrentstandardmethodssuchasMATLABNelder-Meadbasedorgradient-based
estimators(seeGhoshA.etal.(2014);LandmanB.etal.(2007)).
4.2 RealData
Thedataconsistof4596diffusionMR-imagesofthebrainofanhealthyhumanvolunteer,takenfromfour
5mm-thickconsecutiveaxialslices,andmeasuredusingaPhilipsAchieva3.0TeslaMR-scanner.Theimage
resolutionis128×128pixelsofsize1.875×1.875mm2. Aftermaskingouttheskullandtheventricles,we
remainwitharegionofinterest(ROI)containing18764voxels.Intheprotocol,weusedallthecombinations
ofthe32gradientdirectionswiththeb-valuesvaryingintherange0−14000s/mm2,with2−3repetitions,
for a total of 23323644 data points. The average computational cost per voxel by our method the 4th-
6
Figure1:FittedRiciansignalsbytheproposedMLEmethod.Thebluecurvesdepictthesignalintensitiesof
theGT.Theredandgreencurvesshowthefittedsignalsunderthe2nd-ordertensormodelfromthehigh-
and low- noise level datasets. Correspondingly, the black crossings and the cayn stars are the empirical
valuesunderthe4th-ordertensormodel.
Figure 2: Empirical SNR as functions of b values. The GT of the upper curve is from the high SNR cor-
responding the lower noise level with σ2 = 12,8821. The bottom one has the high noise level with
σ2 = 93,0405. The red curves are fitted SNR under the 2nd-order tensor model. While the green stars
representtheempiricalSNRunderthe4th-ordertensormodel.
ordertensormodelfromthisdatasetis1.8331seconds. Weillustratetheresultsmainlyunderthe4th-order
tensormodel. Fig.6showsthemeandiffusivity(MD)andthefractionalanisotropy(FA)ofdiffusionfrom
two consecutive slices, where FA is computed from the results under 2nd-order tensor model, which is
7
Figure3:SamplemeanofSNR.Thesamplemeansarecalculatedfrom100simulateddatasets. TheSNRare
estimatedbydifferentmethods. ThebluecurverepresentstheGT.Thecyancurveandthegreenstarsare
theestimatesbytheLSandtheWLSwiththetruncateddatasets,respectively. Theredcurveistheresults
throughtheWLS,andtheblackcrossingsaretheempiricalvaluesbyourMLEmethod.
givenby
(cid:112)
3((λ −E[λ])2+(λ −E[λ])2+(λ −E[λ])2)
1 2 3
FA= (cid:112) . (4.13)
2(λ2+λ2+λ2)
1 2 3
TheaveragevaluesofFAfromthesetwoROIare0.2769mm2/sand0.2861mm2/s,respectively. Thecolor
inFArepresentstheorientationsofthefibers. Underthe4th-ordertensormodel,MDisexpressedas
1
MD= (D +D +D +2D +2D
5 1111 1122 1133 2222 3333
1
+2D )= trace(D). (4.14)
2233 5
TheaveragevaluesofMDfromSlice3and4are6.248e-03mm2/s,6.045e-03mm2/s,respectively,andwe
havethesameestimatedvaluesofMDunder2nd-ordertensormodel.
We also plot the Rician noise map of σ from the two consecutive slices shown in Fig. 7, where the
artefactsareclearlydepictedbywhitecolorrepresentingveryhighnoise, whichrevealsthetruescenario
fromtherawMRimages.
VisualizationofangularresolutionofDTIdataunderdifferenttensormodelsfromtheregionofinterest
(ROI) of two consecutive slices are displayed in Fig. 8, where the ROI is near the hippocampus and the
empty spaces inside of left parts of the diffusion profiles (DP) are the masked ventricle. DP depict under
the 4th-order tensors providing much angular information of diffusion, where the colors represents the
principleorientationsofdiffusionateachvoxel. ThesetensorprofilesareplottedbyMATLABfanDTasia
toolboxBarmpoutisA.etal.(2007).
8
(a)lowbvalue,b=62s/mm2
(b)highbvalue,b=14000s/mm2
Figure4: SamplemeanoftheRiciansignalsfromlowandhighbvalues. Theplotsillustratethemeansof
thesignalintensitiesatb=62and14000s/mm2,respectively,estimatedbybydifferentmethods.
Figure5: MSEontensorcoefficients
9
(a)Meandiffusivity(MD)
(b)Fractionalanisotropy(FA)
Figure6: MDandFAmapsfromtwoconsecutiveslices,wheretheestimatedFAarecomputedunderthe
2nd-ordertensormodel. ThecolorinFArepresentstheorientationsofthefibers. ThecolorcodedFAmaps
are drawn by using the software ExploreDTI Leemans A et al. (2009). The corresponding MD maps are
fromtheresultsunderthe4th-ordertensormodel, wherethewhitespotscorrespondingtothecorrupted
data(artefacts)withmeasuredmagnitudesincreasingathighb-value.
Figure 7: Rician noise map from two consecutive slices. The white curves in the left bottom of the slices
depicttheartefactscorrespondingtoveryhighnoise.
5 Discussion
Ourmethodsubstantiallydiffersfromthepreviousonesintheliteratureandtheadvantagesaresumma-
rized by the following points: 1) We introduce a novel data augmentation, which allows the non-linear
10