Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Neuroimage. Author manuscript; available in PMC 2014 January 1.
Published in final edited form as:
PMCID: PMC3508305

Bessel Fourier Orientation Reconstruction (BFOR): An Analytical Diffusion Propagator Reconstruction for Hybrid Diffusion Imaging and Computation of q-Space Indices


The ensemble average propagator (EAP) describes the 3D average diffusion process of water molecules, capturing both its radial and angular contents. The EAP can thus provide richer information about complex tissue microstructure properties than the orientation distribution function (ODF), an angular feature of the EAP. Recently, several analytical EAP reconstruction schemes for multiple q-shell acquisitions have been proposed, such as diffusion propagator imaging (DPI) and spherical polar Fourier imaging (SPFI). In this study, a new analytical EAP reconstruction method is proposed, called Bessel Fourier orientation reconstruction (BFOR), whose solution is based on heat equation estimation of the diffusion signal for each shell acquisition, and is validated on both synthetic and real datasets. A significant portion of the paper is dedicated to comparing BFOR, SPFI, and DPI using hybrid, non-Cartesian sampling for multiple b-value acquisitions. Ways to mitigate the effects of Gibbs ringing on EAP reconstruction are also explored. In addition to analytical EAP reconstruction, the aforementioned modeling bases can be used to obtain rotationally invariant q-space indices of potential clinical value, an avenue which has not yet been thoroughly explored. Three such measures are computed: zero-displacement probability (Po), mean squared displacement (MSD), and generalized fractional anisotropy (GFA).

Keywords: Ensemble Average Propagator, Bi-Gaussian Model, Extrapolation, q-Space Indices, High Angular Resolution Diffusion Imaging, Heat Diffusion Smoothing

1. Introduction

The aim of diffusion-weighted imaging (DWI) is to non-invasively recover information about the diffusion of water molecules in biological tissues. The most common form of DWI is diffusion tensor imaging (DTI) (Basser et al., 1994), which is a good model of diffusion-weighted signal behavior at low levels of diffusion weighting. However, DTI is limited by the Gaussian assumption, which is invalid at higher levels of diffusion weighting (b > 2000 s/mm2) and its inability to resolve multiple fiber orientations within a voxel (Wiegell et al., 2000; Alexander et al., 2001; Frank, 2001). In order to recover complex white matter (WM) geometry, high angular resolution diffusion imaging (HARDI) (Tuch et al., 2002), which reduces the diffusion signal sampling to a single sphere (i.e. single level of diffusion weighting) within q-space, was proposed. Many HARDI techniques (Tuch, 2004; Hess et al., 2006; Descoteaux et al., 2007; Canales-Rodriguez et al., 2009; Aganj et al., 2010) seek to extract the orientation distribution function (ODF), a probability density function describing the angular distribution of water molecules during diffusion. Unlike apparent diffusion coefficient (ADC) profiles, the maxima of the ODF are aligned with the fiber directions, making it useful in fiber tractography applications. However, the ODF only retrieves the angular content of the diffusion process.

The ensemble average propagator (EAP) provides more information about tissue microstructure than the ODF because it captures both the radial and angular information contained in the diffusion signal. The ODF, mathematically defined as the radial projection of the EAP, is simply an angular feature of the EAP. Unlike the diffusion tensor, the EAP profiles illustrate and recover crossing fibers. The authors in (Cohen and Assaf, 2002) have suggested that the radial part of the diffusion signal may be sensitive to WM disorders caused by demyelination and could be used to infer the axonal diameter.

The estimation of the EAP requires combination of high angular sampling at multiple levels of diffusion weighting. Under the narrow pulse assumption (Stejskal and Tanner, 1965), the diffusion signal attenuation, E(q) in q-space and the EAP, P(p), are Fourier Transform (FT) pairs (Callaghan, 1991):


where E(q) = S(q)/So is the normalized q-space diffusion signal, S(q) is the diffusion signal measured at position q in q-space, and So is the baseline image acquired without any diffusion gradients (q = 0). We denote q = q u(θ, [var phi]) and p = p r(θ′, [var phi]′), where u and r are 3D unit vectors. The wave vector q is q = γδG/2π, where γ is the nuclear gyromagnetic ratio and G = gu is the applied diffusion gradient direction. The norm of the wave vector, q, is related to the diffusion weighting level (b-value) via b = 4π2q2(Δ − δ/3) (Basser, 2002), where δ is the duration of the applied diffusion gradients and Δ the time between the two pulses. Eq. (1) is valid only if the narrow pulse condition is met, which is rarely the case for q-space diffusion MRI performed under experimental conditions. Several studies (Mair et al., 2002; Weeden et al., 2005; Bar-Shir et al., 2008) however, have shown that even when these assumptions do not hold, the Fourier relationship in Eq. (3) is still a reasonable approximation of the microstructural features. The diffusion displacements, however, will be consistently underestimated (Weeden et al., 2005).

Various methods already exist to reconstruct the EAP. Using the diffusion tensor frame-work, the EAP is simply described by a multivariate Gaussian function (Basser et al., 1994). The authors in (Ghosh and Deriche, 2010) presented a closed-form approximation of the EAP using higher order tensors, specifically the 4th order diffusion tensor. The diffusion orientation transform (DOT) (Ozarslan et al., 2006) is a HARDI technique that computes the iso-radius of the EAP. DOT assumes the radial diffusion follows a mono-exponential decay, which allows the radial integration in Eq. (1) to be solved analytically. The spherical integration is then solved numerically. The application of this technique, however, is limited by its mono-Gaussian assumption of the radial diffusion decay. In addition, the single shell approach of DOT gives an incomplete picture of the EAP, whose estimation requires signal measurements along all of q-space.

EAP reconstruction techniques using multiple diffusion weighting acquisitions can be divided into two strategies: Fast Fourier Transform (FFT) based and analytical. FFT based methods include diffusion spectrum imaging (DSI) (Weeden et al., 2005; Canales-Rodriguez et al., 2010) and hybrid diffusion imaging (HYDI) (Wu and Alexander, 2007; Rathi et al., 2011). DSI is based on direct sampling of the diffusion signal on a Cartesian q-space lattice. The FT in Eq. (1) is then numerically evaluated via FFT to obtain the EAP. A major advantage of DSI is that the EAP is estimated without any prior assumptions of behavior of the diffusion signal. However, DSI requires dense sampling of the Cartesian lattice, resulting in very long acquisition times. HYDI samples the diffusion signal along concentric spherical shells in q-space, with the measurements then being interpolated and regridded onto a 9 × 9 × 9 Cartesian lattice so that the EAP can be similarly reconstructed as in DSI. HYDI uses much fewer samples than DSI, making it more clinically feasible. However, the HYDI propagator reconstruction may suffer from the ad hoc signal interpolation and regridding.

The FFT is impractical for methods employing spherical q-space sampling schemes, such as HYDI, since the FFT requires data to lie on a Cartesian grid. It is also quite computationally expensive. Solving the FT in spherical coordinates (i.e. spherical Fourier Transform) instead, obviates the need for FFT and ad hoc processing. Analytic methods, seeking to obtain a closed-form solution of the EAP, pursue such a route. Currently, the two main analytical EAP reconstruction schemes are diffusion propagator imaging (DPI) (Descoteaux et al., 2011) and spherical polar Fourier imaging (SPFI) (Assemlal et al., 2009a; Cheng et al., 2010a,b).

DPI assumes that E(q) is a solution to the 3D Laplace’s equation [nabla]2E = 0, which results in the signal basis being composed of the regular and irregular solid harmonics. It is fast, and seems to work well with only a small number of samples. However, the DPI signal basis is an unrealistic model of E(q) because Laplacian modeling of diffusion signal entails that (1) E(0) does not exist, which arises from the irregular solid harmonic term, and (2) MSD of water molecules is zero, which will be proved in the Theory section. In addition, the DPI signal basis lacks orthonormality, and hence does not possess the robust numerical stability that would otherwise feature in an orthonormal basis.

SPFI models the diffusion signal in terms of an orthonormal basis comprising the spherical harmonics (SH) and Gaussian-Laguerre polynomials. The SPFI signal basis is, in fact, a modified solution of the 3D quantum mechanical simple harmonic oscillator problem. It is robust to noise and low anisotropy, and works well with just a few number of samples. However, SPFI has not been tested at b > 3000 s/mm2. A slightly modified version of the SPFI signal basis was proposed just recently by the authors in (Caruyer and Deriche, 2012). This paper, however, will only be concerned with the original SPFI basis.

A closely related basis to SPFI was proposed in (Ozarslan et al., 2008, 2009), which use the Hermite polynomials to estimate the 1D q-space diffusion signal. In addition to forming a complete orthogonal basis, the Hermite polynomials are also eigenfuctions of the Fourier transform. However, the 3D EAP solution has yet to be derived using this basis.

