Minimal Energy Surfaces using

Parametric Splines

by

Gregory E. Fasshauer and Larry L. Schumaker

Abstract. We explore the construction of parametric surfaces which interpolate prescribed 3D scattered data using spaces of parametric splines defined on a 2D triangulation. The method is based on minimizing certain natural energy expressions. Several examples involving filling holes and crowning surfaces are presented, and the role of the triangulation as a parameter is explored. The problem of creating closed surfaces is also addressed. This requires introducing spaces of splines on certain generalized triangulations.

1. Introduction

This paper deals with the problem of constructing smooth parametric surfaces which interpolate to given scattered data in IR^{3}. This problem has been heavily studied in the literature. There are two main approaches:

1) construct a wireframe of curves that interpolate to the data, and then fill in the faces with a series of individual parametric patches by choosing their parameters so that they join together smoothly; see e.g. [Coons ’64], [Gordon ’69], [Peters ’91a], [Sarraga ’87,’89,’90], [Shirman & Séquin ’87,’91] and references therein. This method is common in the CAGD world;

2) construct the surface as a patch complex defined parametrically in terms of coordinate functions on some convenient parameter domain, and choose their parameters to interpolate and to minimize some energy functional; see e.g. [Gmelig Meyling ’87], [Celniker & Gossard ’91] and references therein. This is the usual approach of Finite Element Theory.

Here we consider parametric patch complexes whose components are drawn from linear spaces of polynomial splines defined on a triangulation of a parameter domain in IR^{2}. This leads to surfaces which exhibit C^{1} smoothness as compared to the G^{1} continuity achieved by most patch methods. The method involves minimizing a natural quadratic energy functional, and allows the inclusion of a penalty term to achieve approximate C^{2} continuity of the interpolating surfaces. To illustrate the method, we present some numerical examples using cubic triangular patches for interpolating scattered data in IR^{3}, filling holes, and crowning surfaces. We also explore the possibility of adapting the method for closed surfaces.

The paper is organized as follows. In Section 2 we review the background material for surfaces defined by scalar-valued functions, and in Section 3 we describe the corresponding minimal energy approach. The definition of the linear space of smooth parametric splines, SP_{d}^{r}(△), and a discussion of existing parametric surface fitting methods is presented in Section 4, while Section 5 deals with the energy of parametric surfaces. Interpolation of scattered points in IR^{3} is the subject of Section 6, while in Section 7 we discuss filling n-sided holes and local editing of surfaces via crowning. In Section 8 we discuss the use of a C^{2} penalty term. Finally, the problem of interpolating with closed surfaces is explored in Section 9. This requires the introduction of generalized triangulations and associated spline spaces.

2. Polynomial Splines on Triangulations

We begin by recalling the definition of polynomial splines on triangulations of the plane IR^{2}. Let V := {(x_{i}, y_{i})}_{i=1}^{V} be a set of points in IR^{2}, and let Ω be its convex hull. A triangulation of Ω is defined as follows.

Definition 1. A collection △ = {T^{ℓ}}_{ℓ=1}^{N} of triangles in the plane is called a triangulation (of Ω) provided that

1) each triangle T^{ℓ} contains no points from V other than its vertices,

2) any two triangles either do not intersect, or share an edge,

3) the union Ω of the triangles {T^{ℓ}}_{ℓ=1}^{N} is a connected set in IR^{2}.

For later use we introduce the following notation:

V_{B} := number of boundary vertices,
V_{I} := number of interior vertices,
V := total number of vertices (data points),
E_{B} := number of boundary edges,
E_{I} := number of interior edges,
E := total number of edges,
N := number of triangles.

If we denote the \binom{d+2}{2} dimensional linear space of bivariate polynomials of total degree d by P_{d}^{(2)}, then given an arbitrary triangulation △, we can define an associated linear space of smooth piecewise polynomials.

Definition 2. Given two integers 0 ≤ r < d and a triangulation △, the space of polynomial splines of degree d and smoothness r is defined by

S_{d}^{r}(△) := {g ∈ C^{r}(Ω) : g|_{T^{ℓ}} ∈ P_{d}^{(2)}, ℓ = 1, …, N} .

Clearly S_{d}^{r}(△) is a linear space. For the applications we have in mind, it usually suffices to work with C^{1} continuity. It is known, see e.g. [Alfeld, Piper & Schumaker ’87] and references therein, that for d > 3

dim S_{d}^{1}(△) = n_{d} := D + σ,  (2.1)

where

D := \binom{d + 2}{2} + \binom{d}{2}E_{I} − \frac{d^{2} + 3d − 4}{2}V_{I},

and

σ := ∑_{i=1}^{V_{I}} ∑_{j=1}^{d−1}(j + 2 − je_{i})_{+}.

Here e_{i} is the number of edges with different slopes attached to the i-th vertex. It is a widely accepted conjecture that this formula also holds for d = 3. Exact formulae for dimensions are known for many (but not all) other choices of r and d, see e.g. [Alfeld & Schumaker ’90], [Schumaker ’91].

For numerical purposes it is useful to think of S_{d}^{1}(△) as the linear subspace of S_{d}^{0}(△) defined by

S_{d}^{1}(△) := S_{d}^{0}(△) ∩ C^{1}(Ω).

Moreover, it is convenient to represent the elements in S_{d}^{0}(△) in the well-known triangular Bernstein-Bézier form, cf. [Farin ’86], [de Boor ’87]. In particular, given g ∈ S_{d}^{0}(△), we suppose that for each 1 ≤ ℓ ≤ N, the polynomial g^{ℓ}(x, y) obtained by restricting g to triangle T^{ℓ} is written in the form

g^{ℓ}(r, s, t) := ∑_{i+j+k=d} c_{i,j,k}^{ℓ}B_{i,j,k}^{d}(r, s, t),

where

B_{i,j,k}^{d}(r, s, t) := \frac{d!}{i!j!k!}r^{i}s^{j}t^{k}, i + j + k = d,

are the classical bivariate Bernstein polynomials, and where (r, s, t) are the barycentric coordinates of the point (x, y) with respect to the vertices v_{1}^{ℓ}, v_{2}^{ℓ}, v_{3}^{ℓ} of T^{ℓ}. As usual, we associate the coefficients c_{i,j,k}^{ℓ} with the equally spaced points ξ_{i,j,k}^{ℓ} := (iv_{1}^{ℓ} + jv_{2}^{ℓ} + kv_{3}^{ℓ})/d. It will be convenient to introduce the notation

ccc_{i,j,k}^{ℓ} := (ξ_{i,j,k}^{ℓ}, c_{i,j,k}^{ℓ})  (2.2)

for the corresponding control points lying in IR^{3}. The collection of all such control points is called the control net or Bézier net. We say that c_{d,0,0}^{ℓ} is at the vertex v_{1}, while c_{d−1,1,0}^{ℓ} and c_{d−1,0,1}^{ℓ} are in the first ring around v_{1}. Similarly, we say that c_{d−j,j,0}^{ℓ}, j = 0, …, d, are on the edge \overline{v_{1}v_{2}}.

To force C^{0} continuity between adjoining patches, we identify the coefficients on the common edges between neighboring triangles. If we renumber the remaining coefficients, we get a vector ccc := [c_{1}, …, c_{n_{c}}]^{T} with

n_{c} := V + (d − 1)E + \binom{d − 1}{2}N.  (2.3)

Each spline in the space S_{d}^{0}(△) is uniquely defined by its coefficient vector ccc.

It is well-known (cf. [Farin ’86], [de Boor ’87]), that a spline in S_{d}^{0}(△) will belong to C^{1}(Ω) if and only if its coefficients satisfy a homogeneous linear system of equations of the form

Hccc = 00 0.  (2.4)

Here ccc is the vector of coefficients defined above, and H is a sparse matrix of size n_{cont} × n_{c}, where n_{cont} is the total number of continuity conditions which are imposed.

Two patches sharing a common edge as in Fig. 1 will join together with C^{1} continuity if and only if the following set of linear conditions are satisfied:

c_{i} = rb_{i} + sa_{i} + tb_{i+1}, i = 1, …, d,  (2.5)

where (r, s, t) are the barycentric coordinates of v_{4} with respect to the triangle with vertices v_{1}, v_{2}, v_{3}, and the a_{i}, b_{i}, and c_{i} are as shown in the figure. This implies that each row of the matrix H in (2.4) will contain at most four nonzero entries.

Fig. 1. A typical C^{1} continuity condition (d = 3).

The C^{1} continuity conditions (2.5) can also be expressed geometrically: they are satisfied if and only if

aa a_{i}, bbb_{i}, bbb_{i+1}, ccc_{i} are coplanar,  i = 1, …, d.  (2.6)

Since there are d conditions for each edge, to enforce C^{1} continuity, we have a total of

n_{cont} = dE_{I} = 3dV_{I} + dV_{B} − 3d  (2.7)

conditions. It is well-known that not all of these conditions are independent; there are a number of redundancies. In using (2.4) as side conditions for minimization problems, it is important to remove as many of the redundant rows of H as possible.

Lemma 3. The number of rows n_{r} in the matrix H corresponding to C^{1} continuity conditions which are redundant is

