Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2800042

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. RELATED WORK
- III. DEFINITIONS
- IV. INVERTIBILITY CONDITIONS
- V. CONSTRUCTING SELF-INVERTIBLE MULTISCALE FILTER BANKS
- VI. EXPERIMENTS
- VII. DISCUSSION AND CONCLUSION
- REFERENCES

Authors

Related links

IEEE Trans Image Process. Author manuscript; available in PMC 2009 December 30.

Published in final edited form as:

PMCID: PMC2800042

NIHMSID: NIHMS88084

Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA

Boon Thye Thomas Yeo: ude.tim.liasc@samohty; Wanmei Ou: ude.tim.liasc@iemnaw; Polina Golland: ude.tim.liasc@anilop

The publisher's final edited version of this article is available at IEEE Trans Image Process

See the article "Shape Analysis with Overcomplete Spherical Wavelets" in *Med Image Comput Comput Assist Interv* on page 468.

See other articles in PMC that cite the published article.

The theories of signal sampling, filter banks, wavelets, and “overcomplete wavelets” are well established for the Euclidean spaces and are widely used in the processing and analysis of images. While recent advances have extended some filtering methods to spherical images, many key challenges remain. In this paper, we develop theoretical conditions for the invertibility of filter banks under continuous spherical convolution. Furthermore, we present an analogue of the Papoulis generalized sampling theorem on the 2-Sphere. We use the theoretical results to establish a general framework for the design of invertible filter banks on the sphere and demonstrate the approach with examples of self-invertible spherical wavelets and steerable pyramids. We conclude by examining the use of a self-invertible spherical steerable pyramid in a denoising experiment and discussing the computational complexity of the filtering framework.

Multiscale filtering methods, such as wavelets [7] and “overcomplete wavelets” [9], [23], [25] have many applications in feature detection, compression, and denoising of planar images. Extending the theories and the methods of filtering to spherical images promises similar benefits in the fields that give rise to such images, including shape analysis in computer vision [5], illumination computation in computer graphics [21], [22], cosmic background radiation analysis in astrophysics [26], [29] and brain cortical surface analysis in medical imaging [32], [33].

We consider a two-stage filtering framework [Fig. 1(a)], conceptually equivalent to the usual Euclidean filtering framework except the planar images and filters are replaced by spherical ones. We can think of the first set of filters as analysis filters which project the input image onto the space spanned by the analysis filters. A reconstructed image is then obtained by passing the intermediate outputs through the second layer of filters. Fig. 1(b) shows a modification of the framework in Fig. 1(a) that introduces sampling between the first and the second layers of filters. The sampling is useful if one is interested in processing the outputs of the analysis filters before passing them through the synthesis filters.

In this work, we analyze the relationship between the reconstructed image and the original image, and establish conditions under which the reconstructed image is the same as the original one, i.e., conditions for invertibility. We demonstrate the use of our results on continuous invertibility and generalized sampling for designing filter banks that enable explicit control of both analysis and synthesis filters. We illustrate the framework by creating examples of self-invertible spherical wavelets and steerable pyramids.

Self-invertible filter banks employ identical filters for analysis and synthesis. Self-invertibility is desirable for image manipulation in the wavelet domain, leading to an intuitive notion that a convolution coefficient corresponds to the contribution of the corresponding filter to the reconstructed signal. Without self-invertibility, the effects of nonlinear processing of wavelet coefficients could propagate to spatial locations and frequencies other than those which were used to compute the coefficients [23]. To the best of our knowledge, this is the first approach demonstrated on a sphere that enables the design of self-invertible filter banks.

In the Euclidean domain, convolution is computed efficiently using the fast Fourier transform (FFT) [6]. Once we move to the sphere, FFT must be replaced with an alternative efficient method for computing convolutions. An original algorithm for axisymmetric convolution kernels on the sphere was derived in [11] and was recently extended to arbitrary functions [28], [30]. These results allow us to efficiently compute the outputs of the first set of filters. Unfortunately, they do not apply to the inverse convolution with the second layer of filters because the outputs of the first layer of filters are in general not spherical images as we will see in Section III.

In the past decade, there has been much work on extending the general paradigm of linear filtering to the spherical domain [1], [4], [8], [12], [13], [17], [21], [22], [26], [29]. For example, the lifting scheme in [21] and [22] adopts a nonparametric approach to computing the wavelet decomposition of arbitrary meshes by generalizing the standard two-scale relation of Euclidean wavelets. This method enables a multiscale representation of the original mesh (image) with excellent compression and speed performance. However, the lifting wavelets are not overcomplete, i.e., exactly one wavelet coefficient is created per sample point, causing difficulties in designing filters for oriented feature detection and invariance to rotation [33].

A similar problem in the Euclidean domain leads to the invention of “overcomplete wavelets,” such as steerable pyramids [14], [23]. A group theoretic formulation of overcomplete continuous spherical wavelets is proposed in [1]. In particular, it can be shown that the stereographic projection of an admissible planar wavelet to the sphere is also admissible under the group theoretic framework, providing a straightforward framework for the design of analysis filters for specific features of interest, such as oriented edges [29].

In the group theoretic approach, defining the mother wavelet completely determines the analysis and synthesis filters. However, while the analysis filters are related by stereographic dilation, the synthesis filters are in general not related by dilation. In fact, the support of corresponding analysis and synthesis filters is guaranteed to be the same in frequency domain but not in the spatial domain. Bogdanova *et al.* [4] discretize the group theoretic wavelets, providing a sampling guarantee for the framework of Fig. 1(b) for the restricted class of axisymmetric filters. An axisymmetric spherical function is one which is symmetrical about the north pole. This work is, therefore, the most similar to ours. In contrast, we study general filter banks, without any restriction on the relationships among the cascade of filters. We derive the analogue of the Papoulis generalized sampling theorem [18] on the sphere, applicable to both axisymmetric and nonaxisymmetric filters.

Driscoll and Healy [11] provide the equivalent of the Nyquist–Shannon sampling theorem on the sphere. While the Nyquist–Shannon sampling theorem provides reconstruction guarantees for bandlimited signals in Euclidean space under perfect sampling (convolution with a delta function), the Papoulis generalized sampling theorem provides guarantees for bandlimited signals sampled via convolutions with kernels of sufficient bandwidth.

An earlier version of this work was first presented at the International Conference on Image Processing [31]. In this paper, we include proofs of the invertibility conditions and demonstrate the generation of self-invertible spherical steerable pyramids. In Section III, we introduce the notation used throughout the paper. In Section IV, we present the main theoretical contributions of this paper: continuous invertibility and the generalized sampling theorem. We propose a procedure for generating self-invertible multiscale filter banks on the sphere in Section V. In Section VI, we illustrate the procedure to design wavelets and steerable pyramids and employ a steerable pyramid in denoising. We conclude with the discussion of future research and outstanding challenges in the proposed framework.

To summarize, our contributions are as follows.

- We present theoretical conditions for the invertibility of axisymmetric and nonaxisymmetric filter banks under continuous spherical convolution.
- We present a generalized sampling theorem of signals on the 2-Sphere for both axisymmetric and nonaxisymmetric filter banks. This generalizes the works of Bogdanova
*et al.*[4] and Starck*et al.*[26] to nonaxisymmetric filters and opens a way for nonlinear processing of the wavelet coefficients generated from general filter banks. - We present a mechanism for generating invertible, as well as self-invertible, wavelets, and steerable pyramids. We provide an analysis of the computational complexity of the filtering framework.

Let *x*(θ,ϕ) *L*^{2}(*S*^{2}) be a square-integrable function on the 2-D unit sphere, where (θ,ϕ) are the spherical coordinates. Suppose *P* = (θ,ϕ) is a point on the sphere. Then, θ [0, π] is the co-latitude, which is the angle between the positive *z*-axis (north pole) and the vector corresponding to *P*. ϕ [0, π] is the longitude and is taken to be the angle between the positive *x*-axis and the projection of *P* onto the *x* – *y* plane. ϕ is undefined on the north and south poles.

The spherical harmonics ${Y}_{l}^{m}\left(\theta ,\varphi \right)$ [20] form an orthonormal set of basis functions for *L*^{2}(*S*^{2}): i.e.,

$$x\left(\theta ,\varphi \right)={\displaystyle \sum _{l=0}^{\infty}{\displaystyle \sum _{\left|m\right|\le l}{x}^{l,m}{Y}_{l}^{m}}}\left(\theta ,\varphi \right)$$

(1)

where is *x ^{l,m}* the spherical harmonic coefficient of degree

$${x}^{l,m}={\displaystyle {\int}_{{S}^{2}}\phantom{\rule{thinmathspace}{0ex}}x\left(\theta ,\varphi \right){Y}_{l}^{m*}\left(\theta ,\varphi \right)}d\mathrm{\Omega}$$

(2)

where *d*Ω = sin θ*d*θ*d*ϕ and * denotes complex conjugation. We call ${Y}_{l}^{m}\left(\theta ,\varphi \right)$ a spherical harmonic of degree *l* and order *m*. We note that for axisymmetric functions (independent of ϕ), only the order 0 harmonics are nonzero. A more detailed background of spherical harmonics is found in Appendix A.

We choose to parameterize rotations on the sphere by the Euler angles, α,β,γ (α [0,2π],β [0,π],γ [0,2π]). The rotation operator *D*(α,β,γ) first rotates a function by γ about the *z*-axis [Fig. 2(b)], then by β about the *y*-axis [Fig. 2(c)], and finally by α about the *z*-axis [Fig. 2(d)]. The direction of positive rotation follows the right-hand screw rule. The three angles specify an element of the rotation group *SO*(3) and provide a natural parametrization of convolution on the sphere. The effects of rotation on the spherical harmonic coefficients of a function is expressible in terms of the so-called Wigner-D functions. The Wigner-D functions form an irreducible representation of the rotation group [20]. Appendix A provides the explicit expressions for the Wigner-D functions.

Rotation via euler angles (α, β, γ). (a) Original spherical image. (b) Rotation by γ about *z*-axis. (c) Rotation by β about *y*-axis. (d) Rotation by α about *z*-axis.

On the plane, convolution is defined in terms of the inner product between two functions translated relative to each other, and is parameterized by the amount of translation. On the sphere, it is more natural to talk about rotation rather than translation, and, therefore, spherical convolution is parameterized by rotation. Given a spherical image *x*(θ,ϕ) and a spherical filter (θ,ϕ), their spherical convolution

$$Y(\alpha ,\beta ,\gamma )={\displaystyle {\int}_{{S}^{2}}[D(\alpha ,\beta ,\gamma )\tilde{h}]*\left(\theta ,\varphi \right)x\left(\theta ,\varphi \right)d\mathrm{\Omega}}$$

(3)

is a function of *L*^{2}(*SO*(3)) rather than *L*^{2}(*S*^{2}). This definition of convolution is identical to that in [28], [30], although [30] calls it directional correlation.

By convention, we shall consider the center (origin) of a spherical filter to be at the north pole (θ = 0). Then intuitively,*y*(α, β, γ) is the inner product between the re-oriented filter *D*(α, β, γ) [e.g., Fig. 2(d)] and the spherical image. In other words, we obtain *y*(α, β, γ) by first re-orienting the spherical filter by a rotation of γ about the *z*-axis (center still at north pole), then bringing the center of the filter to the point (β, α) of the spherical image, and performing an inner product between the image and filter. Therefore, *y*(α, β, γ) is the correlation of the rotated version of with *x*, or the projection coefficient of *x* onto [*D*(α, β, γ)]. In the case of the filter shown in Fig. 2, a high value of *y*(α, β, γ) would imply a presence of an oriented edge at spherical coordinate (β, α) with orientation γ.