With respect to analytical EAP reconstruction methods, one valuable though overlooked use is in extracting rotationally invariant quantitative measures from them. Recently, the authors in (Assemlal et al., 2011) used the SPFI signal basis to compute the novel fiber population dispersion (FPD), an index which assess the presence of crossing fibers within a voxel. The FPD, however, is a relatively new measure that has not yet been computed for an actual human brain. More well-established q-space metrics include generalized fractional anisotropy (GFA) (Tuch, 2004), mean squared displacement (MSD) (Assaf et al., 2000; Wu and Alexander, 2007), and zero-displacement probability (Po) (Assaf et al., 2000). All three are simply scalar features of the EAP, and the GFA & MSD can be viewed as high angular resolution analogues of the DTI indices fractional anisotropy (FA) and mean diffusivity (MD) (Basser and Pierpaoli, 1996), respectively. An analytical representation of the EAP (and hence diffusion signal) can facilitate either analytic computation of such features or numerical efficiency in estimating them.

In this paper, we present Bessel Fourier Orientation Reconstruction (BFOR) (Hosseinbor et al., 2011). Rather than assuming the signal satisfies Laplace’s equation, we reformulate the problem into a Cauchy problem and assume E(q) satisfies the heat equation. The heat equation is a generalization of Laplace’s equation, which the latter approaches at the steady state (i.e. t → ∞). BFOR provides an analytical reconstruction of the EAP profile from diffusion signal and models the diffusion signal in terms of an orthonormal basis. In addition, it contains an intrinsic exponential smoothing term that allows one to control the amount of smoothing in the EAP estimation. The last point is significant because, although the Laplacian modeling intrinsically smoothes the diffusion signal, the amount of smoothing can not be controlled, and hence it may oversmooth the signal. In addition to heat diffusion smoothing, we also look at linear signal extrapolation as a potential means to mitigate the effects of common artifacts afflicting the reconstructed EAP profile, such as Gibbs ringing and signal truncation. Employing a hybrid, non-Cartesian encoding scheme in both synthetic and in vivo datasets, we reconstruct the EAP using BFOR, SFPI, and DPI and assess their performances. Lastly, we use BFOR to compute GFA, Po, and MSD, and compare BFOR’s accuracy in estimating such indices to that of DPI and SPFI.

The paper is organized as follows: in Theory Section, we develop BFOR, first by describing how to estimate the diffusion signal, and then deriving the analytical solution for the EAP using Eq. (1). Scalar features of the EAP are also introduced in this section. The Appendix carefully details the derivations of the BFOR signal basis, EAP, and q-space indices. In Material and Methods Section, we describe the implementation details of BFOR and present the synthetic and in vivo human brain datasets that will be used to validate and illustate BFOR and compare it to SPFI and DPI in Results Section. Lastly, we discuss our results and future applications of analytical EAP methods in Discussion Section.

2. Theory

The heat equation is ubiquitous in the natural sciences, arising in Fick’s law of diffusion, Brownian motion, Fourier’s law of thermal conduction, and the price variation over time of stocks. A nice feature of any solution to the heat equation is its temporal dependence, which can be viewed as an inherent smoothing control mechanism. DPI, by solving the heat equation at the steady state, has no smoothing control mechanism, and so can potentially oversmooth the signal. Smoothing can be useful in situations where reconstrutions greatly suffer from noise. Although our method is similar in spirit to DPI and SPFI, it significantly differs from them due to its inclusion of a smoothing term.

Consider the eigenvalue/boundary condition problem


which we use to solve the Cauchy problem


where f(q) is simply the acquired signal and I is some self-adjoint linear operator. We require λ > 0. Chung et al. in (Chung et al., 2007) derived a unique solution for (3):


where eλit is a smoothing term controlled by parameter t ≥ 0 and the coefficients are given by ai = left angle bracketf, ψiright angle bracket. The implication of Eq. (4) is that the solution decreases exponentially as t increases and smoothes out high spatial frequency noise much faster than low-frequency noise. In DPI, however, the steady state assumption permanently removes any temporal term, which governs the extent of smoothing, so there is no smoothing control mechanism. Note that t = 0 corresponds to no smoothing being applied.

Assuming that I = [nabla]2, where [nabla]2 is the 3D Laplacian operator in spherical coordinates, Eq. (2) becomes


Eq. (5) can be solved via separation of variables to obtain an orthonormal basis, which we show in Appendix A:


where αnl(j) is nth root of lth order spherical Bessel function of first kind jl and τ is the radial distance in q-space at which the Bessel function goes to zero. Yj are a modified real and symmetric SH basis proposed in (Descoteaux et al., 2011) to reflect the symmetry and realness of the diffusion signal. The index j: = j(l, m) = (l2 + l + 2)/2 + m is defined for l = 0, 2, 4, … and m = −l, …, 0, …, l. Hence, for j = {1, 2, 3, 4, 5, 6, 7, …}, l(j) = {0, 2, 2, 2, 2, 2, 4, …}. The eigenvalues are -λnl(j)=-αnl(j)2τ2. The SH are also used to model the angular profile of the diffusion signal in DPI and SPFI.

Within the context of our problem, g(q, t) is the diffusion signal. The assumption of a Laplacian operator results in Eq. (3) becoming the heat equation: 2E(q,t)=E(q,t)t . From Eq. (4) then, the diffusion signal can be expanded in terms of the spherical orthonormal basis ψnj given in Eq. (6):


where Cnj are the expansion coefficients, R=(L+1)(L+2)2 is the number of terms in the modified SH basis of truncation order L, and N is the number of roots for any spherical Bessel function of order l. The total number of coefficients in the expansion is W=N(L+1)(L+2)2. Note that the actual acquired signal from scanner is given at t = 0. In DWI, E(0) = 1, and so for our basis, we obtain the following identity (derived in Appendix B):


which holds for any u within the unit sphere S2 (i.e. u [set membership] S2).

An important property of the diffusion signal is that it asymptotically approaches zero as q → ∞. However, the spherical Bessel functions infinitely oscillate about zero, as shown in Fig. 1, so a finite upper bound τ is needed at which the BFOR signal model becomes zero. The fact that the radial basis in BFOR does not radially decay to zero but becomes zero at some point in q-space is the main limitation of the BFOR algorithm.

Figure 1
Plots of spherical Bessel functions of first kind, which form the radial basis of the BFOR signal solution, for different orders l. As q approaches infinity, they infinitely oscillate about zero.

In deriving the EAP, the spherical integration of Eq. (1) is made easier by expressing the Fourier kernel as a plane wave expansion:


Substituting Eq. (7) and (9) into Eq. (1), we obtain


where we use the orthonormal property of SH, i.e. ∫Yj(u)Yj(u)d2u = δjj and define


The integral in (11) is solved in Appendix C, and we can write the EAP as


The BFOR theoretical solution can be summarized in five steps:

  1. [nabla]2ψi(q) = −λiψi(q), ψi(q = τ, u) = 0
  2. E(q,t)=i=0aihi(t)ψi(q)
  3. 2E(q,t)=E(q,t)t, E(q, t = 0) = f(q)
  4. hi(t) = eλit; ai = left angle bracketf(q), ψi­(q)right angle bracket
  5. P(p,t)=E(q,t)e-2πiq·pd3q=4πj=1(-i)l(j)Yj(r){E(q,u,t)jl(j)(2πqp)Yj(u)d3q}

2.1. Rotationally Invariant q-Space Indices

In addition to analytical EAP reconstruction, the DPI, BFOR, and SPFI modeling bases can be used to obtain rotationally invariant q-space indices. Here, we look at three such measures: Po, MSD, and GFA.

Po = P(p = 0) is the probability density that a water molecule returns back to its initial position within the diffusion time (Assaf et al., 2000; Wu and Alexander, 2007). In a healthy adult brain, Po is greater in white matter (WM) than gray matter (GM) because WM has more restricting barriers including multi-layer myelin sheaths, axonal membranes, and microtubules. Such restrictivity increases the likelihood of a water molecule returning to its original position, whereas it has a very low probability of returning to its original position in areas of unrestricted isotropic diffusion (CSF). Hence, Po can be viewed as a measure of restricted diffusion. Several studies have shown Po to be sensitive to brain pathology, and suggesting that changes in myelin are the primary mechanism for differences in Po (Assaf et al., 2002; Bar-Shir et al., 2009; Wu et al., 2011).

Po can be evaluated either numerically or analytically. The authors in (Wu et al., 2008) computed Po by numerically summing the normalized diffusion signal E(q) over all diffusion measurements in q-space, and then correcting the sum by the sampling density. Analytical formulations of Po were derived for the SPFI and DPI signal bases (Cheng et al., 2010b; Descoteaux et al., 2011). Similarly, an analytical Po expression can be obtained using the BFOR basis, which is derived in Appendix D:


