# Locally Supported Orthogonal Wavelet Bases on the Sphere via Stereographic Projection

Authors: Daniela Roşca; Jean-Pierre Antoine
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

The stereographic projection determines a bijection between the two-sphere, minus the
North Pole, and the tangent plane at the South Pole. This correspondence induces a unitary
map between the correspondingL2
spaces. This map in turn leads to  equivalence between
the continuous wavelet transform formalisms on the plane and on the sphere. More precisely,
any plane wavelet may be lifted, by inverse stereographic projection, to a wavelet on the sphere.
In this work we apply this procedure to orthogonal compactly supported wavelet bases in the
plane, and we get continuous, locally supported orthogonal wavelet bases on the sphere. As
applications, we give three examples. In the first two examples, we perform a singularity
detection, including one where other existing constructions of spherical wavelet bases fail. In
the third example, we show the importance of the local support, by comparing our construction
with the one based on kernels of spherical harmonics.


---

## Body

## 1. Introduction

Two-dimensional wavelets are by now a standard tool in image processing, under the two concurrent approaches, the Discrete Wavelet Transform (DWT), based on the concept of multiresolution analysis, and the Continuous Wavelet Transform (CWT). While the former usually leads to wavelet bases or frames, the CWT has to be discretized for numerical implementation and produces in general only frames [1, 2].Nowadays, many situations yield data onspherical surfaces or closed sphere-like surfaces, that is, surfaces obtained from a sphere by a smooth deformation, for instance, in Earth and Space sciences (geography, geodesy, meteorology, astronomy, cosmology, etc), in crystallography (texture analysis of crystals), in medicine (some organs are regarded as sphere-like surfaces), or in computer graphics (modelling of closed surfaces as the graph of a function defined on the sphere). So one needs a suitable analysis tool for such data. As in the flat case, the Fourier transform is a standard tool, but it amounts to an expansion in spherical harmonics, whose support is the whole sphere. Fourier analysis on the sphere is thus global and cumbersome. Therefore many different methods have been proposed to replace Fourier analysis with some sort of wavelet analysis.In addition, some data may live on more complicated manifolds, such as atwo-sheeted hyperboloid, in cosmology, for instance (an open expanding model of the universe), or a paraboloid. In optics also, data on such manifolds are essential for the treatment of omnidirectional images via the catadioptric procedure, for instance, in robotic vision. This last topic is particularly relevant for engineering purposes, because of the many applications in navigation, surveillance, and visualization. In the catadioptric image processing, a sensor overlooks a mirror, whose shape may be spherical, hyperbolic, or parabolic. However, instead of projecting the data from that mirror onto a plane, an interesting alternative consists in processing them directly on the mirror, and thus wavelets on such manifolds are needed [3]. Among the three shapes, the parabolic one is the most common (think of the headlights of a car). Now this case brings us back to the topic of this paper. Indeed it has been shown in [4] that the reconstruction of the orthographic (i.e., vertical) projection from a parabolic mirror can be computed as the inverse stereographic projection from the image plane onto the unit sphere—which is precisely the tool we are going to use in the sequel for designing wavelets on the sphere.For an efficient wavelet analysis of signals or images, including on the sphere, the following properties are desirable (a thorough discussion of this topic may be found in [5]).(i)
Basis. The redundancy of the frames leads to nonunique expansions. Moreover, the existing constructions of spherical frames (see e.g., [6, 7]) are computationally heavy and often applicable only to band-limited functions. In particular, for the compression of large data sets, it is essential to have wavelet bases, not frames.(ii)
Orthogonality. This is the most economical method since it leads to orthogonal reconstruction matrices. The inversion of such matrices, needed, for example, in data compression, is trivial (the inverse equals the adjoint). Thus, orthogonal bases are ideal for compression, but this is not always sufficient: sparsity of reconstruction matrices is still needed in the case of large data sets. We should also mention that orthogonality of a wavelet basis can be easily achieved by a Gram-Schmidt procedure, but the locality of the support is usually lost.(iii)
Local support. A wavelet has local support if it vanishes identically outside a small region. It is localized if it is negligible outside a small region, so that it may have (small, but nonzero) “tails" there. Since these tails may spread in the process of approximation of data and spoil their good localization properties, local support is definitely preferred (see the example in Figure 4). More important, local support is crucial when working with large data sets, since it yields sparse matrices.(iv)
Continuity, smoothness. These properties are always desirable in approximation, but not easily achieved.There are many constructions in literature that fulfil at least one of the properties above, so we have to restrict ourselves by mentioning only some of them. In [8–11], smooth wavelet bases were obtained by different approaches. Besides smoothness, the wavelets constructed in [11] have local support. Orthogonal, smooth and localized (but not locally supported) wavelets were constructed in [12]. In [13–15], by making use of a radial projection, one has obtained either continuous, semiorthogonal, locally supported wavelet bases, or piecewise constant orthogonal locally supported wavelet bases.However, none of these methods so far has led to wavelet bases on the sphere which aresimultaneously continuous (or smoother), orthogonal and locally supported. The aim of the present paper is precisely to fill this gap. The method we propose consists in lifting wavelets from the tangent plane to the sphere by inverse stereographic projection. It yields simultaneously smoothness, orthogonality, local support, vanishing moments. The disadvantage is that it could entail distortions around the North Pole N. Of course, mathematically the construction presented here applies to the pointed sphere 𝕊̇2:=𝕊2∖{N} (and therefore we have to avoid the data around the North Pole N), but in practice this situation is harmless. Indeed, we can choose N in a region where there are no relevant data. To give an example, European climatologists routinely put the North Pole of their spherical grid in the middle of the Pacific Ocean. Actually, most of the spherical data are given in polar coordinates not containing information at the poles.To summarize our work, we aim at detecting and quantizinglocal singularities in data, at small scales. This is, from the practical point of view, the main purpose of wavelet analysis. We will come back to this point in Section 5.


