An Examination of New Paradigms for Spline Approximations

Lavery splines are examined in the univariate and bivariate cases. In both instances relaxation based algorithms for approximate calculation of Lavery splines are proposed. Following previous work Gilsinn, et al. [7] addressing the bivariate case, a rotationally invariant functional is assumed. The version of bivariate splines proposed in this paper also aims at irregularly spaced data and uses Hseih-Clough-Tocher elements based on the triangulated irregular network (TIN) concept. In this paper, the univariate case, however, is investigated in greater detail so as to further the understanding of the bivariate case.

ring to approximations by regular rectangular grids. They adapt naturally to disparate data densities. Importantly, they support strategic selection of data points to be used as the support points for a TIN surface.
Conceptually, TINs define which points may be considered "neighbors" of each other. This proximity concept based on a direct neighbor criterion differs from concepts based on a fixed distance cutoff in that it adjusts automatically to differing data densities such as in data collected by ground based LADAR. Figure 3 exemplifies the dramatic differences in densities encountered in ground-based LADAR scans. Also, because of the neighbor relation provided by a TIN, the utility provided by that construct is not restricted to surface generation, but extends to data editing and analysis.
The TIN surface, as defined according to Fig. 1, is C 0 , that is, it is continuous but usually not smooth: it is typically not differentiable along triangle edges, since the plane of any of the spatial triangles constituting the surface generally differs from the planes of its adjacent triangles. The TIN concept of triangulation-based surface generation is frequently understood to imply the use of such planar "elements" resulting in a nonsmooth, piecewise linear TIN surface. It should be emphasized, however, that the TIN concept also supports the use of non planar, that is, curved elements for a smooth or C 1 surface. Planar elements offer some advantages besides simplicity and a certain robustness to be discussed below. Visualization still requires piecewise linear representations in order to delineate hidden surfaces. Volume computations similarly tend to be either grid-based or based on piecewise linear surfaces for ease of computation. In those cases, a smooth surface would have to be discretized or approximated by a piecewise linear surface. This then begs the question, why not use planar elements in the first place?
Without any doubt, however, there are many instances, in which a smooth surface representation would provide definite advantages: • More accuracy could be achieved without increasing memory requirements because curved elements encode more geometric information than planar elements. • Clusters of coplanar points are avoided when sampling the surface at a discrete set of foot prints such as regular grid points. • When simulating movement over a terrain surface, a smooth ride is preferable. A ride over a piecewise linear surface is necessarily bumpy at transitions from one triangle to another. • In many applications, several point clouds which were collected with reference to different coordinate systems need to be "registered", that is, combined within a common coordinate system. Some common registration methods, such as point-to-surface iterative closest point (ICP) proposed by Besl [2], require minimizing the deviations of points in one point cloud from the surface representation of another point cloud. This minimization process should work much better if that surface were smooth.
In computer aided design (CAD), design-defined objects are routinely represented by surfaces that are C 1 or even C 2 . Terrain representation, however, poses a  different challenge: The location of "crease lines" or "break lines", along which the actual terrain surfaces are not differentiable, are usually not known ahead of time, whereas in a CAD environment break lines tend to be specified as part of the design. Unspecified break lines, on the other hand, along with actual verticalities, tend to give rise to spurious oscillation and, what may be called, "Gibbs phenomena" in analogy to the phenomenon well known from the theory of Fourier series. TIN surfaces based on planar elements, on the other hand, are more robust, and can be constructed so as to automatically represent and, if necessary, report break lines. Susceptibility to spurious oscillations and Gibbs phenomena are one of the reasons why the terrain modeling community has been slow to accept smooth surfaces. There has thus been a long quest for "non-oscillatory splines" which would obviate this pesky conundrum.
In 1994, Lavery [11,12,13] (see also Gilsinn and Lavery [8,9,10]) proposed successful paradigms for univariate as well as bivariate non-oscillatory splines, which could be used for representing 2D or 3D data sets, respectively. Lavery introduced the term "L 1 splines" for his brand of nonoscillatory splines. The term L 1 splines is, however, frequently misinterpreted as minimizing an L 1 measure-of-fit when approximating a point set by, say, classical splines. For this reason, we prefer the term "Lavery splines".
Classical splines are characterized by their minimizing energy functionals. Lavery splines, on the other hand, minimize different functionals. In the bivariate case, especially, the computational effort of minimizing these functionals, however, exceeds the effort required by the classical approach by an order of magnitude. For the bivariate case we have, therefore, considered an approximation to the calculation of Lavery splines, along with a prior modification of the functional proposed by Lavery for the bivariate case. This modified bivariate functional is still an extension of Lavery's functional for univariate functions, but it extends a different aspect of the latter. It also is invariant under planar rotations of the coordinate system, which Lavery's bivariate functional is not. We will present preliminary results for both univariate and bivariate nonoscillatory splines.
The main thrust of this paper, however, remains to illustrate the utility and performance of bivariate nonoscillatory splines, building on the previous study by Gilsinn et al. [7] , and also on an early terrain modeling study by Mandel et al. [18]. Both studies employ TIN techniques in conjunction with the "reduced Hsieh-Clough-Tocher (rHCT)" element, the approach was pioneered by Lawson [14,15] and is followed in this work. It involves specifying elevations z i and two slopes z ix , z iy at each triangle corner or vertex v i = (x i , y i ) in a triangulation, and filling in, at each triangle, the thereby defined rHCT elements, results in a smooth surface over an entire TIN.
The classical spline approach to, say, interpolating the elevations at vertices of a triangulation would be to prescribe the elevations z i , and select the remaining parameters, namely the partial slopes z ix , z iy , so as to minimize some energy functional. This was indeed the approach taken in the work by Mandel et al. [18] The task to be accomplished there was to represent terrain given by digitized contour lines. The elevations at the data points were thus determined by the elevation of the contour line on which the data point was located. In the end, elevations were to be evaluated at the points of a square 900 × 900 grid. At the time, -1985 -visualization procedures that are now routine were not yet commonly available. A display of the results had to wait another year. So assessing the quality of the representation was restricted to manual spot checking. It thus took until a day before a scheduled presentation that, to our dismay, huge oscillations were detected. Fortunately, it turned out that it was not due to a problem with the method but was caused by an error in the data that had been provided: the elevation of a single contour line had been tagged 100 feet too high, causing the disruption. In this instance, the oscillations proved to be actually beneficial in that they uncovered data errors.
In Sec. 2 we discuss aspects of univariate Lavery splines in order to shed light on related issues for bivariate splines, which are addressed in Sec. 3. The algorithm proposed here for bivariate non-oscillatory splines requires solving large sparse systems of linear equations. Experience was gathered concerning the performance of the well known Gauss-Seidel algorithm.

Univariate Spline Interpolation
While the emphasis of this work is on bivariate splines for surface generation, an examination of univariate splines is taken up in order to highlight some of the issues pertaining to splines in general. In the univariate case, cubic "spline functions" are most commonly used and are considered here. They form a linear space of piecewise cubic C 1 functions f (x) defined locally over intervals between "knots" that is, they consist of cubic polynomials Adjacent cubic polynomials are required to assume the same values y i at common interior knots, This ensures continuity of the complete spline function f (x) over the entire interval In addition, the polynomials are to assume the same slopes The spline functions f (x) are thus continuously differentiable, that is, they belong to class C 1 . In what follows, the linear spaces F ′, F ″ of first and second derivatives of spline functions are also considered, in spite of the fact that, at common knots, the second derivatives of adjacent cubic polynomials may not agree, so that the spline functions f (x) ∈ F are generally not twice differentiable at such knots. However, they are twice differentiable everywhere but on this set of measure zero. For the purposes of integration below, it does not matter that the function f ″(x) may not be defined for those arguments. Each of the constituent cubic polynomials f i (x) is uniquely determined by the values y i-1 , y i at the knots x i-1 , x i and the slopes m i-1 , m i at those locations (Fig. 4), in fact, the polynomial is linear in these parameters.
The entire function f (x) is thus uniquely determined by its values and slopes at the knots, and it, too, depends linearly on these parameters, so that the space F is isomorphic to the 2(n+1)-dimensional vector space of values y i and slopes m i , i = 0, ..., n.
We now turn to the task of interpolation. Here the values y i at the knots x i are fixed and specified. Given a particular specification of values y i , a corresponding "interpolating spline function" depends only on the parameters m i : Collectively, these functions define affine manifolds or flats, S, S ′, S ″, within the linear spaces F, F ′, F ″, respectively. That is, if then it follows that correspondingly The question then becomes, how to select slopes m i so as to achieve a "satisfactory" interpolation. That selection is generally made by minimizing a functional on the affine space S defined by an integral. In this work, we reserve the term "spline" -as opposed to "spline function", for the results of such a minimization.

Paradigms for Univariate Splines
"Classical splines" are uniquely defined as those interpolating cubic spline functions which are (i) C 2 , that is, twice differentiable, such that holds at interior knots, and for which (ii) the second derivative vanishes (1) , 0,1,..., , 2 at the two exterior knots. Holladay proved early on (see Ahlberg [1]) that classical splines are also defined as the unique minimizers of the "thin beam energy" (3) over the affine space S of all C 1 interpolating spline functions.
Condition (2) is familiar to structural engineers as the vanishing of the second derivative at the "free end" of a beam. It is, however, remarkable that this energy minimization enforces also a higher level of compatibility across knots (1) so that the minimizing C 1 spline functions are, in fact, in class C 2 . On the other hand, this stiffness contributes to the tendency of classical splines to produce spurious oscillations, Gibbs phenomena and undesirable inflections.
Several attempts have been made to avoid these problems. Taking a clue from mechanics, Schweikert [20] and Cline [3] have introduced "tension splines", where the arclength of the spline function is made part of the defining minimization. Reinsch (see Stoer and Bulirsch [21]) moved to "exponential splines" by adding a multiple of the square of first derivatives to the integrand in (3). Those efforts were only partially successful. Drawbacks include the need to specify an additional parameter in order to balance conflicting minimization requirements, and the fact that these techniques are not readily generalized to the bivariate case.
Lavery splines, on the other hand, appear to avoid such shortcomings. In Figs. 5 and 6 we compare Lavery splines against classic splines. We are particularly impressed with the performance of Lavery splines for the example in Fig. 8, as compared to the example in Fig. 7. What is the secret of Lavery splines? How are they defined as opposed to classical splines?
Lavery defines his nonoscillatory splines, in essence, as minimizing the integral of the absolute value rather than the square of the second derivative of a spline function:

Expressing Holladay and Lavery Integrals
For the purposes of this paper, we will refer to the integrals (3) and (4) as the "Holladay integral" and the "Lavery integral", respectively. The goal of this section is to derive expressions for the values of these integrals in terms of the slopes m i , i = 1, ..., n, at the knots of the spline function f (x) under consideration. Both are the sums of the corresponding integrals of the second derivatives of the individual cubic polynomials f i (x), which constitute the spline function f (x).
As pointed out in Sec. 2.1, each such cubic polynomial is uniquely defined by its end points (x i-1 , y i-1 ), (x i , y i ) and its end slopes m i-1 , m i . The various formulas describing this polynomial are commonly referred to as "Hermite" formulas. For the versions used here, we introduce the quantities where M i represents the slope of the straight line between end points. Instead of referring to the variable x directly, the following formulas are in terms of the weights (5) where λ i + µ i = 1, and λ i , µ i ≥ 0, for x in the interval [x i-1 , x i ] . Such weights are often referred to as "barycentric coordinates" or, in the bivariate case, as "triangle coordinates". With these conventions, we find, for instance, An alternate expression for the first derivative is readily derived: Here the quantity (8) vanishes if and only if the polynomial f (x) has degree less than three. Because by (5) - . D i < 0 indicates that the function is concave up to its inflection point and convex thereafter. Conversely, D i > 0 indicates that convexity precedes concavity in the direction of the x-axis. The quantity D i will play a major role in what follows. The same is true for the quantities U i , V i in the expression (9) where (10) The coefficients U i , V i thus relate to the two-sided second derivatives of f i (x) at knots x i . Note that and also that, in view of (5), An expression for the thin beam energy of the individual polynomial f i (x) is now readily derived: The full Holladay integral is the inhomogeneous quadratic function of the variables m i , arrived at by adding the energies E(f i ) of all partial functions f i (x). Introducing the vector of slopes where and Minimizing this expression for the Holladay energy integral is to solve the linear system of equations (12) 2Hm  The reader should note that the system (12) for classical splines differs from the linear system mostly offered in the current literature. There the second order differentiability of the classical splines is already assumed and the spline functions are formulated in terms of those second order derivatives, n i = f ″(x i ), i = 0, ..., n. However for weighted classical splines, to be encountered in Sec. 2.4.1, second order differentiability no longer holds, and a weighted version of the linear system (12) needs to be considered.

Expressing the Lavery Integral
As to the integral of the absolute value of the second derivative, it is readily available in the case that the derivative does not change sign inside the subinterval: This includes the case in which the polynomial is of lesser degree than cubic, and thus has constant second derivatives. This case is signaled by the vanishing of the quantity D i introduced earlier in (8).
The function f ″(x), however, is a linear function in x and, unless constant, changes sign at some location x, which also marks the location of the inflection point of f i (x). Suppose this location falls into the interior of the subinterval: Then the integral has to be calculated as the sum of two integrals of linear functions: The integrals between the absolute value bars are of opposite signs, so that the total integral can be written as the absolute value of a difference of integrals These integrals can be separately evaluated as differences of slopes. Let denote the inflection slope, then -in the case of an interior inflection - The quantities U i , V i depend linearly on the two slopes m i-1 , m i . Conversely, those slopes can be expressed in The next step is to express the inflection slope in terms of U i , V i . To this end, we express the inflection argument by its barycentric weights: parabolic or linear. Thus D i = 0 leads back to the previous case of no sign changes by the second derivative.
The weights λˆi, µˆi may now be inserted into the expression (7) for the derivative f ′ i (x). This gives Substituting for m i-1 , m i according to (14) yields as well as (16) in view of (19). In terms of the slopes m i-1 , m i : The above expressions for e(f i ) are also valid if the inflections occur at the ends We are now ready to examine the full Lavery integral. At first blush, all that remains to be done is to sum over the partial integrals e(f i ) in their various forms. We will show, however, that many terms of the expressions (13) cancel each other out as these partial integrals are added together. To this end, we distinguish five separate kinds of polynomials f i (x) depending on their behavior in the interior of the interval between its knots, The last two categories are the ones with an interior inflection point x and inflection slope m .
The interior inflection points considered so far may not be the only inflection points of the spline function f (x). Inflections occur also at knots x i , 0 < i < n if a concave or convex-concave polynomial is followed by a convex or convex-concave polynomial and, analogously, if a convex or concave-convex polynomial is followed by a concave or concave-convex polynomial. We refer to such knots as "inflection knots". To make matters more complicated, however, inflection may occur along an entire stretch of consecutive linear polynomials of equal slope, the inflection slope in this case, provided there are adjacent nonlinear polynomials at both ends of such a stretch exhibiting the same convexity/concavity pattern that would cause an inflection at a knot. In this case, we choose an arbitrary knot, say, the leftmost one in the linear stretch, as the inflection knot representing the inflection.
Clearly, slopes at interior knots that are not inflections cancel out as the expressions (13) (Note that the indices l of m l do not refer to the interval in which they are located).

Properties of Lavery Splines
In this section, we gather some information about Lavery Splines, namely, interpolating spline functions in the affine space S which minimize their respective Lavery integrals. A first general observation concerns the convexity of the Lavery integral.
The Holladay and Lavery integrals (3) and (4)  The quadratic function E( f ) representing the Holladay integral can be shown to be positive definite and, therefore, strictly convex. Its minimum is unique on the affine manifold S interpolating spline functions. The restriction to S is, of course, necessary as the value of E( f ) would not change if the spline function f (x) were modified by adding a linear function in x. Adding a non-zero linear function to the function f (x) would not preserve its interpolation property.

Convexity
Next, we establish the convexity of the Lavery integral. It may be viewed as the extension of the L 1 vector norm to a norm on the linear space F″ of second derivatives of spline functions: The following generic seminorm properties are easily verified for the piecewise linear functions in F″, along with the triangle inequality, Suppose the two spline functions f (1) , f (2) are actually two interpolating spline functions and, therefore, in the affine space S. Then their mean is again in S and, from the triangle inequality, In terms of Lavery integrals, which establishes convexity. This leads immediately to Observation B: Any positive linear combination of Lavery splines for the same interpolation problem, -in particular, their mean -, is again a Lavery spline for this problem.

Uniqueness of Inflections
Contrary to the Holladay functional, which is strictly convex, the Lavery integral is not. As a result, uniqueness does not follow and, in fact, does not hold, as an example in Sec. 2.3.4 will show. Such minima of a convex function, however, must form a convex set.
Note that, in general, If both f (1) and f (2) are Lavery splines for the same interpolation problem, then so is their mean, and all three functions return the same optimal value for their Lavery integrals. Thus equality holds in the above relation. This implies that both f 1 ″(x) and f 2 ″(x) have the same sign pattern: This can be rephrased as , , 2 Observation C: Two Lavery splines for the same optimization problem share essentially the same inflections: if one of them has an inflection point at x = x, then so has the other unless it is linear at this point.

Free Ends of Lavery Splines
Here we will examine the free ends of Lavery splines, in particular, the cubic polynomial f 1 (x) and the corresponding first summand e(f 1 ) of the Lavery integral e(f ) . As the end slope m 0 may vary freely, it must optimize e(f 1 ) while keeping the slope m 1 fixed. This fact determines the behavior of Lavery splines at free ends. As seen in the previous section, the first summand e(f 1 ) is a convex function in the variables m 0 , m 1 . If the slope m 1 is held fixed, e(f 1 )(m 0 ) is convex as a function in m 0 alone. If the fixed slope equals the straight-line slope, then m 0 = m 1 = M 1 obviously represents the optimal value for m 0 , since f 1 (x) in that case is a straight line with e(f 1 ) = 0. We suppose, therefore, that m 1 ≠ M 1 , and we examine the case, that f 1 (x) has an inflection x in the interval between the first two knots, Consider the function where m as well as U 1 , V 1 also depend on the variable m 0 . Clearly, the absolute value of e ∼ (f 1 ) is given by e(f 1 ). Note that in view of (10), and Solving the quadratic equation for λˆ, and taking into account that 0 ≤ λˆ≤ 1, yields The value of m 0 for which this value for λˆ1 is realized can be inferred from the general definition (15) of an inflection, giving the barycentric coordinate λˆ1 in terms of the slopes m 0 , m 1 , as follows: which -for the particular value (17)  This implies that the value (18) for m 0 is also a locally unique stationary value of e(f 1 )(m 0 ). In view of the convexity of this function, it is also its minimizer. At the last end we differentiate (19) and find Symmetrically, we thus have µˆn = λˆ1, λˆn = µˆ1. Observation D enables us to determine universal values for partial Lavery integrals at the end-intervals of Lavery splines. By (15), (16), (18), and again by (20), which implies we find Substituting for m 0 and m n , respectively in D 1 = 6M 1 -3m 0 -3m 1 and D n = 6M n -3m n-1 -3m n , giving rise to Observation E: A necessary, but far from sufficient, condition for a spline function f to be a Lavery spline is that

Examples of Non-unique Lavery Splines
In this section, we present an example in which the Lavery splines are not unique. Consider the three points   Note that both partial Lavery integrals are end-integrals. For an interpolating spline function f (x) to be a Lavery spline it will be necessary by Observation E that so that whereas either m 1 < -1 or m 1 > 1 would imply |1 -m 1 | + |-1 -m 1 | = 2 ≥ 2 |m 1 |.
Thus condition (22) characterizes all Lavery splines for the example. Consider now any slope m 1 from -1 through +1. For m 1 = -1, we have by (18), (20), and in view of M 1 = 1, M 2 = -1: Symmetrically, we find for m 1 = +1 that These slopes, respectively, determine the two extreme Lavery splines, because each partial function is at a free end, and is optimized according to Observation D in the previous section. The two resulting Lavery splines are extreme in that each has a straight line segment as a partial function.
Using Hermite's formula (6) and substituting for x, we find for the choice m 1 = -1, For m 1 = +1, the resulting Lavery spline is the symmetric image of the previous one: Both splines are shown in Figs. 10, 11. The self-symmetric spline from the choice m 1 = 0 is shown in Fig.  12. In the latter case, and This Lavery spline is the mean of the two splines with linear free ends. This is an instance of observation B about positive linear combinations of Lavery splines in Sec. 2.3.1.

