Scattered Data Fitting on the Sphere

Gregory E. Fasshauer and Larry L. Schumaker

Abstract. We discuss several approaches to the problem of interpolating or approximating data given at scattered points lying on the surface of the sphere. These include methods based on spherical harmonics, tensorproduct spaces on a rectangular map of the sphere, functions defined over spherical triangulations, spherical splines, spherical radial basis functions, and some associated multi-resolution methods. In addition, we briefly discuss sphere-like surfaces, visualization, and methods for more general surfaces. The paper includes a total of 206 references.

§1. Introduction

Let S be the unit sphere in IR^{3}, and suppose that {v_{i}}_{i=1}^{n} is a set of scattered points lying on S. In this paper we are interested in the following problem:

Problem 1. Given real numbers {r_{i}}_{i=1}^{n}, find a (smooth) function s defined on S which interpolates the data in the sense that

s(v_{i}) = r_{i},  i = 1, …, n,  (1)

or approximates it in the sense that

s(v_{i}) ≈ r_{i},  i = 1, …, n.  (2)

Data fitting problems where the underlying domain is the sphere arise in many areas, including e.g. geophysics and meteorology where the sphere is taken as a model of the earth. The question of whether interpolation or approximation should be carried out depends on the setting, although in practice measured data are almost always noisy, in which case approximation is probably more appropriate.

In most applications, we will want s to be at least continuous. In some cases we may want it to be C^{1} so that the associated surface F := {s(v)v : v ∈ S} is tangent plane continuous.

Mathematical Methods for Curves and Surfaces II 117 Morten Dæhlen, Tom Lyche, and Larry L. Schumaker (eds.), pp. 117–166. Copyright o^{c} 1998 by Vanderbilt University Press, Nashville, TN. ISBN 0-8265-1315-8. All rights of reproduction in any form reserved.

The aim of this paper is to survey the spectrum of methods which have been developed (mostly in the past ten years) for solving Problem 1. The paper is divided into sections as follows:

1) Introduction,

2) Spherical harmonics,

3) Methods based on mapping the sphere to a rectangle,

4) Methods based on dividing the sphere into subsets consisting of spherical triangles,

5) Spherical splines (piecewise spherical harmonics),

6) Methods based on linear combinations of radial basis functions,

7) Multiresolution methods,

8) Additional topics, including sphere-like surfaces, general surfaces, visualizing surfaces on surfaces, and numerical quadrature on the sphere.

We conclude the paper with a bibliography containing 206 references.

§2. Spherical Harmonics

Many classical interpolation and approximation methods are based on polynomials. The appropriate analog of polynomials on the sphere are called spherical harmonics. They can be defined in several different (equivalent ways) which we now discuss. For details, see e.g. [24,80,105,118,165].

Let P_{d} be the space of trivariate polynomials of total degree at most d, and let H_{d} := P_{d}|_{S} be its restriction to the sphere. A trivariate polynomial p is called homogeneous of degree d provided p(λx, λy, λz) = λ^{d}p(x, y, z) for all λ ∈ IR. It is called harmonic provided Δp ≡ 0, where Δ is the Laplace operator defined by Δf := (D_{x}^{2} + D_{y}^{2} + D_{z}^{2})f.

Definition 2. The linear space

H_{d} := {p|_{S} : p ∈ P_{d} and p is homogeneous of degree d and harmonic}

is called the space of spherical harmonics of exact degree d.

It is well known (see e.g. [24], page 314) that the dimension of H_{d} is 2d+1, and that it is the eigenspace corresponding to the eigenvalue λ_{d} = −d(d + 1) of the Laplace-Beltrami operator Δ^{∗} on S given by

Δ^{∗}f = \frac{1}{sin^{2} θ} D_{ϕ}^{2}f + \frac{1}{sin θ} D_{θ} (sin θ D_{θ}f),  (3)

where θ ∈ [0, π] and ϕ ∈ [0, 2π] are the spherical coordinates of a point on S. It is also known that H_{d} is the orthogonal complement of H_{d−1} in the space H_{d} with respect to the L_{2}-inner product on S. Using this fact repeatedly, it follows that

H_{d} = H_{d} ⊕ H_{d−1} ⊕ ⋯ ⊕ H_{0},  (4)

and thus dim(H_{d}) = (d + 1)^{2}.

For applications, it is important to have an explicit basis for H_{d}. The classical construction (cf. [24,80,105,118,165]) is based on the spherical coordinates θ and ϕ. Let P_{d} be the Legendre polynomial of degree d normalized such that P_{d}(1) = 1, and let P_{d}^{ℓ} be the associated Legendre function of degree d and order ℓ defined by

P_{d}^{ℓ}(x) = (1 − x^{2})^{ℓ/2}D_{x}^{ℓ}P_{d}(x), −1 ≤ x ≤ 1.

Then the functions

Y_{d,2ℓ+1}(θ, ϕ) := cos(ℓϕ) P_{d}^{ℓ}(cos θ),  ℓ = 0, …, d,
Y_{d,2ℓ}(θ, ϕ) := sin(ℓϕ) P_{d}^{ℓ}(cos θ),  ℓ = 1, …, d,

form an orthogonal (but not orthonormal) basis for H_{d}. Each of the Y_{d,ℓ} can be expanded in terms of sine and cosine functions. The formulae are simple for d = 0, 1. Indeed, Y_{0,1}(θ, ϕ) = 1 and

Y_{1,1}(θ, ϕ) = cos(θ),
Y_{1,2}(θ, ϕ) = sin(ϕ) sin(θ),
Y_{1,3}(θ, ϕ) = cos(ϕ) sin(θ).

The formulae become increasingly complicated for larger values of d.

As shown in [70], it is also possible to construct a basis for H_{d} directly in terms of Cartesian coordinates. Every trivariate polynomial p in P_{d} which is homogeneous of degree d has the form

g(x, y, z) = ∑_{i+j+k=d} a_{ijk}x^{i}y^{j}z^{k}.  (5)

Now g will be harmonic if and only if the coefficients of all powers of x, y, z which remain in the equation Δg = 0 are zero. Then a basis for H_{d} can be constructed by finding linearly independent vectors a satisfying the corresponding homogeneous system of equations Ca = 0, where a is the coefficient vector in (5). Applying this process for d = 0, 1, 2, one finds the bases

H_{0} = span{1},
H_{1} = span{x, y, z},
H_{2} = span{xy, xz, yz, x^{2} − y^{2}, x^{2} − z^{2}}.

Both of the above basis constructions are somewhat cumbersome. Some extremely convenient basis functions (which depend on certain rotation invariant spherical barycentric coordinates) will be constructed in Sect. 5.2 for the spaces

B_{d} := \{ \begin{matrix}H_{0} ⊕ H_{2} ⊕ ⋯ ⊕ H_{2k}, & d = 2k, \\ H_{1} ⊕ H_{3} ⊕ ⋯ ⊕ H_{2k+1}, & d = 2k + 1. \end{matrix}  (6)

One of the things that make spherical harmonics interesting is the fact that smooth functions on S can be approximated well by combinations of spherical harmonics (of sufficiently high degree). One way to construct such approximations is via the spherical harmonic expansion (also called the Laplace expansion)

f_{d} = ∑_{k=0}^{d} [ ∑_{ℓ=0}^{k} a_{k,ℓ}Y_{k,2ℓ+1} + ∑_{ℓ=1}^{k} b_{k,ℓ}Y_{k,2ℓ}],  (7)

where

\begin{aligned}a_{k,ℓ} & = \frac{(2k + 1)(k − ℓ)!}{2π(k + ℓ)!} ∫_{0}^{2π} ∫_{0}^{π} f(θ, ϕ)Y_{k,2ℓ+1}(θ, ϕ) sin θ dθ dϕ, \\ b_{k,ℓ} & = \frac{(2k + 1)(k − ℓ)!}{2π(k + ℓ)!} ∫_{0}^{2π} ∫_{0}^{π} f(θ, ϕ)Y_{k,2ℓ}(θ, ϕ) sin θ dθ dϕ. \end{aligned} (8)

This is the analog of the usual Fourier series expansion on the unit circle. It can be shown that for any function f ∈ L_{2}(S), ∥f − f_{d}∥_{2} → 0 as d → ∞, while for functions f ∈ C^{2}(S), ∥f − f_{d}∥_{∞} → 0, see [24], page 513, and also [145,194–196].

There is an extensive literature on the general question of how well smooth functions defined on the sphere can be approximated by spherical harmonics, including both direct and inverse theorems [5,13,44,50,51,71,83,94,98– 101,118,125,128–132,144,145,148–152,157–162,178,186–188,190–192,202,203]. We have space here for just one result which can be regarded as an analog of the classical Jackson’s theorem for polynomial approximation. The statement concerns functions lying in a spherical analog of the classical Sobolev space. Following [100], let W_{p}^{r}(S) be the space of all functions f ∈ L_{p}(S) such that

∥δ_{γ}^{κ}(Δ^{∗})^{ρ}f∥_{L_{p}(S)} ≤ Mγ^{r−2ρ}, 0 < γ < π,

where δ_{γ}^{κ} is a recursively defined spherical difference (computed as the difference between f at some point v and its average on some circle of geodesic radius γ around v), Δ^{∗} is the Laplace-Beltrami operator (3), M is some constant independent of γ, and the integers ρ and κ satisfy 2κ > r −2ρ > 0. Note that the definition of W_{p}^{r}(S) depends on the choice of κ and γ.

Theorem 3. There exists a constant C > 0 such that for every f ∈ W_{p}^{r}(S), there is a function g ∈ H_{d} such that

∥f − g∥_{L_{p}(S)} ≤ Cd^{−r}∥f∥_{W_{p}^{r}(S)}.

Theorem 3 is a direct theorem. For some related inverse theorems, see [100,118]. For other results involving moduli of smoothness on the sphere, see [50,51,132,144].

Formula (7) is not designed for fitting data, but it can be used if there is sufficient data available to create numerical approximations to the integrals in (8). For results based on gridded data, see [7,29].

Since the dimension of H_{d} is n := (d + 1)^{2}, it is natural to use it to interpolate data at a set of n scattered points on S. Writing f in the form (7), this leads to the linear system

∑_{k=0}^{d} [ ∑_{ℓ=0}^{k} a_{k,ℓ}Y_{k,2ℓ+1}(θ_{i}, ϕ_{i}) + ∑_{ℓ=1}^{k} b_{k,ℓ}Y_{k,2ℓ}(θ_{i}, ϕ_{i})] = r_{i},  i = 1, …, n, (9)

for the (d + 1)^{2} coefficients. It is now of interest to identify those point sets for which this system is nonsingular.

Definition 4. A set of points {(θ_{i}, ϕ_{i})}_{i=1}^{n} for which the system in (9) is nonsingular is called H_{d}-unisolvent.

While it is known [118] that unisolvent sets exist, we have not been able to find a general characterization of them in the literature. A special result for points on parallel circles on the sphere can be found in [73]. The papers [84,135,149,150,152] also deal with interpolation.

In contrast, the analogous question of B_{d}-unisolvency for the spaces B_{d} appearing in (6) has a very satisfactory answer.

Theorem 5. Given d, let m := \binom{d+2}{2}. Let v_{1}, …, v_{m} be a set of points on the sphere. Suppose that for each 1 ≤ i ≤ m, there exists a set of d distinct great circle arcs such that

1) v_{i} does not lie on any of the arcs,

2) all of the other v_{j} lie on at least one of the arcs. Then the set {v_{i}}_{i=1}^{m} is a B_{d}-unisolvent set.

This result is an analog of a classical result of Chung and Yao [23] (see also [22]) for interpolation by multivariate polynomials. For a proof and related results based on a study of homogeneous polynomials, see [89–91].

Spherical harmonics have been heavily used for fitting data, particularly in geophysics and meteorology. We make no attempt to list application papers here, but for one interesting application, see [134].

We conclude this section with one final remark. Since spherical harmonics are the direct analog of polynomials, one can expect that they suffer from the same problem as univariate or bivariate polynomials — a tendency to oscillate due to their global nature and lack of flexibility. One way to alleviate this would be to work with piecewise spherical harmonics, which is exactly what is done with the spaces of spherical splines discussed in Sect. 5.

§3. Methods Based on Mapping S to a Rectangle

The idea of fitting scattered data on the sphere by converting the problem to one defined on a rectangle has been exploited in several papers, see e.g. [28– 30,72,177,193]. Early methods suffered from problems at the poles, but this can be overcome with the use of trigonometric splines as discussed in Sect. 3.3 below.

3.1. Mapping S to a Rectangle

We denote the north pole of S by v_{N}, and the south pole by v_{S}. The inverse of the mapping

Φ(θ, ϕ) := ( \begin{matrix}sin(θ) cos(ϕ) \\ sin(θ) sin(ϕ) \\ cos(π − θ) \end{matrix} )  (10)

maps points on S \ {v_{N} , v_{S}} onto

R_{I} := {(θ, ϕ) : 0 < θ < π and 0 ≤ ϕ ≤ 2π}.

A function s defined on the rectangle

R := {(θ, ϕ) : 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ 2π}  (11)

is well-defined on the entire sphere S if and only if it is 2π-periodic in ϕ and has constant values at both the south and north poles, i.e.,

s(0, ϕ) = r_{S} and s(π, ϕ) = r_{N} , 0 ≤ ϕ ≤ 2π, (12)

for some r_{S} and r_{N}.

This identification of the sphere S with the rectangle R makes it possible to recast data fitting problems on S as data fitting problems on R, whereby points v_{i} ∈ S are associated with points (θ_{i}, ϕ_{i}) ∈ R. For example, to solve the interpolation problem (1), it suffices to construct a function defined on R satisfying (12) and

s(θ_{i}, ϕ_{i}) = r_{i}, i = 1, …, n. (13)

Suppose s is such a function. In order to make the corresponding surface F := {s(θ, ϕ)Φ(θ, ϕ) : (θ, ϕ) ∈ R} be smooth, it is necessary to impose some smoothness on s. For example, to get a continuous surface, it is enough to require that s ∈ C(R) and that it be 2π periodic in ϕ. To get a continuously varying tangent plane (everywhere except at the poles), s ∈ C^{1}(R_{I}) and D_{θ}s must be 2π periodic in ϕ.

3.2. Tensor-Product Polynomial Splines

Since the fitting problem has been transformed to a rectangle, it is natural to solve it using tensor product functions of the form

s(θ, ϕ) := ∑_{i=1}^{M} ∑_{j=1}^{N} c_{i,j}B_{i}(θ)T_{j}(ϕ),  (14)

where B_{1}(θ), …, B_{M}(θ) are functions defined on [0, π], and T_{1}(ϕ), …, T_{N}(ϕ) are 2π-periodic functions defined on [0, 2π].

The obvious choice for both sets of functions would be polynomial Bsplines. For the B_{i} one can take B-splines of order m with knots

0 = x_{1} = ⋯ x_{m} < x_{m+1} < ⋯ < x_{M} < x_{M+1} = ⋯ = x_{M+m}. (15)

For the T_{j} one can take periodic B-splines of order n, where the knots are chosen periodically (cf. [175]):

\begin{matrix}0 = y_{n} < ⋯ < y_{N+n} = 2π, \\ y_{j} = y_{j+N} − 2π,  j = 1, …, n − 1. \end{matrix}  (16)

3.3. Trigonometric Splines and Tangent Plane Continuity

In many applications one wants the surface F associated with a function s to be tangent-plane continuous everywhere, including the poles. It was shown in [28] (see also [72]) that this holds if and only if s ∈ C^{1}(R_{I}), s, D_{θ}s and D_{ϕ}s are 2π periodic in ϕ, and

D_{θ}s(π, ϕ) = A_{N} cos(ϕ) + B_{N} sin(ϕ), 0 ≤ ϕ ≤ 2π,  (17)
D_{θ}s(0, ϕ) = A_{S} cos(ϕ) + B_{S} sin(ϕ), 0 ≤ ϕ ≤ 2π,  (18)

at the north and south poles, respectively, where A_{S}, B_{S}, A_{N}, and B_{N} are constants.