## 2. Preliminaries

We consider the unit sphere𝕊2={(η1,η2,η3)∈ℝ3,η12+η22+(η3-1)2=1} and the pointed sphere 𝕊̇2=𝕊2∖{(0,0,2)}, with the parametrization(2.1)η1=cosφsinθ,η2=sinφsinθ,θ∈(0,π],φ∈[0,2π),η3=1+cosθ,
Let p:𝕊̇2→ℝ2 denote the stereographic projection from the North Pole N(0,0,2) and let ν:𝕊̇2→ℝ, ϱ:ℝ2→ℝ be defined as(2.2)ν(η)=22-η3,∀η=(η1,η2,η3)∈𝕊̇2,ϱ(x)=4x2+y2+4,∀x=(x,y)∈ℝ2,
respectively. The relations between  η  and x=p(η) are
(2.3)x=2cotθ2cosφ=η1ν(η),y=2cotθ2sinφ=η2ν(η),
and, conversely,(2.4)η1=xϱ(x),η2=yϱ(x),η3=12(x2+y2)ϱ(x).
Then it is easy to prove the following relations between dx=dxdy, the area element of  ℝ2, and dω(η), the area element of 𝕊2:(2.5)dx=ν(η)2dω(η),dω(η)=ϱ(x)2dx.
For simplicity, we write L2(𝕊̇2)=L2(𝕊̇2,dω(η))(=L2(𝕊2), since the set {N} is of measure zero) and L2(ℝ2)=L2(ℝ2,dx). Then we consider the map π:L2(𝕊̇2)→L2(ℝ2) induced by the stereographic projection, namely,(2.6)[πℱ](x)=ϱ(x)ℱ(p-1(x)),∀ℱ∈L2(𝕊̇2)
and its inverse π-1:L2(ℝ2)→L2(𝕊̇2),(2.7)[π-1F](η)=ν(η)F(p(η)),∀F∈L2(ℝ2).
It is well known that π is a unitary map, hence(2.8)〈πℱ,π𝒢〉L2(ℝ2)=〈ℱ,𝒢〉L2(𝕊̇2),〈π-1F,π-1G〉L2(𝕊̇2)=〈F,G〉L2(ℝ2),
or, equivalently,(2.9)〈ℱ,𝒢〉L2(𝕊̇2)=〈ϱ·(ℱ∘p-1),ϱ·(𝒢∘p-1)〉L2(ℝ2),(2.10)〈F,G〉L2(ℝ2)=〈ν·(F∘p),ν·(G∘p)〉L2(𝕊̇2)
for all ℱ,𝒢∈L2(𝕊̇2),F,G∈L2(ℝ2). The last equality allows us to construct orthogonal bases on L2(𝕊̇2) starting from orthogonal bases in L2(ℝ2). More precisely, we will use the fact that, if the functions F,G∈L2(ℝ2) are orthogonal, then the functions ℱ=ν·(F∘p) and 𝒢=ν·(G∘p) will be orthogonal in L2(𝕊̇2).