Computing Univariate Lavery Splines
We now turn our attention to the computation of Lavery splines. The commonly used approach (see Lavery [11,12]) is to minimize the Lavery integral in discretized form, say, where the integrand is sampled in each subinterval [x i-1 , . 5 5 solving it are available, such as the Simplex Method or Interior Point Methods. The accuracy of the discretization increases with the number of sample points, but computational effort increases accordingly.
For that reason, and also to motivate an analogous approach in the bivariate case, we are proposing to minimize a different approximation to the Lavery integral, one that takes advantage of the ease of computation offered by energy minimization.
By the mean value theorem of integral calculus, there exist arguments u i such that we then propose to approximate the Lavery integral by the following Riemann sum: In contrast to the approximation (23) by discretization, this approximation does not offer the option of further refinement, unless the interpolation problem itself is changed by adding additional knots and ordinates. In a sense, it approximates the Lavery paradigm itself rather than the Lavery integral. We still use, however, the term "approximate Lavery integral" for our proposed approximation (24).

An Iterative Algorithm for Approximate Lavery Splines
For the purpose of computation, we rewrite the nonzero terms in (46) This expression suggests an iterative approach. Starting with the classic spline f (0) (x), a sequence of interpolating spline functions is generated in hopes to converge towards the approximate Lavery integral (24) ) vanishes as the corresponding weight is then not defined?
Simply ignoring such terms may prematurely lock in straight line segments between knots. What first comes to mind is to specify a limit ε > 0 and boost lower values of E(f i (l) ) to this level. A more diligent procedure might be to start with all weights at value 1, -the weight setting that yields the initial classical spline according to Holliday's observation -, and then progressively increase the use of weights. Such strategies remain to be explored.
In general, adding up partial Holliday integrals, each with weight, say, w i leads again to an expression of the energy of a physical structure: a collection of thin beams of different thicknesses given by w i , respectively, and welded together at knots. Minimizing this energy expression requires an adjustment to the linear system (12).
The weighted energy of the partial spline functions f i (x) is just the product of the straight Holliday integral and the respective weight: The total weighted energy is thus given by This means that in the expression (11) for E(f i ), the factor is affected as it is multiplied by w i . In other words, 1the substitution, transforms the linear system (12) for minimizing E(f ) into the linear system for minimizing E W (f ): ( ) ( ) d , 1, ..., .
Note that the weighted classic splines are, in general, not twice differentiable at the knots. However, due to the fact that the first and the last equations have common factors w o , w n , respectively, they are equivalent to the corresponding equations in the Holladay system (12): free ends thus have zero second derivatives in the weighted case, too. Again, this is to be expected from Physics. Note also that some commonly used methods for determining classic splines such as B-splines or linear systems formulated in terms of second derivatives at knots do not carry over to the weighted case. However, the above linear system is still "banded", and many excellent methods are known for its solution. To solve this system we used the venerable Gauss-Seidel method, not just for ease of programming, but also because it seems to work for bivariate weighted splines. Its advantage lies in the fact that the matrix of the linear system need not be changed and can be read, so to speak, in sequence. This is important for the very large systems likely to arise in the bivariate case. In the univariate case, the convergence behavior is well understood (See Varga [22]). Using an iterative method for solving the class of linear systems above will result in a two-tiered iteration procedure: an "outer" iteration, developing new sets of weights, and an "inner" iteration, solving the resulting linear system. Such procedures can be "balanced", that is, the inner iteration may be terminated at a lower accuracy level during the early stages of the outer iteration and may be carried to a higher level of accuracy as the outer method approaches convergence. This is an added advantage of an iterative method for solving the linear systems at hand.
Note finally, that the approximate Lavery method proposed here will definitely not converge to the optimal Lavery integral. This is because the approximate solutions are based on weighted classic splines, and therefore have vanishing second derivatives at the free ends. The second derivatives of Lavery splines, on the other hand, assume their inflections in interiors of the end interval (See Observation D).
Observation F: Unless a free end function of an approximate Lavery spline is linear, it disagrees with the corresponding end function of the true Lavery spline in that the latter has an inflection in the interior of at least one end interval, whereas the former assumes its inflection at the end knot.
The justification for introducing the approximate Lavery splines concept lies in its computational ease and the fact that it appears to retain the anti-oscillatory dynamic of the original Lavery concept. In fact, the examples in Figs. 5, 6, and 8 were calculated using the approximate algorithm outlined in this section.

Bivariate Spline Interpolation
The development of bivariate splines parallels that of univariate ones to a large extent. For both, the first step is to specify a linear space F of bivariate spline functions f (x, y), defined on a polygonally partitioned region Ω In view of our emphasis on the TIN framework of triangulations, however, we choose for our discussion a rectangular region Ω and an irregular set of "knots" in Ω, v i = (x i , y i ) ∈ Ω, i = 1, ..., n.
An interpolation problem results, if "elevations" z i are specified at the vertices v i = (x i , y i ).
The region Ω is then partitioned into triangles The vertices of any such triangle t k are among the specified knots, and those are the only knots in this triangle. In this fashion, a "triangulation" of the region Ω is obtained. Spline functions are then defined by installing partial functions  Institute of Standards and Technology   71   1  1  1  0  1  1  1  1  1   1  1  2  2  1  2  0  1  2  1  2  1  1  2  2  1  2   3  3  3  2  2  2  1  2  3  2  3  2  2  3  3 in their respective triangles t k of the triangulation. Following C. Lawson [14,15], we install above each triangle t k a reduced Hsieh-Clough-Tocher (rHTC) element. The reduced Hsieh-Clough-Tocher (rHCT) element, displayed in Fig. 14, provides a means for representing smooth TIN surfaces (see Lawson [14,15]). The rHCT element is defined as a function over a triangle in the x, y-plane, and is fully determined by the elevations and the partial slopes or derivatives at triangle corners. The triangle is divided into three subtriangles with their common vertex at the centroid of the three corners of the full triangle. In each of the three subtriangles the rHCT function is a bivariate cubic function. Together, these functions form a C 1 function on the full triangle (Fig. 13). This smooth rHCT function is furthermore constrained by the linear perpendicularity condition on the derivatives orthogonal to the outside boundary edges of the triangular element: along each such edge, these derivatives are required to be linear functions. It follows that the orthogonal derivatives along an edge are fully determined by their values at the ends of the edges. These values are inherent to both triangles adjacent to the edge so that the orthogonal derivatives coincide when calculated within each of the two triangles independently. The above condition thus provides the key property ensuring smoothness of a rHTC surface.
As will be seen in Sec. 3.1.2, the restriction inherent in the linear perpendicularity condition may be more severe than commonly expected.
To this end, "elevations" z i and "partial slopes" are prescribed at the knots v i . As the vertices of each triangle t k are among the specified knots, the corresponding elevations and partial slopes furnish the parameters necessary for defining the rHTC element f k (x, y). As pointed out above in the current section, each specification of elevations and partial slopes at the knots yields a C 1 function f (x, y), which also represents a smooth surface over the region of definition Ω. All such functions, -relating to the given triangulation of Ω -, constitute a linear space of "spline functions" and their respective derivative linear spaces If a specific set of elevations z i is to be interpolated, then those spline functions which meet these elevations constitute an affine space of "interpolating spline functions" and its derivative affine spaces S, S x , S y , S xx , S xy , S yy .
Note, again, that second derivatives are not fully defined for rHTC elements. Indeed, along edges separating two subtriangles of an rHTC element as well as along edges separating two triangles of the triangulation there may be two different values in contention depending from which side a point on such an edge is approached. For purposes of integration, this is not an issue since the discrepancies are restricted to a set of measure zero, so that the integrals considered below are still well defined.

Paradigms for Non-Oscillatory Bivariate Splines
For the purpose of interpolation, the elevations z i at knots are specified, and the task at hand is to determine the partial slopes zx i , zy i so as to arrive at a "best" interpolating spline function. Such an optimal spline function is then referred to as a "(bivariate) spline". Obviously, there are various kinds of splines depending on the choice of the linear space of spline functions and the choice of the optimization criterion.
Paralleling the classical univariate approach, "thin plate" minimization exemplifies the "classical" approach, and the resulting spline will be referred to as the "classical spline", the functional requiring minimization being (28) An approach to computing the surface energy formula has been announced by McClain and Witzgall [16] for rHCT element. A newer version of the report is currently being prepared by McClain, Witzgall, and Gilsinn [17].
The authors believe that, contrary to univariate classical splines, bivariate classical splines as defined above are, in general, not in class C 2 , that is, twice differentiable. Nevertheless they are "stiff" and thus subject to the dreaded spurious oscillations and Gibbs phenomena. This can be seen along edges of buildings in the urban data set of Baltimore, MD, given in Fig. 14. Lavery and Gilsinn [13] were able to address this problem by extending his paradigm for univariate nonoscil-latory splines to the bivariate case by minimizing the functional (29) dxdy.
Again, the squares of second derivatives in an energy integral, -here E ∼ (f ) as defined in (28) -, are replaced by absolute values. We will refer to the so defined splines also a "(bivariate) Lavery splines".
In view of the non-uniqueness result for univariate Lavery splines, the bivariate Lavery splines are expected not to be unique either, as will be discussed below. Therefore, the Lavery paradigm typically includes a regulatory term such as adding a small multiple of f ′(x, y) 2 . The resulting bivariate Lavery splines have produced oscillation-free interpolations and approximations in a large array of applications. They are computed by discretizations leading to very large linear programming problems.
For large surface generation tasks based on interpolation, the resulting computational effort may be prohibitive. Moreover, the bivariate Lavery integral (29) functional is not invariant under rotation of the coordinate system. The authors have therefore proposed an alternate extension of the univariate Lavery paradigm to the bivariate case: (30) Note, that this paradigm is indeed an extension of the univariate Lavery paradigm. Indeed, so that the integrand of the univariate Lavery integral may also be interpreted as the squareroot of an energy, -in this case, the elastic energy of a thin beam. Replacing that energy expression by that for a thin plate, yields the alternate extension.

Properties of the Alternative Bivariate Paradigm
The differential expression is readily seen to be rotation invariant, as is to be expected, given its physical meaning as the integrand for the representation (28) of the thin plate energy. It follows that every differential expression of the form Let f = f (x, y) , g = g(x, y) ∈ S be two interpolating spline functions. Invoking the general vector inequality ||u|| + ||v|| ≥ ||u + v||, we find so that and in terms of the spline functionals, Since S is an affine space, the mean (f + g)/2 is again in S. Based on the above inequality, we thus have Observation H: The alternate Bivariate functional (30) is a convex functional. In particular, as expressed in terms of the partial slopes, e ∼ (f ) = e ∼ (zx 1 , zy 1 , ..., zx n , zy n ), it is a convex function in terms of these partial slopes.

A Related Non-Uniqueness Result
We do not know at this point, whether the convexity is moreover strict. We expect that it is not. Strict convexity would imply uniqueness of minima and, again, we do not believe that this is the case. However, an example for such non-uniqueness is eluding us for our specific space of rHTC spline functions.
For different bivariate functions, however, it is possible to establish non-uniqueness of a minimum. Consider a function s(x, y) , -1 ≤ x ≤ +1, 0 ≤ y 1 representing a rectangular strip bent only in the direction of the x-axis. In the direction of the y-axis, the function s(x, y) is constant (Fig. 15), in particular, for all y, s(-1, y) = -1, s(0, y) = 0, s(+1, y) = -1.
The function s(x, y) thus interpolates the values -1 , 0, -1 along the contour lines in the direction of the y-axis at the arguments x = -1, 0, +1. Given any Lavery spline f (x) considered in Sec. 2.3.3 as interpolating the example (21). Then it is readily seen, that any function minimizes both bivariate functionals: e(s), e ∼ (s). For functions corresponding to strips bent in only the strip direction, these bivariate functionals have thus multiple minima. This makes a strong case for non-uniqueness bivariate splines minimizing the functionals (29) and (30). However, we were not able to pinpoint nonuniqueness for the rHTC based spline functions considered in this section. In particular, we were surprised to realize that the function s(x, y) cannot be generated using rHTC elements. As this fact points to an important limitation of using rHTC elements, we formulate the Volume 111, Number 2, March-April 2006 Observation J: The surface function of a rectangular strip bent only in the strip direction cannot be represented by means of rHTC elements.

Computing Bivariate Splines Based on an Alternate Paradigm
The proposed alternate functional (53) is still expensive to minimize. In analogy to the procedure described in Sec. 2.4.1, we propose to approximate that functional as follows. Let For the exploratory implementation reported here, the above minimization of a quadratic function in the partial slopes zx i k,1 , zy i k,1 , zx i k,2 , zy i k,2 , zx i k, 3 , zy i k, 3 .
is again formulated as solving a linear system of equations in these variables, and Gauss-Seidel iteration provided a workable option. Just as discussed in the univariate case, there are two iterative processes, one "on top of" of the other. The "inner iteration" aims to solve the linear system of equations in order to solve the intermediate minimization problem (31), while the "outer" iteration sequentially creates spline functions expected to converge toward a limit that approximates the alternate Lavery spline defined by the minimization of the alternate Lavery functional (30).
This approach raises the same issues as its parallel in the univariate case. The above iterative approximation based on consecutive minimization of weighted sums of integrals (31) cannot be expected to converge exactly to the alternate bivariate Lavery spline. The argument relies on Observation F pertaining to the univariate case. Each iterate f (l) (x, y) minimizes the energy of a physical structure consisting of thin plates of various thicknesses, namely the weights. Any minimizing shape of that structure, however, is known to have some vanishing second derivatives at its boundary. Univariate Lavery splines do not have this property, and we do not expect bivariate Lavery splines to have this property, either.
Provisions need to be made should one of the integrals E ∼ (f k (l-1) ) vanish. This is a serious problem because very large or "infinite" weights will cause certain rHTC elements to remain stuck in linear form. In addition, limits of bivariate Lavery splines, in the original as well as the alternate formulation, may not be unique (see Sec. 3.1). These issues call for regularization procedures which are adequate for particular kinds of applications.
Indeed, spurious spikes were observed for the applications demonstrated here. The problem was solved by interspersing an averaging step, where the partial slopes were replaced by averages of neighboring slopes.
While our experiences with approximating alternate bivariate Lavery splines (Fig. 16) were encouraging, in particular, for very large data sets, much work needs to be done.
The accuracy performance of the Gauss-Seidel iteration is a major issue. Simple examples show that initial steps are not even contracting. When does contraction start? Does it keep contracting after that? Much research has been aimed at these questions. These analyses of the convergence properties of the Gauss- principle based on the integral of the square root of the sum of squares of the second partial derivatives of the bivariate spline. Whereas the Lavery integral in the bivariate case is not rotationally invariant, the alternative principle is. The alternative bivariate functional is shown to be a convex functional and, when expressed in terms of the partial slopes at the vertices of an underlying element, it is a convex function of these slopes. We have also shown that a rectangular strip bent only in the strip direction cannot be represented by means of a triangular rHTC element.
Analogous to the alternative weighted least squares algorithm introduced to compute the approximate univariate splines we introduce an alternative weighted least squares algorithm in the bivariate case. Although we have not introduced a convergence analysis of the alternative algorithms in this paper, computational experience has shown that the iterative, Gauss-Seidel based, algorithm used provides qualitatively satisfactory results in only a few iterations of the main loop.