Estimation and Processing of Ensemble Average Propagator and Its Features in Diffusion MRI

My PhD defense.

Transcript

1. Estimation and Processing of Ensemble Average Propagator and Its Features in Diffusion MRI Jian Cheng Center for Computational Medicine, LIAMA Institute of Automation, Chinese Academy of Sciences Athena Project-Team, INRIA Sophia Antipolis - Méditerranée PhD Defense, May 30, 2012 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 1 / 81
2. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 2 / 81
3. Water diffusion in human brain Water diffusion in human brain is hindered by surrounding tissues. [Johansen-Berg and Behrens 2009] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 3 / 81
4. Diffusion Data in 6D space 3D x-space (spatial space) and 3D R-space (diffusion displacement space) Ensemble Average Propagator (EAP) P(R) EAP field: P(x, R); ODF field: Φk (x, R) [Descoteaux 2008, Hagmann 2006] EAP profile P(R0r) Orientation Distribution Function (ODF) Φk(r) Φk(r) def = 1 Z ∞ 0 P(Rr)RkdR Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 4 / 81
5. Diffusion Weighted Imaging (DWI) The Pulse Gradient Spin-Echo (PGSE) sequence [Stejskal and Tanner 1965] Narrow pulse assumption (δ ∆) E(q) F3D ⇐⇒ P(R) [Callaghan, 1991] q = qu F3D ⇐⇒ R = Rr E(q) = S(q) S(0) q = γδG/2π Fourier relation P(R) = E(q)e−2πiq·Rdq E(q) = P(R)e−2πiq·RdR Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 5 / 81
6. Research Content in Computational dMRI Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 6 / 81
7. Thesis overview Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 7 / 81
8. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 8 / 81
9. Diffusion Tensor Imaging (DTI) Free diffusion, Gaussian diffusion assumption [Basser 1994] P(R) = 1 (4πτ)3|D| exp − 1 4τ RTD−1R D = 3 i=1 λi vi vT i Diffusion tensor model [Basser 1994] E(q) = S(q)/S(0) = exp(−4π2τqTDq) = exp(−buTDu) b = 4π2τ q 2 Tensor estimation (>6 DWIs) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 9 / 81
10. From tensors to Tissue Properties From tensor valued images to scalar images tenor field FA map MD map RGB map FA = √ 3||D − 1 3 Trace(D)I|| √ 2||D|| = 3 2 (λ1 − ¯ λ)2 + (λ2 − ¯ λ)2 + (λ3 − ¯ λ)2 λ2 1 + λ2 2 + λ2 3 MD = 1 3 Trace(D) = λ1 + λ2 + λ3 3 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 10 / 81
11. Riemannian framework for tensors and Gaussian distributions Metric selection for processing tensor field. Parametric family PF (a set of PDFs) parametrized by tensor Σ PF = P(R|Σ) = 1 (2π)3|Σ| exp − 1 2 RTΣ−1R : Σ ∈ Sym+ 3 Fisher metric in information geometry theory: [Rao, 1945, Atkinson 1981, Skovgaard 1984, Lenglet 2006, Pennec 2006, Fletcher 2004] gij = R3 ∂ √ P(R|Σ) ∂Di ∂ √ P(R|Σ) ∂Dj dR = 1 2 Tr(Σ−1Ei Σ−1Ej) Exponential map, Logarithmic map, geodesic mean of tensors (Gaussian distributions) interpolation, filtering, atlas estimation definite positive tensor estimation revisited later Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 11 / 81
12. Advantages and limitation of tensor model Advantages and limitations of DTI  >6 DWIs, normally 20-30 DWIs  useful scalar indices, FA, MD, etc.  DTI works well in the area with isotropic diffusion and single directional diffusion  DTI can not handle non-Gaussian diffusion and complex fiber configurations From Gaussian distributions to general non-Gaussian distributions Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 12 / 81
13. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 13 / 81
14. Modeling Beyond DTI: HARDI Sampling in DTI, in Cartesian grid, single shell and multiple shells sampling in DTI Diffusion Tensor Imaging (DTI) [Basser1994] Dense Cartesian grid Diffusion Spectrum Imaging (DSI) [Wedeen 200, 2005] Single shell sampling sHARDI: Q-Ball Imaging (QBI), exact QBI [Tuch 2004,Descoteaux 2007,Aganj 2010] Diffusion Orientation Transform (DOT) [Özarslan 2006] ... ... Multiple shell or arbitrarily sampling mHARDI: Diffusion Propagator Imaging (DPI) [descoteaux 2010] Spherical Polar Fourier Imaging (SPFI) [Assemlal 2009, Cheng 2010] Simple Harmonic Oscillator based reconstruction and estimation (SHORE) [Özarslan 2009, Cheng 2011a] ... P(R) = E(q)e−2πiq·Rdq Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 14 / 81
15. Diffusion Spectrum Imaging (DSI) Measure the signal in a 3D Cartesian grid in q-space Perform a numerical Fourier transform to obtain EAP [Wedeen 2000, 2005] Advantages and limitations of DSI  model-free, no Gaussian assumption  too many samples, large acquisition time  numerical Fourier transform, re-griding (interpolation, extrapolation)  Numerical integration for ODFs Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 15 / 81
16. Diffusion Orientation Transform (DOT) for P(R0r) Goal: simplify the FT based on mono-exponential decay assumption, [Ozarslan 2006] . E(qu) = E(q0u)q2/q2 0 3D Fourier transform to obtain EAP profile P(R0r) Advantages and limitations of DOT  single shell data, analytical solution via spherical harmonic basis  mono-exponential decay assumption is unrealistic.  the estimated EAP is the true EAP convolved by F3D{E(qu)q2/q2 0 E(qu)−1}  it has difficulty to handle multi-shell data. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 16 / 81
17. Exact QBI for ODFs, mono-exponential decay assumption Two types of ODFs: [Tuch 2004, Wedeen 2005] Φt(r) = 1 Z ∞ 0 P(Rr)dR Φw(r) = ∞ 0 P(Rr)R2dR Projection slice theorem: [Aganj 2010, Canales-Rodriguez 2009, Tristan-Vega 2010] Φt(r) = 1 Z Πr E(qu)qdqdu Φw(r) = 1 4π − 1 8π2 Πr ∆bE(qu) q dqdu Πr = {qu : uTr = 0} Goal: simplify the plane integral based on mono-exponential decay assumption, E(qu) = E(q0u)q2/q2 0 . Φt(r) = − q2 0 2Z S2 1 ln(E(q0u)) δ(rTu)du Φw(r) = 1 4π + 1 16π2 S2 ∆b ln(− ln(E(qu)))δ(rTu)du Funk-Radon Transform (FRT): FRT{f(u)}(r) = S2 f(u)δ(rTu)du SH is the eigenfunction of FRT and Laplace-Beltrami operator S2 Ym l (u)δ(rTu)du = 2πPl(0)Ym l (r) ∆bYm l (u) = −l(l + 1)Ym l (u) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 17 / 81
18. Exact QBI for ODFs, mono-exponential decay assumption Advantages and limitations of exact QBI  single shell data, analytical solution via spherical harmonic basis  mono-exponential decay assumption is unrealistic.  it has difficulty to handle multi-shell data. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 18 / 81
19. Q-ball Imaging for Φt(r), delta function decay assumption Goal: simplify the plane integral based on delta function decay assumption, E(qu) = E(q0u)δ(q − q0). Q-ball Imaging (QBI) for Φt(r), Funk-Radon Transform [Tuch 2004, Descoteaux 2007, Hess 2006] E(q0u) = lm clmYm l (u) Φt(r) ≈ 1 Z FRT{E(q0u)}(r) = 1 Z lm 2πPl(0)clmYm l (r) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 19 / 81
20. Q-ball Imaging for Φt(r), delta function decay assumption Advantages and limitations of QBI  single shell data, analytical solution is fast  delta function decay assumption is unrealistic.  it has difficulty to handle multi-shell data.  the estimated Φt(r) is too smooth, spherical deconvolution Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 19 / 81
21. The state-of-the-art in local estimation sampling Φt(r) = 1 Z ∞ 0 P(Rr)dR Φw(r) = ∞ 0 P(Rr)R2dR P(R) DTI       exp(−4π2τqTDq) DSI       [Wedeen 2000, 2005] original QBI [Tuch 2004] [Tristán-Vega 2009] many negative values δ(q − q0) [Descoteaux 2007] not converge exact QBI [Canales-Rodriguez 2009] [Aganj 2010] DOT exp(−4π2τD(u)) [Tristán-Vega 2010] [Özarslan 2006] SPFI [Assemlal 2009, Cheng 2010a] [Cheng 2010b] {Gn(q)Ym l (u)}  proposition  proposition  proposition SHORE  proposition  proposition {Gnl(q)Ym l (u)} [Özarslan 2009, Cheng 2011a] DPI       {qlYm l (u), q−l−1Ym l (u)} [Descoteaux 2010] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 20 / 81
22. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 21 / 81
23. numerical Spherical Polar Fourier Imaging (nSPFI) Spherical Polar Fourier expression of signals [Assemlal 2009] E(q) = N n=0 L l=0 l m=−l anlmGn(q|ζ)Ym l (u) BSPF nlm (q) = Gn(q|ζ)Ym l (u) Gaussian Laguerre polynomial Gn(q) = κn(ζ) exp −q2 2ζ L1/2 n (q2 ζ ), κn(ζ) = 2 ζ3/2 n! Γ(n+3/2) 1/2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 22 / 81
24. Analytical Spherical Polar Fourier Imaging (SPFI) P(R) is analytically obtained from Fourier dual SPF (dSPF) basis [Cheng 2010a] P(Rr) = N n=0 L l=0 l m=−l anlmFnl(R)Ym l (r) BdSPF nlm (R) = Fnl(R)Ym l (r) Fourier dual SPF (dSPF) basis, orthonormal basis in R-space BdSPF nlm (R) = R3 BSPF nlm (q)e−2πiqTRdq = Fnl(R)Ym l (r) Fnl(R) = 4(−1)l/2 ζ0.5l+1.5πl+1.5Rl 0 Γ(l + 1.5) κn(ζ) n i=0 (−1)i n + 0.5 n − i 1 i! 20.5l+i−0.5Γ(0.5l + i + 1.5)1F1( 2i + l + 3 2 ; l + 3 2 ; −2π2R2ζ) A linear transform from {anlm} to EAP profile represented by SH basis For a given R0 , EAP profile P(R0r) = L l=0 l m=−l cP lm (R0)Ym l (r), cP lm (R0) = N n=0 anlmFnl(R0) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 23 / 81
25. Analytical Spherical Polar Fourier Imaging (SPFI) Linear transform for ODF by Tuch Φt(r) def = 1 Z ∞ 0 P(Rr)dR [Cheng 2010b] Φt(r) = L l=0 l m=−l cΦt lm Ym l (r) cΦt lm = 2πζ Z N n=0 n i=0 κn(ζ) i − 0.5 i (−1)n−iPl(0)anlm Linear transform for ODF by Wedeen Φw(r) def = ∞ 0 P(Rr)R2dR [Cheng 2010b] Φw(r) = L l=0 l m=−l cΦw lm Ym l (r) cΦw lm = 1 √ 4π δ(l)δ(m) − 1 8π N n=1 n i=1 (−1)iκn(ζ) n + 0.5 n − i 2i i Pl(0)(−l)(l + 1)anlm Analytical transform avoids the numerical error {E(qi)}→{anlm }→{cP lm (R0)}, {cΦt lm }, {cΦw lm } Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 24 / 81
26. Scalar measurements in SPFI Return to origin probability (RTO): Po def = P(0) Po = R3 E(q)dq = 4 √ πζ 3 4 N n=0 (−1)n Γ(n + 1.5) n! an00 Generalized Fractional Anisotropy (GFA) for EAPs: GFA def = 1 − S2 P(Rr)Y0 0 (r)dr Y0 0 (r) 2 P(R) 2 = 1 − N n=0 a2 n00 N n=0 L l=0 l m=−l a2 nlm Mean-Squared Displacement (MSD): MSD def = R3 P(R)R2dR MSD 3 8π2.5 N n=0 an00 κn(ζ) ζ (2L3/2 n−1 (0) + L1/2 n (0)) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 25 / 81
27. Scale Estimation based on typical value: ζ = 1 8π2τD0 , D0 = 0.7 × 10−3mm2/s BSPF 000 ∝ exp(−4π2τq2D0) [Cheng 2010a, 2010b] optimization on both A and ζ (slow, local optimal) minA,ζ E − MSPFA 2 based on generalized High Order Tensor (GHOT) model fitting Theorem (Set the scale in SPF basis via GHOT model) When using SPF basis to fit the diffusion signal attenuation represented by GHOT model as E(q) = exp − ∞ n=1 ∞ l=0 l m=−l bnlm q2 ζ1 n Ym l (u) , if b100 > 0, the optimal scale parameter ζ of SPF basis for any given L and large enough N is ζ = ζ1 √ π b100 = 1 8π2τDp Pseudo-ADC in GHOT model, coefficient of 4π2τq2. Dp = b100 8π5/2τζ Spherical Polar Non-Polynomial (SPNP) basis: BSPNP nlm (q) = q2 ζ1 n Ym l (u) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 26 / 81
28. Consider the prior E(0) = 1 in estimation Two ways to consider the prior E(0) = R3 P(R)dR = 1 in estimation. Consider artificial shell E(0u) = 1 [Cheng 2010a, 2010b] min A E − MSPFA 2 + w2 E(0) − M(0) SPF A 2 + ATΛA E(0) = 1 =⇒ N n=0 anlmGn(0) = √ 4πδ0 l , 0 ≤ l ≤ L, −l ≤ m ≤ l a0lm = 1 G0(0)         √ 4πδ0 l − N n=1 anlmGn(0)         , 0 ≤ l ≤ L, −l ≤ m ≤ l N n=1 L l=0 l m=−l anlmYm l (u) Gn(q) − Gn(0) G0(0) G0(q) = E(q) − G0(q) G0(0) A = (M T SPF M SPF + Λ )−1M T SPF E Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 27 / 81
29. Whole estimation process 1. Scale ζ Estimation based on typical value or GHOT model: a) based on typical value: ζ = 1 8π2τD0 , D0 = 0.7 × 10−3mm2/s b) M: basis matrix generated from basis q2 0.5q2 max n Ym l (u) b = (MTM)−1MT ln E, b = (b100 , · · · , bN L L )T is the coefficient vector ζ = 0.5q2 max √ π b100 2. Least Square Estimation of coefficient vector A with E(0) = 1 consideration: M SPF = [Gn(q) − Gn(0) G0(0) G0(q)]: Ns x NA matrix from SPF basis, NA = N(L + 1)(L + 2)/2 E =               E(q1) − G0(q1) G0(0) . . . E(qNS ) − G0(qNs ) G0(0)               : Ns x 1 dimensional data vector Λ =             Λ 100 · · · 0 . . . ... . . . 0 · · · Λ NLL             : NA x NA regularization matrix Λ nlm = λll2(l + 1)2 + λnn2(n + 1)2 A = (M T SPF M SPF + Λ )−1M T SPF E A0 = (a000 , · · · , a0LL)T, a0lm = 1 G0(0)         √ 4πδ0 l − N n=1 anlmGn(0)         A = A0 A 3. Analytical Calculation of EAP and its features from A: P(Rr) = N n=0 L l=0 l m=−l anlmFnl(R)Ym l (r), P(R0r) = L l=0 l m=−l cP lm (R0)Ym l (r), cP lm (R0) = N n=0 anlmFnl(R0) Φt(r) = L l=0 l m=−l cΦt lm Ym l (r), Φw(r) = L l=0 l m=−l cΦw lm Ym l (r), Po, GFA, MSD Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 28 / 81
30. Invariance of diffusion time τ We do not know τ, only know b values and gradients {ui} b = 4π2τq, q = qu Theorem (Invariance for different diffusion time τ) For a given data set with known b values and gradients {uj}Ns j=1 and unknown diffusion time τ, if ζ is set by fitting GHOT model, then the estimated Φt(r), Φt(r, C), Φw(r), Φw(r, C), RTO, MSD, GFA are invariant for different diffusion time τ. If we denote the EAP estimated from τ by P(Rr, τ), then P(Rr, τ2) = τ1 τ2 3/2 P( τ1 τ2 Rr, τ1) The theorem is independent of the estimation (L1, L2, etc) of A. It is also for typical scale setting. N(R|2τ2D) = τ1 τ2 3/2 N( τ1 τ2 R|2τ1D) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 29 / 81
31. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 30 / 81
32. The state-of-the-art in local estimation sampling Φt(r) = 1 Z ∞ 0 P(Rr)dR Φw(r) = ∞ 0 P(Rr)R2dR P(R) DTI       exp(−4π2τqTDq) DSI       [Wedeen 2000, 2005] original QBI [Tuch 2004] [Tristán-Vega 2009] many negative values δ(q − q0) [Descoteaux 2007] not converge exact QBI [Canales-Rodriguez 2009] [Aganj 2010] DOT exp(−4π2τD(u)) [Tristán-Vega 2010] [Özarslan 2006] SPFI [Assemlal 2009, Cheng 2010a] [Cheng 2010b] {Gn(q)Ym l (u)}  proposition  proposition  proposition SHORE  proposition  proposition {Gnl(q)Ym l (u)} [Özarslan 2009, Cheng 2011a] DPI       {qlYm l (u), q−l−1Ym l (u)} [Descoteaux 2010] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 31 / 81
33. Analytical Fourier Transform in Spherical Coordinate (AFT-SC) E(q) is represented by {Bk(qu)}, where Bk(qu) = Rk(q)Qk(u). E(q) = K k=0 ck Bk(qu) = K k=1 ck Rk(q)Qk(u) Theorem (Plane Wave Expansion) e±2πiq·R = 4π ∞ l=0 l m=−l (±i)ljl(2πqR)Ym l (u)Ym l (r) Analytical Fourier Transform P(R) = K k=0 ck Dk(R) = K k=0 ∞ l=0 l m=−l ck Fkl(R)Tklm(r) Radial Integration: Fkl(R) = 4π(−i)l ∞ 0 Rk(q)jl(2πqR)q2dq Spherical Integration: Tklm(r) = S2 Qk(u)Ym l (u)du Ym l (r) If Qk(u) = Ym l (u), then Tklm(r) = Ym l (r) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 32 / 81
34. QBI in AFT-SC framework Delta function decay assumption: E(q) = L l=0 l m=−l alm δ(q − q0)Ym l (u)  P(R) = L l=0 l m=−l alm4π(−1)l/2jl(2πq0R)q2 0 Ym l (r)  Φt(r) = 1 Z ∞ 0 P(Rr)dR = q0 Z L l=0 l m=−l alm2πPl(0)Ym l (r)  Φw(r) = ∞ 0 P(Rr)R2dR = 1 2π2q0 L l=0 l m=−l alm(−1)l/2Ym l (r) ∞ 0 jl(x)x2dx P(R) has many negative values when b is big. Φw(r) does not exist. Addition theorem, Funk-Hecke theorem, plane wave expansion theorem. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 33 / 81
35. SHORE in AFT-SC framework SHO-3D basis: BSHO3 nlm (q|ζ) = Gnl(q|ζ)Ym l (u) [Özarslan 2009, Cheng 2011a] Gnl (q|ζ) = κnl (ζ)( q2 ζ )l/2 exp − q2 2ζ Ll+1/2 n−l/2 ( q2 ζ ) κnl (ζ) = 2 ζ3/2 (n − l/2)! Γ(n + l/2 + 3/2) 1/2 E(q) = N n=0 n l=0 l m=−l anlm Gnl (q|ζ)Ym l (u) ⇐⇒ P(R) = N n=0 n l=0 l m=−l anlm (−1)nGnl (R| 1 4π2ζ )Ym l (r) Theorem (Fourier transform of SHO-3D basis) F3D{BSHO3 nlm (q|ζ)} = (−1)nBSHO3 nlm (R| 1 4π2ζ ) Φk (r) = 1 Z 2N l=0 l m=−l clm Ym l (r) clm = (4π2ζ)− k+1 2 N n=l/2 n−l/2 j=0 anlm (−1)n+j j! n + l/2 + 1/2 n − l/2 − j 2k+l−1 2 +jΓ( k + l − 1 2 + j) MSD, RTO, GFA Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 34 / 81
36. SHORE in AFT-SC framework E(0) = 1 in estimation E(0) = 1 =⇒ N n=0 an00Gn0(0) = √ 4π a000 = 1 G00(0)         √ 4π − N n=1 an00Gn0(0)         N n=1 2n l=0 l m=−l anlm        Gnl(q) − Gn0(0)δ0 l G00(0) G00(q)        Ym l (u) = E(q) − G00(q) G00(0) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 35 / 81
37. DPI in AFT-SC framework DPI0 : solution of Laplace’s equation [Descoteaux, 2010] E(q) = 1 n=0 L l=0 l m=−l cnlm Rnl (q)Ym l (u) Fkl (R) = 4π(−i)l qmax 0 Rk (q)jl (2πqR)q2dq R0l (q) = ( q √ ζ )l R1l (q) = ( q √ ζ )−l−1 P(R0 r) ≈ 1 Z qmax 0 S2 E(qu)q2 exp −2πiqR0 uTr dqdu Φt (r) ≈ 1 Z Rmax 0 P(Rr)dR, Φw (r) ≈ 1 Z Rmax 0 P(Rr)R2dR Advantages and limitations of DPI  analytical solution is fast, generalization of QBI  model-based, ∆E(q) = 0  large modeling error when q is small and q is large, E(0) does not exist  approximation based qmax and Rmax Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 36 / 81
38. DPI in AFT-SC framework DPI0 : solution of Laplace’s equation [Descoteaux, 2010] E(q) = 1 n=0 L l=0 l m=−l cnlm Rnl (q)Ym l (u) Fkl (R) = 4π(−i)l qmax 0 Rk (q)jl (2πqR)q2dq DPI1 : a specific case of SHO-3D basis ( q √ ζ )lYm l (u), ( q √ ζ )−l−1Ym l (u) =⇒ exp − q2 2ζ ( q √ ζ )lYm l (u) = Gnl (q)Ym l (u)|n=l/2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 36 / 81
39. DOT in AFT-SC framework DOT0 : E(q) = exp(−4π2q2D(u)), mono-exponential decay [Özarslan 2009] Il (R, u) = 4π(−1)l/2 ∞ 0 E(q)jl (2πqR)q2dq = Γ(0.5l + 1.5)1 F1 (0.5l + 1.5, l + 1.5, − R2 4τD(u) ) R−l(−1)l/22l+1π0.5(D(u)τ)0.5l+1.5Γ(l + 1.5) P(Rr) = S2 Il(R, u)Ym l (u)du Ym l (r) DOT1 : E(q) = nlm anlm(q2 ζ )nYm l (u) D(u) = lm bm l Ym l (u) E(q) = ∞ n=1 (−4π2q2)n n! ( lm bm l Ym l (u))n = nlm anlmq2nYm l (u) DOT2 : E(q) = nlm anlm exp −q2 2ζ q2 ζ n Ym l (u) (SPNP basis) E(q) = exp(−4π2τb0 0 Y0 0 q2) exp        lm,l 0 4π2τq2bm l Ym l (u)         = exp − q2 2ζ         ∞ n=0 (−4π2τq2)n n! ( lm,l 0 bm l Ym l (u))n         = exp − q2 2ζ nlm anlm q2 ζ n Ym l (u) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 37 / 81
40. Summary of Analytical EAP Estimation Methods It is very easy to propose an analytical EAP estimation methods E(q) = K k=0 ck Bk(qu) = K k=1 ck Rk(q)Qk(u) P(R) = K k=0 ck Dk(R) = K k=0 ∞ l=0 l m=−l ck Fkl(R)Tklm(r) method Rk(q) Qk(u) Fkl(R) = 4π(−i)l ∞ 0 Rk(q)jl(2πqR)q2dq Tklm(r) QBI R(q) = δ(q − q0) Ym l (u) Fl(R) = 4π(−1)l/2jl(2πq0R)q2 0 Ym l (r) SHORE Rnl(q) = Gnl(q|ζ) Ym l (u) Fnl(R|ζ) = (−1)nGnl(R| 1 4π2ζ ) Ym l (r) SPFI Rn(q) = Gn(q|ζ) Ym l (u) Fnl(R) = ∗ ∗ ∗ Ym l (r) DPI0 R0l(q) = ( q √ ζ )l Ym l (u) F0l = (−1)l/2ql+1.5 max ζ−0.5lR−1.5Jl+1.5 (2πqmax R) Ym l (r) R1l(q) = ( q √ ζ )−l−1 F1l = (−1)l/2R−1.5ζ0.5l+0.5((πR)l−0.5 Γ(l+0.5) − Jl−0.5(2πqmaxR) ql−0.5 max ) DPI1 Rl(q) = (q2 ζ )l/2 exp(−q2 2ζ ) Ym l (u) Fl(R) = 2l+1.5ζ0.5l+1.5πl+1.5Rl exp(−2ζπ2R2) Ym l (r) DOT1 Rn(q) = ( q √ ζ )2n Ym l (u) q2n+l+3 max πl+1.5RlΓ(1.5+0.5l+n)1F2(1.5+0.5l+n;1.5+l,2.5+0.5l+n;−π2qmaxR) (−1)l/2ζnΓ(1.5+l)Γ(2.5+0.5l+n) Ym l (r) DOT2 Rn(q) = (q2 ζ )n exp(−q2 2ζ ) Ym l (u) 2n+0.5l+1.5ζ0.5l+1.5πl+1.5RlΓ(n+0.5l+1.5)1F1(n+0.5l+1.5,l+1.5,−2ζπ2R2) (−1)l/2Γ(l+1.5) Ym l (r) (SPNP) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 38 / 81
41. Theoretical Analysis and Comparisons Completeness. {Bk(q)} is complete if it can represent any square integrable E(q) in 3D space. More samples =⇒ less MSE. No modeling error if all samples are known. Representability. P(R) is globally affected by E(q) in R3. Fitting the given {E(qi)} is not enough for a good EAP estimation. Prior 1 (P1): E(0) = 1 because R3 P(R)dR = 1 Prior 2 (P2): E(q) tends to 0 as q increases Prior 3 (P3): E(q) radially decays like (but NOT) a Gaussian Orthogonality. Numerical stability Coefficients under orthonormal basis are independent with the basis order chosen in Least Square approximation if all samples are given. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 39 / 81
42. Theoretical Analysis and Comparisons mHARDI methods for single shell data. Unstable, Rn(q0)Ym l (u) Rn (q0)Ym l (u) is a constant independent with u Considering E(0) = 1, SPF basis with N = 1 is stable Separation information between spherical and radial parts method Rk(q) Qk(u)Completeness P1 P2 P3 orthogonal single shell Separation QBI R(q) = δ(q − q0) Ym l (u) in S2 No No No in S2 Yes Yes SHORE Rnl(q) = Gnl(q|ζ) Ym l (u) Yes YesYesYes Yes No No SPFI Rn(q) = Gn(q|ζ) Ym l (u) Yes YesYesYes Yes Yes, if N = 1 Yes DPI0 R0l(q) = ( q √ ζ )l Ym l (u) No No No No No No No R1l(q) = ( q √ ζ )−l−1 DPI1 Rl(q) = (q2 ζ )l/2 exp(−q2 2ζ )Ym l (u) No YesYesYes No Yes No DOT1 Rn(q) = (q2 ζ )n Ym l (u) in a ball Yes No No No No Yes DOT2 Rn(q) = (q2 ζ )n exp(−q2 2ζ ) Ym l (u) Yes YesYesYes No No Yes (SPNP) DOT0 — — in S2 YesYes No in S2 Yes Yes Prior 1 (P1): E(0) = 1 because R3 P(R)dR = 1 Prior 2 (P2): E(q) tends to 0 as q increases Prior 3 (P3): E(q) radially decays like (but NOT) a Gaussian Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 40 / 81
43. More comparison between SHO-3D basis and SPF basis Spherical Polar Polynomial (SPP) basis: BSPP nlm (q|ζ) = exp − q2 2ζ q2 ζ n Ym l (u), n ≥ l/2 Spherical Polar Non-Polynomial (SPNP) basis: BSPNP nlm (q|ζ) = exp − q2 2ζ q2 ζ n Ym l (u) Theorem (Equivalent bases form the same function space) {BSPP nlm (q|ζ)}n≤N ⇐⇒ {BSHO3 nlm (q|ζ)}n≤N {BSPNP nlm (q|ζ)}n≤N,l≤L ⇐⇒ {BSPF nlm (q|ζ)}n≤N,l≤L Theorem (Linear span of SHO-3D basis and SPF basis) Span{BSHO3 nlm (q|ζ)}n≤N ⊂ Span{BSPF nlm (q|ζ)}n≤N,l≤2N Span{BSHO3 nlm (q|ζ)}n≤N Span{BSPF nlm (q|ζ)}n≤N,l≤2N Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 41 / 81
44. Experimental comparisons: noise free synthetic data Cylinder model [Söderman 1995] mHARDO (SPFI, DPI, SHORE), b=500/1500/3000 s/mm2, 60 × 3 samples, without or with little regularization sHARDI (DOT, QBI, exact QBI), b=1500 s/mm2, 60 samples Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 42 / 81
45. Experimental comparisons: EAPs from noisy data mHARDI, b=500/1500/3000 s/mm2, 60 × 3 samples, SNR=10, λl = λn = 10−8 sHARDI, b=1500 s/mm2, 60 samples, λ = 0.006 EAP results: (A), (B): mHARDI methods with three shell data (C), (D): sHARDI and mHARDI methods with single shell data Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 43 / 81
46. Experimental comparisons: EAPs from noisy data mHARDI, b=500/1500/3000 s/mm2, 60 × 3 samples, SNR=10, λl = λn = 10−8 sHARDI, b=1500 s/mm2, 60 samples, λ = 0.006 ODF results: (A), (B): mHARDI methods with three shell data (C), (D): sHARDI and mHARDI methods with single shell data Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 43 / 81
47. Experimental comparisons: EAPs from noisy data mHARDI, b=500/1500/3000 s/mm2, 60 × 3 samples, SNR=10, λl = λn = 10−8 sHARDI, b=1500 s/mm2, 60 samples, λ = 0.006 Performance summary mHARDI methods with multiple shell data have better results than with single shell data EAP results are generally better than ODF results. EAPs (SPFI, 3 shells) > EAPs (SPFI, SHORE, 1 shell) > EAPs (DOT, 1 shell) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 43 / 81
48. Phantom data Public phantom data with 3 shells (b=500/1500/2000 s/mm2). It is more isotropic than normal data, thus the results needs min-max normalization. mHARDI, 3 shells, 64 × 3 samples, without or with little regularization sHARDI, b=1500 s/mm2, 64 samples Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 44 / 81
49. EAP profiles in phantom data at 25µm SPFI SHORE DPI SPFI SHORE DPI three shell data, b=650,1500,3000 s/mm2 single shell data, b=1500 s/mm2 SPFI SHORE DPI DOT DOT single shell data, b=2000 s/mm2 b=1500 s/mm2 b=2000 s/mm2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 45 / 81
50. Monkey data: scale selection 3 shell, b=500/1500/3000, 30 × 3 samples. scale selection (N , L ) based on pseudo-ADC map. pseudo-ADC, (N , L ) = (1, 4) GFA, (N, L) = (1, 4) MSD, (N, L) = (1, 4) pseudo-ADC, (N , L ) = (2, 4) GFA, (N, L) = (2, 4) RTO, (N, L) = (1, 4) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 46 / 81
51. Monkey data: EAP profile esimation scale selection (N, L) based on GFA map (without regularization) SPFI SHORE DPI 3 shell data, b = 500, 1500, 3000s/mm2 SPFI SHORE DPI DOT single shell data, b = 1500s/mm2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 47 / 81
52. Monkey data: EAP profile esimation scale selection (N, L) based on GFA map (without regularization) SPFI SHORE DPI 3 shell data, b = 500, 1500, 3000s/mm2 SPFI SHORE DPI DOT single shell data, b = 3000s/mm2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 47 / 81
53. Monkey data: ODF esimation No min-max normalization. SPFI SHORE ODF by Wedeen Φw(r), 3 shell data SPFI SHORE exact QBI ODF by Wedeen Φw(r), b=1500 s/mm2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 48 / 81
54. Monkey data: ODF esimation No min-max normalization. SPFI SHORE ODF by Tuch Φt(r), 3 shell data SPFI SHORE exact QBI ODF by Tuch Φt(r), b=1500 s/mm2 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 48 / 81
55. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 49 / 81
56. The basic of information geometry theory Parametric family (PF), parametric space (PS) PF = p(x|c) : χ p(x|c)dx = 1, p(x|c) ≥ 0, c ∈ PS ⊂ RK Fisher information metric gij(c) = χ ∂ log p(x|c) ∂ci ∂ log p(x|c) ∂cj p(x|c)dx = 4 χ ∂ p(x|c) ∂ci ∂ p(x|c) ∂cj dx Exponential map, logarithmic map, geodesic Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 50 / 81
57. Multinomial distribution family Multinomial distribution family PF = p(x|p) = K i=1 pi δ(x = i) : p ∈ RK, pi ≥ 0, K i=1 pi = 1, original parametrization PS = p = (p1 , p2 , · · · , pK)T ∈ RK : K i=1 pi = 1, pi ≥ 0 square root parametrization, gij = 4δij PS = c = (c1 , c2 , · · · , cK)T ∈ RK : K i=1 c2 i = 1, ci ≥ 0 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 51 / 81
58. Riemannian framework for tensors and Gaussian distributions Parametric family PF (a set of PDFs) parametrized by tensor Σ PF = P(R|Σ)= 1 √ (2π)3|Σ| exp − 1 2 RTΣ−1R : Σ is positive definite matrix Fisher metric: [Rao, 1945, Atkinson 1981, Skovgaard 1984, Lenglet 2006, Pennec 2006, Fletcher 2004] gij = R3 ∂ √ P(R|Σ) ∂Di ∂ √ P(R|Σ) ∂Dj dR = 1 2 Tr(Σ−1Ei Σ−1Ej) Exponential map: ExpΣ(W) = Σ1 2 exp(Σ−1 2 WΣ− 1 2 )Σ1 2 Logarithmic map: LogΣ(M) = Σ1 2 log(Σ−1 2 MΣ−1 2 )Σ1 2 Geometric Anisotropy: GA(Σ) = d(Σ, u), where u is the nearest isotropic tensor [¯ λ, ¯ λ, ¯ λ] Log-Euclidean framework: F(Σ) = Logu(Σ), where u is the identity matrix [1, 1, 1] Affine-invariant metric: d(Σ, M) = d(AΣAT, AMAT), d(P(R|Σ), P(R|M)) = d(P(R|AΣAT), P(R|AMAT)) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 52 / 81
59. PDF family PDF family PFK represented by orthonormal basis {Bi(x)} [Cheng 2009] PFK = p(x|c)= K i=1 ciBi(x) 2 : χ p(x|c)dx= K i=1 c2 i =1, K i=1 ciBi(x)≥0 Wavefunction (from Quantum Mechanics) ψ(x|c) = √ P(x|c) = K i=1 ciBi(x) χ can be any field, {Bi(x)} can be any orthonormal basis For ODF, χ is S2 and {Bi} is chosen as Spherical Harmonics {Ym l } For EAP, χ is R3 and {Bi} could be chosen as SPF basis {BSPF nlm } Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 53 / 81
60. Riemannian Framework Parameter space {c}: a closed convex subset of SK−1, contained in a convex cone with 90o Fisher Information metric: gij = 4δij Riemannian distance: dgij (p(·|c), p(·|c )) = dδij (c, c ) = arccos(cTc ) Exponential map: Expc(vc) = c = c cos ϕ + vc vc sin ϕ, where ϕ = vc Logarithmic map: Logc(c ) = vc = c −c cos ϕ c −c cos ϕ ϕ, where ϕ = arccos(cTc ) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 54 / 81
61. Diffeomorphism Invariance Theorem (Diffeomorphism Invariance) Fisher information metric is diffeomorphism invariant by definition. Mixture of Gaussian model, P(R) = w1 P(R|D1) + w2 P(R|D2) Transformed model, (A ◦ P)(R) = w1 P(R|ATD1 A) + w2 P(R|ATD2 A) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 55 / 81
62. Weighted Riemannian Mean/Median & Geometric Anisotropy Weighted mean: µw = arg minc∈PSK N i=1 wid(c, c(i))2 [Cheng 2009] Weighted median: mw = arg minc∈PSK N i=1 wid(c, c(i)) [Cheng 2010c] Theorem (Existance and Uniqueness) Weighted Riemannian mean and median uniquely exist in PS. They can be estimated via a simple gradient descent method. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 56 / 81
63. Weighted Riemannian Mean/Median & Geometric Anisotropy Weighted mean: µw = arg minc∈PSK N i=1 wid(c, c(i))2 [Cheng 2009] Weighted median: mw = arg minc∈PSK N i=1 wid(c, c(i)) [Cheng 2010c] Theorem (Reasonable Riemannian mean and median of PDFs) Weighted Riemannian mean p(x|µw ) and median p(x|mw) satisfies p(x|µw ) ≥ min{p(x|c(i))}N i=1 , p(x|mw) ≥ min{p(x|c(i))}N i=1 , ∀x ∈ χ. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 56 / 81
64. Isotropic Distribution, Geodesic Anisotropy (GA) Isotropic ODF is unique, ψ(r) = 1 4π , ∀r ∈ S2 Isotropic tensor is not unique, D = [λ, λ, λ] Nearest isotropic tensor is [¯ λ, ¯ λ, ¯ λ] Isotropic EAP is not unique, P(Rr) = F(R) Theorem (The nearest isotropic EAP [Cheng 2011b] ) For any EAP P(R) = ψ(R)2, its nearest isotropic EAP is Piso(R) = (ψSiso (R))2, where ψSiso (R) is the normalized version of the isotropic part of ψ(R), i.e. ψSiso (R) = ψLiso (R) ψLiso (R) , ψLiso (R) = R3 ψ(R)Y0 0 (r)dr Y0 0 (r) With SPF basis representation, ψ(R) = nlm cnlm BSPF nlm (R), the nearest isotropic EAP is ψSiso (R) = nlm sc,nlm BSPF nlm (R), sc,nlm = cnlm δ00 lm √ N n=0 c2 n00 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 57 / 81
65. Isotropic Distribution, Geodesic Anisotropy (GA) Geodesic Anisotropy (GA): Riemannian distance from the nearest isotropic EAP/ODF GA(p(x|c)) = d(c, sc) = arccos(cTsc) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 58 / 81
66. Log-Euclidean Framework & Affine-Euclidean Framework project all EAPs/ODFs onto the tangent space of one typical isotropic EAP/ODF with coordinate s = (1, 0, 0, ..., 0)T Projection diffeomorphism: F(c) = Logs(c) (Log-Euc) F(c) = c sTc − s (Affine-Euc) Closed form of mean µw = F−1 N i=1 wiF(c(i)) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 59 / 81
67. Comparison between Euclidean metric and Riemannian metric Different results if the PDFs are far from the uniform PDF. The uniform ODF is the unique isotropic ODF. Similar results. Uniform EAP does not exist. Isotropic EAP is not a unform PDF. Different results. Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 60 / 81
68. Outline 1 Thesis Overview 2 Background and Motivation Diffusion Tensor Imaging (DTI) Modeling Beyond DTI: High Angular Resolution Diffusion Imaging (HARDI) 3 Spherical Polar Fourier Imaging (SPFI) SPF basis and analytical forms Analytical Fourier Transform in Spherical Coordinate (AFT-SC) 4 Riemannian Framework for ODF and EAP Computing Riemannian Framework for general PDF Computing Applications of the Riemannian framework 5 Summary and conclusion Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 61 / 81
69. Wavefunction estimation Represent E(q) using wavefunction ψ(R) = nlm cnlmBnlm(R) E(q|c) = R3 N n=0 L l=0 l m=−l cnlmBnlm(R) 2 e−i2πqTRdR Theorem (Plane Wave Expansion) e±2πiq·R = 4π ∞ l=0 l m=−l (±i)ljl(2πqR)Ym l (u)Ym l (r) E(q|c) = nlm n l m αβ 4π(−1)α 2 cnlmcn l m Inn α(q)Qmm β ll α Yβ α(u) = cTK(q|ζ)c Kn l m nlm (q|ζ) = 2L α=0 α β=−α 4π(−1)α 2 Inn α(q)Qmm β ll α Yβ α(u) Inn α(q) = ∞ 0 Rn(R)Rn (R)jα(2πqR)R2dR Qmm β ll α = S2 Ym l (r)Ym l (r)Yβ α(r)dr Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 62 / 81
70. Wavefunction estimation least square cost function with radial and spherical regularizations c = arg min c =1 M(c), M(c) = 1 2 Ns i=1 cTK(qi|ζ)c − Ei 2 + 1 2 cTΛc gradient in tangent space ∇M(c) = ∂M(c) ∂c − cT ∂M(c) ∂c c ∂M(c) ∂c = Ns i=1 2 cTK(qi |ζ)c − Ei K(qi |ζ)c + Λc gradient descent method c(k+1) = Expc(k) −dt ∇M(c(k)) ∇M(c(k)) , where Expc(v) = c cos v + v v sin v initialization: c(0) = (1, 0, ..., 0)T, typical isotropic Gaussian EAP Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 63 / 81
71. Nonnegative definite ODF/EAP estimation Square Root Parameterized Estimation (SRPE): [Cheng 2012] nonnegative definite ODF/EAP in whole domain the integration of ODF/EAP is 1 From wavefunction to P(R0r) P(R0 r) = (ψ(R))2 = 2L α=0 α β=−α       nlm n l m cnlm cn l m Gn (R0 )Gn (R0 )Qmm β nn α        Yβ α (r) From wavefunction to Φw(r) Φw (r) = ∞ 0 (ψ(R))2 R2dR = nlm n l m ∞ 0 Gn (R)Gn (R)R2dR cnlm cn l m Ym l (r)Ym l (r) = nlm l m cnlm cnl m Ym l (r)Ym l (r) = 2L α=0 α β=−α nlm l m cnlm cnl m Qmm β ll α Yβ α (r) From wavefunction to Φt(r) Φt (r) = 1 Z ∞ 0 (ψ(R))2 dR = 1 Z nlm n l m ∞ 0 Gn (R)Gn (R)dR cnlm cn l m Ym l (r)Ym l (r) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 64 / 81
72. Comparison between SRPE and SPFI [Cheng MICCAI 2012] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 65 / 81
73. Interpolation: swelling effect, reasonable mean value interpolation Mean value interpolation (tensors, ODFs, EAPs) swelling effect: Tensor interpolation with Euclidean metric. Larger determinant. |Σ| > max{Σi } =⇒ N(0|Σ) < min{N(0|Σi)} Definition: Reasonable Mean Value Interpolation For a PDF valued function p(s, x), a mean value interpolation based on Ns samples {ps(i) (x)}Ns i=1 at position {s(i)}Ns i=1 is reasonable if the interpolated PDF ps(x) at position s satisfies ps(x) ≥ min{ps(i) (x)}N i=1 , ∀x ∈ χ Theorem (Reasonable Interpolations) Riemannian and Euclidean mean value interpolations on general EAPs/ODFs are both reasonable. Riemannian mean value interpolation on tensors is reasonable, and more precisely max{Ns(i) (R|Σi)}Ns i=1 ≥ Ns(R|Σ) ≥ min{Ns(i) (R|Σi)}Ns i=1 , ∀R ∈ R3 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 66 / 81
74. Interpolation: swelling effect Two ways in the interpolation of Gaussian EAPs: 1. interpolate tensors =⇒ from tensors to EAPs 2. interpolate EAPs Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 67 / 81
75. Interpolation: 1D Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 68 / 81
76. Interpolation in 2D space and PGA Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 69 / 81
77. Interpolation in 2D space and PGA Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 69 / 81
78. Filtering Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 70 / 81
79. Robustness in atlas estimation Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 71 / 81
80. Monkey data experiments EAP field at 15 µm ODF field GA of EAPs Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 72 / 81
81. Monkey data: EAP and ODF estimation EAPs at 15 µm EAPs at 25 µm ODFs proposed SRPE EAPs at 15 µm EAPs at 25 µm ODFs SPFI, scale is set by typical ADC Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 73 / 81
82. Monkey data: EAP atlas estimation RMM median RMM mean EUM mean original EAP atlas from five subjects re-estimated EAP atlas with one noisy subject Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 74 / 81
83. Summary Local estimation (SPFI, SRPE) AFT-SC framework, comparison of methods Riemannian framework for ODFs/EAPs Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 75 / 81
84. Summary Mean Shift Functional Detection, Riemannian Framework on SK, 2007-2008 [Cheng 2009a] Riemannian framework for ODFs, Spherical Harmonics in S2, 2008-2009 [Cheng 2009b, Cheng 2010c] Similar results for ODFs. Probably different results for EAPs. 2009 3D orthonormal basis for EAPs. 2009-2010 MICCAI 2009, Dr. Assemlal, 2009-2010 SPFI, 2009-2010 [Cheng 2010a, 2010b] Riemannian framework for EAPs, 2010 [Cheng 2011b] L1-SPFI, 2010-2011 comparisons between EAP estimation methods, 2010-2011 [Cheng 2011a] nonnegative definite ODF/EAP estimation 2011-2012 [Cheng 2012] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 76 / 81
85. Publications Jian Cheng, Aurobrata Ghosh, Rachid Deriche, Tianzi Jiang, An Intrinsic Diffeomorphism Invariant Riemannian Framework for Probability Density Function Computing in diffusion MRI, IEEE Transaction on Medical Imaging (in revision) Jian Cheng, Aurobrata Ghosh, Tianzi Jiang, Rachid Deriche, Spherical Polar Fourier Imaging, IEEE transaction on Medical Imaging (to be submitted) Jian Cheng, Tianzi Jiang, Rachid Deriche, Nonnegative Definite EAP and ODF Estimation via a Unified Multi-Shell HARDI Reconstruction, MICCAI 2012 Jian Cheng, Aurobrata Ghosh, Tianzi Jiang, and Rachid Deriche, Diffeomorphism Invariant Riemannian Framework for Ensemble Average Propagator Computing, MICCAI 2011 Jian Cheng, Tianzi Jiang, and Rachid Deriche, Theoretical Analysis and Practical Insights on EAP Estimation via a Unified HARDI Framework, CDMRI 2011 Jian Cheng, Sylvain Merlet, Aurobrata Ghosh, Emmanuel Caruyer, Tianzi Jiang, and Rachid Deriche, Compressive Sensing Ensemble Average Propagator Estimation via L1 Spherical Polar Fourier Imaging, CDMRI 2011 Jian Cheng, Aurobrata Ghosh, Tianzi Jiang, Rachid Deriche, Model-free and Analytical EAP Reconstruction via Spherical Polar Fourier Diffusion MRI, MICCAI 2010 Jian Cheng, Aurobrata Ghosh, Rachid Driche, Tianzi Jiang, Model-Free, Regularized, Fast, and Robust Analytical Orientation Distribution Function Estimation, MICCAI 2010 Jian Cheng, Aurobrata Ghosh, Tianzi Jiang, Rachid Deriche, A Riemannian Framework for Orientation Distribution Function Computing, MICCAI 2009 Jian Cheng, Feng Shi, Kun Wang, Ming Song, Jiefeng Jiang, Lijuan Xu, Tianzi Jiang, Nonparametric Mean Shift Functional Detection on Functional Space for Task and Resting-state fMRI, MICCAI 2009 Workshop Jian Cheng, Aurobrata Ghosh, Tianzi Jiang, Rachid Deriche, Riemannian Median and Its Applications for Orientation Distribution Function Computing, ISMRM 2010 Nianming Zuo, Jian Cheng, Tianzi Jiang, Diffusion Magnetic Resonance Imaging for Brainnetome: A Critical Review, Neuroscience Bulletin (accepted) Emmanuel Caruyer, Jian Cheng, Christophe Lenglet, Guillermo Sapiro, Tianzi Jiang, Rachid Deriche, Optimal Design of Multiple Q-shells experiments for Diffusion MRI, CDMRI 2011 Sylvain Merlet, Jian Cheng, Aurobrata Ghosh, Rachid Deriche, Spherical Polar Fourier EAP and ODF Reconstruction via Compressed Sensing in Diffusion MRI, ISBI 2011 Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 77 / 81
86. Perspective Compressive sensing EAP estimation. [Cheng 2011] Sampling scheme. Registration with Riemannian framework. Better Riemannian framework and SRPE Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 78 / 81
87. Role in Athena and CCM Christophe Lenglet: Riemannian framework for tensors from tensors to ODFs and EAPs Maxime Descoteaux: QBI with single shell data, DPI with multiple shell data from single shell to multiple shell, from ODF to EAP, theoretical comparison, compressive sensing) Aurobrata Ghosh: High Order Tensor, Cartesian Coordinate, positive definite ADC estimation Spherical Coordinate, coding, etc. Emmanuel Caruyer: Sampling Sampling in multiple shell, optimal sampling Sylvain Merlet: Compressive Sensing estimation SPF dual basis, L1-SPFI Yonghui Li, Nianming Zuo, Songma Xie, etc.: diffusion networks ODF/EAP estimation ⇒ fiber tracking ⇒ network analysis Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 79 / 81
88. Acknowledgements Rachid Deriche (Athena, INRIA) Tianzi Jiang (CCM, LIAMA, CASIA) Aurobrata Ghosh (Athena, INRIA) Emmanuel Caruyer (Athena, INRIA) Yonghui Li (CCM, LIAMA, CASIA) Chuishui Yu (Tianjin Medical University General Hospital) Wen Qin (Xuanwu Hospital) Maxime Descoteaux (Sherbrooke University) Gaolang Gong (Beijing Normal University) Théodore Papadopoulo (Athena, INRIA) Hongtu Zhu (UNC-CH) Rong Xue (Institute of Biophysics, CAS) Shan Yu (INRIA) Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 80 / 81
89. Thank you for your interests! Questions? Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 81 / 81