These conditions cannot be satisfied using tensor-product polynomial splines, and so the methods in [28,72] are based on satisfying them only approximately. However, it was shown in [177] that they can be satisfied exactly if the T_{j} are chosen to be periodic trigonometric B-splines (cf. [175]) of odd order m defined on the periodic knot sequence in [0, 2π] described above. The reason for the restriction to odd order is that (12) can only be satisfied if the trigonometric spline space contains constants, which is the case only when m is odd.

Theorem 6. [177] Let s be defined as in (14) where the B_{i} are polynomial splines of order m with knots (15) and the T_{j} are periodic trigonometric splines of order 3 with knots (16). Then the resulting surface F is continuous and tangent plane continuous at all points of S if and only if

\begin{aligned}c_{1,j} & = r_{S} cos ((y_{j+2} − y_{j+1})/2) , \\ c_{M,j} & = r_{N} cos ((y_{j+2} − y_{j+1})/2) , \\ c_{2,j} & = c_{1,j} + (x_{m+1} − x_{m})(A_{S}α_{j} + B_{S}β_{j})/(m − 1), \\ c_{M−1,j} & = c_{M,j} − (x_{M+1} − x_{M})(A_{N} α_{j} + B_{N} β_{j})/(m − 1), \end{aligned}

where α_{j} := cos((y_{j+1} + y_{j+2})/2) and β_{j} := sin((y_{j+1} + y_{j+2})/2), for j = 1, …, N.

3.4. Two-Stage Processes

Tensor-product splines are ideal for fitting gridded data, but not so well-suited to fitting scattered data. One way to handle scattered data is to perform a two-stage process, where the first stage is to use any convenient (local) method to compute values at points on a rectangular grid. This is most conveniently done with a local trigonometric spline quasi-interpolant. An appropriate such quasi-interpolant was given in [177] for m = 3. For a general treatment of trigonometric spline quasi-interpolants, see [104].

The fitting methods discussed in this section lend themselves to use in a multiresolution analysis. We discuss this in detail in Sect. 7.2.

3.5. Map and Blend Methods

Another way to avoid the problems at the poles is to work with two different rectangles, one of which leaves out regions surrounding the north and south pole, and the other of which leaves out regions surrounding the east pole and west pole (which correspond to the points (0, 1, 0) and (0, −1, 0) in Cartesian coordinates). Then the two fits can be blended together, and if both interpolate (assuming there is no data at the poles), it is even possible to arrange for the blended function to also interpolate. For a discussion of some methods based on this idea, see [54].

3.6. Tensor Trigonometric Splines

It is also possible to choose both of the B_{i} and T_{j} in (14) to be trigonometric splines. Since using trigonometric splines on circular arcs is equivalent to using the circular Bernstein-Bézier (CBB) polynomials discussed in [1], this amounts to working with tensor-product CBB-polynomials.

§4. Methods Based on Spherical Triangulations

In this section we discuss several interpolation methods based on spherical triangulations.

4.1. Spherical Triangulations

Any two points v, w ∈ S determine a unique great circle which splits into two pieces. If v and w are not antipodal, then one of these pieces is shorter – we call it the great circle arc connecting v and w, and denote it by ⟨v, w⟩. The length of ⟨v, w⟩ is called the geodesic distance from v to w.

Definition 7. Given three unit vectors v_{1}, v_{2}, v_{3} which span IR^{3}, the associated spherical triangle T = ⟨v_{1}, v_{2}, v_{3}⟩ is defined to be the set

T = {v ∈ S : v = b_{1}v_{1} + b_{2}v_{2} + b_{3}v_{3}, b_{i} ≥ 0}.  (19)

The boundary of T consists of the three great circle arcs ⟨v_{1}, v_{2}⟩, ⟨v_{2}, v_{3}⟩, ⟨v_{3}, v_{1}⟩. A set of triangles Δ := {T_{i}}_{1}^{N} lying on a sphere S is called a spherical triangulation provided that S = ∪T_{i}, and any two triangles intersect only at a common vertex or along a common edge.

As in the planar case, in general there are many different triangulations associated with a given set of vertices V := {v_{i}}_{i=1}^{n}. However, the well-known Euler formulae N = 2n−4 and E = 3n−6 hold for any spherical triangulation, where N is the number of triangles and E is the number of edges of Δ.

By choosing appropriate criteria for comparing triangulations, one can define various types of optimal spherical triangulations, including an analog of the classical Delaunay (Thiessen) triangulation [88,127,153]. Fortran code for constructing the Delaunay triangulation on a sphere is discussed in [154] and is available from netlib.

4.2. The Hermite Interpolation Problem

Several of the methods to be discussed below actually solve a kind of Hermite interpolation problem where in addition to the data (1), derivative information at each of the data points is also given. To describe this in more detail, we need to discuss the derivative of a function defined on the sphere.

Given a point v ∈ S, let g be a unit vector in IR^{3} which is tangent to S at v. Then g together with the center of the sphere defines a plane which cuts the sphere along a great circle C passing through v. Now given a function f defined on S, the directional derivative in the direction g at v can be defined as

D_{g}f(v) := lim_{{w→v \atop w∈C}} \frac{f(w) − f(v)}{d(v, w)},

where d(v, w) is the geodesic distance from v to w. For an alternative definition involving the gradient of a homogeneous extension of f, see Sect. 5.3.

Definition 8. Let v_{1}, …, v_{n} be scattered points on the sphere. Suppose that for each i = 1, …, n, we are given two linearly independent tangent directions g_{i}, h_{i} associated with v_{i}. Then the Hermite interpolation problem is to find a function s defined on the sphere so that

s(v_{i}) = r_{i}, D_{g_{i}}s(v_{i}) = r_{i}^{g}, D_{h_{i}}s(v_{i}) = r_{i}^{h},  i = 1, …, n,  (20)

where {r_{i}, r_{i}^{g}, r_{i}^{h}}_{i=1}^{n} are given real numbers.

In many data fitting situations we will not be given derivative information. However, Hermite interpolation methods can still be applied if we can estimate the derivative information in some reasonable way, see Sect. 5.5.5.

4.3. Transfinite Interpolation Methods

In [87,88] the Hermite interpolation problem (20) is solved by adapting two planar transfinite interpolation methods to the sphere. In particular, a version of the BBG method (see [8]) and of the side vertex method (see [126]) are given. These methods work on one triangle at a time, and are based on first interpolating between its vertices to create a function on the edges of T, and then extending this function to the interior of T by interpolating between a vertex and the opposite side. A variant of these methods which is based on mapping spherical triangles to flat ones instead of working with geodesic distances can be found in [153,154].

4.4. Minimum Norm Networks

The interpolation methods described in the above subsection are based on creating a curve network by simple interpolation. As an alternative, one can build the curve network by minimizing the energy functional (see e.g. [127,147])

σ(f) := ∑_{e_{ij}} \frac{1}{∥v_{i} − v_{j}∥} ∫ [D_{ij}f(t)]^{2}dt,

where e_{ij} := ⟨v_{i}, v_{j}⟩ is the great circle arc connecting v_{i} and v_{j}, t denotes arc length, and D_{ij} is the derivative in the direction described by the arc e_{ij}.

In [137] a similar energy functional based on chord length instead of geodesic distance was used to construct minimum norm networks in a more general setting.

§5. Spherical Splines

In this section we discuss certain linear spaces of splines defined on spherical triangulations which are natural analogs of the classical polynomial splines on planar triangulations. Our treatment is based on [2,3,4], and follows the notation introduced there. We begin by discussing how to define spherical barycentric coordinates on a spherical triangle and how to use them to define spherical functions which are the analogs of Bernstein-Bézier polynomials. Spaces of spherical splines associated with a spherical triangulation are described in Sect. 5.4.

5.1. Spherical Barycentric Coordinates

For many years people in the CAGD community believed it to be impossible to define barycentric coordinates on a spherical triangle. And indeed, it is impossible (cf. the discussion in [16]) if one insists that they sum to 1. However, it was recognized in [2] that a nice theory can be developed without this condition, and that in fact there is a very natural way to define barycentric coordinates with respect to spherical triangles. (It was later discovered that the same coordinates had been introduced and studied more than 100 years ago by Möbius [116]).

Given a nondegenerate spherical triangle T := ⟨v_{1}, v_{2}, v_{3}⟩, every point v ∈ S has a unique representation of the form v = b_{1}v_{1} + b_{2}v_{2} + b_{3}v_{3}. The b_{1}(v), b_{2}(v), b_{3}(v) are called the spherical barycentric coordinates of v relative to T. They are infinitely differentiable functions of v, are nonnegative for all v ∈ T, and satisfy

b_{i}(v_{j}) = δ_{ij}, i, j = 1, 2, 3.

However, in contrast to the usual barycentric coordinates associated with a planar triangle, they do not add up to 1. Instead,

b_{1}(v) + b_{2}(v) + b_{3}(v) > 0, all v ∈ T .

As shown in [2], spherical barycentric coordinates are rotation invariant. Several equivalent formulae for computing them in terms of angles and arc lengths can also be found there.

5.2. Spherical Bernstein-Bézier Polynomials

Given a spherical triangle T and an integer d, the associated spherical Bernstein basis functions of degree d are defined to be the functions

B_{ijk}^{d} := \frac{d!}{i!j!k!}b_{1}^{i}b_{2}^{j}b_{3}^{k}, i + j + k = d.

These \binom{d+2}{2} functions are linearly independent [2]. We write B_{d} for their span.

It was observed in [2] that the B_{ijk}^{d} are actually linear combinations of spherical harmonics. The following proposition shows the precise connection with the spherical harmonic spaces H_{m} introduced in Sect. 2.

Proposition 9. B_{0} = H_{0}. Moreover, H_{d} = B_{d} ⊕ B_{d−1} and

B_{d} = \{ \begin{matrix}H_{0} ⊕ H_{2} ⊕ ⋯ ⊕ H_{2k}, & d = 2k, \\ H_{1} ⊕ H_{3} ⊕ ⋯ ⊕ H_{2k+1}, & d = 2k + 1, \end{matrix}  (21)

for all d ≥ 1.

Proof: The assertion for d = 0 is obvious since both B_{0} and H_{0} consist only of the constant functions. Since B_{ijk}^{d} are restrictions of trivariate polynomials of degree at most d to the sphere, it is clear that all of the functions B_{ijk}^{d} and B_{ijk}^{d−1} lie in H_{d}. We claim they are linearly independent. Suppose d is even and that

∑_{i+j+k=d} a_{ijk}B_{ijk}^{d}(v) + ∑_{i+j+k=d−1} b_{ijk}B_{ijk}^{d−1}(v) ≡ 0

for all v ∈ S. Then examining this expression at −v, and using the fact that B_{ijk}^{d}(−v) = B_{ijk}^{d}(v) while B_{ijk}^{d−1}(−v) = −B_{ijk}^{d−1}(v), we get

∑_{i+j+k=d} a_{ijk}B_{ijk}^{d}(v) − ∑_{i+j+k=d−1} b_{ijk}B_{ijk}^{d−1}(v) ≡ 0.

This implies that each sum is separately identically zero, and the desired linear independence follows. The case d odd is similar. Now since dim B_{d} + dim B_{d−1} = \binom{d+2}{2} + \binom{d+1}{2} = (d + 1)^{2}, we conclude that H_{d} splits into the two subspaces B_{d} and B_{d−1}. Then (21) follows from (4) and the fact that the functions in B_{d} are homogeneous of even (odd) degree when d is even (odd).

Combining this result with (4), it follows that

H_{d} = B_{d} ⊕ B_{d−1}  (22)

and thus

B_{d−1} ⊄ B_{d}  but  B_{d−2} ⊂ B_{d}.  (23)

In [2] an expression of the form

p = ∑_{i+j+k=d} c_{ijk}B_{ijk}^{d}  (24)

is called a spherical Bernstein-Bézier (SBB) polynomial. In view of their definition, it is no surprise that SBB-polynomials can be evaluated efficiently and stably with the usual de Casteljau algorithm. The same algorithm can also be used to perform subdivision, i.e., to find the coefficients of the three SBBpolynomials which represent p on the three subtriangles formed by splitting T about some point w ∈ T — see [2].

As in the planar case, it is also possible to perform degree raising on SBB polynomials, except that now the degree has to be raised by two rather than just one. In particular, if p is as in (24), then

p = ∑_{i+j+k=d+2} c̄_{ijk}B_{ijk}^{d+2},

where

\begin{aligned}c̄_{ijk} = & \frac{1}{(d + 1)(d + 2)}[i(i − 1)c_{i−2,j,k} + β_{110}ijc_{i−1,j−1,k} + j(j − 1)c_{i,j−2,k} \\ & + β_{101}ikc_{i−1,j,k−1} + k(k − 1)c_{i,j,k−2} + β_{011}jkc_{i,j−1,k−1}]. \end{aligned}

Here

β_{011} = \frac{sin^{2} t_{1}}{sin^{2} \frac{t_{1}}{2}} − 2, β_{101} = \frac{sin^{2} t_{2}}{sin^{2} \frac{t_{2}}{2}} − 2, and β_{110} = \frac{sin^{2} t_{3}}{sin^{2} \frac{t_{3}}{2}} − 2,

where t_{i} is the arc length of the edge opposite vertex v_{i}, i = 1, 2, 3.

The restriction of an SBB polynomial to an edge e of a spherical triangle results in a univariate function defined on the circular arc e. Such functions are called circular Bernstein-Bézier (CBB) polynomials [1], and in fact are trigonometric polynomials in arc length. They are also of interest for CAGD purposes, see e.g. [76,85,163].

5.3. Joining SBB Polynomials Smoothly

As a first step towards defining a space of splines on a spherical triangulation, one needs to describe how to make two SBB-polynomials defined on adjoining spherical triangles join together smoothly. Following [2], we now describe how this can be done.

In order to talk about smooth joins of SBB polynomials, we need to work with derivatives. First order directional derivatives of spherical functions were introduced above in Sect. 4.2. To generalize this to higher order (mixed) directional derivatives, it is more convenient to define directional derivatives in terms of the Cartesian coordinate system for IR^{3}. It was shown in [3] that if g is a unit tangent vector as in Sect. 4.2, then

D_{g}f(v) = g^{T} ∇F (v),

where F is a homogeneous extension of f to IR^{3} and ∇F is its gradient. While polynomials have a natural homogeneous extension, a general function f has many homogeneous extensions. However, as shown in [2], if a tangent direction g is chosen, then the value of D_{g}f(v) does not depend on the extension used.

For general functions f, higher order (order 2 or greater) directional derivatives are in general dependent on the way in which f is homogeneously extended, see [2]. However, for an SBB-polynomial, there is a convenient formula for its (unique) derivative. Given vectors g_{1}, …, g_{m} in IR^{3}, it was shown in [2] that the associated m-th order directional derivative is given by

D_{g_{1},...,g_{m}}p := D_{g_{1}} ⋯ D_{g_{m}}p = \frac{d!}{(d − m)!} ∑_{i+j+k=d−m} c_{ijk}^{m}B_{ijk}^{d−m}, (25)

where c_{ijk}^{m} are the coefficients obtained by applying the de Casteljau algorithm m times, starting with c_{ijk} and using the spherical barycentric coordinates b_{1}(g_{ν}), b_{2}(g_{ν}), b_{3}(g_{ν}) at the ν-th step.

Theorem 10. [2] Let T = ⟨v_{1}, v_{2}, v_{3}⟩ and T̃ = ⟨v_{4}, v_{2}, v_{3}⟩ be two spherical triangles sharing the edge ⟨v_{2}, v_{3}⟩, and let {B_{ijk}^{d}} and {B̃_{ijk}^{d}} be the associated Bernstein-Bézier basis functions. Suppose p and p̃ are SBB-polynomials on T and T̃ with coefficients {c_{ijk}} and {c̃_{ijk}}, respectively. Then p and p̃ and all of their directional derivatives up to order m agree on the edge shared by T and T̃ if and only if

c̃_{ijk} = ∑_{r+s+t=i} c_{r,j+s,k+t}B_{rst}^{i}(v_{4})  (26)

for all i = 0, …, m and all j, k such that i + j + k = d.

For later use, we note the following formulae [3] for first and second derivatives at the first vertex of T:

\begin{aligned}D_{g}p(v_{1}) & = d[b_{1}(g)c_{d,0,0} + b_{2}(g)c_{d−1,1,0}], \\ D_{g}^{2}p(v_{1}) & = d(d − 1)[b_{1}^{2}(g)c_{d,0,0} + 2b_{1}(g)b_{2}(g)c_{d−1,1,0} + b_{2}^{2}(g)c_{d−2,2,0}]. \end{aligned}

It is also important to note that the restrictions of derivatives of SBB polynomials to edges of T are CBB polynomials (see the end of Sect. 5.2). For example, if p is cubic and h is a direction which does not lie in the plane through the center of the sphere which contains the edge e := ⟨v_{1}, v_{2}⟩, then for any point v on the edge e, the cross-boundary derivative is the quadratic CBB-polynomial

D_{h}p(v) = 3[c_{200}^{1}b_{1}^{2}(v) + 2c_{110}^{1}b_{1}(v)b_{2}(v) + c_{020}^{1}b_{2}(v)^{2}],

where the c_{ijk}^{1} are computed using one step of the de Casteljau algorithm based on the spherical barycentric coordinates b_{1}(h), b_{2}(h), b_{3}(h) of h relative to T.

5.4. Spherical Splines

Given nonnegative integers r and d, the linear space

S_{d}^{r}(Δ) := {s ∈ C^{r}(S) : s|_{T_{i}} ∈ B_{d}, i = 1, …, N}

is called the space of spherical splines of smoothness r and degree d. Note that this space is defined in terms of B_{d} and not H_{d}. Thus, if s ∈ S_{d}^{r}(Δ), then its pieces are SBB-polynomials whose exact degree is even if d is even, and odd if d is odd, cf. (22)–(23).

The space S_{d}^{r}(Δ) is the analog of the space of polynomial splines defined on a planar triangulation. It seems particularly appropriate for use on the sphere since in view of Proposition 9, it consists of functions whose pieces are spherical harmonics joined together with global smoothness C^{r}, and thus has both the smoothness and high degree of flexibility which makes splines the powerful tools they are.

As in the planar case, it is possible to identify the dimension of S_{d}^{r}(Δ) and construct locally supported bases for them for all values of d ≥ 3r + 2, see [4]. We do not have space to review this theory here. However, we do discuss several scattered data interpolation and approximation methods based on spherical splines in the following sections.

It should be noted that the term spherical spline has several different meanings in the literature. It was used in [64,66,71] for certain radial basis type functions which are spherical analogs of the classical thin-plate splines – see also [197–200]. Unfortunately, recently the term has also been used for curves defined on the surface of the sphere [81].

5.5. Local Spline Interpolation Methods

As shown in [3], the basic interpolation problem (1) can now be solved as follows:

1) construct a triangulation Δ with vertices at the given data points,

