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/mm

^{2}) 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**)/

*S*_{o} is the normalized

*q*-space diffusion signal,

*S*(

**q**) is the diffusion signal measured at position

**q** in

*q*-space, and

*S*_{o} is the baseline image acquired without any diffusion gradients (

*q* = 0). We denote

**q** =

*q*
**u**(

*θ*,

) and

**p** =

*p*
**r**(

*θ*′,

′), where

**u** and

**r** are 3D unit vectors. The wave vector

**q** is

**q** =

*γδ***G**/2

*π*, where

*γ* is the nuclear gyromagnetic ratio and

**G** =

*g***u** 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

*π*^{2}*q*^{2}(Δ −

*δ*/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 4

^{th} 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

^{2}*E* = 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/mm

^{2}. 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**.