We note that the notion of orientation is inherently local here because a continuous unit-norm vector field does not exist on the sphere. Therefore, it is not meaningful to claim an existence of an edge of orientation γ at location (β, α) without specifying a local coordinate system. In our case, we can define such a local coordinate system by first specifying one at the north pole. Our choice of parameterizing rotation via the Euler angles then induces a local coordinate system everywhere, except the south pole.

For axisymmetric filters, (θ,ϕ) = (θ), the rotation by γ about the *y*-axis has no effect, i.e., *y*(α, β, γ) = *y*(α, β) is a spherical image parametrized by θ = β,ϕ = α.

To project a function in *L*^{2}(*SO*(3))onto *L*^{2}(*S*^{2}), we define the inverse convolution of a spherical filter *h*(θ,ϕ) with *y*(α, β, γ) *L*^{2}(*SO*(3)) to be

$${\widehat{x}}_{h}\left(\theta ,\varphi \right)={\displaystyle {\int}_{SO(3)}[D(\alpha ,\beta ,\gamma )h]\left(\theta ,\varphi \right)y\left(\alpha ,\beta ,\gamma \right)d\rho}$$

(4)

where the integration is over the Euler angles: *d*ρ = sin β*d*α*d*β*d*γ. Our definition of inverse convolution generalizes the spherical convolution between two functions on *S*^{2} defined by Driscoll and Healy [11], and is identical to theirs if *y*(α, β, γ) = *y*(β, α), i.e., a spherical image. Our inverse convolution operation is similar to the transpose convolution operation defined in [2], [28]. However, the transpose convolution uses the uniform measure *d*α*d*β*d*γ, ignoring the intrinsic nonuniform component sin β of the measure on *SO*(3).

We can think of the inverse convolution in the following intuitive way. The reconstructed value at a given (θ,ϕ) is obtained by summing (integrating) the contributions of the rotated reconstruction filters *h*, centered at all possible positions (β, α) and oriented by all possible angles γ, where the weights of the contributions are given by the convolution outputs (projection coefficients on the corresponding input filters).

We distinguish between the forward convolution (3) and the inverse convolution (4) because the forward convolution combines two spherical images to give a function on *SO*(3), while the inverse convolution combines a spherical image and a function on *SO*(3) to give a spherical image. Both definitions of convolutions generalize the concept of convolution in Euclidean spaces.

When using a filter bank of *N* analysis-synthesis filter pairs [Fig. 1(a)], the reconstructed signal is obtained by summing the response of all filter pairs

$$\widehat{x}\left(\theta ,\varphi \right)={\displaystyle \sum _{n=1}^{N}{\displaystyle {\int}_{SO(3)}[D(\alpha ,\beta ,\gamma ){h}_{n}]\left(\theta ,\varphi \right){y}_{n}\left(\alpha ,\beta ,\gamma \right)d\rho}}$$

(5)

which is analogous to the definition in [1], with integration over scale replaced by summation over the filter index.

In the Euclidean case, we typically discretize both the input images and the convolution outputs. When working on the sphere, we choose to keep the image domain continuous by working with spherical harmonic coefficients rather than sample values, because this allows us to exploit efficient algorithms for spherical convolution [28], [30]. Since no uniform sampling grid exists on the sphere, performing convolution completely by quadrature would be slow. This is because under each rotation of the filter relative to the spherical image, we would need to re-sample (or re-interpolate) the filter or the image.

For *L*^{2}(*SO*(3))(or equivalently, in the spherical wavelet domain), continuous representation is possible through series of complex exponentials [28] or Wigner-D functions [30]. However, both the complex exponentials and the Wigner-D functions have global support. Therefore, in applications where we want to modify the image in the wavelet domain, manipulating the series coefficients would be tantamount to simultaneously altering all the wavelet coefficients, defeating the purpose of the wavelet decomposition, which is to provide localized control in both spatial and frequency domain. To avoid this, we sample the output of the continuous convolution *y*(α,β,γ) to create its discrete counterpart *y*(α* _{j}*, β

$$\begin{array}{l}{\widehat{x}}_{h}\left(\theta ,\varphi \right)={\displaystyle \sum _{j=0}^{J-1}{\displaystyle \sum _{s=0}^{S-1}{\displaystyle \sum _{k=0}^{K-1}{w}_{j,s,k}[D({\alpha}_{j},{\beta}_{s},{\gamma}_{k})h]\left(\theta ,\varphi \right)}}}\\ \hfill \times y\left({\alpha}_{j},{\beta}_{s},{\gamma}_{k}\right)\end{array}$$

(6)

which includes sampling-dependent quadrature weights *w _{j,s,k}*, introduced so that the discrete case converges to the continuous case as the number of samples increases. This definition allows for an easy transfer of continuous filtering theory to its discrete analogue. In Section IV, we show that “good” choices of

Similar to the continuous case [cf. (4)], the signal reconstructed through *N* analysis-synthesis filter pairs is defined as a sum of contributions of all filter pairs

$$\begin{array}{l}\widehat{x}\left(\theta ,\varphi \right)={\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{w}_{j,s,k,n}}}}}\\ \hfill \text{\hspace{1em}\hspace{1em}}\times [D({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}]\left(\theta ,\varphi \right){y}_{n}\left({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}\right).\end{array}$$

(7)

The sampling grid and the quadrature weights now depend on *n* since different filters in the filter bank might use different sampling schemes.

In this section, we present the main theoretical contributions of our work.

Let ${\{{\tilde{h}}_{n},{h}_{n}\}}_{n=1}^{N}$ be an analysis-synthesis filter bank. Then for any spherical image *x* *L*^{2}(*S*^{2}) and its corresponding reconstructed image *$\widehat{x}$*

$${\widehat{x}}^{l,m}={x}^{l,m}\frac{8{\pi}^{2}}{2l+1}\left\{{\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\prime =-l}^{l}\left[{h}_{n}^{l,m\prime}\right]}}\phantom{\rule{thinmathspace}{0ex}}{\left[{\tilde{h}}_{n}^{l,m\prime}\right]}^{*}\right\}$$

(8)

where *x ^{l,m}* and

Appendix B presents the proof of Theorem 4.1. To draw an analogy with the Euclidean case, we call

$${H}_{\tilde{h},h}(l)=\frac{8{\pi}^{2}}{2l+1}{\displaystyle \sum _{n-1}^{N}{\displaystyle \sum _{m\prime =-l}^{l}\left[{h}_{n}^{l,m\prime}\right]}}\phantom{\rule{thinmathspace}{0ex}}{\left[{\tilde{h}}_{n}^{l,m\prime}\right]}^{*}$$

(9)

the frequency response of the analysis-synthesis filter bank. Note that the degree *l* spherical harmonics coefficients of the reconstructed image are affected only by the degree *l* spherical harmonic coefficients of the filters. However, the degree *l* order *m* spherical harmonic coefficient of the reconstructed signal is affected by all the orders of degree *l* spherical harmonic coefficients of the filters. In contrast, on the plane, the frequency response is simply the sum of products of the Fourier coefficients of the analysis and the synthesis filters

$$\begin{array}{cc}\mathcal{F}\left\{\widehat{x}\right\}({s}_{1},{s}_{2})=& \mathcal{F}\left\{x\right\}({s}_{1},{s}_{2})\hfill \\ \hfill & \times {\displaystyle \sum _{n=1}^{N}F\left\{{h}_{n}\right\}({s}_{1},{s}_{2}{)}^{*}\mathcal{F}\{{\tilde{h}}_{n}\}({s}_{1},{s}_{2}{)}^{*}}\end{array}$$

(10)

where {*$\widehat{x}$*}(*s*_{1},*s*_{2}), {*x*}(*s*_{1},*s*_{2}), {}(*s*_{1},*s*_{2}) and {*h _{n}*}(

Let ${\{{\tilde{h}}_{n},{h}_{n}\}}_{n=1}^{N}$ be an analysis-synthesis filter bank. Then for any spherical image *x* *L*^{2}(*S*^{2}) and its corresponding reconstructed image *$\widehat{x}$*

$$\begin{array}{l}{\widehat{x}}^{l,m}={x}^{l,m}\phantom{\rule{thinmathspace}{0ex}}\text{for all}\phantom{\rule{thinmathspace}{0ex}}(l,m)\phantom{\rule{thinmathspace}{0ex}}\text{iff}{\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\prime =-l}^{l}\left[{h}_{n}^{l,m\prime}\right]}}\phantom{\rule{thinmathspace}{0ex}}{\left[{\tilde{h}}_{n}^{l,m\prime}\right]}^{*}=\frac{2l+1}{8{\pi}^{2}}\\ \hfill \text{for all}\phantom{\rule{thinmathspace}{0ex}}l\phantom{\rule{thinmathspace}{0ex}}s.t.\phantom{\rule{thinmathspace}{0ex}}{x}^{l,m}\ne 0.\hfill \end{array}$$

(11)

We note that the corollary is easily satisfied if there is no constraint on the relationships among the cascade of filters: given a set of analysis filters * _{n}*, there are in general multiple sets of synthesis filters that can achieve invertibility. For example, we can define the synthesis filters to be

$${[{L}_{\psi}{\tilde{h}}_{n}]}^{l,m}=\{\begin{array}{cc}\frac{1}{{H}_{\tilde{h},\tilde{h}}(l)}{\tilde{h}}_{n}^{l,m},& \text{for}\phantom{\rule{thinmathspace}{0ex}}{H}_{\tilde{h},\tilde{h}}(l)>0\\ 0,\hfill & \text{otherwise}\hfill \end{array}$$

(12)

and *H _{,}* (

We now define * _{n}* (

Let ${\{{\tilde{h}}_{n},{h}_{n}\}}_{n=1}^{N}$ be a filter bank with * _{n}* (and, thus,

- α
_{j,n}= 2π*j*/(+_{n}*L*+ 1) for_{n}*j*= 0, 1, … ,(+_{n}*L*);_{n} - γ
_{k,n}= 2π*k*/(*Õ*+_{n}*O*+ 1) for_{n}*k*= 0, 1, …, (*Õ*+_{n}*O*);_{n} *w*and β_{s,n}_{s,n}and are the quadrature weights and knots such thatfor$$\begin{array}{l}{\displaystyle {\int}_{0}^{\pi}{d}_{mm\prime}^{l}}(\beta ){d}_{mm\prime}^{l\prime}(\beta )\text{sin}(\beta )d(\beta )\\ \hfill \text{\hspace{1em}\hspace{1em}}={\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{mm\prime}^{l}({\beta}_{s,n}){d}_{mm\prime}^{l\prime}({\beta}_{s,n})}\end{array}$$(13)*l*≤*L*≤_{n},l′, where ${d}_{\mathit{\text{mm}}\prime}^{l}(\beta )\phantom{\rule{thickmathspace}{0ex}}\text{and}\phantom{\rule{thickmathspace}{0ex}}{d}_{\mathit{\text{mm}}\prime}^{l\prime}(\beta )$ are the Wigner-d functions;_{n}*w*= 4π_{j,s,k,n}^{2}*w*/((_{s,n}+_{n}*L*+1)(_{n}*Õ*+_{n}*O*+ 1)). Then_{n}$${\widehat{x}}^{l,m}={x}^{l,m}\frac{8{\pi}^{2}}{2l+1}\left\{{{\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\prime =-l}^{l}\left[{h}_{n}^{l,m\prime}\right]\phantom{\rule{thinmathspace}{0ex}}\left[{\tilde{h}}_{n}^{l,m\prime}\right]}}}^{*}\right\}.$$(14)

Appendix A provides the definitions and the explicit expressions for the Wigner-d functions. Note that there is an implicit link between the number of samples *S _{n}* and the bandwidths

The measures corresponding to α and γ are constants in *SO*(3), just like in the Euclidean space. We, therefore, use uniform sampling for these parameters in our work. For discrete planar convolution, it is customary to have no weights (or rather, unit weights). On the sphere, however, the nonuniform measure on β, sin β*d*β, presents challenges for sampling. If we are simply interested in convergence, then setting *w _{j,s,k}* = (2π/

The continuous invertibility corollary and the generalized sampling theorem imply the following:

Consider a filter bank of filters with finite maximum spherical harmonic degrees. There exists quadrature schemes such that the filter bank is invertible over a particular frequency range under the discrete spherical convolution if the filter bank is also invertible over the same frequency range under the continuous spherical convolution.

Because functions in *L*^{2}(*S*^{2}) have finite energy, their spherical harmonic coefficients must necessarily decay to zero. Therefore, we can reasonably assume that the filters of the filter bank are of finite bandwidth as required by Theorem 4.3 by representing the filters with a finite number of coefficients up to an arbitrary prespecified precision. We will, therefore, focus on constructing invertible filter banks for continuous convolution.

The discrete invertibility corollary opens a way for nonlinear processing of the wavelet coefficients generated from general filter banks. Applications may include compression, denoising and image enhancement. The corollary also enables the perfect reconstruction of an original signal sampled with equipment that introduces blurring during the acquisition process. For example, the first layer of filters * _{n}* could be the blurring kernels of a set of

In this section, we show how to use the continuous invertibility corollary to generate self-invertible multiscale filter banks. The optimization framework presented here can be easily adapted to design other types of filter banks by altering the structure of the optimization problem according to an application’s needs.

For self-invertible filter banks, * _{n}* (θ,ϕ) is constrained to be the same as

$${\tilde{h}}_{n}(\theta ,\varphi )=\left({\displaystyle \prod _{k=1}^{n}{b}_{k}}\right){D}_{{a}_{n}}\tilde{h}(\theta ,\varphi )$$

(15)

where *b _{k}* ≥ 1 and

In this paper, we adopt the stereographic dilation operator introduced in [1] (Fig. 3), which involves stereographically projecting the function from the sphere onto the tangent plane of the north pole, performing the usual dilation operation on the plane and then projecting the resulting function back onto the sphere^{2}

$$\begin{array}{l}[{D}_{a}f](\theta ,\varphi )=\left(\frac{1+{\text{tan}}^{2}\frac{\theta}{2}}{1+{\left(\frac{1}{a}+\text{tan}\frac{\theta}{2}\right)}^{2}}\right)\\ \hfill \text{\hspace{1em}\hspace{1em}\hspace{1em}}\times \frac{1}{a}f\left(2{\text{tan}}^{-1}\left(\frac{1}{a}\text{tan}\frac{\theta}{2}\right),\varphi \right).\end{array}$$

(16)

The normalization factor ( 1 + tan^{2}(θ/2)) / (1 + ((1 /*a*)tan(θ/2))^{2} ensures the inner product between functions is conserved under stereographic dilation. Stereographic dilation allows for an explicit control of the spatial localization of the wavelets in contrast with previous approaches that define dilation in the frequency domain [4]. Because of the nonlinear nature of stereographic dilation, extreme dilation of a spherical function will eventually lead to high frequencies. Initial dilation of a function localized at the north pole increases the support of the function and results in lower frequencies. As dilation continues, oscillations will accumulate near the south pole, resulting in high frequencies. A useful analysis of this phenomenon can be found in Bogdanova *et al.* [4]. In practice, we will avoid working in that region, since the dilated filter no longer looks like the original filter.

Stereographic dilation. Correspondence between the spherical surface and tangent plane is established via stereographic projection: point *P* on the sphere is mapped to the point *P′* on the tangent plane of the north pole. This allows a spherical **...**

The *b _{k}*’s in (15) are the amplitude scaling parameters that control the tradeoff between self-invertibility and norm-preserving dilation. Corollary 4.2 implies that the sum of squares of the spherical harmonic coefficients of a bank of self-invertible filters must increase linearly with degree. But stretching a function while preserving its norm shifts its spherical harmonic coefficients to the left (spherical harmonic degrees decrease) and magnifies them (Fig. 4).

These extra weights are analogous to the measure of scale 1/*a*^{3}*da* in the group theoretic formulation of spherical wavelets [1] and its discretization in the axisymmetric case [4], resulting in wider filters being assigned smaller weights. On the continuous real line, the measure 1/*a*^{2}*da* nicely cancels out the dilation of the filter (cf. [27, Ch. 5]). On the discrete real line, the convolution outputs of narrower filters are sampled more densely. This suggests two possible approaches: variable sampling of the convolution outputs or variable scaling of the filters. Yet another approach is to sample the scale space unevenly rather than according to a power law. Because the effects of stereographic dilation on the spherical harmonic coefficients of a function are not analytical, none of these approaches leads to a closed-form solution. In this paper, we take the variable scaling approach by finding the appropriate *b _{k}*’s as part of the filter design.

Fortunately, stereographic dilation is distributive over addition. Suppose the template is expressible as a linear combination of the basis functions *B _{i}*(θ,ϕ) i.e., $\tilde{h}(\theta ,\varphi )={\displaystyle {\sum}_{i=1}^{M}{c}_{i}{B}^{i}(\theta ,\varphi )}$. Here, we assume that

$${[{D}_{a}\tilde{h}]}^{l,m}={\left[{D}_{a}{\displaystyle \sum _{i=1}^{M}{c}_{i}{B}^{i}}\right]}^{l,m}={\displaystyle \sum _{i=1}^{M}{c}_{i}{[{D}_{a}{B}^{i}]}^{l,m}}$$

(17)

yields the spherical harmonic coefficients of the analysis filter at another scale. This is useful since the invertibility condition in Corollary 4.2 is expressed in terms of the spherical harmonic coefficients of the filters. We can, therefore, decide on a set of scales ${\left\{{a}_{n}\right\}}_{n=1}^{N}$ and precompute a table of spherical harmonic coefficients of the dilated basis functions [*D _{a}B^{i}*]

After fixing the set of basis functions {*B ^{i}*} and the set of scales {

The quadratic penalty method [3] is effective in solving this optimization problem with nonconvex constraints by incorporating the constraints into the objective function and solving the resulting unconstrained optimization problem using nonlinear least squares minimization

$$E(\overrightarrow{c},\overrightarrow{b})=\lambda {\displaystyle \sum _{l}{\Vert {H}_{\overrightarrow{c},\overrightarrow{b}}(l)-1\Vert}^{2}}+{\displaystyle \sum _{n=1}^{N}{\Vert {\overrightarrow{\mathrm{\Gamma}}}_{n}\left({\tilde{h}}_{{a}_{n}}(\overrightarrow{c},\overrightarrow{b})\right)\Vert}^{2}.}$$

(18)

The first term corresponds to the invertibility conditions by constraining the frequency response *H _{$\stackrel{\u20d7}{c}$,$\stackrel{\u20d7}{b}$}*(

In this section, we demonstrate the optimization procedure formulated in the previous section. We demonstrate the construction of both self-invertible spherical wavelets and spherical steerable pyramids. Similar to the Euclidean domain, we define a spherical wavelet transform to be the decomposition of a spherical signal into component signals at different scales, i.e., employing axisymmetric filter kernels. On the other hand, we reserve the term spherical steerable pyramid transform for the decomposition of a spherical signal into component signals at different scales and orientations, i.e., using nonaxisymmetric filter kernels. We note that in some literature [1], [4], [29], the term “spherical wavelets” includes spherical steerable pyramids.

In designing axisymmetric wavelets, we limit our set of basis functions {*B ^{i}*} to be the first hundred spherical harmonics of order 0, since the spherical harmonic coefficients of axisymmetric functions are zero for orders other than 0.

We define the set of scales to be *a* = {2^{−n/3}}, *n* = −6, −5,…,2, 3,with *a* = 1 (*n* = 0) corresponding to the undilated template. We use S2kit [15]^{3} to create a table of the spherical harmonic coefficients of ${D}_{a}{Y}_{l}^{0}$ for *l* = 0,…99. As mentioned earlier, extreme stereographic dilation and shrinking of spherical harmonics can result in high frequencies. We find the first 600 order 0 spherical harmonic coefficients of each dilated spherical harmonic (a dilated axisymmetric function remains axisymmetric). We verify that for *a* = 4(*n* = −6) and, $a=0.5\phantom{\rule{thinmathspace}{0ex}}\left(n=3\right),\phantom{\rule{thinmathspace}{0ex}}\left|{[{D}_{a}{Y}_{99}^{0}]}^{599,0}\right|/{\text{max}}_{l}\left|{[{D}_{a}{Y}_{99}^{0}]}^{l,0}\right|\phantom{\rule{thinmathspace}{0ex}}<{10}^{-12}$.

For axisymmetric filters, we can use the fast spherical convolution [11] to compute the forward convolution (3). We quote the results here for completeness

$${[y\left(\beta ,\alpha \right)]}^{l,m}=\sqrt{\frac{4\pi}{\left(2l+1\right)}}{x}^{l,m}{\tilde{h}}^{l,0*}.$$

(19)

Furthermore, as noted before, the inverse convolution (4) between two spherical images is the same as the definition of convolution by Driscoll and Healy [11] (except for a conjugation)

$${[\widehat{x}\left(\theta ,\varphi \right)]}^{l,m}=2\pi \sqrt{\frac{4\pi}{\left(2l+1\right)}}{y}^{l,m}{h}^{l,0}.$$

(20)

Because we seek a multiscale decomposition of the original spherical signal, we would like the filter at each scale to act as a bandpass filter. Similar to the Euclidean domain, we require a residual lowpass filter to ensure that the combined wavelet and the lowpass filter bank is invertible up to a particular degree. Since we only use the first 100 spherical harmonics as our basis, the frequency response of _{a=1}(θ,ϕ) will be zero for all degrees higher than 99. If we also penalize the magnitude of the leading spherical harmonic coefficients of _{a=1}(θ,ϕ), the frequency response of _{a=1}(θ,ϕ) will be zeros at both ends, i.e., it will serve as a bandpass filter. To satisfy the self-invertibility conditions, the solution cannot be identically zero, but must rise to a peak somewhere in the middle of the frequency range.

We also penalize the second derivatives of the filters’ frequency responses and spherical harmonic coefficients to force the filters to be relatively smooth and to reduce ringing. Note that the second derivatives are discrete since spherical harmonic degrees are discrete. We can induce a sharper cutoff frequency by penalizing the magnitude of the combined frequency response above a cutoff degree *L _{c}*. In addition, we fix the amplitude scaling factors

Fig. 5(a) illustrates the frequency response of a ten-scale wavelet filter bank ( *a* = {2^{−n/3}}, *n* = −6, −5,…,2, 3) obtained through our optimization procedure. Invertibility is enforced from degree 15 to 79. Furthermore, we impose a quadratic penalty on the magnitude of the combined frequency response for degrees above *L _{c}* = 150. The combined frequency response of the filters is shown in Fig. 5(c).

Ten-scale wavelet filter bank obtained by imposing invertibility from degree 15 to 79 and a combined frequency response cutoff at *L*_{c} = 150. (a) Frequency response of individual filters. (b) Individual filters in the spatial domain (0 ≤ θ **...**

Because the filters are axisymmetric, we can plot the filters in the image domain as a function of θ [Fig. 5(b)]. The existence of a second peak after the peak at θ = 0 (north pole) indicates ringing. When we vary the cutoff frequency penalty, we can trade off the amount of ringing for the slope of the cutoff. For example, Fig. 6(a) shows the combined frequency response of a wavelet filter bank obtained by penalizing the magnitude of the combined frequency response for degrees above *L _{c}* = 100. Notice the combined frequency response drops rapidly after degree 79. However, this results in increased ringing. If we measure ringing by the ratio of the second maxima to the maxima at the north pole, we can measure the tradeoff between ringing and the cutoff frequency, as shown in Fig. 6(a).

(a) Combined frequency response of the filters obtained when the combined frequency response cutoff is set to *L*_{c} = 100. Note the sharper cutoff obtained. However, this is at the expense of ringing. (b) Plot of ringing versus cutoff frequency *L*_{c}. Ringing **...**

As a verification, we apply the resulting filters to a high resolution world elevation map (size > 1000 × 1000). We bilinearly interpolate the map onto the S2kit [15] grid and convert the spatial image into spherical harmonics using S2kit. This is followed by the inverse transform using S2kit to obtain a final spherical image. The spherical harmonics and the final spherical image become the “ground truth” for measuring invertibility. Lower resolution versions of the elevation map are obtained by truncating the higher degrees spherical harmonics. It is necessary to use the spherical image obtained from inverting the spherical harmonics rather than from the linearly interpolated image because the linearly interpolated image might not be samples from a bandlimited signal. In that case, the sampling theorem of Driscoll and Healy [11] does not hold since the quadrature knots and weights were computed assuming samples from a bandlimited signal. This is in contrast to the Euclidean case, where the forward and inverse discrete Fourier transform always return the same results. As a result, the fast spherical harmonic transform of S2kit is not completely invertible for nonbandlimited signals and is instead an approximate fit of the linearly interpolated data. The final spherical image and the spherical harmonics are, however, exact inverses of each other.

We convolve the wavelet filter bank of Fig. 5(a) with the world elevation map [Fig. 7(a)]. The results for four scales are shown in Fig. 7(b)–(e). Upon reconstruction using (20), we find that max_{l,m} | (*$\widehat{x}$ ^{l,m} − x^{l,m}*)/

To demonstrate that the optimization procedure is stable across different settings of parameters, we show a second example where we optimize for a four-scale wavelet filter bank (*a* = {4, 2, 1, 0, 5}. We enforce invertibility from degree 10 to 89 and apply a quadratic penalty on the magnitude of the combined frequency response for degrees above *L _{c}* = 150. The combined frequency response of the resultant filter bank is shown in Fig. 8(a). Once again, by varying the cutoff frequency threshold, we can obtain a tradeoff between ringing and sharpness of the cutoff [see Fig. 8(b)]. We apply the filter bank to the world elevation map [see Fig. 9] and find that invertibility is indeed obtained for degrees between 10 and 89 inclusive. Notice that there is significantly less ringing artifacts than in Fig. 7 as predicted by our measure of ringing at

Just like the Euclidean domain [14], it can be shown that there is a direct tradeoff between angular resolution and steerability of oriented (nonaxisymmetric) filters on the sphere [29], i.e., filters that have higher angular resolving power requires a bigger set of “steering” basis filters. This can also be seen from Theorem 4.3. Filters with higher angular resolving power require spherical harmonics of higher orders, translating to more samples needed to satisfy the conditions of the generalized sampling theorem.

In our experiments, we limit our set of basis functions {*B _{i}*} to be the first two hundred spherical harmonics of order +1and −1. We note that we can increase our angular power by using higher orders, but this decreases the steerability of our filters. By considering only real filters, we can avoid working directly with the order −1 spherical harmonics, since their coefficients are effectively constrained by those of the order +1 spherical harmonics (see Appendix A). For convenience, we further constrain the coefficients of the order +1 spherical harmonics to be real.

We define the set of scales to be *a* = {2^{−n/2}}, *n* = −4,…,4, with *a* = 1 corresponding to the undilated template. Once again, we use S2kit [15] to create a table of the spherical harmonic coefficients of ${D}_{a}{Y}_{l}^{1}$ for *l* = 1,…,200. We find the first 999 order 1 spherical harmonic coefficient of each dilated spherical harmonic (the order of a spherical function does not change under dilation). We verify that for *a* = 4 and *a* = 0.25, ${|[{D}_{a}{Y}_{200}^{1}]}^{999,1}|/{\text{max}}_{l}|\phantom{\rule{thinmathspace}{0ex}}{[{D}_{a}{Y}_{200}^{1}]}^{l,1}|\phantom{\rule{thinmathspace}{0ex}}<{10}^{-12}$.

Similar to the previous subsection, we penalize the magnitude of the leading coefficients of *h*_{a=1}(θ,ϕ). We also penalize the second derivatives of the filters’ frequency responses and spherical harmonic coefficients. Finally, we fix the amplitude scaling factors *b _{k}*’s at all scales to be the same.

Fig. 10(a) and (b) illustrates the frequency response of a ninescale steerable pyramid ( *a* = {2−*n*/2}, *n* = −4,…4) obtained through our optimization procedure. Invertibility is enforced from degree 20 to 170. Note that the frequency response is equal to 0.5 in the invertibility range because we only plot the frequency response contributed by the order +1 harmonics. Fig. 10(c) shows *h*_{a=4}(θ,ϕ) as a spherical image. Note that it looks like a derivative of Gaussian. We can also quantify ringing by plotting *h*_{a=4}(θ,ϕ) as a function of θ while fixing ϕ to correspond to the great circle passing through the maxima and minima of the filter [Fig. 10(d)].

(a)–(d) Nine-scale steerable pyramid (*a* = {2^{−n/2}}, *n* = −4; …, 4) obtained by imposing invertibility from degree 20 to 170. Note that the frequency response is equal to 0.5 in the invertibility range because we only plot **...**

Similarly, Fig. 10(e)–(f) shows the frequency response of a five-scale steerable pyramid (*a* = {4, 2, 1, 0.5, 0.25}) obtained through our optimization procedure. Invertibility is enforced from degree 10 to 180. Fig. 10(g) shows *h*_{a=4}(θ,ϕ) as a spherical image and Fig. 10(h) is a plot of *h*_{a=4}(θ,ϕ) as a function of θ by fixing ϕ.

We now employ the steerable pyramid from Fig. 10(e)–(f) in a denoising experiment. We emphasize that this toy example only serves to illustrate the use of the self-invertible steerable pyramid and is not meant to be the model on how wavelet or steerable pyramid denoising should be done. We first obtain a full pass steerable pyramid filter bank by computing a residual lowpass and highpass filter bank so that the combined filter bank is invertible up to degree 300. By constraining the analysis and synthesis filters to be identical and axisymmetric, the lowpass and highpass filters are uniquely determined.

We use the world elevation map truncated to a maximum degree of 300. A Gaussian noise map is first produced in the spatial domain on the S2kit grid. The forward and inverse spherical harmonic transform is then performed with S2kit [15]. As noted earlier, the resulting noise image map from the inverse spherical harmonic transform will be different from the original noise map, because the fast spherical harmonic transform [11] is not invertible for nonbandlimited signals. Interestingly, noise near the poles tends to be suppressed under the forward and inverse spherical harmonic transform (not shown), possibly related to the concentration of samples near the poles on a latitude- longitude grid. The noise map obtained from the inverse spherical harmonic transform is then added to the world map. Fig. 11(a)–(d) display world elevation maps with different peak signal-to-noise ratios (PSNR).

(a)–(d) Corrupted world elevation maps. (e)–(h) Corresponding denoised world elevation maps, i.e., (e) shows the denoised version of (a), and so on. (a) PSNR = 27.923 dB. (b) PSNR = 21.943 dB. (c) PSNR = 15.9035 dB. (d) PSNR = 9.8804 dB. **...**

We adapt the wavelet shrinkage technique pioneered by Donoho and Johnstone [10] to our denoising problem.

- We perform continuous spherical convolution between the noisy world elevation map and the analysis filters of the steerable pyramid using the algorithm in [28].
- We sample the transform coefficients using the inverse FFT, by using the uniform sampling grid of the second quadrature rule (Appendix D) and the fact that the continuous spherical convolution in step 1 is computed by representing functions
*SO*(3) as sum of complex exponentials. - We threshold the sampled transform coefficientswhere for brevity, we denote the sampled transform coefficients$${y}_{njsk}\leftarrow \text{sign}({y}_{njsk})\phantom{\rule{thinmathspace}{0ex}}\text{max}(|{y}_{njsk}|-{t\sigma}_{n},0)$$(21)
*y*(α_{n}_{j,n},β_{s,n},γ_{k,n}) by*y*. The soft threshold (shrinkage) sets to zero for_{njsk}*y*below_{njsk}*t*σ_{n}in absolute value and pulls other data towards the origin by an amount*t*σ_{n}. The noise standard deviation σ_{n}can be determined empirically or by simulations [26]. In our toy example, we find that σ_{n}consistently doubles from one steerable pyramid scale to the finer one, except for the scale corresponding to the highpass filter.*t*is a user-defined parameter. In our experiments, we vary from*t*0.5 to 5. We note that in the SureShrink [10] algorithm, the parameter*t*is found automatically once the noise standard deviation σ_{n}is known. Unfortunately, naive application of SureShrink gives poor results. One reason is that the coefficients of the overcomplete steerable pyramid is correlated, and, hence, the noise is no longer independently and identically distributed (i.i.d) in the steerable pyramid domain. - We utilize (6) and (7) to reconstruct the denoised image. Details of the computation are found in Appendix E.

Recall that *L _{n}* (

Fig. 12 shows a plot of the PSNR of the denoised images as we vary *t* from 0.5 to 5. In general, the noisier the image before denoising, the greater the improvement. In particular, there is an improvement of more than 10dB for the noisiest map. The best results are obtained for *t* between 1 and 2. The denoised maps corresponding to the best *t*’s are shown in Fig. 11(e)–(f). Artifacts become more noticeable for the noisier images. These artifacts are likely the by-products of ringing.

In this paper, we present theoretical conditions for the invertibility of filter banks under continuous convolution on the 2-Sphere. We discretize the results using quadrature, thus obtaining a generalized sampling theorem. We propose a general procedure for constructing invertible filter banks and demonstrate the procedure by generating self-invertible spherical wavelets and steerable pyramids. Finally, we provide an analysis of the computational complexity of the filtering framework.

Nonlinear dilation of functions on the sphere remains difficult to work with. While we circumvent the problem by using the distributive property of stereographic dilation, the spherical harmonic coefficients table can take up a substantial amount of space. More efficient methods are, therefore, needed. Other definitions of dilation (such as those defined in the frequency domain) might also fit better into the computational framework. More work is needed to understand the space of invertible and self-invertible filter banks. As seen in our experiments, there is an implicit tradeoff between the sharpness of the frequency response of the filters and ringing. It will be useful to formulate an objective function that directly trades off between ringing and the sharpness of the frequency response. The bandlimited spatially-concentrated eigenfunctions proposed by Simons *et al.* [24] could prove to be a useful set of basis functions for the optimization problem in our framework.

This paper introduces theoretical results on invertibility and sampling, and represents a step towards a general framework for filter design on the 2-Sphere. Just as wavelets and steerable pyramids have been useful for the processing and analysis of planar images, we are optimistic that future work will lead to similar applications on the sphere.

The authors would like to thank M. Tappen for discussions on optimization procedures, B. Freeman and T. Adelson for discussions on self-invertibility, V. Goyal for discussions on norm preserving dilations in Euclidean wavelets and wavelet 7 denoising, and F. Durand and B. Horn for reading earlier drafts of this paper. The authors would also like to thank the reviewers for their many useful suggestions and pointers to references. T. Yeo would like to thank C.-L. Cheng for discussions on spherical harmonics and Wigner-D functions.

This work was supported in part by the National Alliance for Medical Image Analysis (NIH NIBIB NAMIC U54-EB005149), in part by the Neuroimaging Analysis Center (NIT CRR NAC P41-RR13218), in part by the Morphometry Biomedical Informatics Research Network (NIH NCRR mBIRN U24-RR021382), in part by the NIH Grant NINDS R01-NS051826, and in part by the National Science Foundation (NSF) under Grant CAREER 0642971. B.T. T. Yeo was supported in part by the Agency for Science, Technology, and Research, Singapore. W. Ou was supported by the NSF graduate fellowship. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Ljubisa Stankovic.

**Boon Thye Thomas Yeo** received the B.S. and M.S. degrees in electrical engineering from Stanford University, Stanford, CA, in 2002. He is currently pursuing the Ph.D. degree in computer science at the Massachusetts Institute of Technology, Cambridge.

His current research focuses on developing and applying techniques in signal processing, statistics, and differential geometry to the problems of registration, segmentation, and shape analysis in neuroimaging.

Mr. Yeo has received several awards, including the MICCAI Young Scientist Award, A * STAR National Science Scholarship, Frederick E. Terman Engineering Scholastic Award, and the Hewlett-Packard and Agilent Technologies Project Award.

**Wanmei Ou** received the B.S. degree in electrical engineering from the City College of the City University of New York (Salutatorian) in 2003, the M.S. degree in electrical engineering and computer science from the Massachusetts Institute of Technology, Cambridge, in 2005, where she is currently pursuing the Ph.D. degree.

She is interested in developing mathematical algorithms to improve brain activity detection for functional brain studies, as well as to quantify difference in anatomical structures. Her research projects include spatial regularization for fMRI, spatio-temporal EEG/MEG inverse problems, neurovascular coupling, wavelets analysis, and joint analysis of fMRI and EEG/MEG data.

**Polina Golland** received the B.S. and M.S. degrees from the Technion—Israel Institute of Technology, Haifa, and the Ph.D. degree in 2001 from the Massachusetts Insitute of Technology (MIT), Cambrige.

She is an Assistant Professor in the Electrical Engineering and Computer Science Department and the Computer Science and Artificial Intelligence Laboratory, MIT. She has worked on various problems in computer vision, including motion and stereo, shape modeling, and representation. Her current research focuses on modeling biological shape and function using images (from MRI to microscopy) as a source of information.

Prof. Golland has served on program committees of several leading conferences and as a referee for international journals in biomedical image analysis and computer vision. She received the NSF CAREER Award in 2007.

Here, we review useful facts on the spherical harmonics.

The spherical harmonics ${Y}_{l}^{m}$ are defined in terms of the associated Legendre polynomials ${P}_{l}^{m}$. For a given degree *l* ≥ 0 and order |*m*| ≤ *l*, 0 ≤ θ ≤ π, and 0 ≤ ϕ ≤ 2π

$${Y}_{l}^{m}(\theta ,\varphi )=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}{P}_{l}^{m}(\text{cos}\theta ){e}^{im\varphi}$$

(22)

where, for *m* ≥ 0 and |*x*| < 1

$${P}_{l}^{m}(x)=\frac{{(-1)}^{m}}{{2}^{l}l!}{(1-{x}^{2})}^{m/2}\frac{{d}^{l+m}}{d{x}^{l+m}}{({x}^{2}-1)}^{l}$$

(23)

$${P}_{l}^{-m}(x)={(-1)}^{m}\frac{(l-m)!}{(l+m)!}{P}_{l}^{m}(x).$$

(24)

Therefore, for *l* ≥ *m* ≥ 0, we have

$$\begin{array}{cc}{Y}_{l}^{m}(\theta ,\varphi )=& \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\frac{{(-1)}^{m}}{{2}^{l}l!}{(1-{\text{cos}}^{2}\theta )}^{m/2}\\ & \times \phantom{\rule{thinmathspace}{0ex}}\frac{{d}^{l+m}}{d{x}^{l+m}}{{({x}^{2}-1)}^{l}|}_{x=\text{cos}\theta}{e}^{im\varphi}\hfill \end{array}$$

(25)

$$\begin{array}{cc}{Y}_{l}^{-m}(\theta ,\varphi )=& \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\frac{1}{{2}^{l}l!}{(1-{\text{cos}}^{2}\theta )}^{m/2}\\ & \times \phantom{\rule{thinmathspace}{0ex}}\frac{{d}^{l+m}}{d{x}^{l+m}}{{({x}^{2}-1)}^{l}|}_{x=\text{cos}\theta}{e}^{-im\varphi}\hfill \end{array}$$

(26)

$$={(-1)}^{m}{Y}_{l}^{m*}(\theta ,\varphi ).$$

(27)

Under rotation, each spherical harmonic of degree *l* is transformed into a linear combination of spherical harmonics of the same degree but possibly different orders. In particular, if we parametrize our rotation by the three Euler angles, α, β, γ, and rotate our original function *f*, it can be shown that

$${[D(\alpha ,\beta ,\gamma )f]}^{l,m}={\displaystyle \sum _{m\prime =-l}^{l}{D}_{mm\prime}^{l}}(\alpha ,\beta ,\gamma ){f}^{l,m\prime}$$

(28)

where ${D}_{mm\prime}^{l}(\alpha ,\beta ,\gamma )$ is the Wigner-D function [20]. We can further decompose ${D}_{mm\prime}^{l}(\alpha ,\beta ,\gamma )$ as follows:

$${D}_{mm\prime}^{l}(\alpha ,\beta ,\gamma )={e}^{-im\alpha}{d}_{mm\prime}^{l}(\beta ){e}^{-im\prime \gamma}$$

(29)

where ${D}_{mm\prime}^{l}(\beta )$ is the Wigner-d function and is real [20]

$$\begin{array}{cc}{d}_{mm\prime}^{l}(\beta )=\hfill & {\displaystyle \sum _{j}{(-1)}^{j-m\prime +m}}\hfill \\ & \times \frac{\sqrt{(l+m\prime )!(l-m\prime )!(l+m)!(l-m)!}}{(l+m\prime -j)!j!(l-j-m)!(j-m\prime +m)!}\hfill \\ & \times {\left(\text{cos}\frac{\beta}{2}\right)}^{2l-2j+m\prime -m}{\left(\text{sin}\frac{\beta}{2}\right)}^{2j-m\prime +m}.\hfill \end{array}$$

(30)

The sum is over all *j* such that none of the denominator terms with factorials is negative. This reflects the fact that only rotations about the *y*-axis mixes orders.

By the Peter–Weyl theorem on compact groups [30]

$$\begin{array}{l}{\displaystyle {\int}_{SO(3)}{D}_{mn}^{l}(\rho ){D}_{m\prime n\prime}^{l\prime *}(\rho )d\rho}\\ \hfill \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}=\frac{8{\pi}^{2}}{2l+1}\delta (l-l\prime ,m-m\prime ,n-n\prime ).\end{array}$$

(31)

By integrating out α and γ, we obtain

$${\int}_{\beta}{d}_{mn}^{l}(\beta ){d}_{mn}^{l\prime}(\beta )\text{sin}\beta d\beta}=\frac{2}{2l+1}\delta (l-l\prime ).$$

(32)

We will also use the identity [28]

$$D\phantom{\rule{thinmathspace}{0ex}}(\alpha ,\beta ,\gamma )=D\phantom{\rule{thinmathspace}{0ex}}\left(\alpha +\frac{\pi}{2},\frac{\pi}{2},0\right)\phantom{\rule{thinmathspace}{0ex}}D\phantom{\rule{thinmathspace}{0ex}}\left(\beta +\pi ,\frac{\pi}{2},\frac{\pi}{2}+\gamma \right).$$

(33)

Here, we prove Theorem 4.1 on continuous frequency response. We first note that by using Parseval’s Theorem and substituting (28), we can re-write the output of the *n*th analysis filter as

$${y}_{n}(\alpha ,\beta ,\gamma )={\displaystyle {\int}_{{S}^{2}}[D(\alpha ,\beta ,\gamma ){\tilde{h}}_{n}{]}^{*}(\theta ,\varphi )x(\theta ,\varphi )d\mathrm{\Omega}}$$

(34)

$$\begin{array}{cc}=& {{\displaystyle \sum _{l\prime =0}^{\infty}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}(\alpha ,\beta ,\gamma ){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\\ & \times \phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}\hfill \end{array}$$

(35)

and the reconstructed image as

$$\widehat{x}(\theta ,\varphi )={\displaystyle \sum _{n=1}^{N}{\displaystyle {\int}_{SO\left(3\right)}[D(\alpha ,\beta ,\gamma ){h}_{n}](\theta ,\varphi ){y}_{n}(\alpha ,\beta ,\gamma )d\rho}}$$

(36)

$$\begin{array}{ll}(\underset{\xaf}{\underset{\xaf}{35}})\hfill & {\displaystyle \sum _{n=1}^{N}{\displaystyle {\int}_{SO(3)}[D(\alpha ,\beta ,\gamma ){h}_{n}](\theta ,\varphi )}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{{\displaystyle \sum _{l\prime =0}^{\infty}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}(\alpha ,\beta ,\gamma ){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}d\mathrm{\rho}.\hfill \end{array}$$

(37)

Projecting *$\widehat{x}$*(θ,ϕ) onto the spherical harmonics basis, we obtain

$${\widehat{x}}^{l,m}={\displaystyle {\int}_{{S}^{2}}\widehat{x}(\theta ,\varphi ){Y}_{l}^{m*}}(\theta ,\varphi )d\mathrm{\Omega}$$

(38)

$$\begin{array}{ll}(\underset{\xaf}{\underset{\xaf}{37}})\hfill & {\displaystyle \sum _{n=1}^{N}{\displaystyle {\int}_{SO(3)}\left[{\displaystyle {\int}_{{S}^{2}}[D(\alpha ,\beta ,\gamma ){h}_{n}](\theta ,\varphi )}{Y}_{l}^{m*}(\theta ,\varphi )d\mathrm{\Omega}\right]}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{{\displaystyle \sum _{l\prime =0}^{\infty}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}(\alpha ,\beta ,\gamma ){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}d\rho .\hfill \end{array}$$

(39)

$$\begin{array}{ll}(\underset{\xaf}{\underset{\xaf}{28}})\hfill & {\displaystyle \sum _{n=1}^{N}{\displaystyle {\int}_{SO\left(3\right)}\left[{\displaystyle \sum _{m\u2034=-l}^{l}{h}_{n}^{l,m\u2034}{D}_{m,m\u2034}^{l}(\alpha ,\beta ,\gamma )}\right]}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{{\displaystyle \sum _{l\prime =0}^{\infty}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}(\alpha ,\beta ,\gamma ){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}d\rho .\hfill \end{array}$$

(40)

$$\begin{array}{ll}=\hfill & {\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\u2034=-l}^{l}\left[{h}_{n}^{l\prime ,m\u2034}\right]{\displaystyle \sum _{l\prime =0}^{\infty}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}{x}^{l\prime ,m\prime}}}}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}[{\tilde{h}}_{n}^{l\prime ,m\u2033}{]}^{*}{\displaystyle {\int}_{SO(3)}{D}_{m,m\u2034}^{l}(\alpha ,\beta ,\gamma )}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{D}_{m\prime ,m\u2033}^{l\prime}{(\alpha ,\beta ,\gamma )}^{*}d\rho \hfill \end{array}$$

(41)

$$(\underset{\xaf}{\underset{\xaf}{31}})\phantom{\rule{thinmathspace}{0ex}}{x}^{l,m}\frac{8{\pi}^{2}}{2l+1}\left\{{{\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\u2034=-l}^{l}\left[{h}_{n}^{l,m\u2034}\right]\phantom{\rule{thinmathspace}{0ex}}\left[{\tilde{h}}_{n}^{l,m\u2034}\right]}}}^{*}\right\}.$$

(42)

Here, we prove the generalized sampling theorem. Recall from Section III that *y _{n}*(α

Since the output of the *n*th synthesis filter * _{n}* (θ, ϕ) is a linear combination of rotated versions of

$${\widehat{x}}_{n}^{l,m}=\frac{8{\pi}^{2}}{2l+1}{x}_{n}^{l,m}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\u2033=-l}^{l}{h}_{n}^{l,m\u2033}{\tilde{h}}_{n}^{l,m\u2033*}=0}.$$

(43)

We will now show that for *l* ≤ *L _{n}*

$${\widehat{x}}_{n}^{l,m}=\frac{8{\pi}^{2}}{2l+1}{x}_{n}^{l,m}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\u2033=-l}^{l}{h}_{n}^{l,m\u2033}{\tilde{h}}_{n}^{l,m\u2033*}.}$$