2) for each triangle T in Δ, use the data at the vertices (along with additional derivative information if necessary) to define a function s_{T} defined on T which is a single SBB polynomial on T or a collection of SBB polynomial pieces on some partition of T. Choose s_{T} to interpolate the data at the vertices of T. With some care one can also make the function s whose restrictions are the s_{T} have some degree of global smoothness.

The fact that the interpolant s is constructed one triangle at a time insures that the method is local in the sense that the restriction of s to a triangle T depends only on the data in that triangle. Methods of this type are called macro-element methods. They have been widely applied in bivariate data fitting, and as observed in [3], every bivariate macro-element method has a natural spherical analog. Here we discuss just three examples:

1) quintic C^{1} macro-elements,

2) cubic C^{1} elements on the Clough-Tocher split,

3) quadratic C^{1} elements on the 6-triangle Powell-Sabin split.

Each of the methods discussed will solve a version of the Hermite interpolation problem (20) associated with a set of vertices {v_{i}}_{i=1}^{n}.

5.5.1. Quintic C^{1} Elements

Theorem 11. Suppose values for

s(v_{i}), D_{g_{i}}s(v_{i}), D_{h_{i}}s(v_{i}), D_{g_{i}}^{2}s(v_{i}), D_{g_{i}}D_{h_{i}}s(v_{i}), D_{h_{i}}^{2}s(v_{i})

are given for i = 1, …, n, where g_{i} and h_{i} are two linearly independent unit vectors which are tangent to S at v_{i}. In addition, suppose a value for a crossboundary derivative at the center of each edge of Δ is given. Then there exists a unique spherical spline s ∈ S_{5}^{1}(Δ) which interpolates these data.

This interpolant can be constructed one triangle at a time. Explicit formulae for the Bernstein-Bézier coefficients of each piece of s can be found in [3]. By construction, s satisfies C^{2} continuity conditions at each of the vertices. In this sense it is a spherical superspline, see [3]. It is the spherical analog of the classical Argyris element.

5.5.2. Clough-Tocher Element

For each triangle T = ⟨v_{1}, v_{2}, v_{3}⟩ in the triangulation Δ, let

v̄ := \frac{v_{1} + v_{2} + v_{3}}{∥v_{1} + v_{2} + v_{3}∥}

be its center. If v̄ is connected to each of the vertices of T with great-circle arcs, T is split into three spherical subtriangles. This is called the Clough-Tocher split of the triangle. Given a triangulation Δ, let Δ_{CT} be the triangulation obtained by applying the Clough-Tocher split to each triangle of Δ.

Theorem 12. Suppose function and first derivative information are given at the vertices v_{i} of the triangulation Δ as in Theorem 11, along with a value for a cross-boundary derivative at the center of each edge of Δ. Then there exists a unique spherical spline s ∈ S_{3}^{1}(Δ_{CT}) which interpolates these data.

This interpolant can also be constructed one triangle at a time. Explicit formulae for the Bernstein-Bézier coefficients of each piece of s can be found in [3].

5.5.3. Powell-Sabin Element

Given a triangulation Δ, for each j = 1, …, N, let v̄_{j} be the incenter of the j-th triangle obtained by radially projecting onto S the incenter of the planar triangle with the same vertices. Suppose the incenters of adjacent triangles are connected with great-circle arcs, and that the incenter of each triangle is also connected to each of its three vertices with great-circle arcs. This splits each triangle of Δ into six spherical subtriangles. The resulting triangulation Δ_{P S} is called the Powell-Sabin refinement of Δ.

Theorem 13. Suppose function and first derivative values are given at each of the vertices of a spherical triangulation Δ as in Theorem 11. Then there exists a unique spline s ∈ S_{2}^{1}(Δ_{P S}) which interpolates these data.

Explicit formulae for the Bernstein-Bézier coefficients of each piece of s can be found in [3].

5.5.4. A Hybrid Rational Element

Let T = ⟨v_{1}, v_{2}, v_{3}⟩ be a spherical triangle. Then a hybrid cubic SBB-element was defined in [97] to be a function of the form

r(v) = ∑_{i+j+k=3} c_{ijk}(v)B_{ijk}^{3}(v),  (27)

where

c_{111}(v) = ∑_{ℓ=1}^{3} α_{ℓ}A_{ℓ}(v),  (28)

and c_{ijk} are constants for the other choices of i, j, k. Here A_{1}, A_{2}, A_{3} are appropriate blending functions. For example, assuming the point v on S has spherical barycentric coordinates b_{1}, b_{2}, b_{3}, one can set

A_{1}(v) := \{ \begin{matrix}0, & v = v_{1}, v_{2}, v_{3}, \\ \frac{b_{2}^{m}b_{3}^{m}}{b_{1}^{m}b_{2}^{m} + b_{2}^{m}b_{3}^{m} + b_{1}^{m}b_{3}^{m}} , & otherwise, \end{matrix}  (29)

with A_{2}(v) and A_{3}(v) defined analogously. Using these elements, the Hermite interpolation problem (20) can be solved as follows:

1) for each triangle T, use the data at the vertices to determine all coefficients c_{ijk} with (i, j, k) ≠ (1, 1, 1).

2) compute α_{1} so that C^{1} continuity between adjoining pieces is guaranteed and ∑_{i=1}^{4} L_{i}^{2} is minimized, where L_{1}, …, L_{4} describe the C^{2} continuity conditions across the edge associated with α_{1}. Repeat to compute α_{2} and α_{3}.

As shown in [97], Step 2 involves solving 3 × 3 linear systems. The method is exact for cubic SBB-polynomials, i.e., if f is such a function, then its piecewise hybrid cubic interpolant s is identically equal to f. Two other methods for computing the α_{ℓ} are also discussed in [97]. The same idea can be used to create a C^{2} quintic hybrid element [18].

5.5.5. Estimating Derivatives on the Sphere

In order to apply the Hermite interpolation methods described above when only the basic data of Problem 1 is given, it is necessary to first estimate the derivatives at each data site v_{i}. The natural approach is to use a numerical differentiation rule based on data at points “near” v_{i}.

One approach is to use trivariate polynomials of low degree, see e.g. [88,153,154]. As a natural alternative, in [97] low degree SBB-polynomials are used instead. This involves choosing a spherical triangle surrounding the point of interest. The effect of varying the size and orientation of this triangle is explored there.

A method for estimating derivatives of functions defined on general surfaces (based on the use of the so-called exponential map) can be found in 5.6. Minimal Energy Interpolants

The idea of creating interpolants which minimize some form of energy has been heavily studied (it has been applied to standard univariate splines, tensorproduct splines, thin-plate splines, and a variety of other situations). The approach can also be carried out on the sphere. As usual, it has the advantage of producing interpolants which are often smoother than those obtained by other methods, but at a higher computational cost due to the global nature of the process.

Following [3], the starting point is the space of spherical splines S_{d}^{0}(Δ). Each spline in this space is uniquely defined by the set of coefficients of its SBB-pieces. By continuity, common coefficients along edges of the spherical triangulation can be identified. This leads to a single coefficient vector ccccccccc whose length is the dimension of S_{d}^{0}(Δ). Assuming that the first n coefficients c_{1}, …, c_{n} correspond to the values of s at the vertices v_{1}, …, v_{n}, it is clear that s will satisfy the interpolation conditions (1) provided c_{i} = r_{i} for i = 1, …, n.

Let Q be a symmetric positive definite matrix, and let E(ccccccccc) := ccccccccc^{T} Qccccccccc. Given r, suppose that

Gccccccccc = 0  (30)

describe the smoothness conditions required for s ∈ S_{d}^{0}(Δ) to belong to C^{r}. Then a minimal energy interpolating spline s is defined to be the function in S_{d}^{r}(Δ) which minimizes E(ccccccccc) subject to (30).

The coefficients ccccccccc of the minimal energy interpolating spline can be computed by solving the linear system

( \begin{matrix}Q & I & G^{T} \\ I & 0 & 0 \\ G & 0 & 0 \end{matrix} ) ( \begin{matrix}c & c & c & c & c & c & c & c & c \\ λ & λ & λ & λ & λ & λ & λ & λ & λ \\ γ & γ & γ & γ & γ & γ & γ & γ & γ \end{matrix} ) = ( \begin{matrix}0 \\ r & r & r & r & r & r & r & r & r \\ 0 \end{matrix} ) ,  (31)

where I is the n × n identity matrix and rrrrrrrrr := (r_{1}, …, r_{n})^{T}. Here γ γ γ γ γ γ γ γ γ, λ λ λ λ λ λ λ λ λ are vectors of Lagrange multipliers.

One way to define the energy functional is to take

E(f) := ∫_{S}(Of)^{2}ds,  (32)

where O is an appropriate differential operator on the sphere, such as

O := (Δ^{∗})^{m/2},  (33)

where m is an even integer and where Δ^{∗} is the Laplace-Beltrami operator (3). The definition of O for the case where m is odd is more complicated and can be found in [198]. The functionals (32) are the same functionals which are minimized in defining spherical thin plate splines, see [198–201].

5.7. Discrete Least Squares

In practice, the data fitting problem usually involves noisy measurements r_{i} = f(v_{i}) + ε_{i} of an unknown function f at n points v_{i} on the sphere. In this case it is better to approximate rather than interpolate. One approach is to create a least squares fit:

1) choose a triangulation Δ with m vertices (m < n),

2) choose a spline space S_{d}^{r}(Δ),

3) find the spline s ∈ S_{d}^{r}(Δ) which minimizes L(s) := ∑_{i=1}^{n}[s(v_{i}) − r_{i}]^{2}.

There does not seem to be a simple automatic way to choose the vertices for the triangulation. In practice it is tempting to choose S_{3}^{1}(Δ) for the spline space, despite the fact that it is not known whether interpolation is even possible with this space, see [3] for a discussion. As usual, the coefficients of s can be computed by solving a linear system.

5.8. Penalized Least Squares

In some fitting problems, particularly when the data are especially noisy, it may be useful to replace the standard discrete least squares problem by a penalized least squares problem. The idea is to minimize a combination

K(ccccccccc) := L(ccccccccc) + λE(ccccccccc),

where E(ccccccccc) is a measure of energy as discussed in Sect. 5.6, and L(ccccccccc) is the sum of squares of the errors as in the previous section. The parameter λ controls the trade-off between these two quantities, and is typically chosen to be a small positive number, see [74].

5.9. Remarks

We do not have space here to give examples of the various spherical spline fits discussed in Sects. 5.5 – 5.8 above. However, a number of numerical examples are discussed in [3] which contains both figures and tables giving a basic comparison of the methods as relates to storage, exactness, and computational time. The effect of thin triangles and the condition numbers of the interpolation matrices are also discussed there, along with appropriate scaling strategies. The effect of near-singular vertices (a vertex is singular when it is formed by two intersecting circular arcs) is also treated there.

In general, the choice of a method will depend on a variety of factors including 1) is the data noisy, 2) how much data is available, 3) what smoothness is required, and 4) what degree of accuracy is needed? Concerning accuracy, we mention that the experiments in [3,97] suggest that the quintic macro element method has accuracy O(h^{6}), the Clough-Tocher and hybrid rational methods have accuracy O(h^{4}), and the Powell-Sabin method has accuracy O(h^{3}), where h is the maximal diameter of the triangles in Δ.

5.10. Other Methods

Virtually any method which works for polynomial splines based on planar triangulations has an analog for the sphere using spherical splines. To see what has been done in the planar case, see the surveys [55,59,63,174]. Here we mention several additional interesting ideas:

1) Data Dependent Triangulations. The essential observation here is that given a fixed set of vertices, there are many possible associated triangulations. Thus, in fitting a set of data using splines based on such triangulations, it may be possible to get much better approximations with certain triangulations rather than others. In particular, the classical Delaunay triangulations may be far from best, and in fact long thin triangles might be much more appropriate. In the classical polynomial spline case on planar triangulations, this idea has been explored in detail in [37–39,142,143]. Algorithmically, the procedure works as follows. One starts with some reasonable initial triangulation, for example the Delaunay triangulation. Then to construct a “best” triangulation, edges are swapped recursively guided by some measure of goodness of fit or smoothness. Since this procedure is essentially trying to solve a very large nonlinear optimization problem, in practice one cannot expect it to always yield a global best approximation. One interesting approach to avoiding getting stuck at local minima is to employ simulated annealing, see [176]. The entire process carries over immediately to spherical splines.

2) Knot Insertion and Deletion. The idea here is to take advantage of the local nature of splines. In particular, if a spline is not doing a good job of fitting in some part of the domain, then one can insert additional knots in that area. Conversely, in regions where the fit is already very good, one can remove knots. Algorithmically, one starts with some initial fit and performs both knot insertion and deletion recursively. For details in the case of classical polynomial splines on planar triangulations, see [92]. Again, the method carries over immediately to spherical splines.

3) Spherical Simplex Splines. In the eighties there was considerable interest in certain spaces of simplex splines. They arose from the process of trying to create a multivariate analog of the univariate polynomial B-splines. For an extensive survey of the theory up to 1983, see [27]. Recently, analogous spherical simplex splines have been introduced and studied in [124,136].

4) Natural Neighbor Methods. Natural neighbor coordinates were introduced by Sibson (see e.g. [179]). In [15] analogous coordinates were defined for S, and used to create a locally supported interpolant to scattered data on S. The resulting surface is C^{0} at the data sites, and C^{1} everywhere else.

§6. Methods Based on Radial Basis Functions

In recent years there has been a great deal of interest in radial basis functions (RBF’s) as a tool for interpolation and data fitting. What is usually referred to as a radial basis function is the composition of a univariate function with some sort of distance function. In Euclidean spaces, an RBF with center v_{i} ∈ IR^{d} is a function of the form