## 3. Multiresolution Analysis (MRA) and Wavelet Bases ofL2(ℝ2)

In order to fix our notations, we will briefly review in this section the standard construction of 2D orthonormal wavelet bases in the flat case, starting from a multiresolution analysis (MRA) [1].LetD be a 2×2 regular matrix with the properties(a)
Dℤ2⊂ℤ2, which is equivalent to the fact that D has integer entries,(b)
λ∈σ(D)⇒|λ|>1, that is, all eigenvalues of D have modulus greater than 1.A multiresolution analysis ofL2(ℝ2) associated to D is an increasing sequence of closed subspaces Vj⊂L2(ℝ2) with ⋂j∈ℤVj={0} and ⋃j∈ℤVj¯=L2(ℝ2), and satisfying the following conditions:(1)
f∈Vj⇔f(D·)∈Vj+1,(2)
there exists a functionΦ∈L2(ℝ2) such that the set {Φ(·-k),k∈ℤ2} is an orthonormal basis of V0.As a consequence,{Φj,k:=|detD|j/2Φ(Dj·-k),k∈ℤ2} is an orthonormal basis for Vj.For eachj∈ℤ, let us define the space Wj as the orthogonal complement of Vj into Vj+1, that is, Vj+1=Vj⊕Wj. The two-dimensional wavelets are those functions which span W0. One can prove (see [16]) that there exist q=|detD|-1 wavelets 1Ψ,2Ψ,…,qΨ∈V1 that generate an orthonormal basis of W0. Therefore, {λΨj,k:=|detD|j2·λΨ(Dj·-k),λ=1,…,q,k∈ℤ2} is an orthonormal basis of Wj for each j, and {λΨj,k,λ=1,…,q,k∈ℤ2,j∈ℤ} is an orthonormal basis of L2(ℝ2).A particular case is that of tensor product wavelets, corresponding to the dilation matrixD=diag[2,2] and a 1D MRA with scaling function and mother wavelet ϕ,ψ. In this case, q=3 and one gets the 2D scaling function Φ(x,y)=ϕ(x)ϕ(y) and the three wavelets(3.1)hΨ(x,y)=ϕ(x)ψ(y),vΨ(x,y)=ψ(x)ϕ(y),dΨ(x,y)=ψ(x)ψ(y).If the one-dimensional functionsϕ and ψ have compact support, then obviously so have Φ and λΨ. This is the case of the well-known Daubechies wavelets dbn that we will use in the sequel.


## 4. Multiresolution Analysis and Orthonormal Wavelet Bases ofL2(𝕊˙2)

