A Multiresolution Tensor Spline Method for Fitting Functions on the Sphere.
We present the details of a multiresolution method which we proposed at the Taormina Wavelet Conference in 1993 (see "L-spline wavelets" in  Wavelets: Theory, Algorithms, and Applications, C. Chui, L. Montefusco, and L. Puccio, eds., Academic Press, New York, pp. 197--212) which is suitable for fitting functions or data on the sphere. The method is based on tensor products of polynomial splines and trigonometric splines and can be adapted to produce surfaces that are either tangent plane continuous or almost tangent plane continuous. The result is a convenient compression algorithm for dealing with large amounts of data on the sphere. We give full details of a computer implementation that is highly efficient with respect to both storage and computational cost. We also demonstrate the performance of the method on several test examples.
Introduction
In many applications (e.g., in geophysics, meteorology, medical modelling, etc.), one
needs to construct smooth functions defined on the unit sphere S which approximate
or interpolate data. As shown in [11], one way to do this is to work with tensor-product
functions of the form
e
defined on the rectangle
where the ' i are quadratic polynomial B-splines on [\Gamma-=2; -=2], and the ~
are
periodic trigonometric splines of order three on [0; 2-]. With some care in the
choice of the coefficients (see Sect. 2), the associated surface
1) Institutt for Informatikk, University of Oslo P.O.Box 1080, Blindern 0316 Oslo, Norway
tom@ifi.uio.no. Supported by NATO Grant CRG951291. Part of the work was
completed during a stay at "Institut National des Sciences Appliqu'ees'' and "Laboratoire
Approximation et Optimisation" of "Universit'e Paul Sabatier'', Toulouse, France.
Department of Mathematics, Vanderbilt University, Nashville, TN 37240,
s@mars.cas.vanderbilt.edu. Supported by the National Science Foundation under
grant DMS-9803340 and by NATO Grant CRG951291.
with will be tangent plane continuous