Φ̂_{i}(v) := φ̂(∥v − v_{i}∥),  v ∈ IR^{d},  (34)

where φ̂ : [0, ∞) → IR and ∥ · ∥ denotes Euclidean distance.

In order to interpolate or approximate given data, it is natural to take linear combinations of RBF’s of the form

s(v) = ∑_{i=1}^{N} c_{i} φ̂(∥v − v_{i}∥),  v ∈ IR^{d}.  (35)

The coefficients c_{i} of this radial basis function expansion are then determined by solving a system of linear equations. The centers are usually chosen to coincide with (some of) the data. This assumption considerably simplifies the analysis of the linear system (see Sect. 6.2 for more details in the spherical case).

An early example of a radial basis function interpolant is provided by Shepard’s method (see [174] for a discussion). Other early uses of RBF’s were in geological modeling [77-79] (as (reciprocal) multiquadrics), and in variational problems, where – in a special case – they are called thin plate splines, see [32,58,60,106,111] and references therein.

Further interest in RBF’s was inspired by the observation [59] that in the bivariate setting, they are among the most effective methods for interpolation and approximation of scattered data. Their use was also encouraged by the fact that radial basis methods are extremely easy to program in any number of variables. We do not attempt to cite the entire classical RBF literature here, but for some typical methods and further references, see [17,33,34,36,40,43,52,115,141,166,180,204].

In analogy with the Euclidean case, cf. (34), it is natural to call

Φ_{i}(v) := φ(d(v, v_{i})), v ∈ S,

a spherical radial basis function (SRBF) provided φ : [0, π] → IR and d(v, v_{i}) is the length of the (shorter) great circle arc connecting v and v_{i}. This is the geodesic distance between v and v_{i}. In the classical literature the function Φ_{i}(v) = φ(d(v, v_{i})) is called a zonal function with pole v_{i}.

Given any Euclidean RBF as in (34), there is a natural way to associate a spherical RBF with it. Indeed, since

∥v − w∥ = \sqrt{2 − 2 v · w} = 2 sin \frac{d(v, w)}{2},

for any v, w ∈ S, it follows that

Φ̂_{i}(v) = φ̂(∥v − v_{i}∥) = Φ_{i}(v) = φ(d(v, v_{i})),

with

φ(t) = φ̂(2 sin(t/2)), 0 ≤ t ≤ π.

See [42,96] for some related discussions.

Recently, basis functions of the form

Φ̃_{i}(v) = φ̃(v · v_{i})

have also been considered, see e.g. [73]. Since v · w = cos(d(v, w)), these can also be interpreted as SRBF’s associated with the function φ(t) = φ̃(cos(t)).

For some recent papers which survey certain aspects of spherical radial basis functions and point to open problems, see [19–21,71].

6.1. Using RBF’s in IR^{3} for the Sphere

One way to use radial basis functions to fit data given at points v_{1}, …, v_{n} on the sphere is to consider these data points to lie in IR^{3}. Then if s is any RBF fit to the data based on Euclidean distance between points, we can simply restrict its evaluation to S.

For example, to interpolate values r_{i} at the points v_{i}, we can solve the system

∑_{j=1}^{n} c_{j}φ̂(∥v_{i} − v_{j}∥) = r_{i}, i = 1, …, n.

There is no complete classification of those functions φ̂ for which this system is nonsingular for all choices of {v_{i}}_{i=1}^{n}. However, wide classes of functions φ̂ have been identified which guarantee that the associated matrix is positive definite, which in turn assures its nonsingularity, see e.g. [115]. In this case, φ̂ is called strictly positive definite.

There are many examples of strictly positive definite RBF’s in the literature, and any of these can be used to solve the interpolation problem (1) on the sphere. However, it is highly advantageous to work with locally supported functions since they lead to sparse linear systems. Wendland [204] found a class of radial basis functions which are smooth, locally supported, and strictly positive definite on IR^{d} for some d. They consist of a product of a truncated power function and a low degree polynomial. For example, one can take

\begin{aligned}φ̂_{h}(t) & = (\frac{h − t}{h})_{+}^{2} , \\ φ̂_{h}(t) & = (\frac{h − t}{h})_{+}^{4} (\frac{4t + h}{h}) , \\ φ̂_{h}(t) & = (\frac{h − t}{h})_{+}^{6} (\frac{35t^{2} + 18rh + 3h^{2}}{h^{2}} ) , \end{aligned}

where h is a (small) positive number. These functions are nonnegative for t ∈ [0, h], are zero for t > h, and belong to C^{0}, C^{2}, and C^{4}, respectively.

Moreover, they are strictly positive definite in IR^{3}. Similar functions with higher-order smoothness, or which are strictly positive definite on IR^{d} for d > 3 can also be constructed. By the remarks at the end of the previous subsection, these functions can be transformed to work directly with geodesic distance on the sphere. Thus, for example, the second function above becomes

φ_{h}(t) = (1 − \frac{2}{h} sin \frac{t}{2})_{+}^{4} (\frac{8}{h} sin \frac{t}{2} + 1) ,

where t now measures geodesic distance. The support of this function is [0, arcsin(h/2)]. In [205] some other locally supported RBF’s in Euclidean spaces are discussed.

6.2. Interpolation with Spherical RBF’s

As observed above, one way to get SRBF’s is to transform Euclidean RBF’s to work with geodesic distance on the sphere. Alternatively, one can study functions φ which work with geodesic distance directly. This approach has spawned a large body of literature which we now briefly review. We begin with interpolation.

For the sphere, the analog of the radial basis function expansion (35) is

s(v) = ∑_{i=1}^{N} c_{i} φ(d(v, v_{i})),  v ∈ S.  (36)

To solve the interpolation problem (1) with linear combinations of spherical RBF’s, it is necessary to find coefficients c_{1}, …, c_{n} such that

∑_{j=1}^{n} c_{j} φ(d(v_{i}, v_{j})) = r_{i}, i = 1, …, n.

As in the Euclidean case, there is no complete characterization of those functions φ for which this n × n linear system is nonsingular. As before, a convenient sufficient condition for nonsingularity is that the corresponding matrix A := (φ(d(v_{i}, v_{j}))) be positive definite.

Definition 14. A continuous function φ : [0, π] → IR is called positive definite of order n on S if

∑_{i=1}^{n} ∑_{j=1}^{n} φ(d(v_{i}, v_{j})) c_{i}c_{j} ≥ 0  (37)

for any n points v_{1}, …, v_{n} ∈ S, and any c = [c_{1}, …, c_{n}]^{T} ∈ IR^{n}. If this inequality holds strictly for any nontrivial c , φ is called strictly positive definite of order n. If φ is (strictly) positive definite of order n for any n, then it is called (strictly) positive definite.

The problem of identifying positive definite functions on the sphere was addressed already by Schoenberg [170]. In particular, he showed that any φ which is positive definite on the sphere has an expansion in terms of Legendre polynomials P_{k} (normalized such that P_{k}(1) = 1) of the form

φ(t) = ∑_{k=0}^{∞} a_{k}P_{k}(cos t),  (38)

with a_{k} ≥ 0, and a_{k} such that ∑ a_{k} converges.

In view of its connection to interpolation on the sphere, the question of which functions φ are strictly positive definite (with respect to geodesic distance on S) has received considerable attention recently. It was shown in [206] that a sufficient condition is that all of the a_{k} in (38) be positive. This condition was recently improved in [171] to allow finitely many zero coefficients; note the comments made by Askey in Math Reviews about a (removable) error in the proof.

As a complement to these sufficient conditions, it has been shown recently that for strict positive definiteness on the circle (and therefore also on the sphere) it is necessary that infinitely many coefficients with even subscripts as well as infinitely many with odd subscripts be positive (see [114] or the survey [20]).

The problem of completely characterizing functions which are strictly positive definite on the sphere remains open. However, C^{∞}-kernels on S × S which are strictly positive definite in a distributional sense were completely characterized in [119]. Also, strictly positive definite functions on the infinitedimensional sphere S^{∞} were completely characterized in [113]. The latter results are of interest since any function which is strictly positive definite on S^{∞} has the same property on S.

For the purpose of interpolation, strict positive definiteness of order n on S suffices. It was shown by [155] and also [146] that a sufficient condition for this is that at least the first (n + 1)/2 coefficients in (38) are positive. To solve an interpolation problem with a (globally supported) SRBF involves evaluating it n^{2} times in order to build the matrix in (37). This means that expansions with (n+1)/2 nonzero coefficients are not very practical for larger values of n, and it would be better to work with a φ which has a simple closed form representation. We conclude this section with such an example.

The function

φ(t) = \frac{1}{\sqrt{1 + γ^{2} − 2γ cos t}},  0 < γ < 1,  (39)

is called the spherical reciprocal multiquadric. This infinitely differentiable function generates the Legendre polynomials, and therefore has the series representation (see e.g. [189])

φ(t) = ∑_{k=0}^{∞} γ^{k}P_{k}(cos t),

which is known to converge for |γ| < 1. By the above, it follows that this φ is strictly positive definite on S.

The use of (39) for interpolation on S was suggested by Hardy and Göpfert [78]. Their motivation came from geophysics: this particular function, which they called the multiquadric biharmonic, arises naturally in the calculation of disturbing potentials of the earth.

6.3. Interpolation with Spherical Harmonic Reproduction

In order to insure that SRBF interpolants of the form (36) do a good job of fitting arbitrary smooth functions on the sphere S, it is important that when used to interpolate a spherical harmonic, they give an exact fit. However, in general, this is not the case, and even constants are not fit exactly. This problem can be addressed by working with functions of the form

s(v) = ∑_{i=1}^{n} c_{i}φ(d(v, v_{i})) + ∑_{k=0}^{m−1} ∑_{ℓ=1}^{2k+1} d_{k,ℓ}Y_{k,ℓ}(v),  v ∈ S,  (40)

where {Y_{k,ℓ}} are the spherical harmonic basis functions for H_{k} described in Sect. 2.

By viewing SRBF interpolation in appropriate reproducing kernel Hilbert spaces (see Sect. 6.4 below), it can be shown that in constructing s as in (40) to satisfy the interpolation conditions (1), it is natural to enforce the additional conditions

∑_{j=1}^{n} c_{j}Y_{k,ℓ}(v_{j}) = 0,  k = 0, …, m − 1, ℓ = 1, …, 2k + 1.  (41)

We now discuss conditions under which this system is nonsingular.

Definition 15. A continuous function φ : [0, π] → IR is called conditionally positive definite of order m on S if

∑_{i=1}^{n} ∑_{j=1}^{n} φ(d(v_{i}, v_{j})) c_{i}c_{j} ≥ 0  (42)

for any n points v_{1}, …, v_{n} ∈ S and any c = [c_{1}, …, c_{n}]^{T} ∈ IR^{n} satisfying (41). If the points v_{1}, …, v_{n} are distinct and (42) is a strict inequality for nonzero c, then φ is called strictly conditionally positive definite of order m on S.

Note that, unfortunately, the term order has a different meaning here than it did in Definition 14. The terminology here seems to be well-established, while its other use for positive definiteness stems from [155].

Theorem 16. Suppose the function φ is strictly conditionally positive definite of order m on S and that the points v_{1}, …, v_{n} contain an H_{m−1}-unisolvent subset. Then there is a unique function of the form (40) satisfying (1) and (41).

Proof: The result follows from the same argument as in the Euclidean case given in e.g. [20].

As pointed out in Sect. 2, there does not seem to be a satisfactory theory of H_{m}-unisolvent sets, but a special result can be found in [73]. A sufficient condition for strict conditional positive definiteness of order m is

Theorem 17. [93] Suppose that

φ(t) = ∑_{k=m}^{∞} a_{k}P_{k}(cos t),

where P_{k} is the Legendre polynomial of degree k, and a_{k} = \frac{2k+1}{b_{k}^{2}4π} for some

real numbers b_{k} > 0 with ∑_{k=m}^{∞} \frac{2k+1}{b_{k}^{2}} < ∞. Then φ is strictly conditionally positive definite of order m on S.

We note that in the statement of this result in [93] there is a restriction that b_{k} ≥ 1 which, according to one of the authors [95], is not needed. A similar result concerning strict conditional positive definiteness in the distributional sense is given in [41]. A complete characterization of strict conditional positive definiteness does not yet exist. Partial results containing various (independent) necessary and sufficient conditions for the case m = 1 can be found in [112–114].

We note that, as in the Euclidean case, strictly conditionally positive definite functions of order one can also be used for interpolation without constant precision (see [20]).

The best-known example of a strictly conditionally positive definite function is the spherical multiquadric (see e.g. [53,77,138])

φ(t) = \sqrt{1 + γ^{2} − 2γ cos t}, γ > 0,

where t measures geodesic distance on the sphere. It was shown in [48] that this function is strictly conditionally positive definite of order one. For γ = 1 it can be seen (using the transformation discussed at the beginning of Sect. 6) that the spherical multiquadric is the restriction to the sphere of the IR^{3}- multiquadric (with multiquadric parameter c = 0). Note that the resulting function is not differentiable at its center in this case. For other values of γ this connection with Euclidean RBF’s no longer holds.

6.4. Variational Interpretation and Spherical Thin Plate Splines

In IR^{2} it was observed by Duchon [32] that the function which minimizes a certain energy functional over all smooth functions satisfying the interpolation conditions (1) is a radial basis function which is the kernel of a corresponding reproducing kernel Hilbert space. Such functions are called thin plate splines, and have been extensively studied — see e.g. [35,58–60,111,117,197] and references therein. Analogous spaces of thin plate splines have also been developed on the sphere, see e.g. [64,66,198,199]. These splines were also called spherical splines in the first two papers, but we shall not use this terminology in order to avoid confusion with the completely different spherical splines discussed in Sect. 5.

A complete variational theory for classical RBF’s was developed by Madych and Nelson [108,109]. Several authors have recently extended this work to the sphere, see [41,71,73,82,93]. The first and fourth of these papers are based on the fact that spherical harmonics are fundamental solutions of the Laplace-Beltrami operator, whereas the other three rely on a certain sequence {b_{k}}_{k=m}^{∞} to define appropriate Hilbert spaces in which to formulate the variational problems. The discussion in [71] and [73,93] is quite similar, the main differences being that in [93] (and its follow-up paper [73]) sharper error bounds are obtained. In [73,93] the case of higher-dimensional spheres is also treated. On the other hand, in [71] the properties of the sequence {b_{k}}_{k=m}^{∞} are given more attention. We will mention some of the ideas of these papers in the next section. The discussion in [41] is the most general since it addresses general Riemannian manifolds.

The main result is that any (conditionally) positive definite function on S can be viewed as the reproducing kernel of an associated Hilbert space, and it can be shown that the interpolant with minimal Hilbert space norm is a spherical radial basis function expansion (with added spherical harmonic terms). Furthermore, this interpolant is unique, and using the variational framework it is possible to give error bounds (see the next section).

6.5. Error Bounds and Stability

Error bounds for interpolation of smooth functions by spherical radial basis functions have recently been announced in [73,82]. While the bounds are very similar, the techniques are quite different. We present the main result of [73] since it is phrased in notation which is similar to what we have used above. It gives error bounds covering conditionally positive definite functions of all orders, and is based on the general theory of Golomb and Weinberger [75].

First we need some additional notation. One can define the associated Hilbert space mentioned in the previous section as the space

X_{m} = {f ∈ C(S) : f = ∑_{k=0}^{∞} ∑_{ℓ=1}^{2k+1} c_{kℓ}Y_{k,ℓ} and |f|_{m} < ∞},

where the semi-norm is given by

|f|_{m} = ( ∑_{k=m}^{∞} b_{k}^{2} ∑_{ℓ=1}^{2k+1} c_{kℓ}^{2})^{1/2},

with b_{k} > 0 such that ∑_{k=m}^{∞} \frac{2k+1}{b_{k}^{2}} < ∞.

Theorem 18. Suppose V_{h} = {v_{1}^{h}, …, v_{n(h)}^{h}} ⊂ S contains an H_{m−1} unisolvent subset and that

max_{v∈S} min_{v_{i}^{h}∈V_{h}} d(v, v_{i}^{h}) ≤ h.

Given f ∈ X_{m}, let s_{h} be the unique interpolant to f at the points in V_{h} satisfying (41). This interpolant has the form

