An Optimal Dimensionality Sampling Scheme on the Sphere for Antipodal Signals In Diffusion Magnetic Resonance Imaging
Abstract
We propose a sampling scheme on the sphere and develop a corresponding spherical harmonic transform (SHT) for the accurate reconstruction of the diffusion signal in diffusion magnetic resonance imaging (dMRI). By exploiting the antipodal symmetry, we design a sampling scheme that requires the optimal number of samples on the sphere, equal to the degrees of freedom required to represent the antipodally symmetric bandlimited diffusion signal in the spectral (spherical harmonic) domain. Compared with existing sampling schemes on the sphere that allow for the accurate reconstruction of the diffusion signal, the proposed sampling scheme reduces the number of samples required by a factor of two or more. We analyse the numerical accuracy of the proposed SHT and show through experiments that the proposed sampling allows for the accurate and rotationally invariant computation of the SHT to near machine precision accuracy.
Alice P. Bates, Zubair Khalid and Rodney A. Kennedy
Email: {alice.bates, zubair.khalid rodney.kennedy}@anu.edu.au
diffusion magnetic resonance imaging; sampling; spherical harmonic transform; antipodal signal; unit sphere.
1 Introduction
Diffusion magnetic resonance imaging (dMRI) allows for the structure and connectivity of white matter in the brain to be determined noninvasively for the detection of neurological diseases and in preoperative planning [1, 2]. White matter has a fibrous structure; the myelin sheath around fibre axons hinders the movement of water molecules, constraining diffusion. The diffusion of water molecules in white matter can therefore be used to determine the location and orientation of fibres.
Acquisition and reconstruction of the diffusion signal from space measurements, where is the diffusion wave vector, is a central problem in dMRI. It is common for samples to be obtained on a single sphere or multiple spheres in space [3, 4, 5]. The number of measurements in space (images) that can be acquired in dMRI is severely limited due to the need for the scan time to be practical in a clinical setting. Hence, a space sampling scheme that allows for the accurate reconstruction of the diffusion signal using the minimum number of samples is desirable [6, 7].
The reconstruction of the diffusion signal can be carried out by expanding the signal in terms of spherical harmonics [8, 9, 10], which serve as natural basis functions for signals defined on the sphere [11]. By choosing a sufficiently large bandlimit, the diffusion signal can be represented in terms of a finite number of spherical harmonic coefficients, which are obtained using the spherical harmonic transform (SHT) calculated numerically from a finite number of measurements of the diffusion signal obtained at a constant space radius (on a sphere) [7, 12]. In order for the accurate reconstruction of the diffusion signal, a sampling scheme on the sphere for obtaining space measurements needs to allow for the accurate computation of the SHT.
1.1 Relation to Prior Work
Many sampling schemes on the sphere in dMRI use an antipodally symmetric sampling grid to reduce the number of measurements; the antipodal symmetry of the diffusion signal enables the signal to be evaluated on the sphere at locations antipodal to where measurements have been obtained [13, 9, 14, 15, 16, 17]. Most of these schemes are designed to obtain uniform (near equal area) sampling on the sphere to ensure that the reconstruction is rotationally invariant, which means that the accuracy does not vary significantly if the sampling grid (or diffusion signal) is rotated [3, 17]. Other sampling scheme designs that focus on the accurate computation of the SHT rather than geometric properties of the diffusion signal have also been used [5].
In order to acquire samples of the diffusion signal bandlimited at (formally defined in Section 2.1) over multiple spheres in space, an equiangular scheme [18] that requires asymptotically samples to accurately compute the SHT has been used [5]. In comparison to the geometric schemes, this scheme [18] has an isolatitude arrangement of samples (where samples along longitude are taken over isolatitude annuli) that enables the fast computation of the SHT. However, the samples in the scheme are nonuniform due to dense sampling near the poles, therefore it is considered that the accuracy of the SHT changes if the signal is rotated [17, 3]. The spherical design with uniform density sampling method in [17] has a uniform and antipodally symmetric arrangement of samples, and allows for the accurate computation of the SHT using least squares. The number of samples required by this scheme does not have a general formula depending on the bandlimit but as reported in [19] is greater than .
For the accurate computation of the SHT, the minimum number of samples, referred as the optimal spatial dimensionality, required by any sampling scheme is equal to the number of spherical harmonic coefficients required to represent a bandlimited signal in spherical harmonics basis [18, 20]. An isolatitude scheme that attains optimal spatial dimensionality ( samples) of an arbitrary bandlimited signal on the sphere has been developed [20], which allows the accurate and fast computation of the SHT. However, the samples in the sampling scheme are neither antipodally symmetric nor uniform (equal areal) by design. Due to antipodal symmetry of the diffusion signal, the bandlimited diffusion signal has degrees of freedom (see Section 3), which is referred to as the optimal spatial dimensionality of the diffusion (antipodal) signal. We have established that the existing sampling schemes require at least samples on the sphere for the accurate computation of the SHT (or the accurate reconstruction) of the diffusion signal.
1.2 Contributions
In this work, we address whether it is possible to develop a sampling scheme on the sphere (at a constant space radius) for measuring the diffusion signal that 1) attains the optimal spatial dimensionality, i.e., requires samples, 2) allows for the accurate computation of the SHT up to and beyond practically relevant bandlimits, 3) permits rotationally invariant reconstruction, and 4) enables fast computation of the SHT with lowest known complexity. In addressing these questions, we develop a sampling scheme on the sphere and corresponding SHT that, in contrast to previous schemes, attains the optimal spatial dimensionality for the accurate reconstruction of the bandlimited diffusion signal at a constant space radius (an antipodally symmetric signal). The proposed scheme is antipodal, but unlike existing schemes that only focus on the geometric properties of the diffusion signal, the SHT associated with the proposed scheme is computationally efficient and accurate, with reconstruction error near machine precision () for the bandlimits of interest in dMRI. Though the scheme is nonuniform by design, we show that the numerical accuracy of the scheme is the same when the diffusion signal on the sphere is rotated, demonstrating that the scheme allows for the rotationally invariant reconstruction of the diffusion signal.
This paper is organised as follows. The mathematical preliminaries necessary to understand this work are contained in Section 2. In Section 3, we state the problem under consideration. The sampling scheme and the SHT for measurement and reconstruction of the diffusion signal are presented in Section 4. This scheme is then evaluated in terms of numerical accuracy and rotational invariance in Section 5. Finally, concluding remarks are made in Section 6.
2 Mathematical Preliminaries
2.1 Signals on the Sphere and Spherical Harmonics
Square integrable complex functions on the unit sphere have the form , where the angles colatitude and longitude parameterise a point on . The inner product of two functions , defined on is given by [11]
(1) 
where denotes the complex conjugate and is the differential area element on the sphere. With the inner product in (1), the set of square integrable complex valued functions on forms a Hilbert space, denoted by . The inner product in (1) induces a norm . We refer to functions with finite induced norm as signals on the sphere.
The spherical harmonic functions (spherical harmonics for short) , defined for integer degree and integer order , form a complete orthonormal set of basis functions [11, 21]. The signal is said to be bandlimited at degree if the signal can be completely represented in spherical harmonics basis functions with and .
2.2 Diffusion Signal on the Sphere
The diffusion signal at a fixed space radius can be represented as a bandlimited signal on the sphere [17, 12]. Furthermore, the diffusion signal, denote by , is antipodally symmetric, with . Since for even and for odd , the expansion of in spherical harmonic basis only includes even degree spherical harmonics, that is,
(2) 
where denotes the spherical harmonic coefficient of degree and order , and is calculated using the SHT, given by
(3) 
The spherical harmonic coefficients form the spectral domain representation of . The reconstruction of the signal from its spherical harmonic coefficients as given in (2) is referred as the inverse SHT. The bandlimit required to accurately represent the diffusion signal depends on the space radius [5], and bandlimits up to are typically used [12, 7].
3 Problem Formulation
The diffusion signal has only even order spherical harmonic coefficients due to its antipodal symmetry. The number of spherical harmonics required to represent in (2) is [7, 22], which also represents the optimal dimensionality attainable by any sampling scheme that allows the accurate computation of the SHT of the diffusion signal. We customise the recently developed sampling scheme [20], which requires the optimal number of samples for the accurate computation of the SHT of an arbitrary bandlimited signal (without antipodal symmetry), for the acquisition and reconstruction of the diffusion signal in dMRI.
3.1 Optimal Dimensionality Sampling Scheme
The optimal dimensionality sampling scheme in [20] has an isolatitude sampling grid, with samples taken over isolatitude rings. The locations along where the isolatitude rings are placed are stored in the vector defined as
(4) 
For the accurate computation of the SHT, there needs to be at least samples along longitude in the ring placed at . The isolatitude rings are placed along such that the SHT can be computed accurately (see [20] for further details). For the accurate computation of the SHT, the number of samples required by the optimal spatial dimensionality scheme is , which is equal to the number of degrees of freedom required to represent an arbitrary bandlimited signal defined on the sphere.
3.2 Research Questions
In this work, we address the following questions:

How can we customise the structure of the optimal dimensionality sampling scheme to make it suitable for the acquisition of measurements of the antipodal diffusion signal?

Can we exploit the antipodal symmetry property of the diffusion signal to develop a scheme that requires only rather than samples while still allowing for the accurate computation of the SHT ?
4 Proposed Sampling Scheme
4.1 Proposed Sampling Scheme — Structure and Design
The diffusion signal on the sphere is antipodally symmetric, with . Hence, by measuring the signal at a location , we also know the value of the signal at a second location . We use this property to design the sampling scheme in order to reduce the number of measurements of the diffusion signal that need to be taken.
With this consideration, we propose to place the isolatitude rings in pairs such that they are antipodal to one another. The vector , given in (4), describing the location of the isolatitude rings is
(5) 
We shortly present the location of these samples along colatitude. Since we require at least samples along in the isolatitude ring placed at for the accurate computation of the SHT, we propose equiangular sampling along longitude, with th sample location, denoted by , in the ring placed at is given by
(6) 
Remark 1.
The samples in the proposed scheme are structured in a way that the samples in the ring are antipodal to the samples in the ring for . Therefore, we only need to take the measurements over the rings for . The values of the signal over the remaining samples can be determined by using the antipodal symmetry of the diffusion signal and exploiting the structure of the proposed sampling scheme. As an example, Fig. 1 shows the proposed sampling scheme for .
Remark 2.
Since the measurements are only required to be taken over rings in the proposed scheme, the total number of measurements required by the proposed sampling scheme is
(7) 
which also represent the degrees of freedom required to represent bandlimited antipodal signal. Thus, the proposed scheme attains the optimal spatial dimensionality.
4.2 Spherical Harmonic Transform (SHT) — Formulation
We here present a variant of the method developed in [20], to compute the SHT of the diffusion signal sampled using the proposed sampling scheme. We define an indexed vector for order , which contains the last points in the vector . Also define a vector
with
(8) 
for each order and each , where denotes scaled associated Legendre functions.
If for each is computed, the spherical harmonic coefficients of order can be recovered from (8) by setting up a system of linear equations [20], given by,
(9) 
where is a matrix of size defined as
and is a vector composed of spherical harmonic coefficients of order , given by
(10) 
4.3 Spherical Harmonic Transform — Computation
The spherical harmonic coefficients of order contained in a vector can be computed accurately by inverting (9), provided is chosen such that is wellconditioned and for each is computed accurately. Since we have considered in the design of the proposed scheme that the ring placed at contains at least samples points along longitude, for each can be exactly computed by evaluating the integral equation in (8) as a summation over samples along longitude [20].
We also require the sampling points along colatitude to be chosen such that the matrix is wellconditioned for each . Since the locations of the isolatitude rings along colatitude in the vector are required to be in pairs, as given in (5), we are only required to obtain measurements of the diffusion signal at positions () along colatitude. We define a set of equiangular samples along colatitude given by
(11) 
and propose the following method to determine the optimal ordering of samples in the vector :