\{ \begin{matrix}n_{r} = 2V_{I} + σ, d ≥ 4, \\ n_{r} ≥ 2V_{I} + σ, d = 3. \end{matrix}

Proof: Using (2.1), for d ≥ 4 we know

\begin{aligned}n_{d} = dim S_{d}^{1}(△) & = \frac{(d + 1)(d + 2)}{2} + \frac{d(d − 1)}{2}E_{I} − \frac{(d^{2} + 3d − 4)}{2}V_{I} + σ \\ & = (d^{2} − 3d + 2)V_{I} + \frac{d(d − 1)}{2}V_{B} − (d^{2} − 3d − 1) + σ. \end{aligned}

The second expression follows immediately from the first and the well-known formulae

\begin{aligned}E_{B} & = V_{B}, \\ E_{I} & = 3V_{I} + V_{B} − 3, \\ N & = 2V_{I} + V_{B} − 2. \end{aligned}

Using these formulae on (2.3) also, we see that the number of coefficients in the vector ccc is

n_{c} = d^{2}V_{I} + \frac{d(d + 1)}{2}V_{B} − d^{2} + 1.

It then follows that the number of linearly independent C^{1} continuity conditions is

n_{ind} := n_{c} − n_{d} = (3d − 2)V_{I} + dV_{B} − 3d − σ.

The number of redundant equations is then given by

n_{r} = n_{cont} − n_{ind},

where n_{cont} is as in (2.7). For the case d = 3, formula (2.1) gives a lower bound for the dimension (see e.g. [Schumaker ’79]), and the proof proceeds as above.

As we mentioned above, it is not only important to know how many redundant rows the matrix H in (2.4) has, but where precisely these redundancies are located. The following lemmas identify exactly where the redundancies are to be found.

Lemma 4. Let v be an interior vertex with m edges attached to it as in Fig. 2, and let λ_{i}, i = 1, …, m, denote the associated C^{1} continuity conditions in the first layer around v. If the first m − 2 of these conditions are satisfied for some spline g in S_{d}^{0}(△), then the two remaining conditions associated with edges m − 1 and m are automatically satisfied.

Proof: Clearly, the C^{1} conditions defined by λ_{1}, …, λ_{m−2} hold if and only if all of the points aa a_{0}, …, aa a_{m} lie in one plane. But then the conditions associated with λ_{m−1} and λ_{m} automatically hold.

Fig. 2. C^{1} conditions at an interior vertex (m = 5, d = 3).

The following lemma shows that at a singular vertex, in addition to the two redundancies mentioned in Lemma 4, one of the continuity conditions in the second ring around v is always redundant (see also [Gmelig Meyling ’87]).

Lemma 5. Let v be a singular interior vertex v as in Fig. 3. Then if any set of three of the continuity conditions in the second ring around v is satisfied for a spline g in S_{d}^{0}(△), then the fourth second-ring condition is also automatically satisfied.

Proof: We show that one of the second-ring conditions depends on the other three second-ring conditions and two of the conditions in the first ring around v. If we label the coefficients as in Fig. 3, we can write these six continuity conditions

Fig. 3. C^{1} conditions at a singular vertex (d = 3).

involved in matrix form as

( \begin{matrix}t_{1} & r_{1} & s_{1} & 0 & −1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ t_{2} & −1 & r_{2} & s_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & t_{1} & 0 & 0 & 0 & s_{1} & 0 & 0 & −1 & r_{1} & 0 & 0 & 0 \\ 0 & 0 & t_{2} & 0 & 0 & −1 & s_{2} & 0 & 0 & 0 & r_{2} & 0 & 0 \\ 0 & 0 & 0 & t_{3} & 0 & 0 & −1 & s_{3} & 0 & 0 & 0 & r_{3} & 0 \\ 0 & 0 & 0 & 0 & t_{4} & 0 & 0 & −1 & s_{4} & 0 & 0 & 0 & r_{4} \end{matrix} ) ( a_{0} ; ⋮ ; a_{4} ; c_{1} ; ⋮ ; c_{4} ; b_{1} ; ⋮ ; b_{4} ) = 00 0,

where (r_{i}, s_{i}, t_{i}) are the barycentric coordinates of v_{i−1} with respect to the triangle with vertices v_{i}, v_{i+1}, and v (and the counting is cyclically modulo 4). ¿From the geometry of the singular vertex we know that all r_{i} are zero. Therefore, after dropping the last four columns and making some column and row interchanges, we obtain the following smaller matrix with the same rank as before:

( \begin{matrix}−1 & 0 & s_{1} & 0 & 0 & 0 & 0 & 0 & t_{1} \\ 0 & −1 & 0 & s_{2} & 0 & 0 & 0 & 0 & t_{2} \\ 0 & 0 & t_{2} & 0 & 0 & −1 & s_{2} & 0 & 0 \\ 0 & 0 & 0 & t_{3} & 0 & 0 & −1 & s_{3} & 0 \\ 0 & t_{1} & 0 & 0 & −1 & s_{1} & 0 & 0 & 0 \\ t_{4} & 0 & 0 & 0 & s_{4} & 0 & 0 & −1 & 0 \end{matrix} ) .

We show this matrix has rank 5. Performing Gaussian elimination, we obtain

( \begin{matrix}−1 & 0 & s_{1} & 0 & 0 & 0 & 0 & 0 & t_{1} \\ 0 & −1 & 0 & s_{2} & 0 & 0 & 0 & 0 & t_{2} \\ 0 & 0 & t_{2} & 0 & 0 & −1 & s_{2} & 0 & 0 \\ 0 & 0 & 0 & t_{3} & 0 & 0 & −1 & s_{3} & 0 \\ 0 & 0 & 0 & 0 & −1 & s_{1} & \frac{t_{1}t_{2}}{t_{3}} & \frac{s_{3}t_{1}t_{2}}{t_{3}} & t_{1}t_{2} \\ 0 & 0 & 0 & 0 & s_{4} & \frac{s_{1}t_{4}}{t_{2}} & \frac{−s_{1}s_{2}t_{4}}{t_{2}} & −1 & t_{1}t_{4} \end{matrix} ) .

Now multiplying the fifth row by s_{4} and adding it to the last row changes all of the entries in the last row to zero. To see this, we use the identities of Lemma 6 below and the fact that s_{i} + t_{i} = 1 (since r_{i} = 0) for i = 1, …, 4.

Lemma 6. At a singular vertex v, the barycentric coordinates satisfy the identities

s_{1}s_{3} = 1, s_{2}s_{4} = 1.

Proof: Given the configuration of Fig. 3, we first write v_{4} in terms of v_{1}, v_{2} and v to obtain

v_{4} = s_{1}v_{2} + t_{1}v.  (2.8)

Then we write v_{2} in terms of v_{3}, v_{4} and v to get

v_{2} = s_{3}v_{4} + t_{3}v.  (2.9)

Now substituting (2.9) into (2.8), we have

v_{4} = s_{1}s_{3}v_{4} + (s_{1}t_{3} + t_{1})v.

Since this equation must hold for arbitrary v_{4} and v, the first of the two identities follows. The second is obtained similarly.

3. Minimal Energy Interpolating Splines

Let S_{d}^{0}(△) be the space of C^{0} continuous splines of degree d defined on a triangulation △ of a domain Ω. Each spline in this space is uniquely defined by the vector ccc of coefficients as described in Sect. 2.

We want to use this space of splines to solve various interpolation problems. We shall be interested in sets of n_{int} interpolation conditions which can be written in the form

Bccc = h h h,  (3.1)

where B is a prescribed n_{int} × n_{c} matrix, and h h h = (h_{1}, …, h_{n_{int}}) is a prescribed vector. Each row of B corresponds to one interpolation condition. If a row of B corresponds to interpolation of function values and derivatives at specified points in Ω, then it contains at most \binom{d+2}{2} entries, see e.g. [Farin ’86]. (For example, the row of B corresponding to interpolation of a value at a vertex of △ has just one element in it, namely a 1 in the position corresponding to the coefficient associated with that vertex in the Bézier net).

To get an interpolating spline in S_{d}^{0}(△) which is C^{1}, we also have to enforce the smoothness conditions (2.4). In view of the discussion in Sect. 2, some of the equations in (2.4) are redundant. Removing the redundant equations, the C^{1} continuity conditions reduce to the homogeneous system

Accc = 00 0,  (3.2)

where A is now an n_{ind}×n_{c} matrix with n_{ind} = n_{cont}−n_{r}. Here n_{cont} is the number of continuity conditions (2.7), and n_{r} is the number of those that are redundant (cf. Lemma 3).

Depending on the interpolation conditions to be enforced, there may not exist any such spline. Thus for the rest of this section we assume that the following hypothesis is satisfied:

U := {g ∈ S_{d}^{0}(△) : Accc = 00 0 and Bccc = h h h} ≠ ∅.  (3.3)

While U itself is not a linear subspace, it is the translate of one, say U^{0}. Assuming (3.3) holds, let n_{f} be the dimension of U^{0}. Then n_{f} is the number of free parameters for elements in U. In order to get visually pleasing smooth surfaces, it has been suggested by several authors (e.g. [Schmidt ’82], [Alfeld ’84b], [Gmelig Meyling ’87]) to use these free parameters to minimize some kind of energy functional for the surface.

As can be found in classical books on mechanics (e.g. [Timoshenko ’59]), the energy of a thin plate of arbitrary shape under small deflections is given by

∫ ∫ [(f_{xx} + f_{yy})^{2} − 2(1 − ν)(f_{xx}f_{yy} − f_{xy}^{2})] dxdy,

where the parameter ν is a constant (Poisson’s ratio) depending on the material at hand. We choose ν = 0 to simplify the integrand as much as possible (ν = 0.3 might be a more realistic setting, since it corresponds to aluminum or steel). This means that on each triangle T^{ℓ}, we should measure the energy by

J_{ℓ} := ∫ ∫_{T_{ℓ}} [(f_{xx})^{2} + 2(f_{xy})^{2} + (f_{yy})^{2}] dxdy,

so that the total energy is given by

J := ∑_{ℓ=1}^{N} J_{ℓ}.  (3.4)

We can now define what we mean by a minimal energy spline interpolant.

Definition 7. A spline g ∈ U is called a minimal energy C^{1} interpolating spline provided

J(g) ≤ J(u) for all u ∈ U.

It turns out that the coefficients (Bézier net) of a minimal energy interpolating spline can be found by solving a sparse system of linear equations. To state this formally, we need some notation. As shown in [Quak & Schumaker ’90], the energy expression (3.4) can be written as a quadratic form

J = ccc^{T} Qccc,  (3.5)

where Q is an n_{c} × n_{c} symmetric positive semi-definite matrix. It was also shown in [Quak & Schumaker ’90] that for cubic splines, the matrix Q has 13 essential entries, while in the quartic case, there are 26. The formulae for these entries are explicitly listed there in terms of inner products of gradient vectors of the barycentric coordinates (r, s, t) of a point (x, y). Using these formulae, the energy matrix Q can be assembled without performing any integrations. This is a major advantage over any method which uses a more complicated energy functional (such as the methods mentioned in Section 5 below), or which uses numerical integration techniques to evaluate the energy integral.

Suppose A is the n_{ind} × n_{c} matrix of (3.2) describing the continuity conditions for a spline g ∈ S_{d}^{0}(△) to lie in S_{d}^{1}(△). Suppose B is an n_{int} × n_{c} matrix describing a set of n_{int} interpolation conditions as in (3.1). Then the following is well-known:

Theorem 8. Assume (3.3) holds. Then there is a unique solution to the following constrained quadratic minimization problem:

min ccc^{T} Qccc, subject to Accc = 00 0, Bccc = h h h.

It can be found by solving the linear system

( \begin{matrix}Q & C^{T} \\ C & 0 \end{matrix} ) ( \begin{matrix}c & c & c \\ λ & λ & λ \end{matrix}) = ( \begin{matrix}0 & 0 & 0 \\ b & b & b \end{matrix}) ,  (3.6)

where

C := ( \begin{matrix}A \\ B \end{matrix} ) , bbb := ( \begin{matrix}0 & 0 & 0 \\ h & h & h \end{matrix}) ,

and λ λ λ is an n_{ind} + n_{int} vector (whose components are called Lagrange multipliers). Moreover, the spline g with Bézier net ccc is the minimal energy interpolating spline.

In the next sections we will show how this approach can be carried over to the construction of parametric surfaces.

4. Parametric Surfaces

First we briefly review the difficulties inherent in the parametric setting. In recent years there has been considerable interest in G^{1} ((oriented) tangent plane or geometric) continuity. One reason for this is that, in general, it is not possible to construct closed surfaces defined over a standard parameter domain in IR^{2} which maintain parametric C^{1} continuity [Herron ’85].

Geometric continuity is weaker than parametric C^{1} continuity, and allows more flexibility in constructing smooth surfaces. However, it has one major drawback: the space of G^{1} splines is a nonlinear space. This is due to the fact that the necessary and sufficient G^{1} continuity conditions are nonlinear (see e.g. [Farin ’86]). One way to formulate the conditions for tangent plane continuity is

det (D_{c}, D_{R}, D_{L}) = 0,

where D_{c} stands for the derivative along the common boundary curve of two adjoining surface patches, and D_{R} and D_{L} denote the cross-boundary derivatives into the right and the left patch respectively.

Many attempts have been made to derive sufficient (linear) conditions for G^{1} continuity (see e.g. [Farin ’86], [Piper ’87], [Sarraga ’87,’89,’90], [Shirman & Séquin ’87,’91], [Hahn ’89a,’89b], [Peters ’90b,’91a], etc.). There has also been a discussion of more theoretical concepts, such as dimensions, for these “linearized” G^{1} spline spaces (cf. [Höllig & Mögerle ’90], [Goodman ’91]). However, the surface interpolation schemes in the literature not only share the common problem of how to linearize, but the need for a sensible way of determining the remaining free coefficients. The various methods use different heuristics to give the user so-called shape parameters to obtain nice looking surfaces. Most of these schemes have default values which produce surfaces with a poor distribution of curvature. This fact was pointed out in the extensive study of many of the afore mentioned schemes in [Lounsbery et al. ’92] and [Mann et al. ’92].

We examine a different approach to smooth parametrically defined surfaces. Instead of working with reparametrizations of one standard parameter domain for all pieces, such as the unit square or unit triangle, we use a different parameter domain for each individual polynomial piece, i.e. in our case the reparametrizations are built into one overall parameter domain. In particular, we use a triangulation as the parameter domain for our surface and define each coordinate as a scalar-valued spline as described in Section 2. We can therefore define a linear space of smooth parametric splines.

Definition 9. Given two integers 0 ≤ r < d and a triangulation △ of the parameter domain in the (u, v) plane, the space of parametric polynomial splines of degree d and smoothness r is defined by

SP_{d}^{r}(△) := {(g_{1}, g_{2}, g_{3})^{T} : g_{j} ∈ S_{d}^{r}(△), j = 1, 2, 3} .  (4.1)

Using a triangulation as the parameter domain for the surface corresponds to affine reparametrizations in the general theory of geometric continuity for parametric surfaces. This means we have a (1,1,1) match in the sense of [Peters ’90b,’91a]. As Peters has shown, this does not allow for the full flexibility of (nonlinear) reparametrizations, but in almost all G^{1} schemes in the literature the authors do not take full advantage of this possibility either.

However, there are great advantages to working with a fixed triangulation of the parameter domain and treating each coordinate of our parametric patch as a spline in S_{d}^{r}(△). For one thing we do not have to worry at all about the vertex enclosure problem which gives so much trouble in the existing literature on smooth parametric surfaces (and which has up to date only been solved for the symmetric case, see e.g. [VanWijk ’86], [Sarraga ’87,’89,’90], [Du ’88], [Hahn ’89b], [Peters ’92]). Our “reparametrizations” are as simple as can be: we reparametrize with the identity and therefore match Hahn’s constraint ([Hahn ’89a], Theorem 7.1) trivially.

¿From a computational point of view it is also very advantageous to treat each coordinate separately, as we will show in the next section.

5. Minimal Energy Interpolating Parametric Surfaces

Some early remarks on minimization of an energy expression in order to obtain smooth parametric surfaces can be found in [Hosaka ’69]. However, no details are given there. Most work on this subject has been done in several very recent papers (see e.g. [Hagen & Schulze ’87], [Celniker & Gossard ’91], [Kallay ’92], [Moreton & Séquin ’92]).

Hagen and Schulze use a variational approach for the nonlinear energy functional

∫_{S}(κ_{1}^{2} + κ_{2}^{2})dS,

where κ_{1} and κ_{2} are the principal curvatures of the surface S. They use biquintic patches to obtain G^{2} smooth surfaces.

Moreton and Séquin choose to minimize the variation in the derivatives of the curvature, i.e.

∫_{S} [(\frac{dκ_{n}}{dê_{1}} )^{2} + (\frac{dκ_{n}}{dê_{2}} )^{2}] dS,

where κ_{n} is the normal curvature, and ê_{1} and ê_{2} are the principal directions. They evaluate this energy integral numerically using Lobatto quadrature. Moreton and Séquin also use biquintic patches but do not even enforce full G^{1} continuity. Instead they only add the necessary and sufficient G^{1} conditions of DeRose [DeRose ’90] as a penalty term (see also Section 7). In this way they achieve approximate G^{1} continuity similarly to [DeRose & Mann ’92]. However, they do have G^{2} continuity at the vertices (a sufficient condition proposed by Peters to avoid the vertex enclosure problem (see [Peters ’91a])).

The method proposed by Celniker and Gossard uses an energy functional very similar to ours, but they also evaluate the energy integrals by Gaussian quadrature. Their method does not allow for arbitrary closed surfaces and uses “12 dof” triangles, corresponding to the well-known Clough-Tocher split (see also Section 6).

In [Kallay ’92] it is not clear whether the author calculates the entries of his matrix A and right hand side bbb in advance, or whether this is done numerically. Also Kallay’s analysis is for tensor-product B-spline surfaces only, which are not as flexible as triangular elements.

Most of the above methods are computationally quite expensive since they either work with nonlinear energy expressions, or use numerical integration techniques to evaluate the quadratic energy integrals.

We now propose a method which avoids all of the difficulties mentioned above, and cuts the amount of work required to solve the interpolation problem to essentially one third by treating each coordinate separately. In particular, given a parametric spline g ∈ SP_{d}^{0}(△), we define its energy by

J̃:= ∑_{ℓ=1}^{N} J̃_{ℓ},  (5.1)

where now

J̃_{T^{ℓ}}(g) := ∫ ∫_{T^{ℓ}} (g_{uu} · g_{uu} + 2g_{uv} · g_{uv} + g_{vv} · g_{vv}) dudv.  (5.2)

As before, g = [g_{1}(u, v), g_{2}(u, v), g_{3}(u, v)]^{T}, and “·” denotes the dot product. Multiplying out the dot products and regrouping we obtain

J̃_{T^{ℓ}}(g) = ∫ ∫_{T^{ℓ}} [(g_{1uu})^{2} + 2(g_{1uv})^{2} + (g_{1vv})^{2}] dudv
+ ∫ ∫_{T^{ℓ}} [(g_{2uu})^{2} + 2(g_{2uv})^{2} + (g_{2vv})^{2}] dudv
+ ∫ ∫_{T^{ℓ}} [(g_{3uu})^{2} + 2(g_{3uv})^{2} + (g_{3vv})^{2}] dudv,

which shows that the energy decouples, and that the total energy of the surface patch is equal to the sum of the coordinate energies.

Suppose now that the Bézier nets associated with the splines g_{1}, g_{2}, g_{3} are ccc_{x}, ccc_{y}, and ccc_{z}, respectively. Then the parametric spline g is C^{1} on △ provided that

Ã( \begin{matrix}ccc & x \\ ccc & y \\ ccc & z \end{matrix} ) = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) ,

where

Ã:= ( \begin{matrix}A & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & A \end{matrix} ) ,

and A is the matrix in (3.2). To force g to interpolate values in IR^{3} with respect to the linear functionals defined by the matrix B in (3.1), we can simply require that each of its components interpolate given values, i.e., if h h h_{x}, h h h_{y}, and h h h_{z} are given vectors of n_{int} real numbers, we require

B̃ ( \begin{matrix}ccc & x \\ ccc & y \\ ccc & z \end{matrix} ) = ( \begin{matrix}h h h & x \\ h h h & y \\ h h h & z \end{matrix} ) ,  (5.3)

where

B̃ := ( \begin{matrix}B & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & B \end{matrix} ) .

Now a minimal energy interpolating parametric spline is defined exactly as in Definition 7, but using the energy expression (5.1) and the set of splines

Ũ := {g ∈ SP_{d}^{0}(△) : Ãccc = 00 0 and B̃ccc = h h h},  (5.4)

where ccc = (ccc_{x}, ccc_{y}, ccc_{z})^{T} and h h h = (h h h_{x}, h h h_{y}, h h h_{z})^{T}.

Definition 10. A spline g ∈ Ũ is called a minimal energy C^{1} interpolating parametric spline provided

J̃(g) ≤ J̃(u) for all u ∈ Ũ.

As in Section 2, the coefficients (Bézier net) of a minimal energy interpolating spline can be found by solving a sparse system of linear equations. Let Q̃ be the 3n_{c} × 3n_{c} matrix

Q̃ := ( \begin{matrix}Q & 0 & 0 \\ 0 & Q & 0 \\ 0 & 0 & Q \end{matrix} ) .

Theorem 11. Assume Ũ ≠ ∅. Then there is a unique solution to the following constrained quadratic minimization problem:

min (ccc_{x}^{T} ccc_{y}^{T} ccc_{z}^{T})Q̃ ( \begin{matrix}ccc & x \\ ccc & y \\ ccc & z \end{matrix} ) ,

subject to

Ã( \begin{matrix}ccc & x \\ ccc & y \\ ccc & z \end{matrix} ) = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) ,  B̃ ( \begin{matrix}ccc & x \\ ccc & y \\ ccc & z \end{matrix} ) = ( \begin{matrix}h h h & x \\ h h h & y \\ h h h & z \end{matrix} ) .

It can be found by solving the linear system with multiple right-hand sides

( \begin{matrix}Q & C^{T} \\ C & 00 0 \end{matrix} ) ( \begin{matrix}ccc_{x} & ccc_{y} & ccc_{z} \\ λ λ λ_{x} & λ λ λ_{y} & λ λ λ_{z} \end{matrix} ) = ( \begin{matrix}00 0 & 00 0 & 00 0 \\ bbb_{x} & bbb_{y} & bbb_{z} \end{matrix} ) ,  (5.5)

where C is the matrix in (3.6) and

bbb_{x} := ( \begin{matrix}0 & 0 & 0 \\ h h h & x \end{matrix} ) , bbb_{y} := ( \begin{matrix}0 & 0 & 0 \\ h h h & y \end{matrix} ) , bbb_{z} := ( \begin{matrix}0 & 0 & 0 \\ h h h & z \end{matrix} ) .

Here λ λ λ_{x}, λ λ λ_{y}, λ λ λ_{z} are n_{ind} + n_{int} vectors (whose components are the Lagrange multipliers). Moreover, the parametric spline g with Bézier net ccc is the minimal energy interpolating spline.

Proof: Because of the nature of the matrices Ã, B̃, and Q̃, this quadratic minimization problem decouples into three problems of the same type as in Theorem 8.

The linear systems that arise from this global surface construction grow quite rapidly with the size of the triangulation. Therefore it is essential to have fast and accurate methods for solving large (partially sparse) linear systems. We have experimented with several methods for solving the problem: direct sparse matrix solvers (cf. [Duff ’77], [Sherman ’78], [Zlatev et al. ’81], [Eisenstat et al. ’82]), iterative sparse matrix solvers (conjugate gradient, generalized minimum residual, iterative refinement) [Seager ’89], three-stage relaxation methods on the energy part [Dyn & Ferguson ’83], quadratic minimization [Powell ’89]. Certainly, those methods which consist of a factorization as well as a solve phase have a considerable advantage in treating problems with multiple right-hand sides. We have found that we can reliably solve problems with system matrices of size up to 3000 × 3000 in only a few seconds on a Silicon Graphics Iris4D workstation using direct methods. Most other methods require setting a number of parameters whose default choice is not obvious for our type of problems, since the structure of the system matrix depends very heavily on the geometry of the domain triangulation.

Due to the global nature of our method, we find that it is particularly well suited for constructing smooth interpolating surfaces with moderately sized constellations, or localized problems such as filling n-sided holes, filleting or local modifications of existing surfaces via crowning. We discuss some of these applications in Sections 6, 7 below.

One possibility for dealing with interpolation problems involving a larger number of points would be to create several surfaces which interpolate locally to smaller subsets of the data, and then blend the constructed parts together to form an overall smooth surface. This would also permit working on smaller pieces of the problem where the surface is relatively “flat”.

6. Interpolation of Points in IR^{3}

In this section we consider the following problem.

Problem 12. Suppose Γ := {(γ_{i}^{1}, γ_{i}^{2}, γ_{i}^{3})}_{i=1}^{n} are given points in IR^{3}, and that F := {F_{i}}_{i=1}^{n} is a set of non-intersecting triangles in IR^{3} with vertices at the points Γ. Let S be the surface in IR^{3} defined by the union of the triangles in F. Find an interpolating parametric surface which is C^{1}, and which has a shape similar to S.

The surface S can be thought of as a kind of “control surface”. It is continuous and passes through the points in the set Γ. For the remainder of this section, we asssume that it is a relatively “flat” open surface (we come back to the general case including closed surfaces in Sect. 9).

Algorithm 13.

1) Choose a plane π in IR^{3} which is as “parallel” to the surface S as possible, and introduce a (u, v) Cartesian coordinate system on π. Choose a projection direction.

2) For each i = 1, …, n, let (u_{i}, v_{i}) be the point on π obtained by parallel projecting the point (γ_{i}^{1}, γ_{i}^{2}, γ_{i}^{3}) onto π.

3) Let △ be the triangulation in π obtained by projecting the edges of the triangles F_{i} onto π.

4) Find a minimal energy parametric spline g = (g_{1}, g_{2}, g_{3})^{T} in the space SP_{d}^{1}(△) which interpolates in the sense that

g_{j}(u_{i}, v_{i}) = γ_{i}^{j},  i = 1, …, n, j = 1, 2, 3.  (6.1)

The projection in step 2) can be in any arbitrary direction, as long as no triangles in IR^{3} are mapped into triangles in π with zero area. To the extent possible, it is desirable to preserve the ratios of the areas of the original triangles. Usually we would project perpendicularly onto the plane π. In step 4) we have required that the spline g interpolate the points γ_{i} at the vertices of the triangulation △. The reason for this is that for splines defined on triangulations, little is known about interpolation at points other than the vertices. For practical problems, we may choose d = 3, 4, i.e., cubic and quartic splines. For the remainder of this section we work with C^{1} cubic splines.

To illustrate this algorithm, we consider the following example. For the sake of reproducibility, we give explicit coordinates for the points in Γ, and lists of the triangles in F forming the control surface S.

Example 14. Let γ_{1} = (0, 0, 1), γ_{2} = (1, 0, 1), γ_{3} = (1, 0, 0), γ_{4} = (1, 1, 0), γ_{5} = (0, 1, 0), γ_{6} = (0, 1, 1), and γ_{7} = (.8, .8, .8). Suppose the indices of the vertices of the triangles F are given by (1, 2, 6), (2, 3, 4), (4, 5, 6), (2, 7, 6), (2, 4, 7), and (4, 6, 7).

Discussion: These 7 points lie on a rounded corner of a cube, i.e., on a suitcase corner. Fig. 4b shows the projection of Γ onto the plane x + y + z = 1, and the associated triangulation △. For each component of the spline, we have

n_{c} = total number of parameters = 37,
n_{d} = conjectured dimension = 21,
n_{int} = number of interpolation conditions = 7,
n_{f} = number of free parameters for the minimization = 14.

The total number of free parameters for the parametric spline is 3 × 14 = 42. The matrix in the linear system of Theorem 11 is of size 60 × 60. Fig. 4a shows the surface corresponding to the cubic C^{1} parametric spline interpolant. It is a reasonable approximation to the suitcase corner.

Clearly, for a given set of points Γ, if we choose to connect them in different ways to get a different control surface, we get a different parametric spline interpolant.

Example 15. Let Γ be as in Example 14, but suppose F now consists of the following set of triangles: (1, 2, 7), (2, 3, 7), (3, 4, 7), (4, 5, 7), (5, 6, 7), and (6, 1, 7).

Discussion: Fig. 4d shows the projection of Γ onto the plane π of Example 14, and the associated triangulation △. Here the control surface S is geometrically identical with the one in Example 14, but consists of different triangles. Fig. 4c shows the surface corresponding to the cubic C^{1} parametric spline interpolant. The change in the parameter domain did not cause a major change in the shape of the surface.

For a given set of points Γ and control surface S, the triangulation △ underlying the parametric spline space depends quite heavily on the choice of the projection plane π in step 1) of our algorithm. The following example illustrates this effect.

Example 16. Let Γ and F be as in Example 14, but suppose the projection plane π is now defined by x+5y +z = 1 and the projection direction is given by (5, 2, 3)^{T}.

Fig. 4. A suitcase corner with three different parametrizations. 18 Discussion: Fig. 4f shows the projection of Γ onto the plane π and the associated triangulation △. It is, of course, just a distortion of the triangulation in Fig. 4b. Fig. 4e shows the associated cubic C^{1} parametric spline interpolant. In this case the suitcase corner has been significantly distorted.

To get a better approximation to the suitcase corner discussed in the above examples, we can interpolate at more points.

Example 17. Let Γ consist of the 7 points of Example 14 plus the following additional 9 points: (1, 0, .5), (1, .5, 0), (1, .5, .5), (.5, 0, 1), (.5, .5, 1), (.5, 1, 0), (.5, 1, .5), (0, .5, 1), and (0, 1, .5). The plane π and projection direction are the same as in Example 14.

Discussion: We do not list the triangles in F since it is clear what they are from Fig. 5b. For each component of the spline, we now have

n_{c} = total number of parameters = 100,
n_{d} = conjectured dimension = 45,
n_{int} = number of interpolation conditions = 16,
n_{f} = number of free parameters for the minimization = 29.

The total number of free parameters for the parametric spline is 3 × 29 = 87. The matrix in the linear system of Theorem 11 is of size 171 × 171. Fig. 5a shows the associated cubic C^{1} parametric spline interpolant, which is now an excellent approximation to the suitcase corner.

To conclude this section, we look at an even more extreme parametrization.

Example 18. Let Γ be as in Example 14, and suppose we choose the triangulation shown in Fig. 6b.

Discussion: This triangulation has only four boundary edges, but since it has 7 vertices, we can still construct a cubic C^{1} parametric spline surface which interpolates the given data. Fig. 6a shows the resulting surface. The reason for its extreme shape is due to the fact that here we have not specified a reasonable control surface, and so the triangular patches connect the given points on our “suitcase corner” in an unnatural way.

Fig. 5. A suitcase corner obtained by interpolating 16 points.

Fig. 6. A bad parametrization for a “suitcase corner”.

7. Filling n-sided Holes and Crowning

The problem of filling an n-sided hole in a surface is one of the central problems of CAGD. For our purposes, we assume that the hole is cut out of an already existing smooth piecewise polynomial surface, and that the boundary of the hole consists of n curves. The problem is to create a surface patch which fills the hole, and which joins smoothly to the outside.

The problem of crowning a surface is closely related. We cut an n-sided hole in a given surface, and then fill it with a patch which is close to the original surface, but which may be raised or lowered somewhat. This can be accomplished by making the patch interpolate some specified point(s) in IR^{3} in addition to matching the boundary curves and the cross derivatives along the boundary curves of the hole.

There have been many suggestions for smoothly filling in an n-sided hole in an already existing smooth surface. The available methods include blending (convex combination) schemes ([Gregory ’86], [Gregory & Hahn ’89]), splitting schemes with triangular ([Shirman & Séquin ’87], [Peters ’90a]) or rectangular pieces ([Hahn ’89a,’89b], [Sarraga ’90]), overlap patches [Varady ’90], subdivision ([Doo & Sabin ’78], [Storry & Ball ’89], [Peters ’93]), and methods that use special nonplanar domains ([Sabin ’83], [Hosaka & Kimura ’84], [Loop & DeRose ’89,’90]). For some good surveys on this subject see [Gregory, Lau & Zhou ’90], [Peters ’90b] or [Varady ’90].

In order to create a surface to fill a hole with a patch maintaining C^{1} parametric continuity, we need to make a compatibility assumption on the initial surface surrounding the hole. We assume that it is given in piecewise polynomial form, and that the derivatives across the boundaries of the hole are such that they maintain C^{1,1} continuity at the corners of the hole (for similar assumptions see [Gregory, Lau & Zhou ’90], [Loop & DeRose ’90]).

In this section we examine the use of our parametric spline method to solve this problem. To apply the method, we have to create a planar triangulation △ to serve as the domain of the spline. In this setting, it seems natural to choose △ so that it has n vertices on the boundary.

To make our surface patch join smoothly with the rest of the surface, we require it to interpolate the surface values at each of the boundary vertices of △, and in addition

1) choose the degree of the patch to match the degree of the outside surface patches

2) construct the interpolant so that the boundary curves match the boundary curves of the hole, and so that the cross boundary derivatives of the patch match the cross-boundary derivatives defined by the outside surface.

One way to create the triangulation △ is to choose some point in IR^{3} which we want our crowned surface to go through, connect it to the n points on the boundary of the hole, and then follow the first three steps in Algorithm 13. This typically gives a triangulation of an n-gon with one interior vertex. Unfortunately, using low degree C^{1} spline spaces we do not get enough free parameters to match all of the required boundary information (in the cubic case), or there are very few free parameters (e.g., in the quartic case). Thus, this triangulation usually has to be further subdivided.

The alternative is to construct the triangulation △ directly in a (u, v)-plane. Given an n-gon which is to be the domain Ω of our spline space, we propose two natural ways to triangulate Ω. The first corresponds to triangulating the n-gon with a single interior vertex, and then splitting each of the resulting triangles into three subtriangles as is done in the well-known Clough-Tocher method of finite element theory.

Triangulation 1. (see Fig. 7a)

1) Find the centroid of the vertices surrounding the hole and connect all the vertices to the centroid.

2) Find the centroid of each of the subtriangles constructed in step 1) and connect the centroids to the vertices of the subtriangles.

Fig. 7. The configurations for (a) Triangulation 1, (b) Triangulation 2 (n = 5).

Our second triangulation involves the same number of interior vertices, but positioned in a somewhat more symmetric way.

Triangulation 2. (see Fig. 7b)

1) Find the centroid of the vertices surrounding the hole.

2) Find the midpoints of the boundary edges (marked with ∗) and then the midpoints (marked with •) of the line segments connecting the boundary midpoints and the centroid obtained in step 1).

3) Connect all of the vertices marked with • and also connect each of them with the centroid of the hole.

4) Connect each of the boundary vertices with the two closest vertices marked with •.

Neither triangulation requires Ω to be convex. However, it should be such that its centroid lies inside the boundary polygon.

The following theorem shows that if we use cubic C^{1} splines matching values and derivatives on the boundary of Triangulation 1, the basic hypothesis (5.4) for the existence of an interpolating spline is satisfied.

Theorem 19. Suppose we are given an n-sided hole in a piecewise polynomial surface with (at most) cubic boundary curves and (at most) quadratic cross-boundary derivatives, and satisfying the basic compatibility condition discussed above. Let △_{1} be a triangulation as shown in Fig. 7a. Then there exists a unique C^{1} cubic parametric spline such that the corresponding surface fills the hole, matches the given surface smoothly (C^{1}) along the boundary, and has a prescribed value and tangent plane at the parameter point corresponding to the center vertex of the triangulation.

Proof: The fact that there is a unique scalar-valued C^{1} cubic spline which takes on given values and gradients at the n boundary vertices and the center vertex of △_{1} is well-known from finite element theory (see e.g., [Alfeld ’84a]). We can think of the resulting surface as a kind of macro element consisting of 3n cubic triangular patches which fill the hole and join with the outside surface with C^{1} smoothness across all patch boundaries.

It is not known if the analogous result holds for Triangulation 2 using cubic splines, but our computational experience suggests that it does. We discuss the use of quartic splines with Triangulation 2 at the end of this section.

Example 20. Consider the function f(x, y) = x^{2}+y^{2} defined on a regular hexagon surrounding the origin with one vertex at (2, 0). Suppose we create a symmetric six-sided hole surrounding the origin as shown in Fig. 8. Fill this hole with a patch which interpolates the value 1.5 at (0, 0) with a horizontal tangent plane.

Discussion: We use Triangulation 2 and cubic C^{1} (scalar) splines. Fig. 8 shows the resulting surface with the fill. The views on the left correspond to a slightly higher viewing angle. The top two images show lines corresponding to a few steps of subdivision of the patches. These clearly reveal the hole which has been filled. The middle two images show a curvature plot of the fill (unfortunately, this is hard to depict in black and white). The bottom images show the surface and the patch using Gouraud shading.

As observed above, if we use C^{1} cubic splines to fill a hole which has been triangulated as in Triangulation 2, there is no guarantee that the resulting system of equations is nonsingular. However, we can establish the non-singularity for C^{1} quartic splines.

Theorem 21. Suppose we are given an n-sided hole surrounded by a piecewise polynomial surface with (at most) quartic boundary curves and (at most) cubic Fig. 8. Filling a symmetric 6-sided hole.

24 cross-boundary derivatives, and satisfying the basic compatibility condition discussed above. Let △_{2} be a triangulation as shown in Fig. 7b. Then there exists a unique C^{1} quartic parametric spline such that the corresponding surface fills the hole and matches the given surface smoothly (C^{1}) along the boundary. Moreover, at the middle splitting point, we can specify the position of the surface and the tangent plane to it, provided n is odd. If n is even, then we can specify the same

information at the middle point provided ∏_{i=1}^{n} r_{i} ≠ ∏_{i=1}^{n} t_{i} (see the proof for details).

Proof: Throughout the proof we will refer to Fig. 9 and the notation used there. We will show that, given positional and cross-boundary information from the outside surface, we have enough free coefficients inside the hole to satisfy a C^{1} join across all patch boundaries while choosing the center point and its tangent plane arbitrarily.

Fig. 9. The notation for a quartic fill in the case n = 5.

The key to the proof is the analysis of a cyclic nonhomogeneous linear system involving the coefficients numbered 1 through n in the figure.

Let us denote all coefficients determined by the outside world by ∗. If we choose those marked with ⋄_{1} arbitrarily, then the coefficients represented by •_{1} are determined by C^{1} continuity conditions across the respective edges. In the next step we choose the ⋄_{2} coefficients freely. Together with the ∗, ⋄_{1} and •_{1} coefficients they uniquely determine those marked with •_{2}. In order to fix the tangent planes at the vertices of the inner ring, we set the coefficients with ⋄_{3} and then also have those with •_{3}. The last set of coefficients we get to choose freely is those denoted by ⋄_{4} in the second layer around the center vertex; they then determine the •_{4} coefficients (in case the •_{4} coefficient cannot be determined because the quadrilateral has a straight edge through the center vertex, we get to choose this coefficient also). Now we are left with the coefficients numbered 1 through n and n + 1 through 3n. They are involved in one continuity condition across each of the edges emanating from the center vertex, and in two C^{1} conditions across the boundaries between the two rings of triangles (the continuity conditions are crosshatched in the figure). If we label the latter edges counterclockwise from e_{1} to e_{n}, then we can write the first set of continuity conditions as follows:

\begin{aligned}c_{2i+n−1} & = t_{i}c_{i} + α_{i} \\ c_{2i+n} & = r_{i}c_{i} + β_{i}, \end{aligned}  (7.1)

where i = 1, …, n, the c’s are the unknown coefficients numbered as in the figure, and (r_{i}, s_{i}, t_{i}) are the barycentric coordinates associated with the quadrilateral to both sides of edge e_{i}, where we always express the center vertex as r_{i} times the “right” vertex of edge e_{i} plus s_{i} times the associated “outside” vertex plus t_{i} times the “left” vertex of edge e_{i}. The quantities α_{i} and β_{i} on the right-hand side contain the already determined coefficients also present in the continuity conditions.