s_{h}(v) = ∑_{i=1}^{n(h)} c_{i}φ(d(v, v_{i}^{h})) + ∑_{k=0}^{m−1} ∑_{ℓ=1}^{2k+1} d_{kℓ}Y_{k,ℓ}(v), v ∈ S,

with φ as in Theorem 17. If φ ∈ C^{λ}[0, ρ] for some 0 < ρ < π/2, then there is a constant C > 0 such that

∥f − s_{h}∥_{∞} ≤ Ch^{λ}|f|_{m},

for h sufficiently small and all f ∈ X_{m}.

For strictly positive definite SRBF’s, analogous error bounds were given in [82]. As mentioned above, the techniques used there are different from those in [73], and have the advantage that they lead to an explicit constant in the error bound. We also mention that [41] contains error bounds on general compact Riemannian manifolds, as well as results for points lying on an equiangular grid on the sphere.

In using SRBF’s for interpolating data on the sphere, another issue which has to be addressed is the stability of the method. This amounts to investigating the condition number of the interpolation matrix, which in turn reduces to estimating the norm of its inverse. This has been done for a certain type of strictly positive definite functions φ in [121].

Theorem 19. Assume φ is strictly positive definite with a_{k} > 0 in (38). Let q ∈ (0, π] be the minimal data separation of the data sites v_{1}, …, v_{n}, and let ν be the smallest integer such that

ν ≥ \frac{24}{1 − π^{2}/80} \frac{log(1/t(q))}{q^{2}},

where

t(q) = \frac{2q}{π^{3}[max\{8(π^{2}/q), 25π/2\}]}.

Then

∥A^{−1}∥_{2} ≤ \frac{C}{min\{a_{0}, min_{1≤k≤⌈ν/2⌉−1} \frac{a_{k}}{\sqrt{k}} \}},

where C is a constant. In particular, if a_{k} ∼ 1/(1 + k^{α}), α > 1, then

∥A^{−1}∥_{2} = O [(\frac{log(1/q^{2})}{q^{2}})^{α+1/2}] , q → 0.

The results in this section show that the error becomes smaller as the data sites move closer together, but at the same time the condition number becomes poorer. This is the “trade-off principle” well known from Euclidean radial basis function theory (see [167]).

Fundamental sets of continuous functions on the sphere were discussed in [183,184].

6.6. Locally Supported SRBF’s

In this section we briefly discuss some locally supported SRBF’s which were constructed directly for the sphere (see [14,71,172]) in contrast to those of Wendland (mentioned in Sect. 6.1) which were originally built for Euclidean space. Given h > 0, let