The construction of multiresolution analysis and wavelet bases inL2(𝕊̇2) is based on the equality (2.10). To every function F∈L2(ℝ2), one can associate the function Fs∈L2(𝕊̇2) as(4.1)Fs=ν·(F∘p).
In particular, if the functions {Fj,k}j,k are orthogonal, so are(4.2)Fj,ks=ν·(Fj,k∘p),forj∈ℤ,k∈ℤ2.
Then, taking F=Φ and F=Ψ, we obtain the spherical functions(4.3)Φj,ks=ν·(Φj,k∘p),λΨj,ks=ν·(λΨj,k∘p).
For j∈ℤ, we define the spaces 𝒱j as(4.4)𝒱j={ν·(F∘p),F∈Vj}.
Using (2.10) and the unitarity of the map π, it is immediate that 𝒱j is a closed subspace of L2(𝕊̇2), thus a Hilbert space. Moreover, these spaces have the following properties:(1)
𝒱j⊂𝒱j+1forj∈ℤ,(2)
⋂j∈ℤ𝒱j={0} and ⋃j∈ℤ𝒱j¯=L2(𝕊̇2),(3)
the set{Φ0,ks,k∈ℤ2} is an orthonormal basis for 𝒱0.We will say that a sequence of subspaces ofL2(𝕊̇2) with the properties above constitutes a multiresolution analysis of L2(𝕊̇2). Unlike the construction in [13], the functions Φ0,ks considered here are scaled versions of the same function Φs. Moreover, in condition (3), we have here a genuine orthogonality of the wavelet basis, instead of the Riesz basis obtained in [13].Remark 4.1.
It should be noticed that, in the spherical MRA defined here, the scale parameterj runs over the whole of ℤ, whereas in the usual constructions on 𝕊2, j runs only from 0 to +∞. This is due to the stereographic projection, which removes one pole, so that effectively 𝕊̇2=𝕊2∖{N} is noncompact and thus arbitrary dilations are now allowed (see, however, Section 5).Once the multiresolution analysis is determined, we construct the wavelet spaces𝒲j in the usual manner. Let 𝒲j denote the orthogonal complement of the coarse space 𝒱j in the fine space 𝒱j+1, so that(4.5)𝒱j+1=𝒱j⊕𝒲j.
One can easily prove that, for each j∈ℤ, {λΨj,ks,k∈ℤ2,λ=1,…,q} is an orthogonal basis for 𝒲j and therefore {λΨj,ks,j∈ℤ,k∈ℤ2,λ=1,…,q} is an orthonormal basis for ⊕j∈ℤ𝒲j¯=L2(𝕊̇2).The conclusion of the analysis may be summarized as follows.(i)
IfΦ has compact support in ℝ2, then Φj,ks has local support on 𝕊̇2 (indeed diam(suppΦj,ks)→0 as j→∞).(ii)
Anorthonormal 2D wavelet basis leads to an orthonormal spherical wavelet basis.(iii)
Smooth 2D wavelets lead to smooth spherical wavelets.(iv)
In particular, plane tensor product Daubechies wavelets lead to locally supported and orthonormal wavelets on𝕊̇2, and so do plane tensor product Haar wavelets.(v)
The decomposition and reconstruction matrices needed in the spherical case are the same as in the plane 2D case, so that the latter can be used (with existing toolboxes).


## 5. Further Comments

The relationdω(η)=ϱ(x)2dx linking the area elements of ℝ2 and of 𝕊2 implies that, near the origin of the plane, corresponding to the South Pole, the areas on the sphere and on the plane are almost the same, but near the North Pole, the ratio of the areas approach infinity. In the standard cases, the elements of a planar orthogonal wavelet basis have supports of the same size all over the plane. But the corresponding spherical wavelets obtained here will have a support that almost vanishes when they are in the neighborhood of the North Pole, which corresponds to regions far away from the origin in the plane. This might clearly lead to analytical as well as numerical problems whenever the information is nearly uniform over the whole sphere, and here we see the limitation of our method. As we have precised already at the end of Section 1, we have to avoid data around the North Pole. More precisely, our method is applicable on the whole sphere, except a spherical cap θ⩽θ0. The value of θ0 may actually depend on the type of data.We emphasize that, in practice, wavelet transforms are useful mainly for alocal analysis. Thus we do not pretend that our method is the best for all kinds of applications. Actually, no construction of spherical wavelets is perfect, as we have mentioned in Section 1. However, we claim that our method is more efficient than most other ones when one deals with large amounts of data, in all applications where one has to make a decomposition and a reconstruction, which implies to invert a large matrix. In the case of Daubechies wavelets, this matrix becomes orthogonal and sparse. Here we see the advantage of being able to use Daubechies wavelets rather than other 2D wavelets.We can be more precise about the latitude effect, that is, the choice of the limiting valueθ0. When we project stereographically onto the plane xOy spherical data situated outside of the spherical cap θ⩽θ0⩽π/2, the projected data will be situated inside the square [-xmax,xmax]×[-xmax,xmax], where xmax=2cot(θ0/2). In the case of tensor product wavelets, the planar grid at level j has the dimension nj×nj and is taken as(5.1){(-xmax+k12j,-xmax+k22j),k1,k2=1,2,…,nj}.
Here nj is taken such that -xmax+(nj/2j)=xmax, that is, nj=⌊2j+1·xmax⌋. Table 1 gives the size of the planar grid for j=4 and j=5 (and therefore of the matrix which is used in the decomposition algorithm) for different values of θ0. For instance, in our first two examples below, we will use the values θ0=8° and j=4. This gives matrices of reasonable size, yet yields good results.Table 1
Dimension of the planar grid at levelj (j=4 and j=5), for a data set localized at a colatitude θ≥θ0.
θ0 (deg.)θ0 (rad.)xmaxdimension of planar gridnjj=4j=510.0175229.177373331466620.0349114.60153666733230.052476.37692444488840.069857.27251832366450.087345.80751466293260.104738.16231221244270.122232.69971046209280.139628.6013915183090.157125.41248131626100.174522.86017311462150.261815.1915486972200.349111.3426363726


