An Optimal Dimensionality Sampling Scheme on the Sphere for Antipodal Signals in diffusion magnetic resonance imaging

# 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 band-limited 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.

\ninept\name

Alice P. Bates, Zubair Khalid and Rodney A. Kennedy 1 \addressResearch School of Engineering, The Australian National University, Canberra, ACT 2601, Australia
Email: {alice.bates, zubair.khalid rodney.kennedy}@anu.edu.au

{keywords}

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 non-invasively 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 band-limit, 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 band-limited 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 iso-latitude arrangement of samples (where samples along longitude are taken over iso-latitude annuli) that enables the fast computation of the SHT. However, the samples in the scheme are non-uniform 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 band-limit 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 band-limited signal in spherical harmonics basis [18, 20]. An iso-latitude scheme that attains optimal spatial dimensionality ( samples) of an arbitrary band-limited 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 band-limited 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 band-limits, 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 band-limited 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 band-limits of interest in dMRI. Though the scheme is non-uniform 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 co-latitude and longitude parameterise a point on . The inner product of two functions , defined on is given by [11]

 ⟨f,g⟩≜∫S2f(θ,ϕ)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯g(θ,ϕ)sinθdθdϕ, (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 band-limited 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 band-limited 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,

 d(θ,ϕ)=L−1∑ℓ=0ℓevenℓ∑m=−ℓ(d)mℓYmℓ(θ,ϕ),L odd, (2)

where denotes the spherical harmonic coefficient of degree and order , and is calculated using the SHT, given by

 (d)mℓ≜∫S2d(θ,ϕ)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ymℓ(θ,ϕ)sinθdθdϕ. (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 band-limit required to accurately represent the diffusion signal depends on the -space radius [5], and band-limits 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 band-limited 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 iso-latitude sampling grid, with samples taken over iso-latitude rings. The locations along where the iso-latitude rings are placed are stored in the vector defined as

 θ≜[θ0,θ1,…,θL−1]T. (4)

For the accurate computation of the SHT, there needs to be at least samples along longitude in the ring placed at . The iso-latitude 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 band-limited 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 iso-latitude rings in pairs such that they are antipodal to one another. The vector , given in (4), describing the location of the iso-latitude rings is

 θ≜[0,…,π−θL−3,θL−3,π−θL−1,θL−1]T,L odd. (5)

We shortly present the location of these samples along co-latitude. Since we require at least samples along in the iso-latitude 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

 ϕnk≜⎧⎨⎩2kπ2n+1,n=0,2,…,L−1,k∈[0,2n],π(2k+1)2n+3,n=1,3,…,L−2,k∈[0,2(n+1)]. (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

 L−1∑n=0neven(2n+1)=(L+1)L2=NO, (7)

which also represent the degrees of freedom required to represent band-limited 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

 gm≡Gm(θm)≜[Gm(θ|m|),Gm(θ|m|+1),…,Gm(θL−1)]T,

with

 Gm(θn)≜∫2π0f(θn,ϕ)e−imϕdϕ=2πL−1∑ℓ=|m|(f)mℓ˜Pmℓ(θn), (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,

 gm=PmLfm,|m|≤L, (9)

where is a matrix of size defined as

 PmL≜2π⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝˜Pm|m|(θ|m|)˜Pm|m|+1(θ|m|)⋯˜PmL−1(θ|m|)˜Pm|m|(θ|m|+1)˜Pm|m|+1(θ|m|+1)⋯˜PmL−1(θ|m|+1)⋮⋮⋱⋮˜Pm|m|(θL−1)˜Pm|m|+1(θL−1)⋯˜PmL−1(θL−1)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

and is a vector composed of spherical harmonic coefficients of order , given by

 fm=[(f)m|m|,(f)m|m|+1,…,(f)mL−1]T. (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 well-conditioned 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 co-latitude to be chosen such that the matrix is well-conditioned for each . Since the locations of the iso-latitude rings along co-latitude 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 co-latitude. We define a set of equiangular samples along co-latitude given by

 Θ={π(2t+1)2L−1},t=0,1,…,L−12, (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 co-latitude ensures that is well-conditioned, 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 co-latitude 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 co-latitudes and/or longitude are not equiangular.

###### Remark 3.

The computational complexity of the SHT for the proposed sampling scheme is ([20], which is much smaller than the complexity of least squares methods () used by many existing sampling schemes [3, 17, 23].

## 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 band-limited 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 band-limited 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

 Emax ≜max|(ft)mℓ−(fr)mℓ|, (12) Emean ≜1L2L−1∑ℓ=0ℓ∑m=−ℓ|(ft)mℓ−(fr)mℓ|, (13)

are recorded and plotted in Fig. 2 for band-limits (the band-limits of interest in dMRI). It is evident that, although the errors and grow as the band-limit 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 band-limited antipodally symmetric signal on the sphere for band-limits .

### 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 band-limits . 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 band-limit , 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

1. thanks: This work is supported by the Australian Research Council’s Discovery Projects funding scheme (Project no. DP150101011).

### References

1. 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.
2. H. Jahanian, A. Yazdan-Shahmorad, and H. Soltanian-Zadeh, “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.
3. J. Cheng, D. Shen, and P.-T. Yap, “Designing single- and multiple-shell 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.
4. W. Ye, S. Portnoy, A. Entezari, S. J. Blackband, and B. C. Vemuri, “An efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators,” IEEE Trans. Med. Imag., vol. 31, no. 5, pp. 1043–1050, May 2012.
5. 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.
6. H. E. Assemlal, D. Tschumperlé, and L. Brun, “Evaluation of q-space 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.
7. J.-D. Tournier, F. Calamante, and A. Connelly, “Determination of the appropriate b value and number of gradient directions for high-angular-resolution diffusion-weighted imaging,” NMR Biomed., vol. 26, no. 12, pp. 1775–1786, Dec. 2013.
8. J. Haldar and R. Leahy, “New linear transforms for data on a fourier 2-sphere with application to diffusion MRI,” in Proc. 9th IEEE Int. Symp. Biomed. Imag., ISBI’2012, Barcelona, Spain, May 2012, pp. 402–405.
9. M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche, “Regularized, fast, and robust analytical q-ball imaging,” Magn. Reson. Med., vol. 58, no. 3, pp. 497–510, Sep. 2007.
10. 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.
11. R. A. Kennedy and P. Sadeghi, Hilbert Space Methods in Signal Processing.   Cambridge, UK: Cambridge University Press, Mar. 2013.
12. C. P. Hess, P. Mukherjee, E. T. Han, D. Xu, and D. B. Vigneron, “Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis,” Magn. Reson. Med., vol. 56, no. 1, pp. 104–117, Jul. 2006.
13. 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.
14. D. S. Tuch, “Q-ball imaging,” Magn. Reson. Med., vol. 52, no. 6, pp. 1358–1372, 2004.
15. L. Fejes-Tóth, “On the densest packing of spherical caps,” Amer. Math. Monthly, vol. 56, pp. 330–331, 1949.
16. 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.
17. 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.
18. 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.
19. E. Caruyer, “Q-space diffusion MRI: acquisition and signal processing,” these, Université Nice Sophia Antipolis, Jul. 2012.
20. Z. Khalid, R. A. Kennedy, and J. D. McEwen, “An optimal-dimensionality sampling scheme on the sphere with fast spherical harmonic transforms,” IEEE Trans. Signal Process., vol. 62, no. 17, pp. 4597–4610, Sep. 2014.
21. J. J. Sakurai, Modern Quantum Mechanics, 2nd ed.   Reading, MA: Addison Wesley Publishing Company, Inc., 1994.
22. 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.
23. 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.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters