Introduction

Interpolation of data points in ℝ 3 by smooth surface is a fundamental problem in applied mathematics which finds applications in a variety of fields such as medicine, architecture, archeology, computer graphics and animation, bioinformatics, scientific visualization, and more. In general the problem can be formulated as follows: Given a set of points (x i , y i , z i ) ∈ ℝ 3 , i = 1, ..., n, find a bivariate function F(x, y) defined in a certain domain D containing points V i = (x i , y i ), such that F possesses continuous partial derivatives up to a given order, and F(x i , y i ) = z i .

Various methods and approaches for solving this problem were proposed and discussed, see, e.g., the surveys [1,2,3], and also [4,5]. A standard approach to solve the problem consists of two steps, see [1]:

1. Construct a triangulation T = T(V 1 , ..., V N ); 2. For every triangle in construct a surface which interpolates the data in the three vertices.

The interpolation surface constructed in Step 2 is usually polynomial or piecewise polynomial. Typically, the patches are computed with a priori prescribed normal vectors at the data points G 1 or G 2 smoothness of the resulting surface is achieved either by increasing the degree of the patches, or by the so called splitting [6] in which for each triangle in T a macro-patch consisting of a fixed number of Bézier sub-patches is constructed. Splitting allows to keep the degree of the Bézier patches low by increasing the degrees of freedom. In practice using patches of least degree and splitting is preferable since it is computationally simple and efficient.

Shirman and Séquin [7,8] construct a G 1 smooth surface consisting of quartic triangular Bézier surfaces. Their method assumes that the normal vectors at the data points are given as part of the input. Shirman and Séquin construct a smooth cubic curve network defined on the edges of T, first, and then degree elevate it to quartic. Next, they apply splitting where for each triangle in T a macro-patch consisting of three quartic Bézier sub-patches is constructed. To compute the inner Bézier control points closest to the boundary of the macro-patch, Shirman and Séquin use a method proposed by Chiyokura and Kimura [9,10]. The interpolation surfaces constructed by Shirman and Séquin's algorithm often suffer from unwanted bulges, tilts, and shears as pointed out by the authors in [11] and more recently by Hettinga and Kosinka in [12].

Nielson [13] proposed a three-steps method for solving the interpolation problem as follows:

Step 1. Triangulation. Construct a triangulation T of V i , i = 1, ..., n. The domain D is the union of all triangles in T.

Step 2. Minimum norm network. The interpolant F and its first order partial derivatives are defined on the edges of T to satisfy an extremal property. The resulting MNN is a cubic curve network, i. e. on every edge of T it is a cubic polynomial.

Step 3. Interpolation surface. The MNN obtained is extended to F by an appropriate blending method based on convex combination schemes. Nielson's interpolant F is a rational function on every triangle in T.

Andersson et al. [14] focused on Step 2 of the above method, namely the construction of the MNN. The authors gave a new proof of Nielson's result by using a different approach. They constructed a system of simple linear curve networks called basic curve networks and then represented the second derivative of the MNN as a linear combination of these basic curve networks. The new approach allows to consider and handle the case where the data are convex and we seek a convex interpolant. Andersson et al. formulate the corresponding extremal constrained interpolation problem of finding a minimum norm network that is convex along the edges of the triangulation. The extremal network is characterized as a solution to a nonlinear system of equations.

Vlachkova and Radev [15] proposed an algorithm for interpolation of data in ℝ 3 which improves on Shirman and Séquin's approach in the following way. First, they use Nielson's MNN which is degree elevated to quartic curve network. Second, they extend it to a smooth surface consisting of quartic triangular Bézier patches by applying different strategy for computation of the control points. A significant advantage of Nielson's method is that the normal vectors at the data points are obtained through the computation of the MNN. Moreover, during the computation of control points closest to the boundary of a macro-patch, Vlachkova and Radev [15] adopt additional criteria so that to avoid unwanted distortions and twists which appear in surfaces constructed by Shirman and Séquin's method. As a result, the quality of the resulting surfaces is improved.

Vlachkova [16] extended Nielson's MNN to minimum L p -norm networks for 1 < p < ∞. The most important conclusions in [16] are as follows:

• the results allow the efficient computation of the minimum L p -norm networks for 1 < p < ∞; • the normal vectors at the data points are obtained simultaneously with the computation of the minimum L p -norm networks;

• the minimum L p -norm networks are obtained through a global optimization which improves their shape;

• the results allows the construction of optimized polynomial minimum norm networks of a priori given degree. In this paper, we consider the following two curve networks: (i) Nielson's MNN which is degree elevated to quartic; (ii) the minimum L p -norm network for p = 3/2 which is slightly modified to quartic. Using the algorithm proposed in [15], we construct the corresponding two interpolation surfaces consisting of quartic triangular Bézier patches. Our goal is to evaluate and compare the quality and the shape of these surfaces. We have chosen the following criteria for comparison:

(i) the highlight-line algorithm [17];

(ii) the color plot of the Gaussian curvature;

(iii) the maximum distance between the function sampled at the data points and the corresponding interpolant.

Our work and contributions presented here are in the field of experimental algorithmics. We share and comment on the observations from the experiments performed, which will help to further optimize the quality and the shape of the interpolation surfaces generated with the algorithm proposed in [15].

Preliminaries and related work

Let n ≥ 3 be an integer and P i : = (x i , y i , z i ), i = 1, ..., n be different points in ℝ 3 . We call this set of points data. We assume that the projections V i = (x i , y i ) of the data Oxy into the plane are different and non-collinear. Stressing the fact that V i are in a general position, although not necessarily irregularly placed, we call such data scattered. Definition 1. A collection of non-overlapping, non-degenerate triangles in is a triangulation of the points V i , i = 1, ..., n, if the set of the vertices of the triangles coincides with the set of the points V i , i = 1, ..., n.

Hereafter we assume that a triangulation T of the points V i , i = 1, ..., n, is given and fixed. Furthermore, for the sake of simplicity, we assume that the domain D formed by the union of the triangles in T is connected. In general D is a collection of polygons with holes. The set of the edges of the triangles in T is denoted by E. If there is an edge between V i and V j in E, it will be referred to by e ij or simply by e if no ambiguity arises. Definition 2. A curve network is a collection of real-valued univariate functions {f e } e∈E defined on the edges in E.

With any real-valued bivariate function defined on we naturally associate the curve network defined as the restriction of on the edges in , i. e. for , The smoothness of the interpolation curve network F ∈ C p (E) geometrically means that at each point P i there is a tangent plane to F, where a plane is tangent to the curve network at a point P i if it contains the tangent vectors at P i of the curves incident to P i .

Inner product and L p -norm are defined in C p (E) by where F ∈ C p (E) and G: = {g e } e∈E ∈ C p (E). We denote the networks of the second derivative of F by F″: = {f e ″} e∈E and consider the following extremal problem:

For i = 1, ..., n let m i denote the degree of the vertex V i , i. e. the number of the edges in E incident to V i . In [14] for any pair of indices is, such that i = 1, ..., n and s = 1, ..., m i -2, a basic curve network B is is defined on E. The basic curve networks are linear curve networks of minimal support.

Let q be the conjugate of p, i.e. 1/p + 1/q = 1. In [16] a full characterization of the solution F * to the extremal problem (P p ) for 1 < p < ∞ was made. It was shown that the curve network F * ∈ C p (E) solves problem (P p ) for 1 < p < ∞ if and only if where

The coefficients α is are obtained as the unique solution to a system of equations which is nonlinear except in the case p = 2 where it is linear. We note that Nielson's MNN is obtained for p = 2. We also note that in the case where q is an even number then the corresponding minimum L p -norm network is a polynomial curve network of degree q + 1. In the case where q is an odd number then on every edge of the triangulation the corresponding minimum L p -norm network is either a polynomial of degree q + 1, or a spline of degree q + 1 with one knot. Now we briefly discuss the constrained interpolation problem of finding a minimum norm network which is convex along the edges of T. We recall that this problem was set and solved in [14].

For a given triangulation, there is a unique continuous function L : D → ℝ 1 that is linear inside each of the triangles of T and interpolates the data. Definition 3. Scattered data D in are convex if there exists a triangulation T of V i , i = 1, ..., n, such that the corresponding function L is convex. The data are strictly convex if they are convex and the gradient of L has a jump discontinuity across each edge inside D.

We introduce the class of smooth interpolation edge convex curve networks and consider the following constrained extremal interpolation problem

In [14] it was shown that in the case of strictly convex data problem has a unique solution (the edge convex MNN) such that where (x) + : = max(x, 0), and α is ∈ ℝ. The coefficients α is are obtained as the unique solution to a nonlinear system of equations. The solution F ̂* on each edge in T is either a convex cubic polynomial, or a convex cubic spline with one knot consisting of a linear function and a convex cubic polynomial.

In this paper we consider the unconstrained problem (P p ) for p = 3/2 and also the constrained problem (P ̂). The solution of (P p ) for p = 3/2 (hence q =3) on each edge in T is a quartic spline with at most one knot. For both problems, we slightly modify their solutions to obtain a polynomial curve networks (quartic for (P p ) for p = 3/2, and cubic for (P ̂)) while preserving the same tangent planes at the data points.

Results from the experiments

Here we present and comment on two examples from the experiments performed and we share our observations and conclusions. We use BezierView [18] to visualize the highlight lines and the Gaussian curvature of the interpolation surfaces.

Example 1. We consider data obtained from a symmetric triangular pyramid. We have n = 4, V 4 = (0,0), and z i = 0, i = 1,2,3, z 4 = -1/2. The triangulation and the corresponding MNN are shown in Figure 1. The highlight lines on the interpolation surface are visualized in Figure 2 (left). The Gaussian curvature is visualized in Figure 2 (right) where the color scale goes from blue (low) to red (large). The minimum value of the Gaussian curvature is -11.8318, the maximum value is 34.3368. The surface that interpolates the edge convex MNN is shown in Figure 3. The highlight lines of this surface are visualized in Figure 3 (left) and its Gaussian curvature is visualized in Figure 3 (right). The minimum value of the Gaussian curvature is -0.6046, and the maximum value is 36.1598.

Table 1

The data for Example 2 Example 2. We consider the convex function f = exp((x -0.5) 2 + (y -0.5) 2 ) which is sampled at 25 points shown in Table 1. The triangulation, which is shown in Figure 4 (left), is the Delaunay triangulation. The corresponding edge convex MNN is shown in Figure 4 (right). The highlight lines on the surface are visualized in Figure 6.

For Example 1, see Figure 2, we can see that our interpolant for the unconstrained case is visually pleasant: the highlight lines are smooth with a small number of inflection points (if any), and the Gaussian curvature is evenly distributed. The highlight lines of the surface interpolating the edge convex MNN in Figure 3 (left) also look nicely as in the unconstrained case but the Gaussian curvature in Figure 3 (right) is not evenly distributed.

For Example 2 where the data are sampled from the convex exponential function, it is clearly seen that the interpolation surface generated from our algorithm presented in [15], see Figure 5 (right), significantly improves on the surface generated from Shirman and Sequin's algorithm [7,8], see Figure 5 (left). Although at first glance the surface in Figure 5 (right) visually even appears convex, when we examine the highlight lines in Figure 6 in detail, it is clearly seen that it is not convex.

The above observations shows that the algorithm presented in [15] can be refined and improved further. Another conclusion obtained from our preliminary numerical experiments is that using the modified quartic minimum L p -norm network for p = 3/2 gives better results with respect to the maximum distance than using the degree elevated MNN.

Conclusion and future work

Given the importance of surface modeling and simulation techniques in practice, it is important to better understand and create interpolation surfaces with good approximation properties based on different criteria. In the future we intend to enlarge our work by analyzing in-depth more data sets, which will support further optimization of our algorithm [15]. More details will be presented in the full version of this paper. The highlight lines of the surface interpolating the edge convex MNN and generated using Algorithm 2 in [15]