Choose farthest from the poles, which is a natural choice for the ring of (the largest number of) samples along . Then, to satisfy the antipodal nature of the scheme, .

For each , choose and from the set , given in (11), that minimises the sum of the condition numbers of the matrices and .

Choose or .
Such placement of samples along colatitude ensures that is wellconditioned, resulting in the accurate computation of the proposed SHT as we demonstrate in the next section.
We have considered the selection of sample positions along colatitude from the set of equiangular samples as it has been shown in [20] that the selection from the equiangular set results in more accurate computation of the SHT. Furthermore, the equiangular placement of samples along longitude, given in (6), permits the use of the fast Fourier transform for the computation of (8). We note that the proposed sampling scheme can be easily customised for the case when the samples along colatitudes and/or longitude are not equiangular.
5 Analysis of Proposed Sampling Scheme
5.1 Numerical Accuracy
For the accurate reconstruction of the diffusion signal, using its expansion in the spherical harmonic basis in (2), the accuracy of the SHT given in (3) is of the upmost importance. A sampling scheme is numerically accurate if the inverse SHT of a bandlimited signal followed by the SHT gives an error between the original and reconstructed signal on the order of the numerical precision [20, 18].
We carry out the following experiment to test the numerical accuracy of the proposed sampling scheme and the corresponding SHT. We first obtain a bandlimited antipodally symmetric test signal in the spectral domain by generating spherical harmonic coefficients for with real and imaginary parts uniformly distributed in the interval ( for ). The inverse SHT is then used to obtain in the spatial domain over the proposed sampling grid (described in Section 4), followed by the forward SHT to compute the spherical harmonic coefficients of the reconstructed signal, . The experiment is repeated 10 times and the average values for the maximum error and the mean error , given by
(12)  
(13) 
are recorded and plotted in Fig. 2 for bandlimits (the bandlimits of interest in dMRI). It is evident that, although the errors and grow as the bandlimit increases, the errors are on the order of numerical precision (less than ), showing that the proposed sampling scheme and SHT allows for the accurate reconstruction of any bandlimited antipodally symmetric signal on the sphere for bandlimits .
5.2 Rotational Invariance
Many existing sampling schemes on the sphere focus on uniform sampling of the diffusion signal on the sphere to ensure that the reconstruction accuracy is rotational invariant i.e. does not change significantly when the signal or sampling grid is rotated [17]. The proposed sampling scheme is not uniform by design, however it does not have dense sampling on any region of the sphere. The rotational invariance of the reconstruction accuracy of the SHT was tested using the numerical accuracy experiment (described in Section 5.1) performed on 5 rotated versions of the antipodally symmetric test signal . For each rotation, we randomly obtain Euler angles from the uniform distributions, where and , and apply the rotation using the convention [11]. Fig. 3 shows the mean reconstruction error , averaged over 10 experiments, for the original test signal and 5 rotated versions of the test signal, for bandlimits . It can be seen that the reconstruction errors are on the same order of magnitude regardless of the rotation angle, demonstrating that the accuracy of the scheme does not depend on the angle of rotation.
6 Conclusions
In this work we have proposed a new sampling scheme on the sphere for the accurate reconstruction of the diffusion signal in dMRI. The proposed scheme exploits the antipodal symmetry of the diffusion signal and attains an optimal number of samples , given by the degrees of freedom required to represent the diffusion signal of bandlimit , in contrast to the existing schemes which require at least samples. The smaller number of samples required by this scheme will enable a reduction in the scan time required to obtain sufficient measurements for the accurate reconstruction of the diffusion signal. The proposed sampling scheme enables accurate and rotational invariant computation of the SHT of the diffusion signal to near machine precision accuracy. Applying the scheme to dMRI data and extending the proposed scheme to multiple shell sampling is the subject of our future work.
Footnotes
 thanks: This work is supported by the Australian Research Council’s Discovery Projects funding scheme (Project no. DP150101011).
