Optimal-Dimensionality Sampling on the Sphere: Improvements and Variations

# Optimal-Dimensionality Sampling on the Sphere: Improvements and Variations

## Abstract

For the accurate representation and reconstruction of band-limited signals on the sphere, an optimal-dimensionality sampling scheme has been recently proposed which requires the optimal number of samples equal to the number of degrees of freedom of the signal in the spectral (harmonic) domain. The computation of the spherical harmonic transform (SHT) associated with the optimal-dimensionality sampling requires the inversion of a series of linear systems in an iterative manner. The stability of the inversion depends on the placement of iso-latitude rings of samples along co-latitude. In this work, we have developed a method to place these iso-latitude rings of samples with the objective of improving the well-conditioning of the linear systems involved in the computation of the SHT. We also propose a multi-pass SHT algorithm to iteratively improve the accuracy of the SHT of band-limited signals. Furthermore, we review the changes in the computational complexity and improvement in accuracy of the SHT with the embedding of the proposed methods. Through numerical experiments, we illustrate that the proposed variations and improvements in the SHT algorithm corresponding to the optimal-dimensionality sampling scheme significantly enhance the accuracy of the SHT.

\IEEEoverridecommandlockouts{IEEEkeywords}

unit sphere, sampling, spherical harmonic transform, optimal-dimensionality, condition number minimization, harmonic analysis

## 1 Introduction