The second set of continuity conditions across the innermost edges, labelled counterclockwise with e_{n+1} through e_{2n}, turns out to be

c_{2i−n−1} = s_{i}c_{2i−n−2} + γ_{i},  i = n + 1, …, 2n,  (7.2)

where now the barycentric coordinates (r_{i}, s_{i}, t_{i}) reflect the shapes of two neighboring triangles when going around the center vertex counterclockwise. Again the γ_{i} contain already set coefficients. This set of equations is of a cyclic nature, so we therefore have to identify the coefficient c_{n} as c_{3n}.

We can now combine (7.1) and (7.2) and obtain a cyclic linear system in the coefficients c_{1} through c_{n} which has the following structure

( \begin{matrix}−s_{n+1}r_{1} & t_{2} & 0 & 0 & … & 0 \\ 0 & −s_{n+2}r_{2} & t_{3} & 0 & … & 0 \\ ⋮ & & ... & ... & & ⋮ \\ 0 & & & ... & ... & 0 \\ 0 & & & & ... & t_{n} \\ t_{1} & 0 & 0 & … & 0 & −s_{2n}r_{n} \end{matrix} ) ( c_{1} ; c_{2} ; ⋮ ; c_{n} ) = ( ρ_{1} ; ρ_{2} ; ⋮ ; ρ_{n} ) ,

where (r_{i}, s_{i}, t_{i}) are as above and the ρ_{i} (i = 1, …, n) are the known right-hand sides. We can easily calculate the determinant of the coefficient matrix as

Det = (−1)^{n} ∏_{i=n+1}^{2n}s_{i} ∏_{i=1}^{n} r_{i} + (−1)^{n−1} ∏_{i=1}^{n} t_{i}.

Since we have ∏_{i=n+1}^{2n}s_{i} = (−1)^{n} due to the fact that the cycle around the center

vertex closes on itself and barycentric coordinates can be viewed as ratios of areas, the expression for the determinant simplifies to

Det = ∏_{i=1}^{n} r_{i} + (−1)^{n−1} ∏_{i=1}^{n} t_{i}.

We now have to distinguish two cases:

Case 1. (n is odd) In this case we have

Det = ∏_{i=1}^{n} r_{i} + ∏_{i=1}^{n} t_{i} ≠ 0,

since in our particular split all relevant quadrilaterals are convex, which means that all of the r_{i} and t_{i} (i = 1, …, n) are positive. So the linear system has a unique solution and we can fill in a hole with an odd number of vertices.

Case 2. (n is even) Now the determinant turns out to be

Det = ∏_{i=1}^{n} r_{i} − ∏_{i=1}^{n} t_{i},

and we cannot positively say whether this expression is nonzero. Therefore a unique solution for the even case is guaranteed to exist for data that satisfy this condition.

It is easy to construct examples to show that the condition in Theorem 21 for the case n even is sufficient, but not necessary (any symmetric setting will do). It is obvious from the proof of this theorem that instead of specifying just the center point of the split with its tangent plane, we could also set all other points marked with a ⋄ in Fig. 9.

Instead of using C^{1} cubic or quartic splines, we could also fill holes using a C^{1} quadratic spline defined on the triangulation obtained by splitting each initial triangle into 6n subtriangles according to the well-known Powell-Sabin split (cf. [Dierckx ’93]). This would assure that values and tangent planes can be specified at all vertices of the original triangulation, and requires that fewer coefficients be determined (and stored). (Here there are 14n + 1 coefficients as compared with 15n + 1 for the cubic scheme on the Clough-Tocher split). However, working with quadratics only also restricts the maximal degree of the surface representation on the outside to quadratics, which might be too low for many applications.

Fig. 10. A typical C^{2} continuity condition (d = 3).

The macro elements discussed here can be viewed as locally supported elements similar to B-splines. However, it is not clear if there is a way of using them to construct a basis for S_{d}^{r}(△) or SP_{d}^{r}(△).

8. Adding a C^{2} Penalty Term

Most surface interpolation methods (cf. those mentioned in Sections 4 and 5) achieve only C^{1} or G^{1} continuity, and some of them have only approximate G^{1} continuity. Using our approach, however, it is not much additional work to incorporate conditions that yield approximate C^{2} continuity, and in fact we can do this with total degree cubic patches as opposed to e.g. the (bi-)quintic patches used by [Hagen & Schulze ’87], or [Moreton & Séquin ’92] who achieve only approximate G^{1} continuity with G^{2} continuity at the vertices!

The idea is to add the C^{2} continuity conditions to the energy functional we are minimizing in the form of a least squares penalty term. We first discuss this idea in the case of surfaces corresponding to scalar-valued splines. It will also suffice to look only at the case of cubic splines; the general case is completely analogous.

As is well-known, cf. [Farin’86], a C^{1} cubic spline is C^{2} continuous across an edge of the triangulation (cf. Fig. 10), if and only if the sum

∑_{i=1}^{2} [c_{i}^{rr} − (r^{2}c_{i+2}^{m} + 2rsc_{i+1}^{l} + s^{2}c_{i}^{ll} + 2rtc_{i+1}^{m} + 2stc_{i}^{l} + t^{2}c_{i}^{m})]^{2}

is zero. Here (r, s, t) are the barycentric coordinates associated with the edge, and the coefficients are as in Fig. 10.

Thus, to make a spline g ∈ S_{3}^{1}(△) be approximately C^{2}, we should minimize the sum of all such terms across all interior edges of △. In view of the quadratic form of the expressions in this sum, it can be written in the form

ccc^{T} Pccc.

Thus, we look for a spline which minimizes a combination of energy and this smoothness term:

J_{ω}(g) = ccc^{T} (Q + ωP)ccc.

Here ω is a parameter controlling the relative weight of energy versus smoothness. For more on penalized least squares in approximation theory and the use of generalized cross validation for the choice of the smoothing parameter, see [von Golitschek & Schumaker ’90].

Interpolating splines which minimize J_{ω} can be found by solving a simple linear system of equations. Indeed, Theorem 8 remains valid if we simply replace Q by Q+ωP. Thus to compute these splines, the only extra work that has to be performed is the creation of the entries of the matrix P. This can be done in the same step in which the C^{1} continuity conditions are being constructed, and therefore does not require much overhead.

The extension of this idea to surfaces defined by parametric splines is obvious. In particular, Theorem 11 also holds if we replace Q by Q + ωP. One could slightly modify the C^{2} penalty term defined above by assigning individual weights ω_{i}, i = 1, …, N, to each triangular patch. This would give an even better control over the smoothing effect.

Fig. 11. The effects of the C^{2} penalty term (a) ω = 0, (b) ω = 100.

To see the effects of the C^{2} penalty term, we solve a simple scalar interpolation problem both with and without the penalty term.

Example 22. Consider the function

f(x, y) = (1 + 2e^{−3[9(x^{2}+y^{2})^{1/2}−6.7]})^{−1/2}

defined on a regular hexagon with one vertex at (2, 0). Let △ be the Delaunay triangulation corresponding to two concentric hexagonal rings and the origin. Find a cubic C^{1} spline on △ which interpolates function values at the 13 vertices of △.

Discussion: The surface on the left in Fig. 11 is constructed with ω = 0, while for the surface on the right ω = 100. Areas of high positive or negative Gaussian curvature are highlighted (unfortunately, in this gray-scale rendering one cannot distinguish between the two extremes). Nevertheless, the C^{2} smoothing effect is clearly demonstrated. This high degree of smoothness is achieved using only cubic splines.

9. Closed Surfaces

In this section we explore the possibility of using the above ideas to create closed surfaces which interpolate given sets of points in IR^{3}. We restrict our attention to closed surfaces which are topologically equivalent to the sphere.

Problem 23. Suppose Γ := {(γ_{i}^{1}, γ_{i}^{2}, γ_{i}^{3})}_{i=1}^{n} are given points in IR^{3}, and that F := {F_{i}}_{i=1}^{n} is a set of non-intersecting triangles in IR^{3} with vertices at the points Γ. Suppose the surface S in IR^{3} defined by the union of the triangles in F is closed. Find an interpolating parametric surface which is smooth, and which has a shape similar to S.

As in Section 6, the surface S can be thought of as a kind of “control surface”. To model S with a parametric spline, we shall use an analog of Algorithm 13. Before describing an algorithm, we need to introduce two new ideas: 1) generalized triangulations, and 2) splines on generalized triangulations. The triangulations introduced in the following definition are generalizations of ordinary planar triangulations in the sense that triangles are allowd to lie on top of each other.

Definition 24. Let Δ = {T^{ℓ}}_{ℓ=1}^{N} be a collection of triangles, and let V = {V_{i}}_{i=1}^{V} be the set of all vertices involved. We say that Δ is a generalized triangulation of

Ω := ∪_{ℓ=1}^{N}T^{ℓ}

provided that

1) each triangle T^{ℓ} contains no points from V other than its vertices,

2) any two triangles either do not intersect at all, or share one, two, or three edges,

3) no more than two triangles share any edge,

4) Ω is connected.

This definition extends the notion of generalized triangulations introduced in [Alfeld, Schumaker & Whiteley’93]. There generalized triangulations were defined by projecting a set of 3D triangles in a special way.

Here is an example showing how such a generalized triangulation can be created by projecting a control surface as described in Problem 23 onto a plane.

Example 25. Find a closed surface which interpolates the six points v_{1} = (1, 0, 0), v_{2} = (0, 1, 0), v_{3} = (−1, 0, 0), v_{4} = (0, −1, 0), v_{5} = (0, 0, 1), and v_{6} = (0, 0, −1).

Discussion: These six points form the vertices of an octahedron in IR^{3}, each of whose faces is a triangle. We can consider the union of these eight triangles to be a closed control surface S. Now if we project S along the z-axis onto the xy-plane, the eight triangular faces are mapped into eight planar triangles, four of which lie on top, and four of which lie on the bottom, see Fig. 12. Note that the images of v_{5} and v_{6} both appear to be at the origin, but v_{5} is on top, and v_{6} is on the bottom.