In practice we often encounter very large data sets, and to get good fits using
tensor product splines (1.1), a large numbers of knots are required, resulting in many
basis functions and many coefficients. Since two spline spaces are nested if their
knot sequences are nested, one way to achieve a more efficient fit without sacrificing
quality is to look for a multiresolution representation of (1.1), i.e., to recursively
decompose it into splines on coarser meshes and corresponding correction (wavelet)
terms. Then compression can be achieved in the standard way by thresholding out
small coefficients.
The paper is organized as follows. In Sect. 2 we introduce notation and give
details on the tensor product splines to be used here. In Sect. 3 we describe the
general decomposition and reconstruction algorithm in matrix form, while in Sect. 4
we present a tensor version of the algorithms. The required matrices corresponding
to the polynomial and trigonometric spline spaces, respectively, are derived in
Sections 5 and 6. Sect. 7 is devoted to details of implementing the algorithm. In
Sect. 8 we present test examples, and in Sect. 9 several concluding remarks.
x2. Tangent plane continuous tensor splines
'm be the standard normalized quadratic B-splines associated with the
knot sequence
Recall that ' i is supported on the interval [x and that the B-splines form
a partition of unity on [\Gamma-=2; -=2]. Let T
m be the classical trigonometric
B-splines of order 3 defined on the knot sequence ~
m+3 , where
and ~
x ~
see Sect. 6. Recall that T j is supported on the
interval
~
m,
be the associated 2-periodic trigonometric B-splines, see [10]. These splines can
be normalized so that for OE 2 [0; 2-]
e
cos
x
e
cos
x
e
sin
x
~
Since the left and right boundaries of H map to the north and south poles,
respectively, a function f of the form (1.1) will be well-defined on S if and only if
x
and
where fS and fN are the values at the poles. Now since f is 2-periodic in the OE
variable and is C 1 continuous in both variables, we might expect that the corresponding
surface S f has a continuous tangent plane at nonpolar points. However,
since we are working in a parametric setting, more is needed. The following theorem
shows that under mild conditions on f which are normally satisfied in practice,
we do get tangent plane continuity except at the poles.
Theorem 2.1. Suppose f is a spline as in (1.1) which satisfies the conditions
and (2.5), and that in addition f('; OE) ? 0 for all ('; OE) 2 H. Then the
corresponding surface S f is tangent plane continuous at all nonpolar points of S.
Proof: Since f is a C 1 spline, the partial deriviatives f ' and f OE are continuous
on H. Now are two
tangents to the surface S f at the point f('; OE)vv v v v v v vv ('; OE). The normal vector to the
surface at this point is given by the cross product nn n n n n n nn := t 1 \Theta t 2 . By the hypotheses,
continuous, and thus to assure a continuous tangent plane, it suffices to show
that nn n n n n n nn has positive length (which insures that the surface does not have singular
points or cusps). Using Mathematica, it is easy to see that
\Theta cos(') 2 f(';

which is clearly positive for all values of ('; OE) 2 H with ' 6= \Sigma-=2.
With some additional side conditions on the coefficients of f , we can make the
surface S f also be tangent plane continuous at the poles. The required conditions
(cf. [3,11]) are that
AS cos
x
sin
and
AN cos
x
+BN sin
m, where AS ,BS ,AN , and BN are constants.
x3. Basic decomposition and reconstruction formulae
Suppose are a nested sequence of finite-dimensional linear subspaces of
an inner-product space X, i.e.
Let
be the corresponding orthogonal decompositions.
For our application, it is convenient to express decomposition and reconstruction
in matrix form, cf. [12]. Let ' k;mk be a basis for V k , and let
be a basis for W Then by the
nestedness, there exists an m k \Theta m k\Gamma1 matrix P k such that
where
The equation (3.1) is the usual refinement relation. Similarly, there exists an m k \Theta
such that
where
Let
be the Gram matrices of size m k \Theta m k and n k\Gamma1 \Theta n k\Gamma1 , respectively. It is easy to
see that
Clearly, the Gram matrices G k and H are symmetric. The linear independence
of the basis functions OE k;i and of / k;i implies that both G k and H are positive
definite, and thus nonsingular.
The following lemma shows how to decompose and reconstruct functions in V k
in terms of functions in V k\Gamma1 and W k\Gamma1 .
Lemma 3.1. Let f
k a k be a function in V k associated with a coefficient
vector a k 2 IR mk , and let
be its orthogonal decomposition, where
Then
a
Moreover,
a
Proof: To find a k\Gamma1 , we take the inner-product of both sides of (3.5) with ' k\Gamma1;i
. Using the refinement relation (3.1) and the orthogonality of
the ' k\Gamma1;i with / k\Gamma1;j , we get
which gives the formula for a k\Gamma1 . If we instead take the inner-products with /
we get the formula for b k\Gamma1 . In view of the linear independence of the functions
, the reconstruction formula (3.6) follows immediately from (3.5)
and the refinement relations.
x4. Tensor-product decomposition and reconstruction
In this section we discuss decomposition and reconstruction of functions in tensor
product spaces V k \Theta e
are as in the previous section, and where e
are similar subspaces of an inner-product space e
X. In particular, suppose
e
and that
e
be as in the previous section, and let e
be the analogous matrices associated with the spaces e
Theorem 4.1. Let f
k A k;' ~
' ' be a function in V k \Theta e
associated with a
coefficient matrix A k;' . Then f k;' has the orthogonal decomposition
with
f
~
(2)
~
where the matrices A
are computed from
the system of equations
G
e
G
e
e
e
e
G
e
e
e
with
e
Moreover,
e
e
e
e
Proof: To find the formula for A k\Gamma1;'\Gamma1 , we take the inner-product of both sides
of (4.1) with ' k\Gamma1;i for
' '\Gamma1;j for
. The
formulae for the B (i)
are obtained in a similar way. The reconstruction formula
follows directly from (4.1) after inserting the refinement relations and using
the linear independence of the components of the vectors ' k and in ~
Note that computing the matrices A k\Gamma1;'\Gamma1 and B (i)
in a decomposition
step can be done quite efficiently since several matrix products occur more than
once, and we need only solve linear systems of equations involving the four matrices
G
G
H . As we shall see below, in our setting the first two
of these are banded matrices, and the second two are periodic versions of banded
matrices. All of them can be precomputed and stored in compact form.
x5. The decomposition matrices for the polynomial splines
In this section we construct the matrices P k , Q k , and G k needed for the decomposition
and reconstruction of quadratic polynomial splines on the closed interval
[\Gamma-=2; -=2]. Consider the nested sequence of knots
where
with
be the associated normalized
quadratic B-splines with supports on the intervals [x k
. For
each k, the span V k of ' k;mk is the m k dimensional linear space of C 1
quadratic splines with knots at the x k
i . These spaces are clearly nested. In addition
to the well-known refinement relations
a simple computation shows that
Equations (5.2), (5.3) provide the entries for the matrix P k . In particular,
the first two and last two columns are determined by (5.3), while for any 3 -
2, the i-th column of P k contains all zeros except for the four rows
which contain the numbers 1=4, 3=4, 3=4, and 1=4. For example,
In general, P k has at most two nonzero entries in each row and and at most four
nonzero entries in each column.
In order to construct the matrices Q k , we now give a basis for the wavelet
space W k\Gamma1 . Here we work with the usual L 2 inner-product on L 2 [\Gamma-=2; -=2]. Let
Theorem 5.1. Given k - 1, let
and for k - 2, let
\Gamma6864\Gamma4967\Gamma4061
In addition, for k - 2, let
form a basis for W k\Gamma1 .
Proof: The wavelets in (5.5) are just the well-known quadratic spline wavelets,
see e.g., [1]. As described in [5], the coefficients of the remaining wavelets can
be computed by forcing orthogonality to V k\Gamma1 . In view of (3.2), the wavelets
are linearly independent if and only if the matrix Q k is
of full rank. This follows since the submatrix of Q k obtained by taking rows
easily be seen to be
diagonally dominant. For an alternate proof of linear independence, see Lemma 11
of [5].
In view of properties of B-splines, it is easy to see that
2:
We now describe the matrices Q k . By Theorem 5.1,
and
For general k - 2, the nonzero elements in the third column of Q k are repeated in
in each successive column they are shifted down
by two rows. The first two and last two columns of Q k contain the same nonzero
elements as Q 2 . Clearly, Q k has at most 4 nonzero entries in each row and at most
8 nonzero entries in each column.
We now describe the Gram matrices G k , which in general are symmetric and
five-banded. To get G k , we start with the matrix with 66h k =120 on the diagonal,
26h k =120 on the first subdiagonal, and h k =120 on the second diagonal. Then replace
the entries in the 3 \Theta 3 submatrices in the upper-left and lower-right corners by
@
For example,
and
x6. The decomposition matrices for the trigonometric splines
In this section we present the matrices e
needed for the decomposition
and reconstruction of periodic trigonometric splines of order 3. Suppose ' - 1, and
that
~
is a nested sequence of knots, where ~ h
where
0; otherwise.
is the usual trigonometric B-spline of order three associated with uniformly spaced
knots (0; h; 2h; 3h). Set
~
ae M ';i (OE);
For later use we define ~
';em '
' ';i for
The span e
';em '
is the space of periodic trigonometric splines
of order three. Clearly, these spaces are nested, and in fact we have the following
refinement relation:
Theorem 6.1. For all
~
where
Proof: By nestedness and the nature of the support of T h ,
for some numbers u; v; w; z. By symmetry, it is enough to compute u and v. To
find u, we note that on [0; h],
Then using (6.2) we can solve for u. To find v we note that
and then solve for v using (6.2).
Theorem 6.1 can now be used to find the entries in the matrix e
needed
in Sect. 2. In particular, each column has exactly the four nonzero elements
starting in the first row in column one, and shifted down
by two rows each time we move one column to the right (where in the last column
the last two elements are moved to the top of the column). For example,
e
Next we describe a basis for the wavelet space f
W '\Gamma1 which has dimension
~
1. In this case we work with the usual L 2 inner-product on
Theorem 6.2. Given ' - 1, let
~
~
where
and
~
with
Then ~
is a basis for the space f
Proof: To construct wavelets in f
we apply Theorem 5.1 of [6] which gives
explicit formulae for the ~
q i in terms of inner-products of ~
' ';i with ~
' '\Gamma1;j . To show
that ~
are linearly independent, it suffices to show that e
is of
full rank. To see this, we construct a ~ n '\Gamma1 \Theta ~
by moving the last
column of e
in front of the first column, and then selecting rows 2;
We
now show that this matrix is strictly diagonally dominant, and thus of full rank.
First, we note that in each row of B ' the element on the diagonal is ~
while the sum of the absolute values of the off diagonal elements is j~q 1 ( ~ h ' )j
simple computation shows that each of the functions D(h)
and r i (h) := ~
has a Taylor expansion which is an alternating series. In
particular, using the first two terms of each series, we get
Now it is easy to see that
\Theta j~q 3
also has an alternating series expansion, and we get
for the same range of h. This shows that B ' is strictly diagonally dominant, and
the proof is complete.
The formulae for the ~
q i in Theorem 6.2 are not appropriate for small values of
~
. In this case we can use the following Taylor expansions:
~
~
~
~
Rather than computing them each time we need them, we can precompute
and store the necessary values of ~
see

Table

1 in Sect. 7. We can now describe the matrix e
needed in Sect. 2 for
decomposing and reconstructing with trigonometric splines. For
e
~
~
~
~
~
~
where all ~
are evaluated at ~ h 1 . For ' - 2, each column of e
contains the 8 entries
~
evaluated at ~ h ' . In particular, these entries start in row
1 in column 1, and are shifted down by two each time we move one column to the
right (where in the last three columns, entries falling below the last row are moved
to the top). Clearly, e
Q ' has exactly four nonzero entries in each row. For example,
e
~
~
~
~
~
~
~
~
where all ~
are evaluated at ~ h 2 .
Finally, we describe the Gram matrices.
Theorem 6.3. For ' - 1, the 3
associated with the
~
' ';i is given by
e
I 00 I 01 I
I 01 I 00 I 01 I
I 02 I 01 I 00 I 01 I
I 01 I 00 I 01 I 02
I I 01 I 00 I 01
I 01 I
where
I
Z ~ x '
~
~
I 01 :=
Z ~ x '
~
~
I 00 :=
Z ~ x '
~
~
with
Moreover,
e
I
I
Proof: Using (6.2), the necessary integrals can be computed directly.
The formulae in Theorem 6.3 are clearly not appropriate for small values of
~
, in which case the following formulae can be used:
I

I
~

I

We can precompute and store the values of I 00 , I 01 , and I 02 for various levels
see

Table

2 in Sect. 7 for the values up to
x7. Implementation
7.1. Decomposition
The decomposition procedure begins with a tensor spline of the form (1.1) based on
polynomial splines ' k;i (') at a given level k - 1 and periodic trigonometric splines
~
';j (OE) at a given level ' - 1 with coefficient matrix C := A k;' of size m k \Theta e
To
carry out one step of the decomposition, we solve the systems (4.2) for A
, and set
To continue the decomposition, we now carry out the same procedure on the matrix
A k\Gamma1;'\Gamma1 . This process can be repeated at most min(k; times, where at each
step the new spline coefficients and wavelet coefficients are stored in C. Thus,
the entire decomposition process requires no additional storage beyond the original
coefficient matrix.
Because of the banded nature of the matrices appearing in (4.2), with careful
programming and the use of appropriate band matrix solvers, the j-th step of the
decomposition can be carried out with O(m k\Gammaj+1 e
To help keep
the number of operations as small as possible, we precompute and store the entries
of the matrices G
. appearing in (4.2). Table 1 gives
the values of ~
needed for the e

Table
gives the values of I
h ' and I needed for the e
G ' . The matrices
H k are symmetric positive definite and seven-banded, while the e
H ' are symmetric
positive definite periodic seven-banded matrices.
To check the robustness of the decomposition process, we computed the exact
condition numbers of the matrices G k , H k , e
H ' for up to eight levels. None
of the condition numbers exceeded 10, and we can conclude that the algorithm is
highly robust.
6 -28.996175484404513950 146.95303891951439472 -302.86111072242944246
9 -28.999940238853933503 146.99926613409589417 -302.99782947935722381
Tab. 1. Trigonometric spline wavelet coefficients for various '.
7.2. Thresholding
Typically, in the j-th step of the decomposition, many of the entries in the matrices
k\Gammaj;'\Gammaj of wavelet coefficients will be quite small. Thus, to achieve compression,
these can be removed by a thresholding process. In view of (5.6), tangent plane
continuity will be maintained at the poles if we retain all coefficients in the first two
and last two rows of these matrices. Given ffl, at the j-th level we remove all other
9 0.5500021215280431720 0.21666764608909212581 0.008333384794548569195
Tab. 2. Inner products of Trigonometric B-splines for various '.
wavelet coefficients in B (1)
k\Gammaj;'\Gammaj whose absolute values are smaller than
ffl=2 j . We do the same for B (3)
k\Gammaj;'\Gammaj using a threshold value of ffl=(300
smaller threshold is applied because of the scaling of the wavelets.
7.3. Reconstruction
In view of (4.3), to carry out one reconstruction step simply involves matrix multiplication
using our stored matrices. Because of the band nature of these matrices,
the computation of A k\Gammaj;'\Gammaj requires O(m k\Gammaj e
operations. At each step of
the reconstruction we can store these coefficients in the same matrix C where the
decomposition was carried out.
x8. Examples
To test the general performance of the algorithms, we begin with the following
simple example.
Example 1. Let s be the tensor spline with coefficients
Discussion: Since the normalized quadratic B-splines form a partition of unity,
it follows from (2.1) that with these coefficients, s j 1 for all ('; OE) 2 H, i.e., the
corresponding surface is exactly the unit sphere. In this case the coefficient matrix
is of size 770 \Theta 1536, and involves 1,182,720 coefficients. To test the algorithms, we
performed decomposition with various values of ffl, including zero. In all cases, after
reconstruction we got coefficients which were correct to machine accuracy (working
in double precision). The run time on a typical workstation is just a few seconds
for a full 7 levels of decomposition and reconstruction.
The illustrate the ability of our multiresolution approach to achieve high levels
of compression while retaining important features of a surface, we now create a
tensor spline fit to a smooth surface with a number of bumps.
Example 2. Let B be the surface shown in the upper left-hand corner of Figure 1.
Discussion: The surface B was created by fitting a spline f 8;8 to data created
by choosing 10 random sized subrectangles at random positions in H, and adding
tensor product quadratic B-splines of maximum height 3/4 with support on each
such rectangle to the constant values corresponding to the unit sphere. For
the coefficient matrix is of size 770 \Theta 768 and involves 591,360 coefficients.
To test the algorithms, we performed decomposition with the thresholding values
9. Table 3 shows the results of a typical run with
step nco
6 9746
Tab. 3. Reduction in coefficients in Example 2 with
Almost 3=4 of the coefficients are removed in the first step of decomposition, and
after 7 steps we end up with only 9745 coefficients (which amounts to a 60:1 compression
ratio). Table 4 shows the differences between the original coefficients and
the coefficients obtained after reconstruction. The table lists both the maximum
norm
and the average ' 1 norm
mem
are the original coefficients, and ~ c ij are the reconstructed ones. Due to the
scaling of the wavelets these numbers are somewhat larger than the corresponding
ffl.
The surfaces corresponding to the values are shown in

Figure

1. At near perfect looking reconstruction, while at
the major features are reproduced with only small wiggles in the surface. At
Fig. 1. Compressed surfaces for Example 2.
we have larger oscillations in the surface. This example shows that there is a critical
value of ffl beyond which the surface exhibits increasing oscillations with very little
additional compression.
x9. Remarks
Remark 9.1. The approach discussed in this paper was first presented at the
Taormina Wavelet Conference in October of 1993, and as far as we know was the
first spherical multiresolution method to be proposed. The corresponding proceedings
paper [6] focuses on the general theory of L-spline wavelets, and due to space
limitations, a full description of the method could not be included. In the meantime
Tab. 4. Coefficient errors in Example 2 for selected ffl.
we have become aware of the recent work [2,4,8,9,13]. In [2] the authors use tensor
splines based on exponential splines in the OE variable. The method in [4] uses
discretizations of certain continuous wavelet transforms based on singular integral
operators, while the method in [8] uses tensor functions based on polynomials and
trigonometric polynomials. Finally, the method in [9] utilizes C 0 piecewise linear
functions defined on spherical triangulations. Except for the last method, we are
not aware of implementations of the other methods.
Remark 9.2. In our original paper [6], an alternative way of making sure that
tangent plane continuity is maintained at the poles was proposed. The idea is to
decompose the original tensor product function s into two parts s H and s P , where
mk \Gamma2
e
and reconstruction can be
performed on s H . After adding s P , the reconstructed spline possesses tangent plane
continuity at the poles. Our implementation of this method exhibits essentially the
same performance in terms of compression and accuracy as the method described
here, but for higher compression ratios produces surfaces which are not quite as
visually pleasing near the poles.
Remark 9.3. The method described here can be extended to the case of nonuniform
knots in both the ' and OE variables. In this case the computational effort
increases considerably since the various matrices can no longer be precomputed
and stored.
Remark 9.4. In Sect. 4 we have presented the details of the tensor-product decomposition
and reconstruction algorithms assuming that the initial function f k;'
lies in the space V k \Phi e
not necessarily the same. Since these spaces
can always be reindexed, this is not strictly necessary in the abstract setting, but
was convenient for our application where there is a natural indexing for our spaces.
Remark 9.5. In computing the coefficients needed in Sections 5 and 6, we found
it convenient to use Mathematica.
Remark 9.6. There are several methods for computing approximations of the
form (1.1). An explicit quasi-interpolation method using data on a regular grid
(along with derivatives at the north and south poles) can be found in [11]. The
same paper also describes a two-stage method which can be used to interpolate
scattered data, and a least squares method which can be used to fit noisy data. A
general theory of quasi-interpolation operators based on trigonometric splines can
be found in [7].
Remark 9.7. A closed, bounded, connected set U in IR 3 which is topologically
equivalent to a sphere is called a sphere-like surface. This means that there exists a
one-to-one mapping of U onto the unit sphere S. Moreover, there exists a point O
inside the volume surrounded by U , such that every point on the surface U can be
seen from O. Such surfaces are also called starlike. For applications, we can focus
on the class of sphere-like surfaces of the form
where ae is a smooth function defined on S. Then each function f defined on U is
just the composition ae of a function g defined on S.
Remark 9.8. As indicated in [9], compression methods on the sphere can be
adapted to the problem of creating multiresolution representations of bidirectional
reflection distribution functions (BRDF's), althought the basic domain for such
functions is actually a hemisphere. We will explore the use of our method for this
purpose in a later paper.
Remark 9.9. It is well-known that the polynomial B-splines are stable. In particular
for quadratic B-splines (' i ) with general knots3 kck1 - k
for all coefficient vectors c. The same bounds hold for trigonometric splines since
the linear functionals
introduced in [11] are dual to the ~
Analogous stability results hold for general
p-norms.



--R

An Introduction to Wavelets
Multiresolution analysis and wavelets on S 2 and S 3
Algorithms for smoothing data on the sphere with tensor product splines
Spherical wavelet transform and its discretiza- tion

in Wavelets: Theory


efficiently representing functions on the sphere
Basic Theory
Fitting scattered data on spherelike surfaces using tensor products of trigonometric and polynomial splines
Wavelets for Computer Graphics
Biorthogonale Wavelets auf der Sph-are
--TR

--CTR
Thomas W. Sederberg , David L. Cardon , G. Thomas Finnigan , Nicholas S. North , Jianmin Zheng , Tom Lyche, T-spline simplification and local refinement, ACM Transactions on Graphics (TOG), v.23 n.3, August 2004
El Bachir Ameur , Driss Sbibih, Quadratic spline wavelets with arbitrary simple knots on the sphere, Journal of Computational and Applied Mathematics, v.162 n.1, p.273-286, 1 January 2004
John E. Lavery, Shape-preserving interpolation of irregular data by bivariate curvature-based cubic L
spherical data compression;tensor splines;multiresolution