## 6. Examples

In order to illustrate the merits of our construction, we present three examples.In the first example, we take the functiongθ0:𝕊2→ℝ,gθ0(θ,φ)=g(θ+θ0,φ), where(6.1)g(θ,φ)={1,θ⩾π2,(1+3cos2θ)-1/2,θ⩽π2.
The function gθ0 and its gradient are continuous, but the second partial derivative with respect to θ has a discontinuity on the latitudinal circle θ=(π/2)-θ0 (see Figure 1(a)). Detecting properly such a discontinuity requires a wavelet with two vanishing moments at least, so that none of the existing constructions of discrete wavelets mentioned in Section 1 could detect this discontinuity.(a) The graph of the functiong0∘(θ,φ). (b) Its analysis with the spherical wavelet associated to the Daubechies wavelet db3; (c) Same analysis by the discretized CWT method with the wavelet ψH2s, at scale a=0.0165. (d)–(f) The same for the function g50∘(θ,φ). (g)–(i) The same for the function g80∘(θ,φ).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Thus, following [6, 17], we take the discretized spherical CWT with the spherical wavelet ΨH2s, defined as in (4.1), associated to the plane wavelet(6.2)ΨH2(x,y)=Δ2[e-(1/2)(x2+y2)]=(x4+y4+2x2y2-8(x2+y2)+8)e-(1/2)(x2+y2).
The plane wavelet (6.2) has vanishing moments of order up to 3 (here too, a simpler wavelet with less vanishing moments could not detect the discontinuity). In order to compare the two methods, we present in Figure 1 the analysis of the function gθ0(θ,φ) for three different values of θ0 (a) θ0=0°; (d) θ0=50°; (g) θ0=80°. The second column shows, in panels (b), (e), (h), the result of the analysis with the spherical wavelet associated to the Daubechies wavelet db3. Finally, the third column, on panels (c), (f), and (i), does the same for the analysis by the discretized CWT method with the wavelet ψH2s.Comparing the second and the third column, it is clear that the Daubechies waveletdb3 lifted on the sphere by (4.1) outperforms the wavelet ΨH2s mentioned above. The precision is much better, in the sense that the width of the detected singular curve is narrower at all three latitudes. With both methods, the detection is still efficient very close to the North Pole, but the computation with the lifted db3 wavelet then requires considering matrices of large dimension, as shown in Table 1.Looking at the third column suggests to increase the resolution of the CWT method, by choosing smaller scales. In order to test this effect, we show in Figure2 a systematic analysis of the function gθ0∘(θ,φ) by the discretized CWT method with the wavelet ψH2s. Panels (a), (b), (c) and (d) present the spherical CWT at smaller  and smaller scales, a=0.04,0.0165,0.012 and 0.0085, respectively. The discretization grids are those considered in [6]. From the panels (a)–(c), it appears that the discontinuity along the equator is detected properly, and the precision increases as the scale decreases. However, there is a limit: when the scale a is taken below a=0.01, the singularity is no more detected. For instance, for a=0.0085 (see panel (d)), the wavelet becomes too narrow and “falls in between" the discretization points, and “ripples" appear. This effect is described in detail in [17]. Here again, the spherical Daubechies wavelet performs better.Analysis of the functiong0∘(θ,φ) by the discretized CWT method with the wavelet ψH2s, at scales: (a) a=0.04; (b) a=0.0165; (c) a=0.012; and (d) a=0.0085.
(a)(b)(c)(d)The next question to examine is the possible distortion in latitude around the North Pole, for a signal which is no longer invariant under rotation around thez-axis. Thus, in the second example, we take the function f:𝕊2→{0,1},(6.3)f(θ,φ)={0,ifθ⩾π2,1,ifθ<π2.
and then fα, the function obtained from f by performing a rotation around the axis Ox by an angle α. We apply a tensor product Haar wavelet transform to fα, for α=70°. Let us proceed in two steps. First, the discontinuity is well detected in the tangent plane (where the wavelet transform is computed), as shown in Figure 3(a). It is indeed a circle, with uniform intensity, as it should, in virtue of the conformal character of the stereographic projection. However, when it is lifted back onto the sphere, some artifacts appear, that can be seen in Figure 3(b). Namely, the part of the discontinuity that lies closer to the North Pole is less clearly visible than the ones that are farther away from it and there are some gaps. This is caused by the unavoidable fact that the grids on the plane and on the sphere never match, the inverse projection does not conserve areas, so that some points are lost in the lifting process.Analysis of the functionf70∘(θ,φ) with the spherical wavelet associated to a Haar wavelet: (a) the wavelet transform in the tangent plane; (b) the same lifted back on the sphere, (c) analysis of the same function by the discretized CWT method with the wavelet ψH2s, at scale a=0.0012.
(a)(b)(c)(a) An example of kernel of spherical harmonics: localized, but not locally supported; (b) the initial data set, in spherical coordinates(θ,φ); (c) Its approximation at level 6 using a spherical harmonics kernel; (d) approximation at first level using the spherical wavelet associated to the Daubechies wavelet db2.
(a)(b)(c)(d)On the other hand, the analysis of the same function by the discretized spherical CWT, at the smallest possible scalea=0.0012, shows no such effect, as seen in Figure 3(c). But here too, the singularity is detected by the discrete WT with a better precision. Note that, in the CWT case, simpler planar wavelets, with less vanishing moments, can be used, but the result is the same. One can also try to analyze the function obtained from g0° by a rotation around the axis Ox, but we do not expect much difference with the present example.In the third example, we consider a data set from the texture analysis of crystals; see Figure4(b). It consists of 180×360 measurements on the sphere and its main characteristic is that the values over the sphere are constant, except for some peaks. Figure 4(c) shows the approximation at level 6 using the spherical wavelet frames constructed in [7], which are localized, but not locally supported. An example of such kernel is given in Figure 4(a). Figure 4(d) shows the approximation at level 1 using the spherical wavelet associated to db2. One can easily see that our wavelets are more efficient in approximating the given data set and also the algorithms are much faster than the ones for spherical harmonics.


## 7. Conclusion

We have presented a new method for constructing spherical wavelets. Its advantage over most existing methods is that it yields locally supported orthonormal bases, which is crucial when one is confronted to very large data sets. The price to pay is that one has to avoid a small region around the North Pole, whose size may depend on the data. Within that region, serious distortions may occur. However, the examples given above show that this “forbidden" region can be very small: although we have usedθ0=70° or 80°, much higher values can be chosen without damage. Hence we believe this limitation does not prevent applying the method to practical situations. In any case, if the data are sufficiently well localized, one can always make a rotation that brings them around the South Pole, as we already mentioned in the Introduction. After all, one should remember that wavelet analysis is primarily a local analysis.We may also mention that the method described in this paper could detect singularities of arbitrarily high order, by choosing appropriate planar wavelets. By contrast, the other existing spherical wavelets do not have enough vanishing moments and thus are not able to detect such discontinuities.As a final remark, we may point out that the method should work for any manifold with a bijective orthogonal projection onto a fixed plane, since the latter induces a unitary map between the respectiveL2 spaces. Such are, for instance, the upper sheet of the two-sheeted hyperboloid z2-x2-y2=1 with projection onto the plane z=0 or the paraboloid z=x2+y2 with the same projection (in addition, we have seen in Section 1 that the latter case may be treated also by the stereographic projection from the sphere). It might also be extended to a local analysis on more general manifolds. A systematic discussion may be found in our recent work [18].



---

*Source: 124904-2009-10-18.xml*