Signal analysis on spherical bodies has widespread applications in the fields of cosmology, geodesy, geomagnetics, acoustics and computer graphics [1, 2, 3, 4, 5, 6]. Data measured over the surface of a spherical object, i.e., in the spatial domain, can be transformed to the harmonic domain using the spherical harmonic transform (SHT) which is the analogue in spherical geometry of the renowned Fourier transform in Euclidean geometry [7]. Sampling schemes utilized for computing SHTs are categorized as theoretically exact, accurate or approximate [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. In this work, we consider those schemes which enable exact or accurate computation of the SHT of band-limited signals. Different sampling schemes have different spatial dimensionality defined as the number of sample points needed to accurately or exactly compute the SHT and thus capture the information content of band-limited signals. For the computation of SHTs of a signal band-limited at  (defined in Section 2.2), the optimal spatial dimensionality attainable by any sampling scheme on the sphere is , which is equal to the degrees of freedom of the band-limited signal in harmonic space.

Driscoll and Healy [8] developed an exact method to compute the SHT of a signal, that is band-limited at , which requires (asymptotically, as ) equiangular samples on the sphere, where the complexity of most stable algorithm to compute SHT is . In comparison, the sampling scheme presented by McEwen and Wiaux [18] requires equiangular samples to exactly compute the SHT with complexity . The Gauss-Legendre sampling scheme [18, 19] also requires for exact computation of the SHT, where the complexity to compute the SHT is . To the best of our knowledge, there does not exist any theoretically exact sampling scheme with dimensionality less than . On the other hand, the SHT can also be computed approximately using the least-squares based method proposed by Sneeuw [9], which, although requiring samples, becomes inaccurate and computationally inefficient scaling as for large band-limits.

Recently, an optimal-dimensionality sampling scheme has been proposed in [20] for the accurate computation of the SHT of band-limited signals using only samples. Optimal-dimensionality sampling has been customized to serve the needs of applications in acoustics [4] and diffusion MRI [21]. Although the SHT associated with this sampling scheme requires the optimal number of samples, it has computational complexity of . The computation of the SHT for optimal-dimensionality sampling involves inversion of a series of systems of linear equations. For accurate inversion of these systems, a condition number minimization method has been proposed in [20] to determine the locations of samples.

This paper aims to improve the accuracy of the SHT associated with the optimal-dimensionality sampling scheme. We serve this objective by developing a new method for the placement of samples and proposing a variation in the computation of the SHT. We develop a method, referred to as the elimination method, for the placement of iso-latitude rings of samples such that the condition number (ratio of the largest to the smallest singular value) of the matrices used in the computation of the SHT is minimized. Due to the iterative nature of the resulting SHT algorithm, the error builds up in the computation of the SHT. To resolve this issue, we also propose a multi-pass SHT algorithm which iteratively reduces the residual between the given signal and the reconstructed signal. We also analyze the changes in the complexity of the SHT with the use of these methods. Through numerical experiments, we demonstrate the improvement in accuracy with the use of the proposed methods. The remainder of the paper is structured as follows. We present the necessary mathematical background in Section 2 before reviewing the optimal-dimensionality sampling scheme in Section 3. Section 4 presents the proposed developments and also contains the accuracy analysis. Concluding remarks are then made in Section 5.

## 2 Mathematical Background

### 2.1 Signals on the Sphere

Let denote a complex-valued, square integrable function on the unit sphere , where and denote the co-latitude and longitude respectively. The space formed by these functions is a Hilbert space, denoted by , equipped with the following inner product given by

 ⟨f,h⟩≜∫S2f(θ,ϕ)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯h(θ,ϕ)sinθdθdϕ, (1)

for any two functions . In (1), denotes the complex conjugate operation, is the differential surface element and is an integral over the whole sphere. The inner product in (1) induces a norm , and signals with finite induced norm are referred to as “signals on the sphere”.

### 2.2 Harmonic Domain Representation

Signals can be transformed to the harmonic domain using the natural basis – spherical harmonic basis functions (or simply spherical harmonics). Spherical harmonics, denoted by for integer degree and integer order , are defined as

 Ymℓ(θ,ϕ)≜√2ℓ+14π(ℓ−m)!(ℓ+m)!Pmℓ(cosθ)eimϕ, (2)

where is the associated Legendre function [7]. Any function can be expanded in terms of spherical harmonics as

 f(θ,ϕ) =∞∑ℓ=0ℓ∑m=−ℓ(f)mℓYmℓ(θ,ϕ). (3)

Here denotes the spherical harmonic coefficient of degree and order and is given by the spherical harmonic transform (SHT) as

 (f)mℓ≜⟨f,Ymℓ⟩=∫S2f(θ,ϕ)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ymℓ(θ,ϕ)sinθdθdϕ. (4)

The synthesis equation, (3), to reconstruct the signal from its spherical harmonic coefficients is referred to as inverse SHT. A signal is said to band-limited if = 0 for , where is the band-limit of the signal, and can be expressed in terms of spherical harmonics as

 f(θ,ϕ) =L−1∑ℓ=0ℓ∑m=−ℓ(f)mℓYmℓ(θ,ϕ). (5)

The signals, band-limited at , form an dimensional subspace of , which we denote by .

## 3 Problem Formulation

### 3.1 Optimal Dimensionality Sampling on the Sphere

The optimal-dimensionality sampling scheme on the sphere requires (optimal number) samples to accurately compute the SHT for a signal with band-limit  [20]. In this scheme, iso-latitude rings are placed on the sphere at locations (to be explained shortly) given in vector , defined as

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

The ring placed at contains equiangular points along longitude .

### 3.2 SHT Formulation

For a signal sampled using optimal-dimensionality sampling scheme, we define a vector , for every as

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

where for each is given as

 Gm(θk)≜∫2π0f(θk,ϕ)e−imϕdϕ =2πL−1∑ℓ=m(f)mℓ~Pmℓ(θk). (8)

Here denotes scaled associated Legendre functions. The second equality in (8) is obtained by using (2) and (5) and employing the orthogonality of complex exponentials. By defining another vector as

 fm=[(f)m|m|,(f)m|m|+1,…,(f)mL−1]T, (9)

containing the spherical harmonic coefficients of order , we formulate a linear system given as

 gm=Pmfm, (10)

where the is an matrix with elements given by

 Pm(i,j)=~Pm|m|+j−1(θ|m|+i−1). (11)

### 3.3 Problem Under Consideration

For each order , the spherical harmonic coefficients contained in can be recovered by solving the linear system given in (10). Computation of the SHT, i.e., the computation of spherical harmonic coefficients of the signal sampled according to the optimal-dimensionality sampling scheme, involves the inversion of a series of linear systems formed by the matrix  (defined in (11)) for  [20]. A condition number minimization method has been proposed in [20] to determine the locations of the iso-latitude rings indexed in (6) such that the matrix for each is well-conditioned and the SHT can be accurately computed. With an objective to improve the accuracy of the SHT, we consider the problem of determining the locations of iso-latitude rings of samples which reduce (improve) the condition number (ratio of the largest to the smallest eigenvalue) of the matrices . The vector containing the locations of iso-latitude rings, initially in the ascending order of co-latitude angle, is re-ordered such that every matrix has minimum condition number. To further increase the accuracy of the SHT, we also propose a multi-pass SHT algorithm which iteratively reduces the error between the given signal (samples in spatial domain) and the signal synthesized using the computed spherical harmonic coefficients.

## 4 Optimized Samples Placement and Multi-Pass SHT

With an objective to improve the accuracy of the SHT using optimal number of samples, we first present our proposed method for the placement of the iso-latitude rings of samples and later we present iterative method for the computation of the SHT.

### 4.1 Condition Number Minimization

The recovery of for each using (10) requires inversion of the matrix for each . For accurate computation of the SHT, it is therefore necessary that the matrix is invertible and well-conditioned. Since is a matrix of associated Legendre polynomials of order and degrees evaluated at , its accurate inversion depends on the locations of the iso-latitude rings indexed in (6). To determine the locations of the iso-latitude rings, we propose a condition number minimization technique, herein referred to as the elimination method, for the construction of the vector .

Let be a set of equiangular co-latitude angles between and defined as

 Ω≜{π(2t+1)2L−1},t=0,1,…,L−1. (12)

For , the matrix is formed by inserting all elements of set in (11) and has dimension . Since , for =1, requires co-latitude angles, we eliminate one element, say , from the set and calculate the condition number, denoted by , of using all possible combinations of residual elements . The element , whose elimination results in the lowest condition number of , is then selected as the first element of the vector. The set is then updated as . The same procedure is carried out for the construction of the vector for which we summarize below in the form of an algorithm.

The vector constructed using the proposed elimination method is optimized in a sense that it generates matrices of lower condition number as compared to the optimal-dimensionality sampling scheme. This improvement in the condition number comes from the fact that the proposed elimination method has choices for such that the condition number of matrix is minimized. In contrast, the method proposed in [20] has choices for the selection of and minimization of the condition number of matrix the . As an illustration, the condition number of the matrix using the proposed optimized placement of iso-latitude rings and the design proposed in [20] is plotted in Fig. 1 for band-limit . We also plot the maximum of the condition number obtained for different band-limits in Fig. 2. It is evident that the proposed elimination method improves the well-conditioning of the systems involved in the computation of the SHT algorithm associated with the optimal-dimensionality sampling on the sphere.

### 4.2 Multi-pass SHT

The spherical harmonic coefficients of a band-limited signal sampled according to the optimal-dimensionality sampling scheme are computed iteratively for each order in a sequence . The SHT is inherently iterative in nature as the spherical harmonic coefficients of order are used in the computation of the SHT of order . Consequently, the error propagates and builds up in the iterative computation of spherical harmonic coefficients. To reduce this error building-up, we propose a multi-pass SHT algorithm which iteratively improves the accuracy of the SHT.

For a signal sampled by the optimal-dimensionality sampling scheme, the spherical harmonic coefficients can be accurately computed by the algorithm presented in [20]1. We define the residual (error) between the signal and the signal synthesized from the recovered spherical harmonic coefficients as

 rk(θ,ϕ)=f(θ,ϕ)−L−1∑ℓ=0ℓ∑m=−ℓ(~fk)mℓYmℓ(θ,ϕ) (13)

where denotes the spherical harmonic coefficient computed using the proposed SHT algorithm and  (indicating the number of times the transform has been carried out). Once residual is computed, we use the SHT algorithm to compute its spherical harmonic coefficients, denoted by , which we use to update as

 (~fk+1)mℓ=(~fk)mℓ+(~rk)mℓ. (14)

We propose to iteratively use (13) and (14) to compute for , until the following stopping criterion is met

 max|rk+1(θ,ϕ)|>max|rk(θ,ϕ)|, (15)

where is taken over the samples of the sampling scheme. Since the proposed method requires to compute the SHT multiple times, we refer to the proposed method for the computation of spherical harmonic coefficients as the multi-pass SHT. Later, we illustrate that the proposed method significantly improves the accuracy of the SHT.

### 4.3 Computational Complexity Analysis

Here we briefly discuss the computational complexity of the proposed elimination method for the placement of iso-latitude rings and the multi-pass SHT algorithm. The elimination method has the computational complexity of . However, it only needs to be run once for the determination of for each . Furthermore, we note that the complexity of the method presented in [20] for the placement of samples is also . For the optimal-dimensionality sampling scheme, the SHT can be computed with complexity of . For the proposed multi-pass SHT algorithm, the complexity scales with the number of iterations, denoted by , needed for the convergence of error. In the next section, we provide examples to illustrate that the proposed multi-pass SHT algorithm converges quickly in number of iterations.

### 4.4 Accuracy Analysis

In this section, we analyse the accuracy of the proposed multi-pass SHT algorithm of a band-limited signal evaluated using the optimal-dimensionality sampling scheme with iso-latitude rings placed using the proposed elimination method. Comparison between the proposed developments and the SHT proposed in [20] has been carried out through numerical experiments. In our experiment, we first take a band-limited signal by randomly generating its spherical harmonic coefficients . The real and imaginary parts of the coefficients are uniformly distributed in . Using inverse SHT, we obtain the signal in the spatial domain, that is, over the samples of the optimal-dimensionality sampling scheme (proposed sampling or [20]). We then apply the SHT presented in [20] and the proposed multi-pass SHT algorithm to recover the spherical harmonic coefficients, denoted by and respectively. We conduct experiments for 10 different signals to obtain the average value of the maximum error between reconstructed and original spherical harmonic coefficients defined as

 Emax ≜max|(~f)mℓ−(f)mℓ|, (16) Ekmax ≜max|(~fk)mℓ−(f)mℓ|, (17)

which we plot for band-limits in Fig. 3, where it can be observed that the proposed multi-pass SHT algorithm and optimized placement of samples results in the more accurate computation of the SHT.

We also analyse the convergence of the multi-pass SHT algorithm and the improvement in the accuracy of the SHT enabled by the proposed multi-pass SHT algorithm. We plot the maximum absolute error for band-limits and in Fig. 4, where it can be observed that the proposed multi-pass SHT improves the accuracy of SHT and converges (quickly) in number of iterations.

## 5 Conclusions

In this work, we have proposed variations in the spherical harmonic transform (SHT) associated with the optimal-dimensionality sampling scheme which consist of iso-latitude rings of samples and enables accurate computation of the SHT of band-limited signals using the optimal number of samples given by the degrees of freedom required to represent a band-limited signal in harmonic space. We have presented the elimination method for the iterative placement of iso-latitude rings of samples. The proposed placement reduces the condition number of matrices involved in the computation of SHT and consequently improves the accuracy of the SHT. We have also proposed the multi-pass SHT algorithm which iteratively reduces the residual between the given signal and the reconstructed signal and therefore improves the overall accuracy of the SHT. We have analyzed the changes in the computational complexity and improvement in accuracy with the use of proposed variations in the computation of the SHT. We have also conducted numerical experiment to illustrate the improvement in accuracy enabled by the proposed methods.

### Footnotes

1. SHT can be computed accurately for band-limited signals sampled over optimal-dimensionality sampling scheme [20] using the MATLAB  based package Novel Spherical Harmonic Transform (NSHT) publicly available at www.zubairkhalid.org/nsht.

### References

1. M. Tegmark, M. A. Strauss, M. R. Blanton, K. Abazajian, S. Dodelson, H. Sandvik, X. Wang, D. H. Weinberg, I. Zehavi, N. A. Bahcall et al., “Cosmological parameters from sdss and wmap,” Physical Review D, vol. 69, no. 10, p. 103501, 2004.
2. M. A. Wieczorek and F. J. Simons, “Minimum variance multitaper spectral estimation on the sphere,” J. Fourier Anal. Appl., vol. 13, no. 6, pp. 665–692, 2007.
3. F. Lowes, “Spatial power spectrum of the main geomagnetic field, and extrapolation to the core,” Geophysical Journal International, vol. 36, no. 3, pp. 717–730, 1974.
4. A. P. Bates, Z. Khalid, and R. A. Kennedy, “Novel sampling scheme on the sphere for head-related transfer function measurements,” IEEE/ACM Trans. Audio Speech Language Process., vol. 23, pp. 1068–1081, June 2015.
5. W. Zhang, M. Zhang, R. A. Kennedy, and T. D. Abhayapala, “On high-resolution head-related transfer function measurements: An efficient sampling scheme,” IEEE Trans. Acoust., Speech, Signal Process., vol. 20, no. 2, pp. 575–584, 2012.
6. R. Ramamoorthi and P. Hanrahan, “A signal-processing framework for reflection,” ACM Transactions on Graphics (TOG), vol. 23, no. 4, pp. 1004–1042, 2004.
7. R. A. Kennedy and P. Sadeghi, Hilbert Space Methods in Signal Processing.   Cambridge, UK: Cambridge University Press, Mar. 2013.
8. J. R. Driscoll and D. M. Healy, Jr., “Computing Fourier transforms and convolutions on the 2-sphere,” Adv. Appl. Math., vol. 15, no. 2, pp. 202–250, Jun. 1994.
9. N. Sneeuw, “Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective,” Geophys. J. Int., vol. 118, no. 3, pp. 707–716, Sep. 1994.
10. M. J. Mohlenkamp, “A fast transform for spherical harmonics,” J. Fourier Anal. and Appl., vol. 5, no. 2-3, pp. 159–184, 1999.
11. K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere,” Astrophys. J., vol. 622, pp. 759–771, Apr. 2005.
12. R. G. Crittenden and N. G. Turok, “Exactly azimuthal pixelizations of the sky,” Arxiv preprint astro-ph/9806374, 1998.
13. D. M. Healy, Jr., D. Rockmore, P. J. Kostelec, and S. S. B. Moore, “FFTs for the 2-sphere - improvements and variations,” J. Fourier Anal. and Appl., vol. 9, pp. 341–385, 2003.
14. J. Keiner, S. Kunis, and D. Potts, “Efficient reconstruction of functions on the sphere from scattered data,” J. Fourier Anal. Appl., vol. 13, no. 4, pp. 435–458, May 2007.
15. P. J. Kostelec and D. N. Rockmore, “FFTs on the rotation group,” J. Fourier Anal. and Appl., vol. 14, pp. 145–179, Apr. 2008.
16. J. D. McEwen, “Fast, exact (but unstable) spin spherical harmonic transforms,” All Res. J. Phys., vol. 1, no. 1, 2011.
17. K. M. Huffenberger and B. D. Wandelt, “Fast and exact spin-s spherical harmonic transforms,” Astrophys. J. Supp., vol. 189, no. 2, pp. 255–260, Aug. 2010.
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. W. Skukowsky, “A quadrature formula over the sphere with application to high resolution spherical harmonic analysis,” J. Geodesy, vol. 60, no. 1, pp. 1–14, Mar. 1986.
20. Z. Khalid, R. A. Kennedy, and J. D. McEwen, “An optimal-dimensionality sampling scheme on the sphere with fast spherical harmonic transforms,” IEEE Transactions on Signal Processing, vol. 62, no. 17, pp. 4597–4610, 2014.
21. A. P. Bates, Z. Khalid, and R. A. Kennedy, “An optimal dimensionality sampling scheme on the sphere with accurate and efficient spherical harmonic transform for diffusion MRI,” IEEE Signal Process. Lett., vol. 23, pp. 15–19, January 2016.
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