The MSD, which we will denote as left angle bracketp2right angle bracket, is simply the second moment of the EAP (Wu and Alexander, 2007): left angle bracketp2right angle bracket = ∫p2P(p)d3p. It is related to the MD, which in the case of Gaussian diffusion is given by the well-known Einstein relation left angle bracketp2right angle bracket = 6(Δ − δ/3)MD. Thus far, an analytical formulation of MSD exists only within the DTI framework. It is calculated numerically in q-space imaging, either by extracting the full width at half maximum of the EAP (Assaf et al., 2000) or taking the geometric mean of the diffusion signal over all directions on a HDYI shell (Wu et al., 2008) In Appendix E, we will show that, in general, the MSD is related to the diffusion signal by the relation


Since DPI assumes [nabla]2E = 0, it predicts the MSD to be zero: left angle bracketp2right angle bracketDP I = 0. The BFOR MSD is (derived in Appendix E)


Tuch in (Tuch, 2004) introduced the concept of GFA and defined it as std(ODF)/rms(ODF). Since ODF is only a feature of the EAP, the subsequent GFA map is derived solely from the angular content of the diffusion profile. Incorporating both the angular and radial contents of the diffusion profile into the definition of GFA will result in a radial dial of GFA maps, illustrating how anisotropy varies with diffusion displacement p. Therefore, we define a new GFA:


Another advantage of Eq. (16) is that it is better suited for multiple diffusion weighted MR experiments, unlike Tuch’s definition, which is single-shell HARDI-based. In order to capture the 3D anisotropy of the EAP-defined GFA maps, 1000 uniformly distributed vertices on a unit sphere in propagator space (i.e. 1000 values of θ′ and [var phi]′) were acquired using the approach described in (Wong and Roos, 1994).

3. Materials and Methods

In general, we are given k HARDI shell datasets. The number of encoding directions in each shell does not have to be the same. Each HARDI dataset corresponds to a different b-value. Across all k shells, we have total of M diffusion measurements (including the b = 0 measurement). Hence, from these multiple shell datasets, we want to reconstruct the EAP, P(p).

3.1. Numerical Implementation of BFOR

The task is to estimate coefficients Cnj in Eq. (12) from the observed signal E(q, u, t = 0). We achieve this by carrying out a linear least square (LLS) fitting with regularization in the radial and angular parts. We let S = [E(q1, t = 0) … E(qM, t = 0)]T be the M × 1 vector representing the M diffusion signal measurements across all k shells. We also let C represent the W × 1 vector of unknown expansion coefficients Cnj, where W=N(L+1)(L+2)2. Defining Znj(q, u) = jl(j)(αnl(j)q/τ)Yj(u), we let Z denote our M × W design matrix:


Thus, we have a simple linear model of the form S = ZC. This system of over-determined equations is solved with a regularized LLS solution yielding vector Ĉ given by


where Lreg is the Laplace-Beltrami regularization diagonal matrix with l2(l + 1)2 entries on the diagonal and Nreg is the regularization diagonal matrix for the radial basis, with entries n2(n + 1)2 on the diagonal. The angular and radial regularization matrices penalize, respectively, high degrees of the angular and radial parts of Eq. (7) in the estimation under the assumption that they are likely to capture noise (Assemlal et al., 2009a). They also serve to reinforce the positivity contraint of the EAP. λl and λn are the regularization terms for angular and radial bases, respectively.

3.2. Visualization of EAP

Lastly, from the estimated vector Ĉ, we can extract the Cnj coefficients needed to compute the EAP, Po, MSD, and GFA. The spherical function P(p, r, t) is the iso-probability profile at some instant of smoothening t for a given p- that is, the probability density that a water molecule, initially at the origin, diffuses a distance p along the direction r. It is computed by generating 800 equidistant points along the equator of a sphere of radius p i.e. the polar angle θ is fixed at π/2 and the azimuthal angle [var phi] is uniformly varied from 0 to 2π. The EAP profile P(p, r, t) is then interpolated along these 800 points. Thus, the resulting profiles are 2D with the equator perpendicular to the z-axis. It is important to note that in this paper smoothening was applied only to the EAP, itself, and not on the diffusion signal.

3.3. Value of τ Parameter

An important point to consider in the implementation is how to determine the parameter τ in the signal basis. In practice, the diffusion signal is bounded by the maximal q-value qmax achievable by the imaging system. The authors in (Assaf and Cohen, 1998) have shown that, depending on the length of diffusion time, the amount of signal present at b-values near 30,000 s/mm2 varies from about half a percent to about five percent, which means that the signal does not approach zero at qmax unless the diffusion weighting/diffusion time are very high/long. Thus, we conclude τqmax. Based on numerical simulations, we find the value of τ that best reconstructs the EAP to be τoptimal = qmax + Δq, where Δq is the (uniform) q-space sampling interval.

3.4. Diffusion MRI Data Acquisitions for Synthetic and Human Brain Data

The synthetic and in vivo datasets use a hybrid, non-Cartesian sampling scheme (Wu and Alexander, 2007), shown in Table 1. Since EAP reconstruction is sensitive to angular resolution, the number of encoding directions is increased with each shell to increase the angular resolution with the level of diffusion weighting. The number of directions in the outer shells were increased to better characterize complex tissue organization. Diffusion tensor elements for measurements in the second shell were calculated using non-linear least squares estimation with the Camino software package (Cook et al., 2006), which were then used to obtain the FA and principal eigenvector.

Table 1
HYDI Encoding Scheme for Synthetic and Human Datasets

3.4.1. Synthetic Data

The mono-exponential (also referred to as mono-Gaussian) mixture model (Tuch et al., 2002) is frequently used to generate synthetic data to validate a given EAP reconstruction, such as in (Assemlal et al., 2009a; Cheng et al., 2010b), where the maximum b-value used was 3000 s/mm2. However, diffusion MR imaging experiments using high b-values (> 2000 s/mm2) have shown that the diffusion signal decay is no longer mono-exponential. Studies in normal human brain, with b-values over an extended range of up to 6000 s/mm2, have shown that the signal decay is better described with a bi-exponential i.e. bi-Gaussian curve (Mulkern et al., 1999; Clark and Le Bihan, 2000). Similar findings were made for rat brain, using multiple b-values of up to 10000 s/mm2 (Niendorf et al., 1996). According to (Assaf and Cohen, 1998), a bi-exponential fit gives very good agreement with the observed water signal attenuation in excised brain tissue from rats for b-values of up to 2 – 3 × 104 s/mm2. Thus, BFOR, SPFI, and DPI were applied to simulations of crossing fiber configurations generated by a bi-Gaussian mixture model. Fig. 2 illustrates mono-exponential and bi-exponential decay curves, where the latter has a pronounced tail at high q values, indicating that it takes longer for the signal to decay to zero than under the mono-exponential assumption. The head and tail of the bi-exponential decay curve can be viewed as the fast and slow diffusion components, respectively (Clark and Le Bihan, 2000; Maier et al., 2004).

Figure 2
Plots of diffuison signal attenuation as a function of b-value illustrating mono-exponential and bi-exponenial decays. The tail of the bi-exponential can be interpreted as the slow diffusion component, while the head the fast diffusion component.

In bi-Gaussian mixture,


where Nb is the total number of simulated fibers, fkf the volume fraction of the fast component of the kth fiber, and fks the volume fraction of the slow component. The summation of all volume fractions is 1, i.e., k=1Nb[fkf+fks]=1. Dkf and Dks describe the diffusion tensor for the fast and slow components, respectively, of the kth fiber assuming no exchange between the fast- and slow-diffusion compartments. It should be noted that there is controversy over the assignment of these components and whether the bi-Gaussian model should take into account exchange between compartments (Mulkern et al., 1999). The ground truth of EAP is then


where ε = Δ − δ/3. For the synthetic data, the diffusion gradient duration is δ = 45 ms and diffusion gradient separation Δ = 56 ms.

In reconstructing the EAP, we look at two equally weighed fibers crossing at 60°, and set eigenvalues of each diffusion tensor to be [1.6,0.4,0.4]e-3, which gives an FA value of 0.7071. The values of the fast and slow Gaussian diffusion functions were taken from (Maier et al., 2004) and are shown in Table 2. Monte Carlo noise simulations were then performed to investiage the effect of SNR on the estimation of Po, MSD, and GFA for a single voxel for each EAP method. Seven SNR levels ([10 20 30 40 50 60 100]) for the b = 0 image were simulated, 1000 times each, by adding Rician noise in a similar manner as in (Descoteaux et al., 2007) for four different scenerios: a fast isotropic component (D = 0.00115 mm2/s); a slow isotropic component (D = 0.00045 mm2/s); fast anisotropic components of a corpus callosum fiber and internal capsule fiber crossing at 60°; and the slow anisotropic components for the previous scenario. The BFOR parameters are {L = 4, N = 6, τ = 91.2 mm−1, λl = 10−6, λn = 10−6}, DPI parameters {L = 4, λl = 0 (no noise)/λl = 0.006 (with noise)}, and SPFI parameters {L = 4, N = 3, ζ = 500, λl = 10−8, λn = 10−8}. For each method, model parameters were chosen based on giving the optimal EAP reconstruction when no noise was present.

Table 2
Fast/Slow Diffusion ADCs & Component Size Fractions (from (Maier et al., 2004))

3.4.2. Human Brain Data

HYDI was performed on a healthy, adult human using a 3.0 T GE-SIGNA scanner with an 8-channel head coil and ASSET parallel imaging. The DW pulse sequence was a single-shot, spin-echo, echo-planar imaging (SS-SE-EPI) with pulse-oximeter gating. The MR parameters were as follows: TE = 122 ms, TR 12 s, FOV = 256 mm, matrix = 128 × 128, voxel size = 2 × 2 mm2, 30 slices with slice thickness = 3 mm, and a total scan time of about 30 min. Diffusion parameters were maximum b-value bmax = 9375 s/mm2, diffusion gradient duration δ = 45 ms, diffusion gradient separation Δ = 56 ms, q-space sampling interval Δq = 15.2 mm−1, maximum length of the q-space wave vector qmax = 76 mm−1, field of view of the diffusion displacement space FOVp = (1/Δq) = 65 μm, and resolution of the diffusion displacement space Δp = (1/2qmax) = 6.6 μm (Callaghan, 1991). The same BFOR, DPI, and SPFI modeling parameters utilized for synthetic data were also used for in vivo data.

4. Results

BFOR, DPI, and SPFI are first applied to the numerical phantom and then on the real dataset. The numerical phantom is used to validate BFOR, compare its performance to those of DPI and SPFI, assess all three methods’ robustness in estimating the scalar measures Po, MSD, and GFA, and answer the following questions: (1) Can these methods properly reconstruct a diffusion signal acquired via hybrid sampling? (2) How does the slow diffusion component affect the EAP reconstruction and the estimations of the scalar quantities? (3) What can be done to reduce the effects of Gibbs ringing on the EAP reconstructions? It is also important to note that the EAP and quantitative scalar measures were reconstructed using only 125 diffusion measurements, while those presented in (Descoteaux et al., 2011) used 256.

4.1. Results of Synthetic Data

Can these methods be used for diffusion signal estimation?

Fig. 3 displays the BFOR, DPI, and SPFI signal fit for each shell and the corresponding ground truth. The BFOR signal basis fits the diffusion signal nearly perfectly for all shells. Both SPFI and DPI reasonably fit the diffusion signal for b ≥ 1500 s/mm2, but poorly reconstruct the inner most shell (b = 375 s/mm2). As expected, DPI also tends to oversmooth the diffusion signal, especially so at b = 3375, 6000 s/mm2. Results of DPI signal fitting were also reported in (Descoteaux et al., 2011), where unform sampling along four spherical shells (b = 2000, 4000, 6000, and 8000 s/mm2) was done. Although a hybrid sampling scheme is used here, our results are consistent with those of (Descoteaux et al., 2011) for b ≥ 1500 s/mm2. SPFI’s and DPI’s poor signal fit of the diffusion signal in inner most shell may be due to the fact that only six measurements were acquired in this shell (see Table 1), which may be inadequate for their respective radial bases.

Figure 3
The ground truth diffusion signal (green) and estimated signal (red) using BFOR, SPFI, and DPI when noise was absent. Two equally weighted WM fibers were simulated crossing at 60°. Measurements from all 5 shells were used.

Performances of BFOR, DPI, and SPFI in reconstructing EAP

Fig. 4 shows the EAP reconstruction for the fast diffusion component across propagator space using each method. Modeling the fast component (i.e. head of the bi-exponential in Fig. 2) is tantamount to fitting a mono-exponential curve, and so the EAP reconstruction for the fast diffusion component can be viewed as if the diffusion signal decay was mono-exponential. Both BFOR and SPFI model the fast component EAP very well, accurately capturing the geometry and orientation of the EAP profile, and the BFOR reconstruction is nearly identical to that of SPFI. DPI performs reasonably well, but tends to overestimate the EAP.

Figure 4
Reconstruction of the fast component EAP (red) using BFOR, SPFI, and DPI compared with the ground truth (green). Two equally weighted WM fibers were simulated crossing at 60°.

Fig. 5 shows the EAP reconstruction for the slow diffusion component, which can be viewed as modeling the tail of the bi-exponential curve in Fig. 2. Note that the BFOR and SPFI reconstructions are quite alike. At p = 1 μm, all three methods capture the correct geometry of the ground truth EAP profile, but underestimate it, DPI more so. At p = 5 μm, all three methods are unsuccessful in capturing the correct geometry, in particular failing to capture the peaks of the ground truth EAP profile. However, they do capture the correct orientation. At p = 10 and especially p = 15 μm, the BFOR and SPFI EAP reconstructions of the slow diffusion component begin to suffer from Gibbs ringing, which arise from the truncation of the signal bases at high q. The DPI reconstruction at p = 15 μm benefits from the inherent smoothening of Laplacian signal modeling, with no spurious peaks present. Although oversmoothened and overestimated, it does a much better job than BFOR and SPFI in resolving the correct fiber orientation at p = 15 μm. The difficulty in reconstructing the EAP for the slow diffusion component is due to the slow diffusion component being sensitive to truncation effects. The reality of finite sampling makes it challenging to capture the tail of the bi-exponential curve. How then should one combat the effects of truncation artifacts?

Figure 5
Reconstruction of the slow component EAP (red) using BFOR, SPFI, and DPI compared with the ground truth (green). Two equally weighted WM fibers were simulated crossing at 60°.

Signal Extrapolation

Extrapolating the diffusion signal to higher q-values so that q-space is more thoroughly explored could mitigate the truncation effects. Signal extrapolation can increase the spatial resolution of the EAP (Cohen and Assaf, 2002), and in the case of DSI, significantly reduce the cumbersome q-space sampling (Yeh et al., 2008). By linearly damping the signal measurments in the outermost shell (b = 9375 shell), we were able to (linearly) extrapolate samples onto three new ‘pseudo-shells.’ Specifically, the outermost signal measurements were attenuated by a factor of 0.7, 0.4, and 0.1 to form the three ‘pseudo-shells.’ The BFOR and SPFI scaling factors, τ and ζ, respectively, were changed for the extrapolation to τ = 136.8 mm−1 and ζ = 1100. Note that the q-space sampling interval Δq was not changed for the extrapolation.

Fig. 6 shows that signal extrapolation improves the reconstruction of the slow EAP component for both BFOR and SPFI. At p = 10 μm, the BFOR and SPFI slow EAP reconstructions with signal extrapolation, although not perfectly capturing the ground truth geometry, better capture the angular features of the ground truth than those without signal extrapolation. The biggest improvement, however, is seen at p = 15 μm, where the pronounced Gibbs ringing is greatly reduced by the signal extrapolation. In particular, the BFOR and SPFI slow EAP reconstructions with signal extrapolation at p = 15 μm are not spiky and much closer to the ground truth, although their orientations are slightly off, than those without extrapolation (Fig. 5). Note that both the BFOR and SPFI slow EAP reconstructions with extrapolation are quite alike. The BFOR fast EAP reconstruction was not affected by extrapolation, being nearly idential to its counterpart without extrapolation. However, the SPFI fast EAP reconstruction with extrapolation was moderately less acurate than that without extrapolation. In general, signal extrapolation can significantly improve EAP reconstructions at larger diffusion displacements (i.e. p = 15 μm).

Figure 6
Extrapolated samples were acquired by linearly damping the signal measurements in the outermost shell. Reconstruction of the EAP (red) using BFOR and SPFI compared with the ground truth (green). Two equally weighted WM fibers were simulated crossing at ...

Estimation of q-space indices

Fig. 7 shows the results for the Po measurements. Without signal extrapolation, BFOR and SPFI asymptotically approach the ground truth fast anisotropic Po, whereas DPI overestimates it. All three methods, without signal extrapolation, severely underestimate slow anisotropic Po, which is due to the truncation of the signal bases at high q. Both BFOR and SPFI asymptotically approach the ground truth fast/slow isotropic Po, while DPI overestimates both. At low levels of SNR (e.g. 10 & 20), which is quite common in diffusion MRI, all three methods (without signal extrap.) have biased estimates of Po, though the variance is fairly small.

Figure 7
Monte Carlo simulation investigating the effect of noise on estimation of Po using fast/slow anisotropic components and fast/slow isotropic components. Error bars denote one standard deviation across 1000 trials.

When signal extrapolation is applied, the estimation of the slow anisotropic Po by BFOR and SPFI significantly improves, as shown in Fig. 7b. According to Fig. 7a, signal extrapolation does not asymptotically affect SPFI’s estimation of fast anisotropic Po, but slightly worsens that of BFOR’s. At low levels of SNR, however, signal extrapolation results in more severe overestimation of fast anisotropic Po than without signal extrapolation for both BFOR and SPFI.