Fig. 12. The generalized triangulation for Example 25.

For later use, we need some additional terminology.

Definition 26. Given a generalized triangulation Δ, an edge which is shared by two triangles is called an interior edge. All other edges are called boundary edges. Vertices which lie on boundary edges are called boundary vertices, and all others are called interior vertices. An interior edge is said to be folded provided that one of the two triangles sharing that edge is on top of the other. A vertex is said to be folded provided it lies on a folded edge. A folded vertex is said to be exceptional provided that its two adjacent folded edges are not collinear.

Clearly, if a generalized triangulation arises by collapsing a closed control surface as above, then all edges and vertices are interior. In the generalized triangulation shown on the left in Fig. 12, the edges drawn with solid lines are folded edges, while those drawn with dashed lines are not. The vertices v_{5} and v_{6} are ordinary vertices (not folded), while v_{i}, i = 1, …, 4, are folded exceptional vertices.

We now introduce a space of splines defined on a generalized triangulation. It should consist of scalar-valued functions defined on the triangles of Δ. Thus, given an integer d > 0, we define the associated space of piecewise polynomials to be

PP_{d}(Δ) := {g : g|_{T} ∈ P_{d}^{(2)}, for all T ∈ Δ},

where as before, P_{d}^{(2)} is the space of bivariate polynomials of total degree d.

Each piece of an element g ∈ PP_{d}(Δ) is a polynomial of total degree d on a planar triangle, and thus can be expressed in Bernstein-Bézier form as described in Sect. 1. This means we can define an associated coefficient vector ccc in the same way as in Sect. 2 for splines on ordinary triangulations (where the coefficients of patches which are associated with adjoining triangles and which lie on the common edge between the two triangles are identified). We denote the length of ccc by n_{c} as before. We can now describe a kind of spline space where adjoining patches are required to “match” along their common edge.

Definition 27. Let Δ be a generalized triangulation, and let 0 ≤ r < d be given integers. We define the space of splines by

S_{d}^{r}(Δ) := {g ∈ PP_{d}(Δ) : H_{e}^{r}ccc = 0 for all interior edges e},

where H_{e}^{r} is the m_{r} × n_{c} matrix with m_{r} := d + (d − 1) + ⋯ + (d − r + 1) which describes the C^{r} continuity across the edge e, except that if e is a folded edge, then in computing the barycentric coordinates, we first fold out the edge. Thus, for example in Fig. 12, since the edge between v_{1} and v_{2} is a folded edge and v_{6} is underneath, we first fold out the triangle on the bottom before computing the barycentric coordinates (r, s, t) of v_{6} in terms of the triangle with vertices v_{2}, v_{5}, and v_{1}. The conditions across non-folded edges are exactly as in the case of ordinary splines.

The space of splines defined on generalized triangulations introduced in [Alfeld, Schumaker & Whiteley’93] is not the same as ours, since there the space was defined using the usual continuity conditions even for folded edges (without folding them out).

It is clear that spaces of parametric splines on a generalized triangulation can be defined in exactly the same way as was done for ordinary triangulations, see Definition 9. Any such parametric spline defines a surface in IR^{3} defined on the parameter domain Ω.

Because of the unusual nature of the space S_{d}^{r}(Δ), we should say something about the smoothness of the resulting parametric surface at a folded vertex v when r = 1. First, if v is not exceptional, then the surface is visually smooth at v in the usual sense, i.e., it has a continuous tangent plane in a neighborhood of v. The reason for this is that in this case folding out the bottom triangles can be regarded as a reparametrization based on an ordinary triangulation.

If v is a folded exceptional vertex, then as the following lemma shows, the Bézier coefficients in the first ring around v must all be equal to the coefficient at the vertex, and thus (as is well-known), the surface may have a cusp at v. This corresponds to a singular parametrization. Such cusps can be removed by appropriately setting Bézier coefficients in the second ring (cf. [Neamtu ’91], [Peters ’91b], [Reif ’93]). We illustrate this in Fig. 14.

Theorem 28. Let v be an exceptional folded vertex with m + k edges attached. Suppose the coefficients of g ∈ S_{d}^{0} satisfy the conditions H_{e}^{1}ccc = 0 for all edges e attached to the vertex v. Then all coefficients in the first ring around v must have the same value as the coefficient at v.

Proof: For the proof it will suffice to consider only those triangles of the generalized triangulation Δ which are attached to the vertex v. Assume the triangles have vertices v_{i}, v_{i+1}, v, and that those on top correspond to i = 1, …, m, while those on the bottom correspond to i = m, …, m + k, where v_{m+k+1} = v_{1}. Now cut along edge vv_{m+k} and fold out all of the bottom triangles. Let ṽ_{i} be the new location of the vertex v_{i} for i = m + 1, …, m + k − 1. By the assumption that v is exceptional, the vertex v_{m+k} moves to two different positions, ṽ_{m+k} and v̂_{m+k} associated with the folded out triangles vṽ_{m+k−1}ṽ_{m+k} and vv̂_{m+k}v_{1}.

Let the Bézier coefficient associated with v be c_{0}. By the C^{1} continuity conditions, all of the points (v_{i}, c_{i}) for i = 1, …, m as well as the points (ṽ_{i}, c_{i}) for i = m + 1, …, m + k and (v̂_{m+k}, c_{m+k}) lie on a plane passing through (v, c_{0}). In addition, the points (v_{i}, c_{i}) for i = m, …, m+k +1 must also lie on a plane passing through (v, c_{0}). Since the points v_{1}, v, v_{m} are noncollinear, the two planes coincide with a common plane. This plane must be the horizontal plane through (v, c_{0}) since it contains the points (v_{m+k}, c_{m+k}), (v̂_{m+k}, c_{m+k}), (ṽ_{m+k}, c_{m+k}). We conclude that c_{i} = c_{0} for all i = 1, …, m + k.

Theorem 28 implies that if gg g = (g_{1}, g_{2}, g_{3})^{T} is a parametric spline, then the corresponding control points (which are then 3-vectors) all collapse to the control point at v.

To solve Problem 23 using the space of splines S_{3}^{1}(Δ) defined on a generalized triangulation, we now introduce the total energy of such a spline surface as the sum of the individual energies on triangles, i.e., J̃(g) can be defined exactly as in (5.1). Now we can define a minimal energy interpolating parametric spline just as in Definition 10: it minimizes the quadratic form J̃(g) subject to smoothness and interpolation conditions. As in Theorem 11, this reduces to solving a linear system of equations similar to (5.5).

We can now present an algorithm for interpolating data in IR^{3} with a closed parametric spline surface.

Algorithm 29.

1) Choose a plane π in IR^{3} and introduce a (u, v) Cartesian coordinate system on it. Choose a projection direction.

2) For each i = 1, …, n, let (u_{i}, v_{i}) be the point on π obtained by projecting the point (γ_{i}^{1}, γ_{i}^{2}, γ_{i}^{3}) onto π.

3) Let Δ be the generalized triangulation in π obtained by parallel projecting the edges of the triangles F_{i} onto π.

4) Find a minimal energy parametric spline g = (g_{1}, g_{2}, g_{3})^{T} in the space S_{3}^{1}(Δ) which interpolates in the sense that

g_{j}(u_{i}, v_{i}) = γ_{i}^{j},  i = 1, …, n, j = 1, 2, 3.  (9.1)

The projection in step 2) can be in any arbitrary direction, as long as no triangles in IR^{3} are mapped into triangles in π with zero area.

Fig. 13. A closed surface with singularities.

Fig. 13a shows the result of using splines in S_{3}^{1}(Δ) with the generalized triangulation shown in Fig. 13b to interpolate the data in Example 25. In this example the singularities in the parametrization occur at the corners of the square, which are folded exceptional vertices in this case. The resulting surface is smooth across the interiors of all edges, but as the figure shows, because of the singularities, there are cusps at the corners.

To smooth out the cusps at the corners, we use the idea mentioned above of setting the Bézier coefficients in the second ring appropriately. To do this, we have to refine the generalized triangulation of Fig. 12 in such a way as to separate the folded exceptional vertices, and to provide enough coefficients to adjust. The top half of our refined generalized triangulation is shown in Fig. 14b.

Fig. 14a shows the minimal energy parametric spline in SP_{3}^{1}(Δ) which interpolates the data of Example 25 based on this triangulation. Here we have set two control points in the second layer around each folded exceptional vertex to form a vertical tangent plane there. It should be noted that although the surface looks smooth, it has infinite curvature at folded corners.

Fig. 14. The singularities of Fig. 13 smoothed out.

A disadvantage when using generalized triangulations and singular parametrizations is that, as soon as we define a common tangent plane at a vertex (e.g. by specifying its normal vector, instead of explicitly fixing certain coefficients as we did for Fig. 14), the decoupling of the three spline coordinates is destroyed and the size of the system matrix grows by a factor of 9. However, it does have the block structure mentioned in Section 5. We have not attempted to exploit this fact.

Another problem associated with the use of generalized triangulations and singular parametrizations is that the redundancies in the C^{1} continuity conditions behave differently than in the regular case. In the regular case one simply ignores the last two conditions in the ring around an interior vertex as well as an additional one at a singular vertex (see also Section 2). In the case of a singular parametrization, this is not necessarily true, and now we have to use some method such as singular value decomposition to numerically find the dependent rows of the matrix A describing the continuity conditions.

References

1. P. Alfeld, A trivariate Clough-Tocher scheme for tetrahedral data, Comput. Aided Geom. Design 1 (1984), 169–181.

2. P. Alfeld, Triangular Extrapolation, MRC Technical Summary Report #2707, University of Wisconsin–Madison, 1984.

3. P. Alfeld, B. Piper, and L. L. Schumaker, An explicit basis for C^{1} quartic bivariate splines, SIAM J. Numer. Anal. 24 (1987), 891–911.

4. P. Alfeld, and L. L. Schumaker, On the Dimension of Bivariate Spline Spaces of Smoothness r and Degree d = 3r + 1, Numer. Math. 57 (1990), 651–661.