(44)

From (35)

$$\begin{array}{l}{y}_{n}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n})\\ \text{\hspace{1em}\hspace{1em}}={{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}.\end{array}$$

(45)

$${\widehat{x}}_{n}(\theta ,\varphi )\phantom{\rule{thinmathspace}{0ex}}(\underset{\xaf}{\underset{\xaf}{6}})\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{w}_{j,s,k,n}}}}[D({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}](\theta ,\varphi ){y}_{n}\left({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}\right)$$

(46)

$$\begin{array}{c}(\underset{\xaf}{\underset{\xaf}{45}}){\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}\frac{4{\pi}^{2}{w}_{s,n}}{{J}_{n}{K}_{n}}}}}[D({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}]\left(\theta ,\varphi \right)\\ \text{\hspace{1em}\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}{{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}\end{array}$$

(47)

$$\begin{array}{ll}{\widehat{x}}^{l,m}& ={\displaystyle {\int}_{{S}^{2}}{\widehat{x}}_{n}(\theta ,\varphi ){Y}_{l}^{m*}(\theta ,\varphi )d\mathrm{\Omega}}\\ & (\underset{\xaf}{\underset{\xaf}{28}})\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}\frac{4{\pi}^{2}{w}_{s,n}}{{J}_{n}{K}_{n}}}}}\phantom{\rule{thinmathspace}{0ex}}\left[{\displaystyle \sum _{m\u2034=-l}^{l}{D}_{m,m\u2034}^{l}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}^{l,m\u2034}}\right]\end{array}$$

(48)

$$\begin{array}{c}\text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}{{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}\left[{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{D}_{m\prime ,m\u2033}^{l\prime}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){\tilde{h}}_{n}^{l\prime ,m\u2033}}\right]}}}^{*}\phantom{\rule{thinmathspace}{0ex}}{x}^{l\prime ,m\prime}\\ =\phantom{\rule{thinmathspace}{0ex}}\frac{4{\pi}^{2}}{{J}_{n}{K}_{n}}{\displaystyle \sum _{m\u2034=-l}^{l}{h}_{n}^{l,m\u2034}}{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}{x}^{l\prime ,m\prime}}{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{\left[{\tilde{h}}_{n}^{l\prime ,m\u2033}\right]}^{*}}}\hfill \\ \text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{w}_{s,n}}}}{D}_{m,m\u2034}^{l}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){D}_{m\prime ,m\u2033}^{l\prime *}\left({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}\right)\hfill \end{array}$$

(49)

For *w j,s,k,n* = 4π

$$\begin{array}{l}{\widehat{x}}_{n}^{l,m}=\frac{4{\pi}^{2}}{{J}_{n}{K}_{n}}{\displaystyle \sum _{m\u2034=-l}^{l}{h}_{n}^{l,m\u2034}}{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}{x}^{l\prime ,m\prime}}}\\ \hfill \times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{\left[{\tilde{h}}_{n}^{l\prime ,m\u2033}\right]}^{*}}\mathrm{\Phi}.\end{array}$$

(50)

To simplify Φ, we write

$$\begin{array}{cc}\mathrm{\Phi}=& {\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{w}_{s,n}}}}{D}_{m,m\u2034}^{l}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n})\\ & \times {D}_{m\prime ,m\u2033}^{l\prime *}({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n})\hfill \end{array}$$

(51)

$$\begin{array}{ll}(\underset{\xaf}{\underset{\xaf}{29}})\hfill & {\displaystyle \sum _{j=0}^{{J}_{n}-1}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{w}_{s,n}{e}^{-im{\alpha}_{j,n}}{d}_{m,m\u2034}^{l}({\beta}_{s,n})}}}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}{e}^{-im\u2034}{\gamma}_{k,n}{e}^{-im\prime {\alpha}_{j,n}}{d}_{m\prime ,m\u2033}^{l\prime}({\beta}_{s,n}){e}^{im\u2033{\gamma}_{k,n}}\hfill \\ =\hfill & {\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{m,m\u2034}^{l}({\beta}_{s,n}){d}_{m\prime ,m\u2033}^{l\prime}({\beta}_{s,n})}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}\left[{\displaystyle \sum _{j=0}^{{J}_{n}-1}{e}^{i(m\prime -m){\alpha}_{j,n}}}\right]\phantom{\rule{thinmathspace}{0ex}}\left[{\displaystyle \sum _{k=0}^{{K}_{n}-1}{e}^{i(m\u2033-m\u2034){\gamma}_{k,n}}}\right].\hfill \end{array}$$

(52)

We note that and |*m′*| ≤ *l′* ≤ * _{n}* |

$$\begin{array}{cc}{\displaystyle \sum _{j=0}^{{J}_{n}-1}{e}^{i(m\prime -m){\alpha}_{j,n}}}& ={\displaystyle \sum _{j=0}^{{\tilde{L}}_{n}+{L}_{n}}{e}^{i(m\prime -m)2\pi j/({\tilde{L}}_{n}+{L}_{n}+1)}}\hfill \\ & =\{\begin{array}{ll}{\tilde{L}}_{n}+{L}_{n}+1,\hfill & \text{if}(m\prime -m)=0\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$

(53)

Similarly, from (48), we observe that |*m′*| ≤ *O _{n}* and |

$$\begin{array}{l}{\displaystyle \sum _{k=0}^{{K}_{n}-1}{e}^{i(m\u2033-m\u2034)}{\gamma}_{k,n}}\\ \text{\hspace{1em}\hspace{1em}}=\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{k=0}^{{O}_{{\tilde{h}}_{n}}+{O}_{{h}_{n}}}{e}^{i(m\u2033-m\u2034)}2\pi k/({O}_{{\tilde{h}}_{n}}+{O}_{{h}_{n}}+1)}\\ \text{\hspace{1em}\hspace{1em}}=\{\begin{array}{ll}{O}_{{\tilde{h}}_{n}}+{O}_{{h}_{n}}+1,\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}(m\u2033-m)=0\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$

(54)

Substituting into (52), we get

$$\begin{array}{ll}\mathrm{\Phi}=\hfill & {\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{m,m\u2034}^{l}({\beta}_{s,n}){d}_{m\prime ,m\u2033}^{l\prime}({\beta}_{s,n})}\hfill \\ \hfill & \times \phantom{\rule{thinmathspace}{0ex}}\left[{\displaystyle \sum _{j=0}^{{J}_{n}-1}{e}^{i(m\prime -m){\alpha}_{j,n}}}\right]\phantom{\rule{thinmathspace}{0ex}}\left[{\displaystyle \sum _{k=0}^{{K}_{n}-1}{e}^{i(m\u2033-m\u2034){\gamma}_{k,n}}}\right]\hfill \end{array}$$

(55)

$$\begin{array}{c}=\phantom{\rule{thinmathspace}{0ex}}{J}_{n}{K}_{n}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{m,m\u2034}^{l}({\beta}_{s,n}){d}_{m\prime ,m\u2033}^{l\prime}({\beta}_{s,n})}\\ \text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}\delta (m-m\prime )\delta (m\u2033-m\u2034)\hfill \end{array}$$

(56)

$$\begin{array}{c}=\phantom{\rule{thinmathspace}{0ex}}{J}_{n}{K}_{n}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{m,m\u2033}^{l}({\beta}_{s,n}){d}_{m,m\u2033}^{l\prime}({\beta}_{s,n})}\\ \text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}\delta (m-m\prime )\delta (m\u2033-m\u2034)\hfill \end{array}$$

(57)

$$(\underset{\xaf}{\underset{\xaf}{32}})\frac{2{J}_{n}{K}_{n}}{2l+1}\delta (l-l\prime )\delta (m-m\prime )\delta (m\u2033-m\u2034)$$

(58)

where in the third equality, *m* = *m″*, *m′* = *m* because of the delta functions, and the last equality was obtained using the assumption that *w _{sn}* and β

$$\begin{array}{cc}{\widehat{x}}_{n}^{l,m}=& \frac{4{\pi}^{2}}{{J}_{n}{K}_{n}}{\displaystyle \sum _{m\u2034=-l}^{l}{h}_{n}^{l,m\u2034}}\hfill \\ & \times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}{x}^{l\prime ,m\prime}}{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{\left[{\tilde{h}}_{n}^{l\prime ,m\u2033}\right]}^{*}\mathrm{\Phi}}}\hfill \end{array}$$

(59)

$$\begin{array}{c}=\frac{4{\pi}^{2}}{{J}_{n}{K}_{n}}{\displaystyle \sum _{m\u2034=-l}^{l}{h}_{n}^{l,m\u2034}}{\displaystyle \sum _{l\prime =0}^{{\tilde{L}}_{n}}}\hfill \\ \text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\prime =-l\prime}^{l\prime}{x}^{l\prime ,m\prime}}{\displaystyle \sum _{m\u2033=-l\prime}^{l\prime}{\left[{\tilde{h}}_{n}^{l\prime ,m\u2033}\right]}^{*}}\hfill \\ \text{\hspace{1em}}\times \phantom{\rule{thinmathspace}{0ex}}\frac{2{J}_{n}{K}_{n}}{2l+1}\delta (l-l\prime )\delta (m-m\prime )\delta (m\u2033-m\u2034)\hfill \end{array}$$

(60)

$$=\frac{8{\pi}^{2}}{2l+1}{x}^{l,m}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{m\u2033=-l}^{l}{h}_{n}^{l,m\u2033}{\tilde{h}}_{n}^{l,m\u2033*}}\text{for all}\phantom{\rule{thinmathspace}{0ex}}l\le {L}_{n}$$

(61)

thus proving (44). Together with (43), we have

$$\begin{array}{ll}{\widehat{x}}^{l,m}\hfill & =\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{n=1}^{N}{\widehat{x}}_{n}^{l,m}}\hfill \\ \hfill & =\phantom{\rule{thinmathspace}{0ex}}\frac{8{\pi}^{2}}{2l+1}{x}^{l,m}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{n=1}^{N}{\displaystyle \sum _{m\u2033=-l}^{l}{h}_{n}^{l,m\u2033}{\tilde{h}}_{n}^{l,m\u2033*}}}\text{for all}\phantom{\rule{thinmathspace}{0ex}}l.\hfill \end{array}$$

(62)

In this appendix, we derive two different quadrature rules that satisfy the conditions of Theorem 4.3, namely, we show quadrature weights *w _{s,n}* and knots β

$$\begin{array}{c}{\displaystyle {\int}_{0}^{\pi}{d}_{mm\prime}^{l}(\beta ){d}_{mm\prime}^{l\prime}(\beta )\text{sin}(\beta )}d\beta \hfill \\ \hfill \text{\hspace{1em}\hspace{1em}}=\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{s=0}^{{S}_{n}-1}{w}_{s,n}{d}_{mm\prime}^{l}({\beta}_{s,n}){d}_{mm\prime}^{l\prime}({\beta}_{s,n}).}\end{array}$$

(63)

We denote $\begin{array}{ccc}f(\beta )& =& {d}_{mm\prime}^{l}(\beta ){d}_{mm\prime}^{l\prime}(\beta )\end{array}$ and observe that *f*(β) consists of a linear combination of even powers of cos(β/2) and sin(β/2) (see (30)). Therefore, if we make the substitution *u* = sin(β/2), and noting that *f*(*u*) (where we are overloading *f*) is now a polynomial with maximum degree, *Q* = 2(*l + l′*) ≤ 2 (*L _{n}* +

$$\begin{array}{ll}{\displaystyle {\int}_{0}^{\pi}f(\beta )\text{sin}\beta}d\beta \hfill & =2\phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\int}_{0}^{\pi}f(\beta )\text{sin}\left(\frac{\beta}{2}\right)}\phantom{\rule{thinmathspace}{0ex}}\text{cos}\left(\frac{\beta}{2}\right)d\beta \hfill \\ \hfill & =4{\displaystyle {\int}_{0}^{1}f(u)udu.}\hfill \end{array}$$