Fig. 8 shows the results for the MSD measurements. BFOR estimates both the fast and slow anisotropic/isotropic MSD very well across SNR levels. However, at low levels of SNR, the variability (given by the standard deviation) of the BFOR estimation of MSD is quite large, indicating strong sensitivity to noise. SPFI without signal extrapolation severely underestimates anisotropic MSD, giving negative values for slow anisotropic MSD, which indicates that it will give inaccurate measurements of the MSD of WM. SPFI also underestimates fast isotropic MSD, but asymptotically approaches ground truth slow isotropic MSD. Interestingly, there is less variability in SPFI’s estimation of MSD than BFOR.

Figure 8
Monte Carlo simulation investigating the effect of noise on estimation of MSD using fast/slow anisotropic components and fast/slow isotropic components. Error bars denote one standard deviation across 1000 trials.

When signal extrapolation is applied to the MSD measurements, SPFI’s estimations of slow anisotropic and slow isotropic MSD significantly improve, well-estimating them across SNR levels, as shown in Figs. 8b and 8d, respectively. Specifically, the SPFI slow anisotropic MSD estimation is no longer negative with signal extrapolation. BFOR’s estimations of slow anisotropic and slow isotropic MSD with extrapolation are nearly identical to those without it. However, the variability of the BFOR and SPFI MSD measurements is much less than those without extrapolation. A negative consequence of signal extrapolation is that it increases the inaccuracy (i.e. asymptotically worsening) of the BFOR/SPFI fast anisotropic and fast isotropic MSD estimations, all of which are greatly underestimated. The simulations indicate signal extrapolation may be more beneficial to SPFI’s estimation of MSD than that of BFOR’s.

Fig. 9 shows the results for the GFA measurements. Asymptotically, BFOR and SPFI without signal extrapolation approach the ground truth fast anisotropic GFA, while DPI underestimates it. Across SNR levels, DPI severely underestimates slow anisotropic GFA, while both BFOR and SPFI greatly overestimate it. For the case of isotropic diffusion, where GFA=0, the GFA estimated by each method approaches zero at high SNR. However, at SNR=10, both SPFI and BFOR severely overestimate the isotropic GFA, giving values comparable to the GFA of WM. DPI’s estimation of the isotropic component is more robust to noise than BFOR and SPFI.

Figure 9
Monte Carlo simulation investigating the effect of noise on estimation of GFA using fast/slow anisotropic components and an isotropic component. GFA was computed at p = 10 μm. Error bars denote one standard deviation across 1000 trials.

SPFI with signal extrapolation still overestimates slow anisotropic GFA, though slightly less so than without it, but the extrapolation increases the estimation’s variability at the same time. Signal extrapolation has negligible effects on BFOR’s estimation of slow anisotropic GFA. Both BFOR’s and SPFI’s estimation of fast anisotropic GFA are not asymptotically affected by signal extrapolation, having similar convergences as those without signal extrapolation, but the extrapolation causes both methods to overestimate fast anisotropic GFA to a larger degree at SNR=10, 20, & 30 than without it. Based on Fig. 9c, both SPFI and BFOR with signal extrapolation overestimate isotropic GFA, across SNR levels, to a larger extent than without extrapolation, implying that extrapolation is quite sensitive to noise in CSF regions.

4.2. Results of Human Brain Data

Resolving single fibers

In Fig. 10, a 4 × 4 ROI was drawn on the splenium of corpus cal-losum. The EAP profiles reconstructed at p = 10 μm by each method have the fundamental peanut shape of a single fiber. Note that the BFOR and SPFI reconstructions in both cases are very similar. We see that application of smoothing t = 550 removes the center peaks of the BFOR EAP profiles. Whether these center peaks in the EAP profiles are the result of Gibbs ringing (i.e. artifical peaks) or describe some underlying biological process is an open question.

Figure 10
Axial slice of FA map of adult human brain at b=1500 s/mm2 (second shell), where 4 × 4 ROI is drawn on splenium of corpus callosum. Plotted is the EAP profile at p = 10 μm overlaid on FA map in 4 × 4 ROI using BFOR at (c) t = 0, ...

Resolving crossing fibers

In Fig. 11, a 4 × 9 ROI was drawn in a region of fiber crossing, where the EAP profiles were reconstructed at p = 10 μm. Although not identical, the BFOR and SPFI reconstructions are quite similar, and they recover and well discriminate crossing fiber configurations in the EAP. DPI, however, tends to oversmooth the EAP profiles of crossing WM fibers, resulting in spherical/oblate shapes that give the impression of isotropic diffusion. Based on Fig. 11d, the application of smoothing t = 60 to BFOR removes the center peaks from several voxels, but at the expense of slight angular smoothening of EAP profiles themselves.

Figure 11
Axial slice of FA map of adult human brain at b=1500 s/mm2 (second shell), where a 4 × 9 ROI is drawn on a region of crossing fibers. Plotted is the EAP profile at p = 10 μm overlaid on FA map in 4 × 9 ROI using (c) BFOR at t = ...

Fig. 12 shows the reconstructed EAP profiles for the same crossing fiber region, but at p = 15 μm. Fiber crossing configurations are recovered and well discriminated in the EAP for each method. Unlike at p = 10 μm, the DPI EAP reconstruction at p = 15 μm is sharper and does not suffer from oversmoothening. In fact, as the propagator radius p increases, the angular resolution improves, at the enpense of the EAP profiles becoming spiky, as is evident in Figs. 12a, 12e, and 12g. When a smoothing of t = 60 is applied to the BFOR EAP reconstruction (without signal extrapolation), the subsequent EAP profiles are still spiky. At t = 350, the spikiness is smoothed out, but many of the WM voxels have EAP reconstructions significantly differing, with respect to orientation and geometry, from those at t = 0. Figs. 12d and 12f show the BFOR and SPFI EAP profiles reconstructed at p = 15 μm with signal extrapolation, respectively, which are much less spiky than the corresponding ones without signal extrapolation, which is consistent with the synthetic results shown in Fig. 6. The signal extrapolation also smoothes the reconstructed EAP profiles, but unlike BFOR at t = 350, none of the WM voxels are oversmoothed to such a degree that their EAP profiles have oblate shapes. Unlike BFOR at t = 350, the underlying EAP geometry and orientation of the BFOR/SPFI reconstructions with signal extrapolation are fairly consistent to those without extrapolation (at t = 0). As observed at p = 10 μm, the BFOR and SPFI EAP reconstructions at p = 15 μm are quite similar, which is consistent with the synthetic data results.

Figure 12
EAP profiles reconstructed at p = 15 μm for same crossing fiber region shown in Fig. 11.

Q-space Indices

Table 3 shows the mean index value and corresponding standard deviation for genu and splenium of corpus callosum (WM) and putamen (GM). Three 4 × 4 ROIs were drawn on both the genu and splenium and one such ROI on both the left and right putamen, across several slices. The table shows that the SPFI MSD (without signal extrapolation) erroneously gives negative values for the MSD of genu and splenium. With extrapolation, the SPFI MSD of genu and splenium are positive.

Table 3
Values of Indices for Various WM & GM Structures

Based on the numerical simulations, signal extrapolation was applied to both BFOR’s and SPFI’s estimation of Po. Fig. 13 displays an axial slice of Po generated by each method. In the first row, we show Po computed from BFOR, SPFI, DPI, and numerically (Wu et al., 2008). The BFOR, SPFI, and numerical Po maps are quite similar, exhibiting rich GM/WM and tissue/CSF contrasts while the DPI Po map has less GM/WM contrast. In particular, based on Table 3, the Po ratio of WM to GM is slightly above 2 for BFOR and SPFI, while less than 2 for DPI.

Figure 13
Axial slice of Po generated by each method.

The MSD maps computed from BFOR, SPFI, DPI, and numerically (Wu et al., 2008) are shown in Fig. 14. Both the BFOR and numerical MSD maps exhibit rich tissue/CSF contrast, but have little WM/GM contrast, which is similar to the DTI MD. Table 3 shows that the BFOR MSD values for WM and GM are quite similar. In the SPFI (without signal extrapolation) MSD map, WM regions are completely dark, having negative MSD values. This is consistent with the results of the noise simulations, which showed that SPFI severely underestimates the MSD of WM. Signal extrapolation has the effect of enforcing the positivity constraint on the MSD for SPFI. However, both the BFOR and SPFI MSD maps with signal extrapolation have poor tissue/CSF contrast because of the noise induced by the signal extrapolation. With regards to BFOR, both the synthetic and in vivo data suggest that it’s best not to use signal extrapolation in estimating MSD. SPFI, however, does not generate reliable MSD maps either with or without signal extrapolation. Although DPI predicts the MSD to be zero, an MSD map was computed for it by numerically estimating the second moment of the diffusion propagator. The constrast of the MSD DPI map is completely inverted, with WM appearing bright and CSF dark.

Figure 14
Axial slice of MSD generated by each method.

Fig. 16 displays axial slices of the GFA computed at p = 5, 10, and 15 μm for each method, illustrating how the anisotropy of different WM regions, such as the corpus callosum and capsules, varies with diffusion displacement p. According to Table 3, at p = 5 μm, the anisotropy of corpus callosum is lower with respect to levels seen in DTI. At p = 10 μm, the corpus callosum is very anisotropic, as can be seen from Table 3, indicating that p = 10 μm is a diffusion displacement worth reconstructing the EAP at. The GFA at p = 15 μm is more noisy, which is due to truncation of signal basis at high q-values and 15 μm being well beyond the resolvable resolution (of diffusion displacement) limit. The BFOR and SPFI GFA maps without signal extrapolation are very similar, while WM regions in the DPI computed GFA maps at p = 5 and 10 μm have lower intensity than those of BFOR and SPFI, which is consistent with the underestimation of GFA(10) by DPI observed in the Monte Carlo noise simulations (Fig. 9). CSF regions in the BFOR GFA(10) map with signal extrapolation are more noisy than without it, which is consistent with the simulation results shown in Fig. 9c. In the case of SPFI, however, the noise level is very severe in CSF regions in the GFA(5) & GFA(10) maps with signal extrapolation than those without it.

Figure 16
GFA maps computed at p = 5, 10, and 15 μm.

5. Discussion

The three analytical EAP reconstruction schemes examined in this paper possess both certain advantages and disadvantages. Among the three, the DPI reconstruction uses the least number of expansion coefficients. According to both synthetic and in vivo data, DPI tends to greatly oversmoothen the EAP, especially p = 10 μm, but performs well at p = 15 μm where it did not make use of signal extrapolation. DPI’s assumption of Laplacian signal modeling, however, entails that the MSD is zero (refer to Eq. (14)). Fig. 9 indicates that DPI greatly underestimates GFA(10), which is reflected in Fig. 16 and Table 3. In addition, the DPI signal basis is only applicable at q > 0, which is unrealistic since the diffusion signal is defined at q = 0.

The SPFI signal basis possesses several advantages in that it radially decays to zero and has no singularity at q = 0. In addition, the EAP is derived via integration over the entire q-space, unlike BFOR and DPI, where the q-space integration is up to a certain bound that is related to qmax. Interestingly, however, the SPFI EAP reconstructions for both synthetic and real datasets are quite similar to those of BFOR, suggesting the two methods may be inherently related. According to the synthetic data, signal extrapolation greatly improves the SPFI EAP reconstruction at p = 15 μm. Although not shown in this paper, heat diffusion smoothening can be applied to SPFI. SPFI’s estimation of Po and GFA, either with or without signal extrapolation, is quite comparable to those of BFOR’s. However, it poorly estimates the MSD, which recalling Eq. (14), may be due to, from a computational standpoint, SPFI’s signal basis not being an eigenfucntion of the Laplacian operator.

The main limitation of the BFOR signal model, as mentioned in the Theory section, is that it infinitely oscillates about zero, which entails a finite integration of the signal over q-space (τ being the upperbound) to retrieve the EAP. However, based on Fig. 3, BFOR outperforms DPI and SPFI in modeling the diffusion signal. Heat diffusion smoothening helps in removing potentially spurious peaks, and signal extrapolation significantly improves the EAP reconstruction at p = 15 μm. According to both the synthetic and real data, BFOR gives reasonable estimates of all three q-space indices.

The slow component of diffusion is the most sensitive to truncation artificats, which can induce severe Gibbs ringing and adversely affect the orientation of reconstructed EAP. In this paper, signal extrapolation was proposed as a means to mitigate the effects of such artifacts, and was oberseved to be most effective at higher radii (i.e. p = 15 μm), where the effects of signal truncation artifacts are most pronounced. The significant improvement in the BFOR/SPFI EAP reconstruction at p = 15 μm via linear signal extrapolation hence suggests that extrapolation may be a useful preprocessing step for EAP reconstruction at large diffusion diplacements. Signal extrapolation also greatly improves the accuracy of the BFOR/SPFI Po estimation, according to the synthetic data. However, signal extrapolation increases the severity of noise in CSF regions in the GFA and MSD maps for both BFOR & SPFI, as evidenced by the synthetic and real datasets, reducing tissue/CSF contrast. Hence, extrapolation may not be desirable for GFA and MSD estimation. Future work includes optimizing the signal extrapolation for a given signal basis.

The degree of heat diffusion smoothing desired depends on the propagator radius and whether the fibers are single or crossing. Based on Fig. 10, a smoothing of t = 550 was applied to splenium at p = 10 μm to remove the center peaks. However, for a crossing fiber region at p = 10 μm, a smoothing of t = 60 was only applied because the EAP profiles of crossing fibers can easily become oversmoothed, resulting in oblate shapes. The fact that the smoothing factor is different for different brain regions poises a problem for whole brain EAP processing. One way to address this issue is to use an optimal bandwidth selection framework from statistics to estimate the optimal t. Specifically, the bandwidth t is selected to minimize a certain cost function. In a spline setting, the cost function will be a generalized cross-validation (GCV) criterion (Katkovnik, 1999). In a more simple setting like ours, we can choose the t that minimizes the sum of the squared residuals, where the residual is simply the difference between the actual data and model fit.

Although the encoding scheme in this study consisted of equally spaced concentric spherical shells, the BFOR framework does not require such a scheme. BFOR only requires a minimum of two diffusion weightings and use of a spherical coordinate system. Random sampling along q-space or even the use of unequally spaced concentric shells is perfectly valid. This, however, leads to the important question of what is the best way to sample N diffusion measurements in q-space, which have started to be addressed (Assemlal et al., 2009b; Merlet et al., 2011; Caruyer et al., 2011). Future work includes optimizing the q-space sampling and applying compressed sensing to BFOR.

Both the ODF and EAP profiles are not sharp enough to extract the true fiber orientation. Rather, the fiber orientation is given by the fiber orientation distribution function (fODF), which can be computed via spherical deconvolution of some assumed kernel (i.e. response function) from q-space diffusion signal (Tournier et al., 2004). Mathematically, the angular convolution is given by


where F(r) is the fODF and K the kernel. The derivation of the fODF using the BFOR, SPFI, and DPI signal bases is worth exploring in the future.

In any future clinical study employing HYDI to examine brain pathology, where rotationally invariant indices like GFA and Po can be used to assess changes between diseased and normal subjects, voxel-wise analysis is desired. However, spatial normalization of multiple b-value datasets is no easy task. Recently, the authors in (Du et al., 2012) proposed a registration algorithm to align HARDI datasets using the ODFs. Specifically, the algorithm seeks an optimal diffeomorphism of large deformation between two ODF fields across a spatial volume domain and at the same time, locally reorients an ODF in a manner consistent with the underlying anatomical structure. HYDI images could be aligned using the same algorithm, except replacing the ODF with the EAP.

The MSD measure is quite sensitive to noise (Assaf and Cohen, 2000; Wu and Alexander, 2007). The authors in (Wu et al., 2008) proposed an alternative measure to MSD called the q-space inverse variance (QIV), which is a pseudo-diffusivity measure. Mathematically, the QIV is defined as


The QIV can thus be interpreted as the inverse of the “variance” of q (i.e. QIV = 1/left angle bracketq2right angle bracket). It is not a real variance in the statistical sense because E(q) does not constitute a probability density function. The QIV is not an arbitrary meausre, but related to the EAP in a manner analogous to which the MSD is related to the diffusion signal - in Appendix F, we will show that QIV-1=-2P(p)p=04π2. The QIV within the BFOR framework is (see Appendix F for derivation)


Fig. 15 displays an axial slice of the BFOR QIV, illusrating rich tissue/CSF contrast. The tissue/CSF contrast in the QIV is more enhanced than that of the MSD, and unlike the MSD, the QIV map also exhibits WM/GM contrast (the right and left putamen are visible in Fig. 15 but not in the MSD maps). According to Table 3, the QIV of the corpus callosum is about a third of that of the putamen.

Figure 15
Axial slice of BFOR QIV. Within the CSF regions, some voxels were zeroed out because they blew up upon the division operation in computing QIV.

6. Conclusion

We have introduced a new orthonormal basis to model the q-space diffusion signal and from which the EAP can be analytically reconstructed using hybrid, non-Cartesian sampling with multiple q-shell measurements. BFOR is a linear and efficient reconstruction based on heat equation estimation of the diffusion signal. Compared to DSI, BFOR employs much fewer diffusion measurements. Rotationally invariant q-space indices such as GFA, Po, and MSD can then be obtained using the derived EAP.

Research Highlights

  • Propose a new analytical diffusion propagator scheme called BFOR.
  • Derive rotationally invariant q-space indices using BFOR, like Po, MSD, & GFA.
  • Validate BFOR using HYDI-acquired datasets, both synthetic and in vivo.
  • Compare BFOR to existing analytical diffusion propagator schemes, like DPI and SPFI.
  • Explore ways to mitigate the effects of Gibbs ringing on propagator reconstruction.