5. P. Alfeld, L. L. Schumaker, and W. Whiteley, The generic dimension of the space of C^{1} splines of degree d ≥ 8 on tetrahedral decompositions, SIAM J. Numer. Anal. 30 (1993), 889–920.

6. C. de Boor, B-Form Basics, in Geometric Modeling: Algorithms and New Trends, G. E. Farin (ed.), SIAM Publications (Philadelphia), (1987), 131– 148.

7. G. Celniker, and D. Gossard, Deformable curve and surface finite-elements for free-form shape design, Computer Graphics (ACM SIGGRAPH) 25 (1991), 257–266.

8. S. Coons, Surfaces for computer aided design, MAC-TR-4, M.I.T., 1964.

9. T. DeRose, Necessary and sufficient conditions for tangent plane continuity of Bézier surfaces, Comput. Aided Geom. Design 7 (1990), 165–179.

10. T. DeRose, and S. Mann, An Approximately G^{1} Cubic Surface Interpolant, in Mathematical Methods in Computer Aided Geometric Design II, T. Lyche and L. Schumaker (eds.), Academic Press (New York), (1992), 185–196.

11. P. Dierckx, Curve and Surface Fitting with Splines, Clarendon Press, Oxford, 1993.

12. D. Doo, and M. A. Sabin, Behaviour of recursive division surfaces near extraordinary points, Computer-Aided Design 10 (1978), 350–355.

13. W. Du, Etude sur la représentation de surfaces complexes : application à la reconstruction de surfaces échantillonnées, dissertation, Télécom Paris, 1988.

14. I. S. Duff, MA28 – a set of Fortran subroutines for sparse unsymmetric linear equations, Report AERE R8730, HMSO, London, 1977.

15. N. Dyn, and W. Ferguson, The numerical solution of equality constrained quadratic programming problems, Math. Comp. 41 (1983), 165–170.

16. S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman, Yale sparse matrix package. 1: The symmetric codes, Internat. J. Numer. Meth. Engr. 18 (1982), 1145–1151.

17. G. Farin, Triangular Bernstein-Bézier patches, Comput. Aided Geom. Design 3 (1986), 83–127.

18. R. H. J. Gmelig Meyling, Approximation by cubic C^{1} splines on arbitrary triangulations, Numer. Math. 51 (1987), 65–85.

19. T. N. T. Goodman, Closed surfaces defined from biquadratic splines, Constr. Approx. 7 (1991), 149–160.

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

21. W. Gordon, Free-form surface interpolation through curve networks, GMR-921, GM Research Labs, 1969.

22. J. A. Gregory, N-sided surface patches, in The Mathematics of Surfaces, J.A. Gregory (ed.), Clarendon Press (Oxford), (1986), 217–232.

23. J. A. Gregory, and J. M. Hahn, A C^{2} polygonal surface patch, Comput. Aided Geom. Design 6 (1989), 69–75.

24. J. Gregory, Vincent K.H. Lau, and Jianwei Zhou, Smooth parametric surfaces and n-sided patches, in Computation of Curves and Surfaces, W. Dahmen, M. Gasca and C. Micchelli (eds.), Kluwer (Dordrecht, Netherlands), (1990), 457–498.

25. H. Hagen, and G. Schulze, Automatic smoothing with geometric surface patches, Comput. Aided Geom. Design 4 (1987), 231–235.

26. J. M. Hahn, Geometric continuous patch complexes, Comput. Aided Geom. Design 6 (1989), 55–67.

27. J. M. Hahn, Filling polygonal holes with rectangular patches, in Theory and Practice of Geometric Modelling, W. Strasser and H. P. Seidel (eds.), Springer, 1989, 81–91.

28. G. Herron, Smooth closed surfaces with discrete triangular interpolants, Comput. Aided Geom. Design 2 (1985), 297–306.

29. K. Höllig, and H. Mögerle, G-splines, Comput. Aided Geom. Design 7 (1990), 197–207.

30. M. Hosaka, Theory of curve and surface synthesis and their smooth fitting, Information Processing in Japan 9 (1969), 60–68.

31. M. Hosaka, and F. Kimura, Non-four-sided patch expressions with control points, Comput. Aided Geom. Design 1 (1984), 75–86.

32. M. Kallay, Constrained optimization in surface design, EDS, 1992.

33. C. Loop, and T. DeRose, A multisided generalization of Bézier surfaces, Transactions on Graphics 8 (1989), 204–234.

34. C. Loop, and T. DeRose, Generalized B-spline surfaces of arbitrary topology, SIGGRAPH’90 24 (1990), 347–356.

35. M. Lounsbery, S. Mann, and T. De Rose, An overview of parametric surface interpolation, Comp. Graphics and Applics. 12(5) (1992), 45–52.

36. S. Mann, C. Loop, M. Lounsbery, D. Meyers, J. Painter, T. DeRose, and K. Sloan, A survey of parametric scattered data fitting using triangular interpolants, in Curve and Surface Design, H. Hagen (ed.), SIAM Publications, SIAM (Philadelphia PA), (1992), 145–172.

37. H. Moreton, and C. Séquin, Functional Optimization for Fair Surface Design, Computer Graphics (ACM SIGGRAPH) 26 (1992), 167–176.

38. M. Neamtu, Multivariate Splines, dissertation, University of Twente, Holland, 1991.

39. J. Peters, Smooth mesh interpolation with cubic patches, Computer-Aided Design 22 (1990), 109–120.

40. J. Peters, Local smooth surface interpolation: a classification, Comput. Aided Geom. Design 7 (1990), 191–195.

41. J. Peters, Smooth interpolation of a mesh of curves, Constr. Approx. 7 (1991), 221–246.

42. J. Peters, Parametrizing Singularly to Enclose Data Points by a Smooth Parametric Surface, in Proceedings of Graphics Interface ’91, Calgary, 1991.

43. J. Peters, Joining smooth patches around a vertex to form a C^{k} surface, Comput. Aided Geom. Design 9 (1992), 387–411.

44. J. Peters, Smooth free-form surfaces over irregular meshes generalizing quadratic splines, CSD-TR-92-063, Purdue University, 1993.

45. B. Piper, Visually smooth interpolation with triangular Bézier patches, in Geometric Modeling: Algorithms and New Trends, G. E. Farin (ed.), SIAM Publications (Philadelphia), (1987), 221–234.

46. M.J.D. Powell, TOLMIN: A Fortran Package for Linearly Constrained Optimization Calculations, Dept. of Applied Mathematics and Theoretical Physics, University of Cambridge, 1989.

47. E. Quak, and L. L. Schumaker, Calculation of the energy of a piecewise polynomial surface, in Algorithms for Approximation II, M. G. Cox and J. C. Mason (eds.), Chapman & Hall (London), (1990), 134–143.

48. U. Reif, Neue Aspekte in der Theorie der Freiformflächen beliebiger Topologie, dissertation, Univ. Stuttgart, 1993.

49. M. Sabin, Non-Rectangular Surface Patches Suitable for Inclusion in a B-Spline Surface, in EUROGRAPHICS’83, 57–69.

50. R. Sarraga, G^{1} interpolation of generally unrestricted cubic Bézier curves, Comput. Aided Geom. Design 4 (1987), 23–39.

51. R. Sarraga, Errata: G^{1} interpolation of generally unrestricted cubic Bézier

curves, Comput. Aided Geom. Design 6 (1989), 167–172.

52. R. Sarraga, Computer modelling of surfaces with arbitrary shapes, Comp. Graphics and Applics. 10(2) (1990), 67–77.

53. R. Schmidt, Eine Methode zur Konstruktion von C^{1}-Flächen zur Interpolation unregelmässig verteilter Daten, in Multivariate Approximation Theory II, W. Schempp and K. Zeller (eds.), Birkhäuser (Basel), (1982), 343–361.

54. L. L. Schumaker, On the dimension of spaces of piecewise polynomials in two variables, in Multivariate Approximation Theory, W. Schempp and K. Zeller (eds.), Birkhäuser (Basel), (1979), 396–412.

55. L. L. Schumaker, Recent progress on multivariate splines, in Mathematics of Finite Elements and Applications VII, J. Whiteman (ed.), Academic Press (London), (1991), 535–562.

56. M. K. Seager, Sparse Linear Algebra Package Version 2.0, Lawrence Livermore National Laboratory, 1989.

57. A. H. Sherman, Algorithm 533. NSPIV, a Fortran subroutine for sparse Gaussian elimination with partial pivoting, ACM Trans. Math. Software 4 (1978), 391–398.

58. L. Shirman, and C. Séquin, Local surface interpolation with Bézier patches, Comput. Aided Geom. Design 4 (1987), 279–295.

59. L. Shirman, and C. Séquin, Local surface interpolation with Bézier patches: errata and improvements, Comput. Aided Geom. Design 8 (1991), 217–221.

60. D. J. T. Storrey, and A. A. Ball, Design of an n-sided patch from Hermite boundary data, Comput. Aided Geom. Design 6 (1989), 111-120.

61. S. Timoshenko, and S. Woinowsky-Krieger, Theory of Plates and Shells, Mc-Graw Hill, 1959.

62. J. J. Van Wijk, Bicubic patches for approximating non-rectangular controlpoint meshes, Computer-Aided Design 3 (1986), 1–13.

63. T. Varady, Survey and new results in n-sided patch generation, in Algorithms for Approximation II, M. G. Cox and J. C. Mason (eds.), Chapman & Hall (London), (1990), 203–235.

64. Z. Zlatev, J. Wasniewski, and K. Schaumburg, Y12M. Solution of large and sparse systems of linear algebraic equations, Lecture Notes in Computer Science 121, Springer Verlag, Berlin, 1981.

Gregory E. Fasshauer and Larry L. Schumaker Department of Mathematics Vanderbilt University, Nashville, TN 37240

fass@math.vanderbilt.edu s@mars.cas.vanderbilt.edu

The second author was supported in part by NSF Grant DMS - 8902331.