References
 P. J. Basser, J. Mattiello, and D. LeBihan, “MR diffusion tensor spectroscopy and imaging,” Biophys. J., vol. 66, no. 1, pp. 259–267, Jan. 1994.
 H. Jahanian, A. YazdanShahmorad, and H. SoltanianZadeh, “4D wavelet noise suppression of MR diffusion tensor data,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., ICASSP’2008, Las Vegas, Nevada, Mar. 2008, pp. 509–512.
 J. Cheng, D. Shen, and P.T. Yap, “Designing single and multipleshell sampling schemes for diffusion MRI using spherical code,” in Med. Image Comput. Comput. Assist. Interv., MICCAI’2014, Boston, Massachusetts, 2014, vol. 8675, pp. 281–288.
 W. Ye, S. Portnoy, A. Entezari, S. J. Blackband, and B. C. Vemuri, “An efficient interlaced multishell sampling scheme for reconstruction of diffusion propagators,” IEEE Trans. Med. Imag., vol. 31, no. 5, pp. 1043–1050, May 2012.
 A. Daducci, J. D. McEwen, D. V. D. Ville, J. P. Thiran, and Y. Wiaux, “Harmonic analysis of spherical sampling in diffusion MRI,” in Proc. 19th Ann. Meet. Int. Soc. Magn. Reson. Med., Jun. 2011.
 H. E. Assemlal, D. Tschumperlé, and L. Brun, “Evaluation of qspace sampling strategies for the diffusion magnetic resonance imaging,” in Med. Image Comput. Comput. Assist. Interv., MICCAI’2009, London, UK, 2009, vol. 12, no. Pt 2, pp. 406–414.
 J.D. Tournier, F. Calamante, and A. Connelly, “Determination of the appropriate b value and number of gradient directions for highangularresolution diffusionweighted imaging,” NMR Biomed., vol. 26, no. 12, pp. 1775–1786, Dec. 2013.
 J. Haldar and R. Leahy, “New linear transforms for data on a fourier 2sphere with application to diffusion MRI,” in Proc. 9th IEEE Int. Symp. Biomed. Imag., ISBI’2012, Barcelona, Spain, May 2012, pp. 402–405.
 M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche, “Regularized, fast, and robust analytical qball imaging,” Magn. Reson. Med., vol. 58, no. 3, pp. 497–510, Sep. 2007.
 A. W. Anderson, “Measurement of fiber orientation distributions using high angular resolution diffusion imaging,” Magn. Reson. Med., vol. 54, no. 5, pp. 1194–1206, Nov. 2005.
 R. A. Kennedy and P. Sadeghi, Hilbert Space Methods in Signal Processing. Cambridge, UK: Cambridge University Press, Mar. 2013.
 C. P. Hess, P. Mukherjee, E. T. Han, D. Xu, and D. B. Vigneron, “Qball reconstruction of multimodal fiber orientations using the spherical harmonic basis,” Magn. Reson. Med., vol. 56, no. 1, pp. 104–117, Jul. 2006.
 D. K. Jones, M. A. Horsfield, and A. Simmons, “Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging,” Magn. Reson. Med., vol. 42, no. 3, pp. 515–525, Sep. 1999.
 D. S. Tuch, “Qball imaging,” Magn. Reson. Med., vol. 52, no. 6, pp. 1358–1372, 2004.
 L. FejesTóth, “On the densest packing of spherical caps,” Amer. Math. Monthly, vol. 56, pp. 330–331, 1949.
 N. G. Papadakis, C. D. Murrills, L. D. Hall, C. L. Huang, and T. Adrian Carpenter, “Minimal gradient encoding for robust estimation of diffusion anisotropy,” Magn. Reson. Med., vol. 18, no. 6, pp. 671–679, Jul. 2000.
 E. Caruyer, C. Lenglet, G. Sapiro, and R. Deriche, “Design of multishell sampling schemes with uniform coverage in diffusion MRI,” Magn. Reson. Med., vol. 69, no. 6, pp. 1534–1540, Jun. 2013.
 J. D. McEwen and Y. Wiaux, “A novel sampling theorem on the sphere,” IEEE Trans. Signal Process., vol. 59, no. 12, pp. 5876–5887, Dec. 2011.
 E. Caruyer, “Qspace diffusion MRI: acquisition and signal processing,” these, Université Nice Sophia Antipolis, Jul. 2012.
 Z. Khalid, R. A. Kennedy, and J. D. McEwen, “An optimaldimensionality sampling scheme on the sphere with fast spherical harmonic transforms,” IEEE Trans. Signal Process., vol. 62, no. 17, pp. 4597–4610, Sep. 2014.
 J. J. Sakurai, Modern Quantum Mechanics, 2nd ed. Reading, MA: Addison Wesley Publishing Company, Inc., 1994.
 E. Caruyer and R. Deriche, “Diffusion MRI signal reconstruction with continuity constraint and optimal regularization,” Med. Image Anal., vol. 16, no. 6, pp. 1113–1120, Aug. 2012.
 C. G. Koay, “A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere,” J. Comput. Sci., vol. 2, no. 4, pp. 377–381, Dec. 2011.