The authors are thankful to Cheng Guan Koay and Steve Kecskemeti of the University of Wisconsin-Madison for insightful discussions on the estimation of the diffusion propagator.

Appendix A. Derivation of BFOR Signal Basis

We want to solve the following boundary value partial differential equation:

(Appendix A.1)

where we require λ > 0. Substituting the separable solution of the form

(Appendix A.2)

into Eq. (Appendix A.1), we obtain

(Appendix A.3)

where ΔLB=1sinθθ(sinθθ)+1sin2θ2φ2 is the Laplace-Bertrami operator and l is some real-valued constant.

We first solve for the second equation in (Appendix A.3):

(Appendix A.4)

The solutions to Eq. (Appendix A.4) are the spherical harmonics (SH)Ylm(θ,φ).

The second equation in (Appendix A.3) can be written as

(Appendix A.5)

Defining a new variable f(q)=π2λqF(q), we can transform Eq. (Appendix A.5) to

(Appendix A.6)

which is simply a sclaed version of the Bessel differential equation. The only bounded solution at the origin is given in terms of the Bessel funciton of the first kind as F(q)=Jl+1/2(λq). The solution to Eq. (Appendix A.5) is then f(q)=π2λqJl+1/2(λq)=jl(λq), where jl is the spherical Bessel function of the first kind and we invoke the relation jl(x)=π2xJl+1/2(x).

Imposing the boundary condition from Eq. (Appendix A.1), we have jl(λτ)=0. Defining αnl as the nth root of the lth order spherical Bessel function of first kind, then the eigenvalues are found to be -λnl=-αnl2τ2. Note that for l = 0, the roots are simply αn0 = .

Multiplying the spherical Bessel functions and the spherical harmonics together, we obtain the eigenfuctions (i.e. our orthonormal basis) to (Appendix A.1): Znlm(q,θ,φ)=jl(λnlq)Ylm(θ,φ). Thus, the complete set of solutions to Eq. (Appendix A.1) is

(Appendix A.7)

where N is the truncation order of the number of roots of spherical Bessel function and L the truncation order of the SH.

Now, to derive the diffusion signal, we make two important assumptions. First, we assume the diffusion signal E(q) is a solution to the heat equation:

(Appendix A.8)

where H(q) is simply the acquired signal. Second, we assume that the diffusion signal can be expressed as a linear combination of the orthonormal basis derived in Eq. (Appendix A.7):

(Appendix A.9)

Substituting (Appendix A.9) back into (Appendix A.8), we obtain

(Appendix A.10)

A unique solution exists if and only if gnl(t)=bnle-αnl2tτ2 (Chung et al., 2007), and so

(Appendix A.11)

Note that all constants are absorbed into Cnlm. In the following sections, we will use the SH basis Yj proposed in (Descoteaux et al., 2011).

Appendix B. Diffusion Signal at Origin

In diffusion weighted imaging (DWI), E(0) = 1. Thus, for our basis, we obtain the following identity:

(Appendix B.1)

which holds for any u within the unit sphere S2 (i.e. u [set membership] S2). In deriving Eq. (Appendix B.1), we invoked a basic property of the spherical Bessel function that


and the identity Y00=14π.

Appendix C. Derivation of Analytical BFOR EAP Solution

In the Theory section, we showed that

(Appendix C.1)

where Inl(j)(p)=0τq2jl(j)(αnl(j)qτ)jl(j)(2πqp)dq. We rewrite In(l(j)(p) in terms of the Bessel function of the first kind:

(Appendix C.2)

Recall the Bessel function of first kind Jk(ax), where k is some real-valued constant, satisfies

(Appendix C.3)

Thus, by definition of the Bessel function, Jl(j)+1/2(αnl(j)qτ) and Jl(j)+1/2(2πqp) satisfy

(Appendix C.4)

(Appendix C.5)

respectively. Multiplying Eq. (Appendix C.4) by Jl(j)+1/2(2πqp) and Eq. (Appendix C.5) by Jl(j)+1/2(αnl(j)qτ) and then subtracting, we obtain


Integrating the above from q = 0 to q = τ via integration by parts and noting that Jl(j)+1/2(αnl(j)) = 0, we have

(Appendix C.6)

The right side of Eq. (Appendix C.6) can be simplified via the Bessel recurrence relations ddxJk(x)=12(Jk-1(x)-Jk+1(x)):

(Appendix C.7)

Using the Bessel recurrence relation Jk+1(x)=2kxJk(x)-Jk-1(x), we obtain Jl(j)+3/2(αnl(j)) = −Jl(j)−1/2(αnl(j)), and so we can rewrite Eq. (Appendix C.7) as

(Appendix C.8)

Thus, substituting Eq. (Appendix C.8) back into Eq. (Appendix C.2), we obtain


and so the diffusion propagator is then

(Appendix C.9)

Appendix D. Derivation of BFOR Zero-Displacement Probability

We can derive Po by evaluating Eq. (Appendix C.9) at p = 0:

(Appendix D.1)

where we used the relation J-1/2(x)=2πxcos(x).

Appendix E. Relationship between MSD and Diffusion Signal in q-Space

We define the wave vector q as q = qxî + qyĵ + qzk and the radius vector in propagator space p as p = pxî + pyĵ + pzk. The norm of p is p=px2+py2+pz2.

Since the diffusion signal and EAP are FT pairs, then the inversion of Eq. (1) gives

(Appendix E.1)

Taking the second derivative of E(q) with respect to qx, qy, and qz gives


The sum of the derivatives is simply the Laplacian operator acting on E(q):

(Appendix E.2)

Note that the Laplacian of E(q) evaluated at q = 0 is related to the second moment of the EAP. Thus, the MSD is

(Appendix E.3)

For the case of DTI, where E(q, u) = e−4π2q2(Δ−δ/3)uTDu, Eq. (Appendix E.3) simplifies to the Einstein relation

(Appendix E.4)

Appendix E.1. Derivation of BFOR MSD

The BFOR signal basis is an eigenfunction of the Laplacian operator, with eigenvalues -αnl(j)2τ2. Hence

(Appendix E.5)

Evaluating the Laplacian of E(q) at q = 0 gives

(Appendix E.6)

Substituting Eq. (Appendix E.6) into (Appendix E.3), we obtain

(Appendix E.7)

where we dropped the smoothing term.

Appendix F. Relationship between QIV and EAP in q-Space

Using the definition of QIV, we have

(Appendix F.1)

The Dirac delta function is defined as

(Appendix F.2)

Taking the second derivative of δ(p) with respect to px, py, and pz gives


The sum of the derivatives is simply the Laplacian operator acting on δ(p):

(Appendix F.3)

Thus, Eq. (Appendix F.1) becomes

(Appendix F.4)

Since the Laplacian operator [nabla]2 is Hermitian and the EAP & δ(p) are real-valued, Eq. (Appendix F.4) can be equivalently stated as


Exploiting the property of the Dirac delta function that -f(x)δ(x)dx=f(0), we have

(Appendix F.5)

which is very similar in form to Eq. Appendix E.3. Thus, whereas the MSD directly varies with the Laplacian of the diffusion signal evaulated at the origin, the QIV inversely varies with the Laplacian of the EAP evaluated at the origin.

Appendix F.1. Derivation of BFOR QIV


The last integral can easily be solved via integration by parts, and so the QIV is

(Appendix F.6)


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Aganj I, Lenglet C, Sapiro G, Yacoub E, Ugurbil K, Harel N. Reconstruction of the orientation distribution function in single- and multiple-shell q-ball imaging within constant solid angle. Magn Reson Med. 2010;64:554–566. [PMC free article] [PubMed]
  • Alexander AL, Hasan KM, Lazar M, Tsuruda JS, Parker DL. Analysis of partial volume effects in diffusion-tensor MRI. Magn Reson Med. 2001;45:770–780. [PubMed]
  • Assaf Y, Ben Bashat D, Chapman J, Peled S, Biton IE, Kafri M, Segev Y, Hendler T, Korczyn AD, Graif M, Cohen Y. High b-value q-space analyzed diffusion-weighted MRI: application to multiple sclerosis. Magn Reson Med. 2002;47:115–126. [PubMed]
  • Assaf Y, Cohen Y. Non-mono-exponential attenuation of water and n-acetyl aspartate signals due to diffusion in brain tissue. J Magn Reson. 1998;131:69–85. [PubMed]
  • Assaf Y, Cohen Y. Assignment of the water slow-diffusing compartment in the central nervous system using q-space diffusion MRS: implications for fiber tract imaging. Magn Reson Med. 2000;43:191–199. [PubMed]
  • Assaf Y, Mayk A, Cohen Y. Displacement imaging of spinal cord using q-space diffusion-weighted MRI. Magn Reson Med. 2000;44:713–722. [PubMed]
  • Assemlal HE, Campbell J, Pike B, Siddiqi K. MICCAI. 2011. Apparent intravoxel fibre population dispersion (FPD) using spherical harmonics; pp. 157–165. [PubMed]
  • Assemlal HE, Tschumperlé D, Brun L. Efficient and robust computation of PDF features from diffusion MR signal. Med Image Anal. 2009a;13:715–729. [PubMed]
  • Assemlal HE, Tschumperlé D, Brun L. Evaluation of q-space sampling strategies for the diffusion magnetic resonance imaging. MICCAI. 2009b:482–489. [PubMed]
  • Bar-Shir A, Avram L, Ozarslan E, Basser PJ, Cohen Y. The effect of the diffusion time and pulse gradient duration ratio on the diffraction pattern and the structural information estimated from q-space diffusion MR. J Magn Reson. 2008;194:230–236. [PubMed]
  • Bar-Shir A, Ducan ID, Cohen Y. QSI and DTI of excised brain of the myelin-deficient rat. NeuroImage. 2009;48:109–116. [PubMed]
  • Basser PJ. Relationships between diffusion tensor and q-space MRI. Magn Reson Med. 2002;47:392–397. [PubMed]
  • Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophysical Journal. 1994;66:259–267. [PubMed]
  • Basser PJ, Pierpaoli C. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. J Magn Reson. 1996;111:209–219. [PubMed]
  • Callaghan PT. Principles of Nuclear Magnetic Resonance Microscopy. Oxford University Press; Oxford: 1991.
  • Canales-Rodriguez EJ, Iturria-Medina Y, Aleman-Gomez Y, Melie-Garcia L. Deconvolution in diffusion spectrum imaging. NeuroImage. 2010;50:136–149. [PubMed]
  • Canales-Rodriguez EJ, Melie-Garcia L, Iturria-Medina Y. Mathematical description of q-space in spherical coordinates: exact q-ball imaging. Magn Reson Med. 2009;61:1350–1367. [PubMed]
  • Caruyer E, Cheng J, Lenglet C, Sapiro G, Jiang T, Deriche R. Optimal design of multiple q-shells experiments for diffusion MRI. MICCAI Workshop on Computational Diffusion MRI-CDMRI’11.2011.
  • Caruyer E, Deriche R. Optimal regularization for MR diffusion signal. IEEE International Symposium on Biomedical Imaging.2012.
  • Cheng J, Ghosh A, Deriche R, Jiang T. Model-free, regularized, fast, and robust analytical orientation distribution function estimation. MICCAI. 2010a:648–656. [PubMed]
  • Cheng J, Ghosh A, Jiang T, Deriche R. Model-free and analytical EAP reconstruction via spherical polar Fourier diffusion MRI. MICCAI. 2010b:590–597. [PubMed]
  • Chung MK, Dalton KM, Shen L, Evans AC, Davidson RJ. Weighted Fourier series representation and its application to quantifying the amount of gray matter. IEEE Transac Med Imaging. 2007;26:566–581. [PubMed]
  • Clark CA, Le Bihan D. Water diffusion compartmentation and anisotropy at high b values in the human brain. Magn Reson Med. 2000;44:852–859. [PubMed]
  • Cohen Y, Assaf Y. High b-value q-space analyzed diffusion-weighted MRS and MRI in neuronal tissues: A technical review. NMR in Biomedicine. 2002;15:516–542. [PubMed]
  • Cook PA, Bai Y, Nedjati-Gilani S, Seunarine KK, Hall MG, Parker GJ, Alexander DC. Camino: open-source diffusion-MRI reconstruction and processing. Proc Intl Soc Mag Reson Med 2006
  • Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Regularized, fast, and robust analytical q-ball imaging. Magn Reson Med. 2007;58:497–510. [PubMed]
  • Descoteaux M, Deriche R, LeBihan D, Mangin JF, Poupon C. Multiple q-shell diffusion propagator imaging. Med Image Anal. 2011;15:603–621. [PubMed]
  • Du J, Goh A, Qiu A. Diffeomorphic metric mapping of high angular resolution diffusion imaging based on riemannian structure of orientation distribution functions. IEEE Transac Med Imaging. 2012;31:1021–1033. [PubMed]
  • Frank LR. Anisotropy in high angular resolution diffusion-weighted MRI. Magn Reson Med. 2001;45:935–939. [PubMed]
  • Ghosh A, Deriche R. Fast and closed-form ensemble-average-propagator approximation from the 4th-order diffusion tensor. IEEE International Symposium on Biomedical Imaging 2010
  • Hess CP, Mukherjee P, Han ET, Xu D, Vigneron DB. Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis. Magn Recon Med. 2006;56:104–117. [PubMed]
  • Hosseinbor AP, Chung MK, Wu YC, Alexander AL. Bessel Fourier orientation reconstruction: an analytical EAP reconstruction using multiple shell acquisitions in diffusion MRI. MICCAI. 2011:217–225. [PubMed]
  • Katkovnik V. A new method for varying adaptive bandwidth selection. IEEE Transac Sig Processing. 1999;47:2567–2571.
  • Maier SE, Vajapeyam S, Mamata H, Westin CF, Jolesz FA, Mulkern RV. Biexponential diffusion tensor analysis of human brain diffusion data. Magn Reson Med. 2004;51:321–330. [PubMed]
  • Mair RW, Sen PN, Hurlimann MD, Patz S, Cory DG, Walsworth RL. The narrow pulse approximation and long length scale determination in xenon gas diffusion NMR studies of model porous media. J Magn Reson. 2002;156:202–212. [PubMed]
  • Merlet S, Caruyer E, Deriche R. Impact of radial and angular sampling on multiple shells acquisition in diffusion MRI. MICCAI. 2011:113–121. [PubMed]
  • Mulkern RV, Gudbjartsson H, Westin CF, Zengingonul HP, Gartner W, RGC, Robertson R, WK, Schwartz R, Holtzman D, Jolesz FA, Maier SE. Multi-component apparent diffusion coefficients in human brain. NMR Biomed. 1999;12:51–62. [PubMed]
  • Niendorf T, Dijkhuizen RM, Norris DG, van Lookeren Campagne M, Nicolay K. Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging. Magn Reson Med. 1996;36:847–857. [PubMed]
  • Ozarslan E, Koay C, Shepherd TM, Blackband SJ, Basser PJ. Simple harmonic oscillator based reconstruction and estimation for three-dimensional q-space MRI. Proc Intl Soc Mag Reson Med 2009
  • Ozarslan E, Koay CG, Basser PJ. Simple harmonic oscillator based estimation and reconstruction for one-dimensional q-space MR. Proc Intl Soc Mag Reson Med 2008
  • Ozarslan E, Shepherd TM, Vemuri BC, Blackband SJ, Mareci TH. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT) NeuroImage. 2006;31:1086–1103. [PubMed]
  • Rathi Y, Michailovic O, Setsompop K, Bouix S, Shenton ME, Westin CF. Sparse multi-shell diffusion imaging. MICCAI. 2011:58–65. [PMC free article] [PubMed]
  • Stejskal E, Tanner J. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J Chem Phys. 1965;42:288–292.
  • Tournier JD, Calamante F, Gadian DG, Connelly A. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage. 2004;23:1176–1185. [PubMed]
  • Tuch DS. Q-ball imaging. Magn Reson Med. 2004;52:1358–1372. [PubMed]
  • Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Weeden VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn Reson Med. 2002;48:577–582. [PubMed]
  • Weeden VJ, Hagmann P, Tseng WYI, Reese TG, Weisskoff RM. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magn Reson Med. 2005;54:1377–1386. [PubMed]
  • Wiegell MR, Larsson HB, Wedeen VJ. Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology. 2000;217:897–903. [PubMed]
  • Wong STS, Roos MS. A strategy for sampling on a sphere applied to 3D selective RF pulse design. Magn Reson Med. 1994;32:778–784. [PubMed]
  • Wu YC, Alexander AL. Hybrid diffusion imaging. NeuroImage. 2007;36:617–629. [PMC free article] [PubMed]
  • Wu YC, Field AS, Alexander AL. Computation of diffusion function measures in q-space using magnetic resonance hybrid diffusion imaging. IEEE Transac Med Imaging. 2008;27:858–865. [PMC free article] [PubMed]
  • Wu YC, Field AS, Duncan ID, Samsonov AA, Kondo Y, Tudorascu D, Alexander AL. High b-value and diffusion tensor imaging in a canine model of dysmyelination and brain maturation. NeuroImage. 2011;58:829–837. [PMC free article] [PubMed]
  • Yeh CH, Cho KH, Lin HC, Wang JJ, Lin CP. Reduced encoding diffusion spectrum imaging implemented with a bi-Gaussian model. IEEE Transac Med Imaging. 2008;27:1415–1424. [PubMed]