φ̃_{h}^{(k)}(t) = \{ \begin{matrix}(\frac{t − h}{1 − h})^{k} , & for t ≥ h, \\ 0, & otherwise, \end{matrix}  (43)

with t = v · w. Although these functions can be transformed to Wendland’s C^{0} function by 1 − t → t and 1 − h → h, it is not clear whether they are strictly positive definite on S. It is more instructive (cf. [172]) to consider

φ_{h}^{(k)}(t) = \{ \begin{matrix}(\frac{cos t − h}{1 − h})^{k} , & for t ≥ arccos h, \\ 0, & otherwise, \end{matrix}

where now t = d(v, w) = arccos(v·w) is the geodesic distance. With h replaced by cos h, the cases k = 1, 2 were already used in [185].

In [172] it was shown that the functions (φ_{h}^{(k)})^{(2)} obtained by spherical convolution of φ_{h}^{(k)} with itself are strictly positive definite on S. In [14] these functions were used to construct approximations via convolutions to an unknown function f assumed to be Lipschitz continuous on S. Since convolution has a smoothing effect, it is clear that in order to obtain a good approximation, the support of the basis functions should not be chosen too large. The following error estimate (see [14]) reflects this observation:

∥f − (φ̃_{h}^{(k)})^{(2)} ∗′ f∥_{∞} ≤ C\sqrt{1 − h},

where ∗′ is used to denote a discrete approximation to the spherical convolution, and (φ̃_{h}^{(k)})^{(2)} is as in (43). The local support of the basis functions is used to devise a hierarchical approximation method similar to the one described in Sect. 7.5. There seems to have been no attempt made to use these locally supported radial basis functions for interpolation purposes directly (although the possibility of doing so is mentioned in [172]). A major disadvantage of using the functions (φ_{h}^{(k)})^{(2)} is that they have to be computed numerically.

Approximation of functions on the sphere via singular integrals has also been studied, see e.g. [13].

Before locally supported basis functions were discovered, local interpolation schemes were constructed by interpolating only to subsets of the data and then blending these partial interpolants together. This idea was outlined in both [117] and [138] for radial basis functions on the sphere.

6.7. Generalized Hermite Interpolation

It is also possible to solve more general interpolation problems with spherical radial basis functions. Suppose we are given data sites v_{1}, …, v_{n}, as before, but that instead of simple function values at these points, we have data which are generated by a linearly independent set of linear functionals {L_{i}}_{i=1}^{n}. Thus, the problem is to find an interpolant s such that

L_{i}s = z_{i}, i = 1, …, n,

where z_{i} = L_{i}f for some unknown function f defined on S. If the functionals are point evaluation functionals, we have the standard interpolation problem (1). If we consider point evaluations of derivatives, then we have an Hermite interpolation problem. Other linear functionals such as local averages are also possible.

To solve this problem, the interpolant is assumed to be of the form

s(v) = ∑_{i=1}^{n} c_{i}L_{i}^{(2)}φ(d(v, v_{i})),  (44)

where the superscript (2) denotes that L_{i} is applied to φ ∘ d as a function of the second variable.

The first to study this problem on the sphere was Narcowich [119]. He showed that if all a_{k} > 0 in (38), then (for a general Riemannian manifold) the generalized Hermite interpolation problem has a unique solution. This class of functions includes e.g. the spherical reciprocal multiquadrics, and functions of the form φ(t) = e^{γ cos t}, γ > 0. In [46] other classes of basis functions were shown to lead to a unique solution for this problem, namely compositions of conditionally positive definite functions of order one with completely monotone functions, or with such functions whose derivative is completely monotone. These results ensure that the spherical multiquadrics can also be used for generalized Hermite interpolation. An expansion such as (44) seems to have been used first on the sphere by Freeden, see [65,67].

6.8. Discrete Least Squares Approximation

For RBF’s in the Euclidean case, discrete least squares approximation was studied in [180]. Essentially, the authors of those papers show that the associated Gramian matrix is nonsingular for some of the most popular radial basis functions (norm, multiquadrics, reciprocal multiquadrics and Gaussians) as long as the data sites are sufficiently well distributed and the centers for the radial basis functions are fairly evenly clustered about the data sites with the diameter of the clusters being relatively small compared to the separation distance of the data. Although there is no literature on least squares approximation with spherical RBF’s, similar results can be expected.

In [47] two algorithmic approaches to adaptive least squares fitting on S were compared. The first method is knot insertion, where starting with a coarse approximation (few centers), additional centers are added iteratively at the data location with the largest error component. This process is repeated until a certain tolerance is satisfied. If the initial knots are also chosen at data locations, then there are no problems with the collocation matrix becoming singular.

The second method is usually referred to as knot removal and proceeds in the reverse direction. It is most valuable if a sparse representation of the data is needed, e.g., for many subsequent evaluations. The algorithm starts with a very accurate initial fit and then tries to delete those centers (or basis functions) whose removal results in the smallest error. The procedure ends when a certain tolerance is reached.

The above algorithms can be made more efficient by dealing with several points at a time, and by optimizing the center locations with the help of nonlinear optimization methods. For the Euclidean setting, adaptive least squares fitting was discussed in [61,62].

§7. Multiresolution Methods

In recent years there has been a great deal of interest in multiresolution methods for compressing and approximating images, signals, and general functions and data. While an extensive theory exists for univariate and bivariate functions, there are only a few results for functions defined on the sphere S.

7.1. The Basic Multilevel Approach

Let

F_{0} ⊂ ⋯ ⊂ F_{k} ⊂ F  (45)

be a nested set of linear spaces, and suppose that for each i, P_{i} is a linear projection mapping F to F_{i}. In general, P_{i}s should provide an approximation in F_{i} to s. Typical choices include L_{2} best approximation and interpolation at appropriate points.

Given an s_{k} at the finest level, i.e., s_{k} ∈ F_{k}, it can be rewritten as s_{k} = s_{k−1} + g_{k−1}, where s_{k−1} := P_{k−1}s_{k} ∈ F_{k−1} and g_{k−1} := s_{k} − s_{k−1} ∈ F_{k}. Repeating this process recursively leads to the multilevel expansion

s_{k} = s_{0} + g_{0} + ⋯ + g_{k−1},  (46)

where

s_{i−1} := P_{i−1}s_{i}, g_{i−1} := s_{i} − s_{i−1},

for i = 1, …, k. The term s_{0} in the expansion lies in the coarsest space F_{0}, while the terms g_{i−1} ∈ F_{i} can be regarded as providing finer and finer levels of detail. Assuming that each of these terms is expanded in terms of appropriate basis functions, compression can be achieved by removing the associated coefficients which fall below some prescribed threshold.

Classical wavelet theory is based on choosing P_{i} to be the L_{2} approximation operator which produces a best approximation in F_{i}. The above decomposition process becomes especially simple if both the spaces F_{i} and their orthogonal complements G_{i} in F_{i+1} have simple bases. Often these are spanned by translates and dilates of a small number of functions. Examples of spaces with simple bases will be discussed in Sects. 7.2 and 7.3 below.

Another very useful way to choose the P_{i} is to take them to be interpolation operators, where P_{i}s is defined by the condition

P_{i}s(v_{j}^{i}) = s(v_{j}^{i}),  j = 1, …, n_{i},  (47)

where V_{i} := {v_{j}^{i}}_{j=1}^{n_{i}} are appropriate points in the domain of the functions in F_{i}, and where

V_{0} ⊂ ⋯ ⊂ V_{k}  (48)

is a prescribed nested sequence of sets of points. A method of this type defined on the sphere will be discussed in Sect. 7.4 below.

It should be pointed out that the above decomposition algorithm is closely related to a classical iterative procedure for creating increasingly accurate approximations to a given function s. Given a nested sequence of spaces as in (45), suppose that Q_{i} are linear projection operators with Q_{j}Q_{i} = Q_{i} for all j ≥ i. Suppose that Q_{i} maps a space F containing F_{k} (and thus all of the spaces in (45)) onto F_{i}. Then, given s ∈ F, Q_{i}s selects an approximation to s from the space F_{i}. Now, starting with a given function s ∈ F, first compute s̃_{0} := Q_{0}s ∈ F_{0}, and then set g̃_{0} := Q_{1}(s − s̃_{0}). This process can be repeated by setting

g̃_{i−1} := Q_{i}(s − s̃_{i−1}),  s̃_{i} := s̃_{i−1} + g̃_{i−1},

for i = 1, …, k, leading to the multilevel approximation

s ≈ s̃_{k} = s̃_{0} + g̃_{0} + ⋯ + g̃_{k−1},  (49)

where s̃_{0} ∈ F_{0} and g̃_{i−1} ∈ F_{i}. As before, s̃_{0} can be regarded as a coarse approximation, and the g̃_{i−1} as providing finer and finer levels of detail. The final s̃_{k} is in general not the best approximation of s from F_{k}, but has the advantage of the convenient multilevel expansion (49). If one starts with s ∈ F_{k}, then (49) becomes an exact multilevel expansion of s.

In general, the two algorithms described above are not equivalent. However, there is an important setting where they are.

Theorem 20. Suppose that a nested sequence of point sets as in (48) is given, and that P_{i} = Q_{i} is the operator which maps a function in F to the unique function in F_{i} which interpolates it at the points of V_{i} as in (47). Then the above two algorithms are equivalent in the sense that for any s_{k} ∈ F_{k}, the expansions (46) and (49) are exactly the same with s_{0} = s̃_{0} and g_{i} = g̃_{i} for all i = 0, …, k − 1.

Proof: Given s_{k} ∈ F_{k}, it is easy to see by induction that s_{i} is the unique function in F_{i} which interpolates to s_{k} on the point set V_{i}. Indeed, assuming this for i + 1, then

s_{i}(v_{j}^{i}) = P_{i}s_{i+1}(v_{j}^{i}) = s_{i+1}(v_{j}^{i}) = ⋯ = s_{k}(v_{j}^{i})

for all j = 1, …, n_{i}. We can prove the same thing for s̃_{i} by induction in the other direction. Indeed,

\begin{aligned}s̃_{i}(v_{j}^{i}) & = s̃_{i−1}(v_{j}^{i}) + Q_{i}(s_{k} − s̃_{i−1})(v_{j}^{i}) \\ & = \{ \begin{matrix}s̃_{i−1}(v_{j}^{i}), & v_{j}^{i} ∈ V_{i−1} \\ s̃_{i−1}(v_{j}^{i}) + s_{k}(v_{j}^{i}) − s̃_{i−1}(v_{j}^{i}), & v_{j}^{i} ∈ V_{i} \ V_{i−1} \end{matrix} = s_{k}(v_{j}^{i}). \end{aligned}

The assertion follows.

The basic multilevel idea has been used in a number of recent papers (see e.g. [14,25,49,52,120,173]). We discuss some of these in the following sections.

7.2. A Multiresolution Method Based on Tensor Splines

In this section we discuss a multilevel scheme for the sphere S based on L_{2} approximation using the tensor-product polynomial–trigonometric splines defined on the rectangle R defined in Sect. 3. Our discussion follows [103].

The process begins with a tensor polynomial-trigonometric spline s satisfying the conditions of Theorem 6. This insures tangent plane continuity of the associated surface defined on S. To simplify matters, equally spaced knots are used.

For each i = 0, …, k, let F_{i} be the space spanned by the normalized polynomial B-splines N_{1}^{3}, …, N_{m_{i}}^{3} of order 3 associated with the knots {x_{j}^{i}}_{j=1}^{m_{i}+3} where x_{1}^{i} = x_{2}^{i} = x_{3}^{i} = 0, x_{m_{i}+1}^{i} = x_{m_{i}+2}^{i} = x_{m_{i}+3}^{i} = π, and

x_{j+3}^{i} = jh_{i},  j = 1, …, m_{i} − 2,  (50)

with h_{i} = π/(m_{i} − 2) and m_{i} = 3 · 2^{i} + 2. Similarly, let F̃_{i} be the space of periodic trigonometric splines of order 3 associated with the knots

x̃_{j}^{i} = (j − 1)h̃_{i},  j = 1, …, m̃_{i},  (51)

where h̃_{i} = 2π/m̃_{i} and m̃_{i} = 3 · 2^{i}. Here the integer i is a nesting parameter, with the knots at level i being obtained from those at level i − 1 by inserting new knots at midpoints of knot intervals. These spline spaces satisfy F_{0} ⊂ ⋯ ⊂ F_{k} and F̃_{0} ⊂ ⋯ ⊂ F̃_{k}. Let

F_{i} = F_{i−1} ⊕ G_{i−1},  F̃_{i} = F̃_{i−1} ⊕ G̃_{i−1}.

Now the idea is to decompose s into a coarse spline and several terms containing details. Starting with s_{k} := s, the main result of [103] is an algorithm for rewriting s_{k} as the orthogonal decomposition

s_{k} = s_{k−1} + g_{k−1}^{(1)} + g_{k−1}^{(2)} + g_{k−1}^{(3)},

where

g_{k−1}^{(1)} ∈ F_{k−1} × G̃_{k−1}, g_{k−1}^{(2)} ∈ G_{k−1} × F̃_{k−1}, g_{k−1}^{(3)} ∈ G_{k−1} × G̃_{k−1}. (52)

The key to creating this decomposition is the fact that each of the spaces in this decomposition has a very simple basis in terms of polynomial and trigonometric B-splines and their associated spline wavelets. For the construction of trigonometric spline wavelets (and more general L-spline wavelets), see [102]. Then the decomposition algorithm (the process of finding the coefficients in the expansions of each of the functions in (52)) is a simple matter of matrix algebra.

The decomposition step (52) can be repeated until one ends up with a coarsest approximation s_{0}. Then standard thresholding methods can be applied to eliminate small coefficients in all expansions, thus achieving compression. Some care must be exercised in the thresholding step in order to maintain smoothness at the poles. In particular, at each level, coefficients of splines or wavelets which contribute to the values of s_{k}(θ, ϕ) or D_{θ}s_{k}(θ, ϕ) for θ = 0, π must be retained.

An associated reconstruction algorithm (which is again nothing but matrix multiplication) is also presented, along with various test examples. A similar method (based on exponential splines) was discussed in [26], see also [168].

Another way to avoid problems at the poles is to first split s into two parts s = f + g, where f is a homogeneous part satisfying f(θ, ϕ) = D_{ϕ}f(θ, ϕ) = 0 for θ = 0, π and all 0 ≤ ϕ ≤ 2π, see [102].

7.3. A Multiresolution Method Based on SRBF’s

In [122,123] wavelets were constructed using spherical radial basis functions. One starts with an infinite sequence of points v_{1}, v_{2}, … which is dense in S. Assuming that V_{i} consists of the first n_{i} of these points, it is clear that the V_{i} are nested. Given a strictly positive definite function φ on S as in Theorem 19, the associated n_{i} dimensional spaces

F_{i} := span{φ(d(v, v_{j}))}_{j=1}^{n_{i}}

are called the sampling spaces. These spaces are nested and their union is dense in L_{2}(S) (see [184]). In particular, any f ∈ F_{i−1} can be represented in F_{i} by simply padding its last n_{i} − n_{i−1} coordinates (with respect to the basis functions listed above) with zeros.

As in the classical wavelet theory, the wavelet spaces G_{i−1} are defined as the orthogonal complement of F_{i−1} in F_{i}, i.e.

F_{i} := F_{i−1} ⊕ G_{i−1}.

It turns out that there is a convenient way to represent functions in G_{i−1}. Let φ ⋆ φ be the SRBF defined via spherical convolution by

φ ⋆ φ(d(v, w)) = ∫_{S} φ(d(v, u))φ(d(u, w))dω(u),

where dω is the surface element for S. It is shown in [123] that φ ⋆ φ is strictly positive definite on S. Let A_{⋆} = (φ ⋆ φ(d(v_{j}, v_{k}))_{j,k=1}^{n_{i}} Then a basis for G_{i−1} is given by the last n_{i} − n_{i−1} columns of A_{⋆}^{−1}, i.e., by those columns corresponding to the centers v_{n_{i−1}+1}, …, v_{n_{i}}.

Decomposition and reconstruction formulae can also be found in [123]. They involve the solution of the interpolation problem involving A_{⋆}. Thus, it is desirable that the sampling spaces be spanned by locally supported functions. In [123] one can find several examples for φ along with explicit formulae for φ ⋆ φ. The locally supported functions of [172] (see also Sect. 6.6) and their iterated versions can also be used.

7.4. A Multiresolution Method Based on Spherical Splines S_{1}^{0}(Δ)

In this section and the next we describe multilevel methods for the sphere in which the projection operators P_{i} of Sect. 7.1 are chosen to be interpolation operators. The method to be described in this section makes use of C^{0} spherical splines of degree one.

Suppose Δ_{0} ⊂ Δ_{1} ⊂ ⋯ is a nested sequence of spherical triangulations which are obtained from a basic triangulation Δ_{0} by the following process: to get Δ_{i}, each triangle in Δ_{i−1} is split into four subtriangles by connecting the midpoints of its adjacent edges. Clearly, the associated spherical spline spaces S_{1}^{0}(Δ_{i}) are nested, i.e., S_{1}^{0}(Δ_{0}) ⊂ S_{1}^{0}(Δ_{1}) ⊂ · · ·. Let V_{i} be the set of vertices of Δ_{i}, and suppose P_{i} denotes the interpolating projector onto S_{1}^{0}(Δ_{i}).

The decomposition process begins with a spline s_{k} ∈ S_{1}^{0}(Δ_{k}) at the finest level k. Then s_{k−1} is the spline in S_{1}^{0}(Δ_{k−1}) which interpolates to s_{k} at the vertices of Δ_{k−1}. This process is almost trivial to implement since on each triangle (cf. Sect. 5), any spline s in S_{1}^{0}(Δ_{k−1}) can be expanded in terms of the first degree SBB-polynomials B_{100}^{1}, B_{010}^{1}, B_{001}^{1} associated with that triangle, and the coefficients are just the values of s at those vertices. Thus, the coefficients of s_{k−1} are obtained from those of s_{k} by simply retaining those associated with vertices of Δ_{k−1} and dropping those associated with midpoints of the edges of triangles in Δ_{k}.

Since g_{k−1} lies in S_{1}^{0}(Δ_{k}), its SBB-coefficients are also easily obtained. Indeed, by linear interpolation, if w is the midpoint of the great circle arc joining two vertices v_{i} and v_{j} of Δ_{k−1}, then the value of g_{k−1}(w) is simply

g_{k−1}(w) = s_{k}(w) − s_{k−1}(w) = s_{k}(w) − \frac{[s_{k}(v_{i}) + s_{k}(v_{j})]}{2}.

This is the key decomposition formula. The associated reconstruction formula is

s_{k}(w) = g_{k−1}(w) + \frac{[s_{k}(v_{i}) + s_{k}(v_{j})]}{2}.

As above, the decomposition process can be repeated until the coarsest level is reached. This gives the decomposition (46). At this point a thresholding process can be applied to remove all small coefficients, thus achieving compression.

This scheme was first presented in [173], along with several other related schemes with different choices of P_{i} (including certain lifting schemes which have additional desirable moment properties). Despite its simplicity, this method is quite effective in compressing data, e.g. see [173] for an example involving compression of a very large set (called ETOPO5) of elevation/bathymetric data taken over the surface of the earth by satellite. As pointed out in [25] (see also [169]) this idea was used already by Faber [45] in 1909. It can be extended to C^{1} splines on appropriate Powell-Sabin triangulations, see [25,169].

7.5. A Multilevel Interpolation Method Based on SRBF’s

In the context of radial basis functions in Euclidean spaces, the basic multilevel interpolation approach was introduced in [52]. It is also the subject of the recent papers [49,120]. Its application to the sphere is straightforward, and our discussion below focusses on this case.

Let v_{1}, …, v_{k} be a sequence of points on S, and let V_{i} consist of the first n_{i} (reasonably uniformly spaced) of these points. Then clearly the V_{i} are nested. Given a function s defined on the sphere, applying the method of Sect. 7.1 (cf. [52]) one starts at the coarsest level with the interpolant s_{0} = P_{0}s of the form

s_{0}(v) = ∑_{j=1}^{n_{0}} c_{j}^{0}φ_{h_{0}}(∥v − v_{j}^{0}∥).

At each step an additional amount of detail is computed via g_{i−1} = P_{i}(s − s_{i−1}), where

g_{i−1}(v) = ∑_{j=1}^{n_{i}} c_{j}^{i}φ_{h_{i}}(∥v − v_{j}^{i}∥).

Here h_{i} is a parameter which is used to scale the support size of the basis functions in accordance with the density of the data sites. The scaling should be done in such a way as to leave the bandwidth of the interpolation matrix (or the number of data sites inside the support of the basis functions) roughly constant throughout the process. This insures that the process will be wellconditioned.

Note that in this implementation the points v_{1}, …, v_{k} can be arbitrarily spaced. This is in contrast to the previous section where they were obtained as vertices of the triangles arising in successive refinement of a basic spherical triangulation. Numerically, one can observe that the algorithm has a rate of convergence which is at least linear in the “meshsize”, where meshsize is understood as in Theorem 18.

A similar algorithm for approximation by singular integrals was proposed in [14]. The authors of that paper obtained their level i approximation by a convolution of the data with an iterate of a locally supported SRBF as described in Sect. 6.6.

The motivation for [120] was to give rates of convergence for the algorithm suggested in [52]. However, the authors of [120] were able to give estimates only for a modified problem which employs convolutions of increasing multiplicity of the basis functions at the various levels, with the most smoothness required at the coarsest level. It is pointed out in [120] that for large problems, one arrives at smaller errors using this version of the multilevel algorithm as compared to solving the problem directly at the finest level with basis functions of relatively little smoothness (in case this is even computationally feasible).

7.6. An Approximate Newton Method

In [49] the decomposition step s_{i} = s_{i−1} + g_{i−1} in the basic multilevel method was interpreted as an instance of an approximate Newton method, i.e., s_{i} − s_{i−1} = g_{i−1}, where g_{i−1} is interpreted as an approximation to the inverse of the derivative of the mapping which computes the residual at level i. It was shown there (for Euclidean spaces) that an additional smoothing step at each level of the iteration can improve the convergence rate from linear to superlinear. With smoothing, the approximation spaces are no longer nested, and (49) is no longer an identity, although the right-hand side still provides an approximation to s. This interpretation of the algorithm is especially suited for the solution of differential equations.

§8. Additional Topics

In this final section we mention several other topics related to fitting scattered data on the sphere which we do not have space to discuss here in detail.

8.1. Sphere-like Surfaces

All of the above discussion has focused on the unit sphere S. As observed in [2,3], much of it is also valid for the following more general surfaces:

Definition 21. Given a smooth positive function ρ defined on the unit sphere S, the set S := {u ∈ IR^{3} : u = ρ(v)v, v ∈ S} is called a sphere-like surface.

When ρ ≡ 1, S reduces to the unit sphere. The earth is an example of a sphere-like surface, although the determination of ρ is no simple matter (and can be regarded as a scattered data fitting problem on the sphere), see [31].

8.2. General Surfaces

The sphere is only one example of an interesting surface for which the problem of fitting scattered data is of importance. The surface could as well be a torus, the wing of an aircraft, or any other physical 3D object. Some work has been done in the more general setting – see e.g. [9–12,56,119,137] and the references therein.

We point out that the simplest approach to fitting scattered data on a general surface — restricting the evaluation of a trivariate method to the surface — is not recommended, since problems are likely to arise as soon as this surface is “thin”, such as the wing of an aircraft (cf. [11]).

The problem of creating a Voronoi diagram associated with a set of points on a general surface (which then leads to an associated Delaunay triangulation) has been addressed in [86].

8.3. Visualizing Surfaces Defined on S

As observed above, the standard way (cf. [2,6,56,57,63,127]) to associate a surface with a function s defined on the sphere is via F := {s(v)v : v ∈ S}. Thus, s(v) can be identified with the distance from the center of the sphere to a point on the associated surface F if one moves out along the normal vector at v. Alternatively, s(v) can be taken as the distance above S, i.e. F := {(s(v) + 1)v : v ∈ S}.

One way to visualize F is to replace it by a surface consisting of planar facets. This can be accomplished by creating a spherical triangulation on S, computing the points s(v_{i})v_{i} ∈ IR^{3} associated with the vertices v_{i} of this triangulation, and then associating a planar facet with each spherical triangle as follows: given T with vertices v_{i}, v_{j}, v_{k}, let F_{T} be the triangle in IR^{3} with vertices at s(v_{i})v_{i}, s(v_{j})v_{j}, s(v_{k})v_{k}. The resulting faceted surface can then be displayed with standard graphical software. The capability to rotate the surface in 3D is very helpful. This method is easily generalized to general (convex) surfaces.

Another simple way to visualize function values defined on a general surface in IR^{3} is to triangulate the domain surface as above, and then assign a color mapped from some color scale to each vertex of the triangulation.

There are several other approaches to visualizing surfaces on surfaces including (color coded) contour regions, isophotes, and projections of 4Dgraphs. For details, see [57,110,139,140,156,182] and references therein.

8.4. Numerical Quadrature on the Sphere

The question of computing approximate values for integrals of functions defined on S is a classical subject in numerical analysis. The standard approach to creating spherical quadrature formulae can be described as follows. Suppose λ is a linear functional defined on C(S), and suppose P is an approximation process mapping C(S) to an approximating space A. Then given any f ∈ C(S), λf ≈ λP f. Such formulae are of particular interest when P f depends only on values of f at some (scattered) points v_{i} on S, and when λP f can be explicitly computed in terms of the values r_{i} = f(v_{i}).

One approach to developing such quadrature formulae is to use spherical harmonics for the approximating functions. If this is done locally on some spherical triangulation, it amounts to using the spherical splines defined in Sect. 5. Useful quadrature formulae have been derived in this way, although in contrast to the case of classical polynomials in Euclidean space, there are no exact formulae for integrals of spherical harmonics. Here we list only a few references [69,107,181].

References

1. Alfeld, P., M. Neamtu, and L. L. Schumaker, Circular Bernstein-Bézier polynomials, in Mathematical Methods for Curves and Surfaces, Morten Dæhlen, Tom Lyche, Larry L. Schumaker (eds), Vanderbilt University Press, Nashville & London, 1995, 11–20.

2. Alfeld, P., M. Neamtu, and L. L. Schumaker, Bernstein–Bézier polynomials on spheres and sphere–like surfaces, Comput. Aided Geom. Design 13 (1996), 333–349.

3. Alfeld, P., M. Neamtu, and L. L. Schumaker, Fitting scattered data on sphere-like surfaces using spherical splines, J. Comput. Appl. Math. 73 (1996), 5–43.

4. Alfeld, P., M. Neamtu, and L. L. Schumaker, Dimension and local bases of homogeneous spline spaces, SIAM J. Math. Anal. 27 (1996), 1482–1501.

5. Babar, C. S. and G. S. Pandey, On approximation to a function by Cesàro-Nörlund means of its ultraspherical series on a sphere, Vikram Math. J. 11 (1991), 29–38.

6. Bajaj, C. and G. Xu, Modeling and visualization of scattered function data on curved surfaces, in Fundamentals of Computer Graphics, J. Chen, N. Thalmann, Z. Tang, and D. Thalmann (eds.), World Scientific Publishing Co., 1994, 19–29.

7. Balmino, G., K. Lambeck, and W. M. Kaula, A spherical harmonic analysis of the earth’s topography, J. Geophys. Res. 78 (1973), 478–481.

8. Barnhill, R. E., G. Birkhoff, and W. J. Gordon, Smooth interpolation in triangles, J. Approx. Theory 8 (1973), 114–128.

9. Barnhill, R. E. and T. A. Foley, Methods for constructing surfaces on surfaces, in Geometric Modeling, Methods and Applications, H. Hagen and D. Roller (eds), Springer Verlag, Berlin, 1991, 1–15.

10. Barnhill, R. E., K. Opitz, and H. Pottmann, Fat surfaces: a trivariate approach to triangle-based interpolation on surfaces, Comput. Aided Geom. Design 9 (1992), 365–378.

11. Barnhill, R. E. and H. S. Ou, Surfaces defined on surfaces, Comput. Aided Geom. Design 7 (1990), 323–336.

12. Barnhill, R. E., B. R. Piper, and K. L. Rescorla, Interpolation to arbitrary data on a surface, in Geometric Modeling: Algorithms and New Trends, G. E. Farin (ed), SIAM Publications, Philadelphia, 1987, 281–290.

13. Berens, H., P. L. Butzer, and S. Pawelke, Limitierungsverfahren von Reihen mehrdimensionaler Kugelfunktionen und deren Saturationsverhalten, Publ. Res. Inst. Math. Sci. (Kyoto), Ser. A. 4 (1968), 201–268.

14. Brand, R., W. Freeden, and J. Fröhlich, An adaptive hierarchical approximation method on the sphere using axisymmetric locally supported basis functions, Computing 57 (1996), 187–212.

15. Brown, J. L., Natural neighbor interpolation on the sphere, in Wavelets, Images, and Surface Fitting, P.-J. Laurent, A. Le Méhauté, and L. L. Schumaker (eds), A. K. Peters, Wellesley MA, 1994, 67–74.

16. Brown, J. L. and A. J. Worsey, Problems with defining barycentric coordinates for the sphere, Math. Modelling Numer. Anal. 26 (1992), 37–49.

17. Buhmann, M., New developments in the theory of radial basis function interpolation, in Multivariate Approximation: From CAGD to Wavelets, Kurt Jetter and Florencio Utreras (eds), World Scientific Publishing, Singapore, 1993, 35–75.

18. Chang, L. H. T., Scattered Data Interpolation Schemes for Curves and Surfaces, dissertation, Univ. Sains Malaysia, 1997.

19. Cheney, E. W., Approximation and interpolation on spheres, in Approximation Theory, Wavelets and Applications, S. P. Singh (ed), Kluwer, Dordrecht, Netherlands, 1995, 47–53.

20. Cheney, E. W., Approximation using positive definite functions, in Approximation Theory VIII, Vol. 1: Approximation and Interpolation, Charles K. Chui and Larry L. Schumaker (eds), World Scientific Publishing Co., Inc., Singapore, 1995, 145–168.

21. Cheney, E. W. and X. Sun, Interpolation on spheres by positive definite functions, to appear in a volume dedicated to Professor Varma and edited by N. K. Govil, 1997.

22. Chui, C. K. and Hang-Chin Lai, Vandermonde determinants and Lagrange interpolation in IR^{s}, in Nonlinear and Convex Analysis, B. L. Lin and S. Simons (eds), Marcel Dekker, New York, 1987, 23–32.

23. Chung, K. C. and T. H. Yao, On lattices admitting unique Lagrange interpolations, SIAM J. Numer. Anal. 14 (1977), 735–743.

24. Courant, D. and D. Hilbert, Methods of Mathematical Physics, Vol. 1, Interscience, New York, 1953.

25. Dæhlen, M., T. Lyche, K. Mørken, R. Schneider, and H.-P. Seidel, Multiresolution analysis based on quadratic Hermite interpolation – Part 2: piecewise polynomial surfaces, preprint, 1997.

26. Dahlke, S., W. Dahmen, I. Weinreich, and E. Schmitt, Multiresolution analysis and wavelets on S^{2} and S^{3}, Numer. Func. Anal. Optim. 16 (1995), 19–41.

27. Dahmen, W. and C. A. Micchelli, Recent progress in multivariate splines, in Approximation Theory IV, C. Chui, L. Schumaker, and J. Ward (eds), Academic Press, New York, 1983, 27–121.

28. Dierckx, P., Algorithms for smoothing data on the sphere with tensor product splines, Computing 32 (1984), 319–342.

29. Dierckx, P., The spectral approximation of bicubic splines on the sphere, SIAM J. Sci. Statist. Comput. 7 (1986), 611–623.

30. Dierckx, P., Fast algorithms for smoothing data over a disk or a sphere using tensor product splines, in Algorithms for the Approximation of Functions and Data, J. C. Mason and M. G. Cox (eds), Oxford Univ. Press, Oxford, 1987, 51–65.

31. Dragomir, V. C. et al., Theory of the Earths Shape, Elsevier Sc. Publ., New York, 1982.

32. Duchon, J., Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces, RAIRO Anal. Numer. 10 (1976), 5–12.

33. Dyn, N., Interpolation of scattered data by radial functions, in Topics in Multivariate Approximation, C. K. Chui, L. L. Schumaker, and F. Utreras (eds), Academic Press, New York, 1987, 47–61.

34. Dyn, N., Interpolation and approximation by radial and related functions, in Approximation Theory VI, C. Chui, L. Schumaker, and J. Ward (eds), Academic Press, New York, 1989, 211–234.

35. Dyn, N., D. Levin, and S. Rippa, Surface interpolation and smoothing by ”thin plate” splines, in Approximation Theory IV, C. Chui, L. Schumaker, and J. Ward (eds), Academic Press, New York, 1983, 445–449.

36. Dyn, N., D. Levin, and S. Rippa, Numerical procedures for surface fitting of scattered data by radial functions, SIAM J. Sci. Statist. Comput. 7 (1986), 639–659.

37. Dyn, N., D. Levin, and S. Rippa, Algorithms for the construction of data dependent triangulations, in Algorithms for Approximation II, J. C. Mason and M. G. Cox (eds), Chapman & Hall, London, 1990, 185–192.

38. Dyn, N., D. Levin, and S. Rippa, Data dependent triangulations for piecewise linear interpolation, IMA J. Numer. Anal. 10 (1990), 137–154.

39. Dyn, N., D. Levin, and S. Rippa, Boundary corrections for data dependent triangulations, J. Comput. Appl. Math. 39 (1992), 179–192.

40. Dyn, N. and C. A. Micchelli, Interpolation by sums of radial functions, Numer. Math. 58 (1990), 1–9.

41. Dyn, N., F. Narcowich, and J. D. Ward, Variational principles and Sobolev-type estimates for generalized interpolation on a Riemannian manifold, Constr. Approx., to appear.

42. Dyn, N., F. Narcowich, and J. D. Ward, A framework for interpolation and approximation on Riemannian manifolds, preprint, 1997.

43. Dyn, N. and A. Ron, Radial basis function approximation: from gridded centers to scattered centers, Proc. London Math. Soc. 71(3) (1995), 76– 108.

44. Džafarov, A. S., The inverse problem of the theory of best approximation of functions on a sphere and on a segment (Russian), Akad. Nauk Azerbaĭdžan. SSR Dokl. 26 (1970), 3–6.

45. Faber, G., Über stetige Funktionen, Math. Ann. 66 (1909), 81–94.

46. Fasshauer, G. E., Hermite interpolation with radial basis functions on spheres, preprint, 1995.

47. Fasshauer, G. E., Adaptive least squares fitting with radial basis functions on the sphere, in Mathematical Methods for Curves and Surfaces, Morten Dæhlen, Tom Lyche, Larry L. Schumaker (eds), Vanderbilt University Press, Nashville & London, 1995, 141–150.

48. Fasshauer, G. E., Radial Basis Functions on Spheres, dissertation, Vanderbilt University, 1995.

49. Fasshauer, G. E. and J. W. Jerome, Multistep approximation algorithms: Improved convergence rates through postconditioning with smoothing kernels, preprint, 1997.

50. Fedorov, V. M., Approximation of functions on a sphere, Moscow Univ. Math. Bull. 45 (1990), 17–23.

51. Fedorov, V. M., Approximation of functions on a sphere. II, Moscow Univ. Math. Bull. 46 (1991), 30–36.

52. Floater, M. S. and A. Iske, Multistep scattered data interpolation using compactly supported radial basis functions, J. Comput. Applied Math. 73 (1996), 65–78.

53. Foley, T. A., Interpolation to scattered data on a spherical domain, in Algorithms for Approximation II, J. C. Mason and M. G. Cox (eds), Chapman & Hall, London, 1990, 303–310.

54. Foley, T. A., The map and blend scattered data interpolant on the sphere, Comput. Math. Appl. 24 (1992), 49–60.

55. Foley, T. A. and H. Hagen, Advances in scattered data interpolation, Surveys Math. Indust. 4 (1994), 71–84.

56. Foley, T. A., D. A. Lane, G. M. Nielson, R. Franke, and H. Hagen, Interpolation of scattered data on closed surfaces, Comput. Aided Geom. Design 7 (1990), 303–312.

57. Foley, T. A., D. A. Lane, G. M. Nielson, and R. Ramaraj, Visualizing functions over a sphere, IEEE Comp. Graphics & Appl. 10 (1990), 32– 40.

58. Franke, R., Smooth interpolation of scattered data by local thin plate splines, Comp. Maths. Appls. 8 (1982), 273–281.

59. Franke, R., Scattered data interpolation: tests of some methods, Math. Comp. 38 (1982), 181–200.

60. Franke, R., Scattered data interpolation using thin plate splines with tension, Comput. Aided Geom. Design 2 (1985), 87–95.

61. Franke, R., H. Hagen, and G. M. Nielson, Least squares surface approximation to scattered data using multiquadric functions, Advances in Comp. Math. 2 (1994), 81–99.

62. Franke, R., H. Hagen, and G. M. Nielson, Repeated knots in least squares multiquadric functions, in Geometric Modelling – Dagstuhl 1993, H. Hagen, G. Farin and H. Noltemeier (eds), Springer Verlag, Berlin, 1995, 177–187.

63. Franke, R. and G. M. Nielson, Scattered data interpolation and applications: A tutorial and survey, in Geometric Modeling, Methods and Applications, H. Hagen and D. Roller (eds), Springer Verlag, Berlin, 1991, 131–160.

64. Freeden, W., On spherical spline interpolation and approximation, Math. Meth. Appl. Sci. 3 (1981), 551–575.

65. Freeden, W., Spline methods in geodetic approximation problems, Math. Meth. Appl. Sci. 4 (1982), 382–396.

66. Freeden, W., Spherical spline interpolation – basic theory and computational aspects, J. Comput. Appl. Math. 11 (1984), 367–375.

67. Freeden, W., A spline interpolation method for solving boundary value problems of potential theory from discretely given data, Num. Meth. Part. Diff. Eq. 3 (1987), 375–398.

68. Freeden, W. and J. C. Mason, Uniform piecewise approximation on the sphere, in Algorithms for Approximation II, J. C. Mason and M. G. Cox (eds), Chapman & Hall, London, 1990, 320–333.

69. Freeden, W. and R. Reuter, Remainder terms in numerical integration formulas of the sphere, in Multivariate Approximation Theory II, W. Schempp and K. Zeller (eds), Birkhäuser, Basel, 1982, 151–170.

70. Freeden, W. and R. Reuter, Exact computation of spherical harmonics, Computing 32 (1984), 365–378.

71. Freeden, W., M. Schreiner, and R. Franke, A survey of spherical spline approximation, Surveys Math. Indust. 7 (1997), 29–85.

72. Gmelig-Meyling, R. H. J. and P. Pfluger, B-spline approximation of a closed surface, IMA J. Numer. Anal. 7 (1987), 73–96.

73. Golitschek, M. von and W. A. Light, Interpolation by polynomials and radial basis functions on spheres, preprint, 1997.

74. Golitschek, M. von and L. L. Schumaker, Data fitting by penalized least squares, in Algorithms for Approximation II, J. C. Mason and M. G. Cox (eds), Chapman & Hall, London, 1990, 210–227.

75. Golomb, M. and H. F. Weinberger, Optimal approximation and error bounds, in On Numerical Approximation, R. E. Langer (ed), University of Wisconsin Press, Madison, 1959, 117-190.

76. Goodman, T. N. T. and S. L. Lee, Interpolatory and variation-diminishing properties of generalized B-splines, Proc. Roy. Soc. Edinburgh Sect. A 96A (1984), 249–259.

77. Hardy, R. L., Theory and applications of the multiquadric-biharmonic method, Comput. Math. Appl. 19 (1990), 163–208.

78. Hardy, R. L. and W. M. Göpfert, Least squares prediction of gravity anomalies, geodial undulations, and deflections of the vertical with multiquadric harmonic functions, Geophys. Res. Letters 2 (1975), 423–426.

79. Hardy, R. L. and S. A. Nelson, A multiquadric-biharmonic representation and approximation of disturbing potential, Geophys. Res. Letters 13 (1986), 18–21.

80. Hobson, E. W., The Theory of Spherical and Ellipsoidal Harmonics, Chelsea, New York, 1955.

81. Hoschek, J. and G. Seemann, Spherical splines, Rev. Francaise Automat. Informat. Rech. Opér., Anal. Numer. 26 (1992), 1–22.

82. Jetter, K., J. Stöckler, and J. D. Ward, Error estimates for scattered data interpolation on spheres, preprint, 1997.

83. Kamzolov, A. I., Approximation of functions on the sphere S^{n} (Russian), Serdica 10 (1984), 3–10.

84. Kel’zon, A. A, Interpolation on the sphere (Russian), Izv. Vysš. Učebn. Zaved. Matematika 1975 (11), 41–46.

85. Koch, P. E., T. Lyche, M. Neamtu, and L. L. Schumaker, Control curves and knot insertion for trigonometric splines, Advances in Comp. Math. 3 (1995), 11–20.

86. Kunze, R., F.-E. Wolter, and T. Rausch, Geodesic Voronoi diagrams on parametric surfaces, Rpt. 2, Informatik, Univ. Hannover, 1997.

87. Lawson, C. L., Subroutines for C^{1} surface interpolation to data defined over the surface of a sphere, Rpt. 487, JPL, Cal. Tech, 1982.

88. Lawson, C. L., C^{1} surface interpolation for scattered data on a sphere, Rocky Mountain J. Math. 14 (1984), 177–202.

89. Lee, S. L., The use of homogeneous coordinates in spline functions and polynomial interpolation, preprint.

90. Lee, S. L. and G. M. Phillips, Interpolation on the simplex by homogeneous polynomials, in Numerical Mathematics, J. Wilson, R. P. Argarwal, and Y. M. Chow (eds), ISNM Vol. 86, Birkhäuser Verlag, Basel, 1988, 295–305.

91. Lee, S. L. and G. M. Phillips, Construction of lattices for Lagrange interpolation in projective space, Constr. Approx. 7 (1991), 283–297.

92. Le Méhauté, A. and Y. Lafranche, A knot removal strategy for scattered data in IR^{2}, in Mathematical Methods in Computer Aided Geometric Design, T. Lyche and L. Schumaker (eds), Academic Press, New York, 1989, 419–426.

93. Levesley, J., W. Light, D. Ragozin, and X. Sun, Variational theory for interpolation on spheres, preprint, 1997.

94. Li, L. Q. and H. Berens, The Peetre K-moduli and best approximation on the sphere (Chinese), Acta Math. Sinica 38 (1995), 589–599.

95. Light, W. A., Private communication, 1998.

96. Light, W. A. and E. W. Cheney, Interpolation by periodic radial basis functions, J. Math. Anal. Appl. 168 (1992), 111–130.

97. Liu, X.-Y. and L. L. Schumaker, Hybrid Bézier patches on sphere-like surfaces, J. Comput. Appl. Math. 73 (1996), 157–172.

98. Lizorkin, P. I., On the approximation of functions on the sphere σ. On the spaces B_{p,q}^{α}(σ). (Russian), Dokl. Akad. Nauk 331 (1993), 555–558.

99. Lizorkin, P. I. and S. M. Nikol’skiĭ, Approximation on the sphere in the metric of continuous functions (Russian), Dokl. Akad. Nauk SSSR 272 (1983), 524–528.

100. Lizorkin, P. I. and S. M. Nikol’skiĭ, A theorem concerning approximation on the sphere, Anal. Math. 9 (1983), 207–221.

101. Lizorkin, P. I. and Kh. P. Rustamov, Nikol’skiĭ-Besov spaces on a sphere that are associated with approximation theory, Proc. Steklov Inst. Math. 204 (1994), 149–172.

102. Lyche, T. and L. L. Schumaker, L–spline wavelets, in Wavelets: Theory, Algorithms, and Applications, C. Chui, L. Montefusco, and L. Puccio (eds), Academic Press, New York, 1994, 197–212.

103. Lyche, T. and L. L. Schumaker, A multiresolution tensor spline method for fitting functions on the sphere, preprint, 1997.

104. Lyche, T., L. L. Schumaker, and S. Stanley, Quasi-interpolants based on trigonometric splines, J. Approx. Theory, to appear.

105. MacRobert, T. M., Spherical Harmonics, Pergamon Press, Oxford, 1967.

106. Madych, W. R. and S. A. Nelson, Multivariate interpolation: a variational theory, manuscript, 1983.

107. Madych, W. R. and S. A. Nelson, Spherical quadrature and inversion of radon transforms, Proc. Amer. Math. Soc. 95 (1985), 453–457.

108. Madych, W. R. and S. A. Nelson, Multivariate interpolation and conditionally positive definite functions, Approx. Theory Appl. 4 (1988), 77–89.

109. Madych, W. R. and S. A. Nelson, Multivariate interpolation and conditionally positive definite functions, II, Math. Comp. 54 (1990), 211–230.

110. Max, N. L. and E. D. Getzoff, Spherical harmonic surfaces, IEEE Comp. Graph. Appl. 8 (1988), 42–50.

111. McMahon, John R. and Richard Franke, Knot selection for least squares thin plate splines, SIAM J. Sci. Statist. Comput. 13(2) (1992), 484–498.

112. Menegatto, V. A., Interpolation on spherical domains, Analysis 14 (1994), 415–424.

113. Menegatto, V. A., Strictly positive definite kernels on the Hilbert sphere, Appl. Anal. 55 (1994), 91–101.

114. Menegatto, V. A., Strictly positive definite kernels on the circle, Rocky Mountain J. Math. 25 (1995), 1149–1163.

115. Micchelli, C. A., Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constr. Approx. 2 (1986), 11– 22.

116. Möbius, A. F., Ueber eine neue Behandlungsweise der analytischen Sphärik, in Abhandlungen bei Begründung der Königl. Sächs. Gesellschaft der Wissenschaften, Jablonowski Gesellschaft, Leipzig, 45-86. (See also A. F. Möbius, Gesammelte Werke, F. Klein (ed.), vol. 2, Leipzig, 1886, 1–54.

117. Montès, P., Local kriging interpolation: Application to scattered data on the sphere, in Curves and Surfaces, P.-J. Laurent, A. Le Méhauté, and L. L. Schumaker (eds), Academic Press, New York, 1991, 325–329.

118. Müller, C., Spherical Harmonics, Springer Lecture Notes in Mathematics, Vol. 17, 1966.

119. Narcowich, F. J., Generalized Hermite interpolation and positive definite kernels on a Riemannian manifold, J. Math. Anal. Appl. 190 (1995), 165–193.

120. Narcowich, F. J., R. Schaback, and J. D. Ward, Multilevel interpolation and approximation, preprint, 1997.

121. Narcowich, F. J., N. Sivakumar, and J. D. Ward, Stability results for scattered data interpolation on Euclidean spheres, Advances in Comp. Math., to appear.

122. Narcowich, F. J. and J. D. Ward, Nonstationary spherical wavelets for scattered data, in Approximation Theory VIII, Vol. 2: Wavelets and Multilevel Approximation, C. Chui and L. Schumaker (eds), World Scientific Publishing, Singapore, 1995, 301–308.

123. Narcowich, F. J. and J. D. Ward, Nonstationary wavelets on the m-sphere for scattered data, Appl. Comput. Harmonic Anal. 3 (1996), 324–336.

124. Neamtu, M., Homogeneous simplex splines, J. Comput. Appl. Math. 73 (1996), 173–189.

125. Newman, D. J. and H. S. Shapiro, Jackson’s theorem in higher dimensions, in Über Approximationstheorie, P. L. Butzer, and J. Korevaar (ed), Birkhäuser, Basel, 1964, 208–219.

126. Nielson, G. M., The side-vertex method for interpolation in triangles, J. Approx. Theory 25 (1979), 318–336.

127. Nielson, G. M. and R. Ramaraj, Interpolation over a sphere based upon a minimum norm network, Comput. Aided Geom. Design 4 (1987), 41–58.

128. Nikol’skiĭ, S. M. and P. I. Lizorkin, Approximation theory on the sphere, Proc. Steklov Inst. Math. 172 (1985), 295–302.

129. Nikol’skiĭ, S. M. and P. I. Lizorkin, On the theory of approximations on a sphere (Russian), Trudy Mat. Inst. Steklov 172 (1985), 272–279, 355.

130. Nikol’skiĭ, S. M. and P. I. Lizorkin, Function spaces on a sphere that are connected with approximation theory (Russian), Mat. Zametki 41 (1987), 509–516, 620.

131. Nikol’skiĭ, S. M. and P. I. Lizorkin, Approximation of functions on the sphere, Math. USSR-Izv. 30 (1988), 599–614.

132. Nikolskiĭ, S. M. and P. I. Lizorkin, Approximation on the sphere – a survey, Banach Center Publ. 22 (1989), 281–292.

133. Opitz, K. and H. Pottmann, Computing shortest paths on polyhedra: applications in geometric modeling and scientific visualization, Int. J. Comput. Geom. Appl. 4 (1994), 165–178.

134. Perrin, F., J. Pernier, O. Bertrand, and J. F. Echallier, Spherical splines for scalp potential and current density mapping, Electroencephalography and clinical neurophsiology 72 (1989), 184–187.

135. Pevnyĭ, A. B., Spherical splines and interpolation on a sphere, Comput. Math. Math. Phys. 35 (1995), 109–112.

136. Pfeifle, R. and H.-P. Seidel, Spherical triangular B-splines with application to data fitting, Comp. Graphics Forum (Proc. Eurographics ’95) 14 (1995), C89–C96.

137. Pottmann, H., Interpolation on surfaces using minimum norm networks, Comput. Aided Geom. Design 9 (1992), 51–67.

138. Pottmann, H. and M. Eck, Modified multiquadric methods for scattered data interpolation over a sphere, Comput. Aided Geom. Design 7 (1990), 313–321.

139. Pottmann, H., H. Hagen, and A. Divivier, Visualizing functions on a surface, J. Visual. Comput. Anim. 2 (1991), 52–58.

140. Pottmann, H. and K. Opitz, Curvature analysis and visualization for functions defined on Euclidean spaces or surfaces, Comput. Aided Geom. Design 11 (1994), 655–674.

141. Powell, M. J. D., The theory of radial basis functions in 1990, in Advances in Numerical Analysis II: Wavelets, Subdivision, and Radial Basis Functions, W. Light (ed), Oxford University Press, Oxford, 1992, 105–210.

142. Quak, E. and L. L. Schumaker, Cubic spline fitting using data dependent triangulations, Comput. Aided Geom. Design 7 (1990), 293–301.

143. Quak, E. and L. L. Schumaker, Least squares fitting by linear splines on data dependent triangulations, in Curves and Surfaces, P.-J. Laurent, A. Le Méhauté, and L. L. Schumaker (eds), Academic Press, New York, 1991, 387–390.

144. Ragozin, D. L., Constructive polynomial approximation on spheres and projective spaces, Trans. Amer. Math. Soc. 162 (1971), 157–170.

145. Ragozin, D. L., Uniform convergence of spherical harmonic expansions, Math. Ann. 195 (1972), 87–94.

146. Ragozin, D. L. and J. Levesley, Zonal kernels, approximations and positive definiteness on spheres and compact homogeneous spaces, in Curves and Surfaces in Geometric Design, A. Le Méhauté, C. Rabut, and L. L. Schumaker (eds), Vanderbilt University Press, Nashville TN, 1997, 371– 147. Ramaraj, R., Interpolation and display of scattered data over a sphere,

M.S. Thesis, Arizona State University, 1986.

148. Reimer, M., Best approximations to polynomials in the mean and norms of coefficient-functionals, in Multivariate Approximation Theory, ISNM 51, W. Schempp and K. Zeller (eds), Birkhäuser Verlag, Basel, 1979, 289– 304.

149. Reimer, M., Interpolation on the sphere and bounds for the Lagrangian square sum, Resultate der Math. 11 (1987), 144–164.

150. Reimer, M., Interpolation mit sphärischen harmonischen Funktionen, in Numerical Methods in Approximation Theory Vol. 8, ISNM 81, L. Collatz, G. Meinardus and G. Nürnberger (eds), Birkhäuser, Basel, 1987, 184–187.

151. Reimer, M. and B. Sündermann, A Remez-type algorithm for the calculation of extremal fundamental systems for polynomial spaces over the sphere, Computing 37 (1986), 43–58.

152. Reimer, M. and B. Sündermann, Günstige Knoten für die Interpolation mit homogenen harmonischen Polynomen, Resultate der Math. 11 (1987), 254–266.

153. Renka, R. J., Interpolation of data on the surface of a sphere, ACM Trans. Math. Software 10 (1984), 417–436.

154. Renka, R. J., Algorithm 623: Interpolation on the surface of a sphere, ACM Trans. Math. Software 10 (1984), 437–439.

155. Ron, A. and X. Sun, Strictly positive definite functions on spheres in Euclidean spaces, Math. Comp. 65 (1996), 1513–1530.

156. Ruhlmann, G. M. and J. C. McKeeman, Local search: a new hidden line elimination algorithm to display spherical coordinate equations, Comput. & Graphics 15 (1991), 535–544.

157. Rustamov, Kh. P., Direct theorems on best L_{p}-approximation on a sphere S^{n−1} (Russian), Izv. Akad. Nauk Azerbaĭdzhan. SSR Ser. Fiz.-Tekhn. Mat. Nauk 5(5) (1984), 3–8.

158. Rustamov, Kh. P., Inverse theorems of best L_{p}-approximation on the sphere S^{n−1} (Russian), Izv. Akad. Nauk Azerbaĭdzhan. SSR Ser. Fiz.- Tekhn. Mat. Nauk 5(6) (1984), 6–12.

159. Rustamov, Kh. P., On direct and inverse theorems for the best L_{p}- approximation on a sphere (Russian), Dokl. Akad. Nauk SSSR 294 (1987), 788–791.

160. Rustamov, Kh. P., On the best approximation of functions on the sphere in the metric of L_{p}(S^{n}), 1 < p < ∞, Anal. Math. 17 (1991), 333–348.

161. Rustamov, Kh. P., Approximation of functions on a sphere (Russian), Dokl. Akad. Nauk SSSR 320 (1991), 1319–1325.

162. Rustamov, Kh. P., On the approximation of functions on a sphere (Russian), Izv. Ross. Akad. Nauk Ser. Mat. 57 (1993), 127–148.

163. Sánchez-Reyes, J., Single valued curves in polar coordinates, Comput. Aided Geom. Design 22 (1990), 19–26.

164. Sánchez-Reyes, J., Single-valued surfaces in spherical coordinates, Comput. Aided Geom. Design 11 (1994), 491–517.

165. Sansone, G., Orthogonal Functions, Interscience, New York (reprinted 1991 by Dover), 1959.

166. Schaback, R., Creating surfaces from scattered data using radial basis functions, in Mathematical Methods for Curves and Surfaces, Morten Dæhlen, Tom Lyche, Larry L. Schumaker (eds), Vanderbilt University Press, Nashville & London, 1995, 477–496.

167. Schaback, R., Error estimates and condition numbers for radial basis function interpolation, Advances in Comp. Math. 3 (1995), 251–264.

168. Schmitt, E., Wavelets and multiresolution analysis on sphere-like surfaces, SPIE Wavelet Applications in Signal and Image Processing, 1995, 92–101.

169. Schneider, R., Multiresolution analysis basierend auf beschränkten Projektionen mit Anwendungen in der Computergraphik, Diplomarbeit, Univ. Erlangen, 1997.

170. Schoenberg, I. J., Positive definite functions on spheres, Duke Math. J. 9 (1942), 96–108.

171. Schreiner, M., On a new condition for strictly positive definite functions on spheres, Proc. Amer. Math. Soc. 125 (1997), 531–539.

172. Schreiner, M., Locally supported kernels for spherical spline interpolation, J. Approx. Theory 89 (1997), 172–194.

173. Schroeder, P. and W. Sweldens, Spherical wavelets: Efficiently representing functions on the sphere, Computer Graphics Proceedings (SIG-GRAPH 95), 161–172.

174. Schumaker, L. L., Fitting surfaces to scattered data, in Approximation Theory, II, G. G. Lorentz, C. K. Chui, and L. L. Schumaker (eds), Academic Press, New York, 1976, 203–268.

175. Schumaker, L. L., Spline Functions: Basic Theory, Wiley, New York, 1981.

176. Schumaker, L. L., Computing optimal triangulations using simulated annealing, Comput. Aided Geom. Design 10 (1993), 329–345.

177. Schumaker, L. L and C. Traas, Fitting scattered data on spherelike surfaces using tensor products of trigonometric and polynomial splines, Numer. Math. 60 (1991), 133–144.

178. Shalaev, V. V., Sharp estimates for the approximation of functions, continuous on a sphere, by linear operators of convolution type, Ukrainian Math. J. 43 (1991), 523–525.

179. Sibson, R., A brief description of natural neighbor interpolation, in Interpreting Multivariate Data, D. V. Barnett (ed), Wiley, New York, 1981, 21–36.

180. Sivakumar, N. and J. D. Ward, On the least squares fit by radial functions to multidimensional scattered data, Numer. Math. 65 (1993), 219–243.

181. Sobolev, S. L., Cubature formulas on the sphere invariant under finite groups of rotations, Soviet Math. 3 (1962), 1307–1310.

182. Suffern, K. G., Perspective views of polar coordinate functions, Computer-Aided Design 24 (1992), 307–315.

183. Sun, X., The fundamentality of translates of a continuous function on spheres, Numerical Algorithms 8 (1994), 131–134.

184. Sun, X. and E. W. Cheney, Fundamental sets of continuous functions on spheres, Constr. Approx. 13 (1997), 245–250.

185. Svensson, S. L., Finite elements on the sphere, J. Approx. Theory 40 (1984), 246–260.

186. Swarztrauber, P. N., On the approximation of discrete scalar and vector functions on the sphere, SIAM J. Numer. Anal. 7 (1970), 934–949.

187. Swarztrauber, P. N., On the spectral approximation of discrete scalar and vector functions on the sphere, SIAM J. Numer. Anal. 16 (1979), 934–949.

188. Swarztrauber, P. N., The approximation of vector functions and their derivatives on the sphere, SIAM J. Numer. Anal. 18 (1981), 191–210.

189. Szegő, G., Orthogonal Polynomials, Amer. Math. Soc. Coll. Publ. Vol. XXIII, Providence, 1959.

190. Terekhin, A. P., Approximation in L_{p} by spherical polynomials, and classes of differentiable functions on a sphere (Russian), Trudy Mat. Inst. Steklov 187 (1989), 216–226.

191. Tissier, G., Approximation d’une fonction définie sur une sphère dans une base de fonctions sphèriques de Legendre (French), Rev. Francaise Automat. Informat. Rech. Opér., Anal. Numer. 5 (1971), 9–28.

192. Tissier, G., Interpolation à plusieurs variables sur la sphère (French), Numer. Math. 19 (1972), 136–145.

193. Traas, C. R., Smooth approximation of data on the sphere with splines, Computing 38 (1987), 177–184.

194. Ugulava, D. K., Approximation of functions on an m-dimensional sphere by the Cesàro means of the Fourier-Laplace series (Russian), Mat. Zametki 9 (1971), 343–353.

195. Ugulava, D. K., On the theory of the approximation of functions on a multidimensional sphere (Russian), Sakharth. SSR Mecn. Akad. Gamothvl. Centr. Šrom. 11 (1972), 139–156.

196. Ugulava, D. K., Approximation of continuous functions on a multidimensional sphere (Russian), Sakharth. SSR Mecn. Akad. Gamothvl. Centr. Šrom. 12 (1973), 36–52.

197. Wahba, G., Convergence rate of thin plate smoothing splines when the data are noisy, in Smoothing Techniques for Curve Estimation, T. G. M. Rosenblatt (ed), Springer Verlag, Berlin, 1979, 232–246.

198. Wahba, G., Spline Interpolation and smoothing on the sphere, SIAM J. Sci. Statist. Comput. 2 (1981), 5–16.

199. Wahba, G., Erratum: Spline interpolation and smoothing on the sphere, SIAM J. Sci. Statist. Comput. 3 (1982), 385–386.

200. Wahba, G., Vector splines on the sphere, with application to the estimation of vorticity and divergence from discrete, noisy data, in Multivariate Approximation Theory II, W. Schempp and K. Zeller (eds), Birkhäuser, Basel, 1982, 407–429.

201. Wahba, G., Surface fitting with scattered noisy data on Euclidean d− space and on the sphere, Rocky Mountain J. Math. 14 (1984), 281–299.

202. Wang, K. Y. and P. Zhang, Strong uniform approximation on sphere, Beijing Shifan Daxue Xuebao 30 (1994), 321–328.

203. Wehrens, M., Best approximation on the unit sphere in R^{k}, in Functional Analysis and Approximation, Birkhäuser, Basel, 1981, 233–245.

204. Wendland, H., Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in Comp. Math. 4 (1995), 389–396.

205. Wu, Z., Multivariate compactly supported positive definite radial functions, Advances in Comp. Math. 4 (1995), 283–292.

206. Xu, Y. and E. W. Cheney, Strictly positive definite functions on spheres, Proc. Amer. Math. Soc. 116 (1992), 977–981.

Acknowledgments. We thank Will Light, Fran Narcowich, Mike Neamtu, and Joe Ward for their helpful comments. The second author was partially supported by NSF Grant 9500643.

Greg Fasshauer  Larry L. Schumaker
CSAM Department  Dept. of Mathematics
Illinois Institute of Technology  Vanderbilt University
Chicago, IL 60616, USA  Nashville, TN 37240, USA

 fass@amadeus.csam.iit.edu       s@mars.cas.vanderbilt.edu