(64)

Now, making the substitution, υ = 2*u* − 1, we have

$$\begin{array}{ll}{\displaystyle {\int}_{0}^{\pi}f(\beta )\text{sin}\beta}d\beta \hfill & =2\phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\int}_{-1}^{1}f\left(\frac{\upsilon +1}{2}\right)\frac{\upsilon +1}{2}}\phantom{\rule{thinmathspace}{0ex}}d\upsilon \hfill \\ \hfill & =2{\displaystyle \sum _{k=0}^{N-1}{r}_{k}\frac{{\upsilon}_{k}+1}{2}f\left(\frac{{\upsilon}_{k}+1}{2}\right)}\hfill \end{array}$$

(65)

where *r _{k}*’s are defined to be the weights of the Gauss–Legendre quadrature on the interval [−1,1], and υ

From the substitution above, we have sin β_{k}/2 = *u _{k}* = (υ

We will derive another rule in this section, using the technique shown in [11]. However, first we need to obtain the Fourier series formula for the square wave, *SQ*(*u*), which is defined to be periodic from −π to π

$$SQ(u)=\{\begin{array}{ll}1,\hfill & -\pi <u<\frac{-\pi}{2}\hfill \\ -1,\hfill & \frac{-\pi}{2}<u<0\hfill \\ 1,\hfill & 0<u<\frac{\pi}{2}\hfill \\ -1,\hfill & \frac{\pi}{2}<u<\pi .\hfill \end{array}$$

(66)

Projecting *SQ*(*u*) onto the Fourier series basis, we get

$$\begin{array}{c}\frac{1}{2\pi}{\displaystyle {\int}_{-\pi}^{\pi}SQ(u){e}^{-iku}du}\hfill \\ \text{\hspace{1em}\hspace{1em}}=\frac{1}{2\pi}[{\displaystyle {\int}_{-\pi}^{-\pi /2}{e}^{-iku}du}-{\displaystyle {\int}_{-\pi /2}^{0}{e}^{-iku}du}\hfill \\ \hfill +{\displaystyle {\int}_{0}^{\pi /2}{e}^{-iku}du}-{\displaystyle {\int}_{\pi /2}^{\pi}{e}^{-iku}du}]\end{array}$$

(67)

$$\begin{array}{c}=\frac{i}{2\pi k}[{e}^{-iku}{|}_{-\pi}^{-\pi /2}-{e}^{-iku}{|}_{-\pi /2}^{0}\\ \hfill +{e}^{-iku}{|}_{0}^{\pi /2}-{e}^{-iku}{|}_{\pi /2}^{\pi}]\end{array}$$

(68)

$$\begin{array}{c}=\frac{i}{2\pi k}[({e}^{-ik\pi /2}-{e}^{ik\pi})-(1-{e}^{-ik\pi /2})\\ \hfill +\phantom{\rule{thinmathspace}{0ex}}({e}^{-ik\pi /2}-1)-({e}^{-ik\pi}-{e}^{-ik\pi /2})]\end{array}$$

(69)

$$=\frac{i}{\pi k}\left[{e}^{ik\pi /2}+{e}^{-ik\pi /2}-1-\frac{1}{2}{e}^{ik\pi}-\frac{1}{2}{e}^{-ik\pi}\right]$$

(70)

$$=\frac{i}{\pi k}\left[2\text{cos}\left(k\frac{\pi}{2}\right)-1-2\text{cos}\left(k\pi \right)\right].$$

(71)

- If
*k*= 4*n*, we get 2 cos(2*n*π) − 1 − cos(4*n*π) = 0 - If
*k*= 4*n*+ 1, we get 2 cos(2*n*π + π/2) − 1 − cos (4*n*π + π) = 0 - If
*k*= 4*n*+ 2, we get 2 cos(2*n*π + π) − 1 −cos(4*n*π + 2π) = −4 - If
*k*= 4*n*+ 3, we get 2 cos(2*n*π + 3π/2) − 1 − cos(4*n*π + 3π) = 0.

Therefore, the Fourier series for *SQ*(*u*) is nonzero for *k* = 4*n*+2, and is equal to −4*i*/π*k*, and we have

$$\begin{array}{cc}SQ(u)& ={\displaystyle \sum _{p=-\infty}^{\infty}-\frac{4i}{\pi (4p+2)}{e}^{i(4p+2)u}}\\ & ={\displaystyle \sum _{p=-\infty}^{\infty}-\frac{2i}{\pi (2p+1)}{e}^{i(4p+2)u}.}\end{array}$$

(72)

Now, we can continue with the derivation of the quadrature. We note that *f*(β) = ∑_{k,k′}*a _{k}* cos(β/2)

$$\begin{array}{c}{\displaystyle {\int}_{0}^{\pi}f(\beta )\text{sin}(\beta )}d\beta \hfill \\ \hfill \text{\hspace{1em}\hspace{1em}}=2{\displaystyle {\int}_{0}^{\pi /2}f(\beta \prime )\text{sin}(\beta \prime )}d\beta \prime \hfill \end{array}$$

(73)

$$\begin{array}{c}=\frac{1}{2}{\displaystyle {\int}_{-\pi}^{-\pi /2}f(\beta \prime )\text{sin}(2\beta \prime )}d\beta \prime \hfill \\ \text{\hspace{1em}\hspace{1em}}-\frac{1}{2}{\displaystyle {\int}_{-\pi /2}^{0}f(\beta \prime )\text{sin}(2\beta \prime )}d\beta \prime \\ \text{\hspace{1em}\hspace{1em}}+\frac{1}{2}{\displaystyle {\int}_{0}^{\pi /2}f(\beta \prime )\text{sin}(2\beta \prime )}d\beta \prime \\ \text{\hspace{1em}\hspace{1em}}-\frac{1}{2}{\displaystyle {\int}_{\pi /2}^{\pi}f(\beta \prime )\text{sin}(2\beta \prime )}d\beta \prime \end{array}$$

(74)

$$=\frac{1}{2}{\displaystyle {\int}_{-\pi}^{\pi}f(\beta \prime )\text{sin}(2\beta \prime )SQ(\beta \prime )}d\beta \prime $$

(75)

where the middle equality was obtained using symmetry arguments since *f* is a linear combination of even positive powers of sin and cos.

Since |*q*| ≤ *Q* implies that the term *f*(β′) sin(2β′) has exponential powers ≤ *Q* + 2, we can eliminate terms in the Fourier series of *SQ*(β′) that falls out of the range (by orthonormality of the exponentials), thus we only require *p* such that

$$|4p+2|\le Q+2$$

(76)

$$\iff -Q-2\le 4p+2\le Q+2$$

(77)

$$\iff -\frac{Q}{4}-1\le p\le \frac{Q}{4}$$

(78)

$$\iff -\lfloor \frac{Q}{4}\rfloor -1\le p\le \lfloor \frac{Q}{4}\rfloor .$$

(79)

Defining $\tilde{SQ}(\beta \prime )={\displaystyle {\sum}_{p=-\lfloor Q/4\rfloor -1}^{\lfloor Q/4\rfloor}-(2i/(\pi (2p+1))){e}^{i(4p+2)(\beta \prime )}}$, we have

$$\begin{array}{c}{\displaystyle {\int}_{0}^{\pi}f(\beta )\text{sin}(\beta )}d\beta \hfill \\ \text{\hspace{1em}\hspace{1em}}=\frac{1}{2}{\displaystyle {\int}_{-\pi}^{\pi}f(\beta \prime )\text{sin}(2\beta \prime )\tilde{SQ}(\beta \prime )}d\beta \prime \end{array}$$

(80)

$$\begin{array}{c}=\frac{1}{2}{\displaystyle \sum _{q=-Q}^{Q}{b}_{q}}{\displaystyle \sum _{p=-\lfloor Q/4\rfloor -1}^{\lfloor Q/4\rfloor}-\frac{2i}{\pi (2p+1)}}\hfill \\ \hfill \times \phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\int}_{-\pi}^{\pi}{e}^{i\beta \prime}q\frac{{e}^{i2\beta \prime}-{e}^{-i2\beta \prime}}{2i}}{e}^{i(4p+2)\beta \prime}d\beta \prime .\hfill \end{array}$$

(81)

Let us define the highest exponential power to be *B* and notice that *B* = *Q* + 2 + 4*Q*/4 + 2 = *Q* + 4*Q*/4 + 4, while lowest exponential power corresponds to −*Q* −2 −4*Q*/4 − 2 = −*Q* − 4*Q*/4 = −*B*. Let *N* be the smallest integer such that *B* < 4*N*, where *N* is an integer. Therefore, *N* = *Q*/4 + *Q*/4 + 1 + 1 = 2*Q*/4 + 2

It is easy to verify the following identity:

$$\begin{array}{cc}\frac{1}{4N}{\displaystyle \sum _{k=-2N}^{2N-1}{e}^{2\pi ikl/(4N)}}& =\frac{1}{2\pi}{\displaystyle {\int}_{-\pi}^{\pi}{e}^{il\beta \prime}}d\beta \prime \hfill \\ & \hfill =\{\begin{array}{ll}1,\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}l=0\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\text{\hspace{1em}\hspace{1em}}\forall \left|l\right|<4N.\end{array}$$

(82)

Substituting the identity into (81), we get β′ → 2π*k*/(4*N*)); see (83)–(87), shown at the bottom of the page, where the second

$$\begin{array}{c}\frac{1}{2}{\displaystyle \sum _{q=-Q}^{Q}{b}_{q}}{\displaystyle \sum _{p=-\lfloor Q/4\rfloor -1}^{\lfloor Q/4\rfloor}-\frac{2i}{\pi (2p+1)}}\frac{\pi}{2N}\hfill \\ \text{\hspace{1em}\hspace{1em}}\times {\displaystyle \sum _{k=-2N}^{2N-1}{e}^{2\pi ikq/(4N)}\frac{{e}^{2\pi ik2/(4N)}-{e}^{-2\pi ik2/(4N)}}{2i}{e}^{2\pi ik(4p+2)/(4N)}}\hfill \\ \text{\hspace{1em}\hspace{1em}}=\frac{\pi}{4N}{\displaystyle \sum _{k=-2N}^{2N-1}{\displaystyle \sum _{q=-Q}^{Q}{b}_{q}{e}^{2\pi ikq/(4N)}\text{sin}\left(\frac{4\pi k}{4N}\right)}}\hfill \\ \text{\hspace{1em}\hspace{1em}}\times {\displaystyle \sum _{p=-\lfloor Q/4\rfloor -1}^{\lfloor Q/4\rfloor}-\frac{2i}{\pi (2p+1)}}{e}^{2\pi ik(4p+2)/(4N)}\hfill \end{array}$$

(83)

$$=\frac{\pi}{4N}{\displaystyle \sum _{k=-2N}^{2N-1}f\phantom{\rule{thinmathspace}{0ex}}\left({\beta}_{k}^{\prime}=\frac{2\pi k}{4N}\right)\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{4\pi k}{4N}\right)}\tilde{\phantom{\rule{thinmathspace}{0ex}}SQ}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{2\pi k}{4N}\right)$$

(84)

$$=\frac{\pi}{4N}{\displaystyle \sum _{k=-2N}^{2N-1}f\phantom{\rule{thinmathspace}{0ex}}\left({\beta}_{k}^{\prime}=\frac{2\pi k}{4N}\right)\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\pi k}{N}\right)}\tilde{\phantom{\rule{thinmathspace}{0ex}}SQ}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\pi k}{2N}\right)$$

