$30 off During Our Annual Pro Sale. View Details »

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

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

My PhD defense.

Jian Cheng

May 30, 2012
Tweet

More Decks by Jian Cheng

Other Decks in Research

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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

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

    View Slide

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

    View Slide

  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

    View Slide

  9. Diffusion Tensor Imaging (DTI)
    Free diffusion, Gaussian diffusion assumption [Basser 1994]
    P(R) =
    1
    (4πτ)3|D|
    exp −
    1

    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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    − 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

    + 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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    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

    View Slide

  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

    View Slide

  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


    δ(l)δ(m) −
    1

    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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    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

    View Slide

  36. SHORE in AFT-SC framework
    E(0) = 1 in estimation
    E(0) = 1 =⇒
    N
    n=0
    an00Gn0(0) =


    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

    View Slide

  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

    View Slide

  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

    (
    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

    View Slide

  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

    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










    n=0
    (−4π2τq2)n
    n!
    (
    lm,l 0
    bm
    l
    Ym
    l
    (u))n








    = exp −
    q2

    nlm
    anlm
    q2
    ζ
    n
    Ym
    l
    (u)
    Jian Cheng (ccm, athena) Estimation and Processing of EAPs/ODFs May 30, 2012 37 / 81

    View Slide

  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

    ) 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

    ) 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

    View Slide

  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

    View Slide

  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

    )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

    ) 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

    View Slide

  43. More comparison between SHO-3D basis and SPF basis
    Spherical Polar Polynomial (SPP) basis:
    BSPP
    nlm
    (q|ζ) = exp −
    q2

    q2
    ζ
    n
    Ym
    l
    (u), n ≥ l/2
    Spherical Polar Non-Polynomial (SPNP) basis:
    BSPNP
    nlm
    (q|ζ) = exp −
    q2

    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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  64. Isotropic Distribution, Geodesic Anisotropy (GA)
    Isotropic ODF is unique, ψ(r) = 1

    , ∀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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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 α

    α(u) = cTK(q|ζ)c
    Kn l m
    nlm
    (q|ζ) =
    2L
    α=0
    α
    β=−α
    4π(−1)α
    2 Inn α(q)Qmm β
    ll α

    α(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

    View Slide

  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

    View Slide

  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 α








    α
    (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 α

    α
    (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

    View Slide

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

    View Slide

  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

    View Slide

  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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

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

    View Slide