Tweet
Tweet

Slide 1

Slide 1 text

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

Slide 2

Slide 2 text

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

Slide 3

Slide 3 text

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

Slide 4

Slide 4 text

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

Slide 5

Slide 5 text

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

Slide 6

Slide 6 text

Research Content in Computational dMRI Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 6 / 81

Slide 7

Slide 7 text

Thesis overview Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 7 / 81

Slide 8

Slide 8 text

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

Slide 9

Slide 9 text

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

Slide 10

Slide 10 text

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

Slide 11

Slide 11 text

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

Slide 12

Slide 12 text

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

Slide 13

Slide 13 text

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

Slide 14

Slide 14 text

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

Slide 15

Slide 15 text

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

Slide 16

Slide 16 text

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

Slide 17

Slide 17 text

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

Slide 18

Slide 18 text

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

Slide 19

Slide 19 text

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

Slide 20

Slide 20 text

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

Slide 21

Slide 21 text

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

Slide 22

Slide 22 text

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

Slide 23

Slide 23 text

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

Slide 24

Slide 24 text

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

Slide 25

Slide 25 text

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

Slide 26

Slide 26 text

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

Slide 27

Slide 27 text

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

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

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

Slide 30

Slide 30 text

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

Slide 31

Slide 31 text

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

Slide 32

Slide 32 text

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

Slide 33

Slide 33 text

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

Slide 34

Slide 34 text

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

Slide 35

Slide 35 text

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

Slide 36

Slide 36 text

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

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

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

Slide 39

Slide 39 text

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

Slide 40

Slide 40 text

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

Slide 41

Slide 41 text

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

Slide 42

Slide 42 text

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

Slide 43

Slide 43 text

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

Slide 44

Slide 44 text

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

Slide 45

Slide 45 text

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

Slide 46

Slide 46 text

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

Slide 47

Slide 47 text

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

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

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

Slide 50

Slide 50 text

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

Slide 51

Slide 51 text

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

Slide 52

Slide 52 text

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

Slide 53

Slide 53 text

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

Slide 54

Slide 54 text

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

Slide 55

Slide 55 text

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

Slide 56

Slide 56 text

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

Slide 57

Slide 57 text

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

Slide 58

Slide 58 text

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

Slide 59

Slide 59 text

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

Slide 60

Slide 60 text

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

Slide 61

Slide 61 text

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

Slide 62

Slide 62 text

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

Slide 63

Slide 63 text

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

Slide 64

Slide 64 text

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

Slide 65

Slide 65 text

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

Slide 66

Slide 66 text

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

Slide 67

Slide 67 text

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

Slide 68

Slide 68 text

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

Slide 69

Slide 69 text

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

Slide 70

Slide 70 text

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

Slide 71

Slide 71 text

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

Slide 72

Slide 72 text

Comparison between SRPE and SPFI [Cheng MICCAI 2012] Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 65 / 81

Slide 73

Slide 73 text

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

Slide 74

Slide 74 text

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

Slide 75

Slide 75 text

Interpolation: 1D Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 68 / 81

Slide 76

Slide 76 text

Interpolation in 2D space and PGA Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 69 / 81

Slide 77

Slide 77 text

Interpolation in 2D space and PGA Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 69 / 81

Slide 78

Slide 78 text

Filtering Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 70 / 81

Slide 79

Slide 79 text

Robustness in atlas estimation Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 71 / 81

Slide 80

Slide 80 text

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

Slide 81

Slide 81 text

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

Slide 82

Slide 82 text

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

Slide 83

Slide 83 text

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

Slide 84

Slide 84 text

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

Slide 85

Slide 85 text

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

Slide 86

Slide 86 text

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

Slide 87

Slide 87 text

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

Slide 88

Slide 88 text

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

Slide 89

Slide 89 text

Thank you for your interests! Questions? Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 81 / 81