(85)

$$=\frac{\pi}{2N}{\displaystyle \sum _{k=0}^{2N-1}f\phantom{\rule{thinmathspace}{0ex}}\left({\beta}_{k}^{\prime}=\frac{2\pi k}{4N}\right)\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\pi k}{N}\right)}\tilde{\phantom{\rule{thinmathspace}{0ex}}SQ}\phantom{\rule{thinmathspace}{0ex}}(\frac{\pi k}{2N})$$

(86)

$$=\frac{\pi}{N}{\displaystyle \sum _{k=0}^{N-1}f\phantom{\rule{thinmathspace}{0ex}}\left({\beta}_{k}^{\prime}=\frac{2\pi k}{4N}\right)\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\pi k}{N}\right)}\tilde{\phantom{\rule{thinmathspace}{0ex}}SQ}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\pi k}{2N}\right)$$

(87)

last equality uses the fact that sin(−2π) = 0 and the function *f* is even, and the last equality uses the fact that sin π = 0 and *f*(β′) is even about β′ = π. Because ${\beta}_{k}^{\prime}={\beta}_{k}/2=2\pi k/(4N)$, hence β_{k} = π*k*/*N* for *k* = 0, 1,…*N* − 1 corresponds to our quadrature knots, with quadrature weights, ${w}_{k}=\pi /N\phantom{\rule{thinmathspace}{0ex}}\text{sin}(\pi k/N)\phantom{\rule{thinmathspace}{0ex}}\tilde{SQ}\phantom{\rule{thinmathspace}{0ex}}(\pi k/2N))$.

We have formulated two possible ways of obtaining an exact quadrature of ${\int}_{0}^{\pi}{d}_{mm\prime}^{l}(\beta ){d}_{mm\prime}^{l\prime}(\beta )}d\beta $ and the integral is equal to ${\sum}_{s=0}^{N-1}{w}_{s}{d}_{mm\prime}^{l}({\beta}_{s}){d}_{mm\prime}^{l\prime}({\beta}_{s})}=(2/(2l+1))\delta (l-l\prime )$ if we select the correct samples and corresponding weights. In particular, this is true for:

- β
_{s}= 2 sin^{−1}((υ_{s}+ 1 and*w*=_{s}*r*(υ_{s}_{s}+ 1),*s*= 0, 1,…*N*−1, where:*N*= [*Q*/2] + 1;*Q*is the highest power of ${d}_{mm\prime}^{l}(\beta ){d}_{mm\prime}^{l\prime}(\beta )$ when viewed as a polynomial in sin(β/2);*r*and υ_{s}_{s}are the weights and nodes of the Gaussian-Lengendre quadrature on the interval [−1, 1].

- β
_{s}= π*s*/*N*and ${w}_{s}=\pi /N\phantom{\rule{thinmathspace}{0ex}}\text{sin}(\pi s/N)\phantom{\rule{thinmathspace}{0ex}}\tilde{SQ}\phantom{\rule{thinmathspace}{0ex}}(\pi s/2N))$,*s*= 0, 1,…*N*− 1, where:*N*= 2*Q*/4 + 2;*Q*is the highest power of ${d}_{mm\prime}^{l}(\beta ){d}_{mm\prime}^{l\prime}(\beta )$ when viewed as a polynomial in*e*^{(iβ/2)}.

Note that *Q* = 2(*l* + *l′*) ≤ 2(*L _{n}* +

For brevity, we denote *y _{n}*(α

$${\widehat{x}}_{n}^{l,m}={\displaystyle \sum _{j,s,k}\frac{4{\pi}^{2}}{JK}}{w}_{s,n}{y}_{njsk}{[D({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}]}^{l,m}.$$

(88)

Since

$$\begin{array}{cc}[D(\alpha ,\beta ,\gamma ){h}_{n}{]}^{l,m}& (33\underset{=}{,}28){\displaystyle \sum _{M,M\prime}{D}_{mM}^{l}\left(\alpha +\frac{\pi}{2},\frac{\pi}{2},0\right)}\\ & \times {D}_{MM\prime}^{l}\left(\beta +\pi ,\frac{\pi}{2},\frac{\pi}{2}+\gamma \right){h}_{n}^{l,M\prime}\end{array}$$

(89)

for a particular (α_{j,n}, β_{s,n}, γ_{k,n}), we can compute [*D*(α_{j,n}, β_{s,n}, γ_{k,n})*h _{n}*]

Recall that *L _{n}*(

However, a faster algorithm can be achieved.

- ${x}_{n}^{l,m}\leftarrow 0$ for all
*l, m*. - For each β
_{s}and γ_{k}:- compute [
*D*(0,β_{s,n}, γ_{k,n})*h*]_{n}^{l,m}for all*l, m*; - for each α
_{j,n}:- compute [
*D*(0, 0,α_{j,n})*D*(0, α_{s,n},γ_{k,n})*h*]_{n}^{l,m}for all*l, m*; - .$${x}_{n}^{l,m}\leftarrow {x}_{n}^{l,m}+\frac{4{\pi}^{2}}{JK}{w}_{s,n}{y}_{njsk}{[D({\alpha}_{j,n},{\beta}_{s,n},{\gamma}_{k,n}){h}_{n}]}^{l,m}.$$

Since step 2*a* dominates step 2*b*, the runtime complexity is $O({L}_{n}^{4}{O}_{n}^{2})$.

^{1}We note that the symbol *D* is overloaded to imply rotation as well as dilation, but the meaning should be clear depending on the context.

^{2}We note that the approach commonly used with planar images of applying a constant filter to a subsampled image fails here because the sphere is periodic and bounded, causing the effective size of the image features (relative to the filter) to stay constant with subsampling. We also note that nonlinear dilation is necessary since the sphere is compact, hence dilating a spherical function by naively scaling the radial component of the spherical function, *f*(θ ϕ) → *f*((θ/*a*);ϕ), leads to undesirable “wrap-around” effects.

^{3}We use the matlab interface from the publicly available Yet Another Wavelet Toolbox (YAWTb): http://rhea.tele.ucl.ac.be/yawtb/

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

1. Antoine J-P, Vandergheynst P. Wavelets on the 2-Sphere: A Group-theoretical approach. Appl. Comput. Harmon. Anal. 1999;vol. 7:262–291.

2. Armitage C, Wandelt BD. Deconvolution map-making for cosmic microwave background observations. Phys. Rev. D 70. 2004;vol. 123007:123013–123013.

3. Bertsekas DP. Nonlinear Programming. Nashua, NH: Athena Scientific; 1998.

4. Bogdanova I, Vandergheynst P, Antoine J-P, Jacques L, Morvidone M. Stereographic wavelet frames on the sphere. Appl. Comput. Harmon. Anal. 2005;vol. 19:223–252.

5. Brechbuhler CH, Gerig G, Kubler O. Parametrization of closed surfaces for 3-D shape description. Comput. Vis. Image Understand. 1994;vol. 61(no 2):154–170.

6. Cooley JW, Tukey JW. An algorithm for the machine calculation of complex Fourier series. Math. Comput. 1965;vol. 19(no 90):297–301.

7. Daubechies I. Ten lectures on wavelets. SIAM. 1992

8. Demanet L, Vandergheynst P. Gabor wavelets on the sphere. Proc. SPIE. 2003;vol. 5207:208–215.

9. Do MN, Vetterli M. The finite ridgelet transform for image representation. IEEE Trans. Image Process. 2003 Jan;vol. 12(no 1):16–28. [PubMed]

10. Donoho DL, Johnstone IM. Adapting to unknown smoothness via wavelet shrinkage. J. Amer. Statist. Assoc. 1995;vol. 90(no 432):1200–1224.

11. Driscoll JR, Healy DM. Computing Fourier transforms and convolutions on the 2-Sphere. Adv. Appl. Math. 1994;vol. 15:202–250.

12. Freeden W, Windheuser U. Spherical wavelet transform and its discretization. Adv. Comput. Math. 1996;vol. 5:51–94.

13. Freeden W, Gervens T, Schneider M. Constructive Approximation on the Sphere: With Applications to Geomathematics. Oxford, U.K.: Clarendon; 1998.

14. Freeman WT, Adelson EW. The design and use of steerable filters. IEEE Trans. Pattern Anal. Mach. Intell. 1991 Sep;vol. 13(no 9):891–906.

15. Kostelec PJ, Rockmore DN. A Lite Version of Spharmonic Kit. [Online]. Available: http://www.cs.dartmouth.edu/geelong/sphere/

16. Kostelec PJ, Rockmore DN. Ffts on the rotation group. Proc. Santa Fe Institute Working Papers Series Paper 03-11-060. 2003. [Online]. Available http://www.cs.dartmouth.edu/geelong/

17. Mhaskar HN, Narcowich FJ, Prestin J, Ward JD. Polynomial frames on the sphere. Adv. Comput. Math. 2000;vol. 13:387–403.

18. Papoulis A. Generalized sampling expansion. IEEE Trans. Circuits Syst. 1977 Nov;vol. 24(no 11):652–654.

19. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C: The Art of Scientific Computing. 2nd ed. Cambridge, U.K.: Cambridge Univ. Press;

20. Sakurai JJ. Modern Quantum Mechanics. 2nd ed. Reading, MA: Addison Wesley; 1994.

21. Schroder P, Sweldens W. Spherical wavelets: Efficiently representing functions on the sphere; Proc. SIGGRAPH; 1995. pp. 161–172.

22. Schroder P, Sweldens W. Spherical wavelets: Texture processing. Rendering Tech. 1995 Aug;:252–263.

23. Simoncelli EP, Freeman WT, Adelson EH, Heeger DJ. Shiftable multi-scale transforms. IEEE Trans. Inf. Theory. 1992 Feb;vol. 38(no 2):587–607.

24. Simons F, Dahlen F, Wieczorek M. Spatiospectral concentration on a sphere. SIAM Rev. 2006;vol. 48(no 3):504–536.

25. Starck J-L, Candes EJ, Donoho DL. The curvelet transform for image denoising. IEEE Trans. Image Process. 2002 Jun;vol. 11(no 6):670–684. [PubMed]

26. Starck J-L, Moudden Y, Abrial P, Nguyen MK. Wavelets, ridgelets and curvelets on the sphere. Astron. Astrophys. 2006;vol. 446:1191–1204.

27. Vetterli M, Kovačevic J. Wavelets and Subband Coding. Englewood Cliffs, NJ: Prentice-Hall; 1995.

28. Wandelt BD, Gorski KM. Fast convolution on the sphere. Phys. Rev. D 63. 2001;vol. 1230002:1–6.

29. Wiaux Y, Jacques L, Vandergheynst P. Correspondence principle between spherical and Euclidean wavelets. Astrophys. J. 2005;vol. 632:15–28.

30. Wiaux Y, Jacques L, Vandergheynst P. Fast directional correlation on the sphere with steerable filters. Astrophys. J. 2006;vol. 652:820–832.

31. Yeo BTT, Ou W, Golland P. Invertible filter banks on the 2-Sphere. Proc. Int. Conf. Image Processing; 2006. pp. 2161–2164.

32. Yu P, Grant PE, Qi Y, Han X, Segonne F, Pienaar R, Busa E, Pacheco J, Makris N, Buckner R, Golland P, Fischl B. Cortical surface shape analysis based on spherical wavelets. IEEE Trans. Med. Imag. 2007 Apr;vol. 26(no 4):582–598. [PMC free article] [PubMed]

33. Yu P, Yeo BTT, Grant E, Fischl B, Golland P. Cortical folding development study based on over-complete spherical wavelets. presented at the IEEE Computer Society Workshop on Mathematical Methods in Biomedical Image Analysis.2007. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |