A common engineering task consists of interpolating a set of discrete points that arise from measurements and experiments. Another traditional requirement implies creating a curve that mimics a given array of points, namely, a polyline. Any of these problems require building an analytical representation of the given discrete set of points. If the geometrical shape represented by the input polyline is complicated, then we may expect that a global interpolant or polynomial will be of a high degree, to honor all imposed constraints, which makes its use prohibited. Indeed, a global interpolant often experiences inflection points and sudden changes in curvature. To avoid these drawbacks, we often seek solving the interpolation/approximation problem using piecewise polynomial functions called “splines.”
- cubic splines
- tension splines
- Bèzier curves
In this chapter, we tackle passing a curve/surface through a given set of data points. We first focus in the case in which the curve to be constructed can be described as . We refer this data set as scalar data. We describe herein cubic and tension splines, which are powerful interpolants suitable to tackle large data sets. We then introduce a parametric case for a vector-valued curve and hence, it is able to represent arbitrary topologies. We explain how to construct piecewise continuous cubic Bèzier curves called “B-splines.” We cover the interpolation and approximation problems with B-splines, this latter denoted as well as “inverse” design. We extend our treatment to tensor product surfaces that are referred as piecewise bicubic B-splines. Applications encompass translational and interpolation surfaces. We briefly introduce nonuniform rational B-spline curves and surfaces (NURBS). We present applications such as approximating conic sections. We finalize the chapter introducing Duchon splines that are radial basis functions to interpolate scattered data sets in two or three dimensions.
2. Scalar splines
We cover herein the scalar case in which a spline function fits a given set of sorted point pairs. We introduce cubic splines and their specialized version that offers a “tension” parameter that allows attracting the interpolant toward the polyline that connects the input points, i.e., linear spline. We refer to this latter as tension splines. The last section presents a couple of numerical examples.
2.1. Cubic splines
A spline is a piecewise continuous function consisting of several polynomials, each specified in a subinterval, bound themselves by certain continuity conditions. Let be sorted points such that whose corresponding values are denoted by . A spline of degree with knots is a function such that:
1. is a polynomial of degree that is continuous up to derivative over .
2. Two adjacent splines need to have continuity at the junction points:
We thus enforce continuity conditions at the junction points which yields to equations to determine unknown spline coefficients. We omit details herein but refer the reader to [1, 2]. We end up with a tridiagonal system for the unknown curvature values , at the junction points:
where and . The last equation provides a system of conditions for . Since both and are arbitrary, a logical choice is choosing , which we refer as “natural spline.” For the latter, we can write in matrix form:
Once we determine the curvature values , by solving Eq. (3), we define the spline function as
where the coefficients are given by
2.2. Tension splines
In some problems of adjusting discrete data, it is useful to have a parameter called “tension, .” When has a small value, the resulting curve approaches a cubic spline. When tends to , the resulting curve approaches a linear spline. For the same sequence of sorted point pairs mentioned above, the tension spline satisfies.
1. and .
2. On every interval .
That is, has continuity globally, interpolates to the data, and satisfies certain ordinary differential equation in each subinterval. It is clear that the prescription leads to cubic polynomials when solving the equation. To determine , we proceed similarly to the case of cubic splines, i.e., :
The solution is given by :
where and , and we compute the curvatures by solving the system:
and , and the arguments are ():
2.3. Numerical examples
2.3.1. Example 1
Fit the following collection of point pairs using a natural cubic spline.
We assume ; thus , and . The tridiagonal system (3) yields
As an illustration, is given by
for instance, and .
2.3.2. Example 2
Figure 1 depicts radial velocity profiles that represent the laminar fluid flow within a pipeline.
These velocity profiles were obtained by solving the Navier-Stokes equations under simplifying assumptions. The symbols represent the discrete point pairs, the abscissas correspond to the normalized radial coordinate from the center, and the -coordinates are the normalized radial velocities. We fit all data sets by using natural splines. To solve the system (3), we recommend the Thomas method, i.e., a direct frontal solver for tridiagonal matrixes [1, 2, 3]. We also recommend employing a quick-search algorithm to evaluate the piecewise function. Indeed, for an arbitrary , we need to determine what is the interval where this abscissa lies, i.e., .
2.3.3. Example 3
We finalize the examples by comparing cubic and tension splines. Figure 2 depicts a car-like profile polygon that we would like to interpolate. We try both splines mentioned above. We highlight in red color the cubic spline (top) interpolant, while the tension spline is black (bottom curve). We observe that the cubic spline experiences inflection points because the car-shaped polygon is challenging. This latter is the kind of application for tension splines where we seek to attract the spline toward the input polyline. We notice that we achieve that goal herein.
3. Bèzier, B-spline, and NURBS curves
The appropriate representation and meshing of the computational domain for the physical problem under study are necessary premises for a satisfactory computer simulation. In fact, one of the most demanding computational tasks in a simulation is defining the geometry because it will impact many aspects of the study such as the grid generation process . Therefore, special methods must be applied to fit discrete data without sudden changes in curvature. The approach should be free of inflection points, and at minimum, it must enforce continuity of the fitted curve. In this chapter, this goal is achieved by using Bèzier, B-spline, and NURBS curves and surfaces [5, 6].
A Bèzier curve (BC), , shown in Figure 3, is obtained by specifying the coordinates of a series of points in space, such that only the first and last ones fall on the originally given curve. All these points are known as control points, and the polyline resulting from connecting them with straight lines is called control polygon, which mimics the original curve, allowing an easy control of its shape. Although inflection points may be present in Bèzier curves, they are less common than in polynomials or other analytical functions [5, 6].
Global Bèzier curves, i.e., only one curve represents the given polyline, provide a powerful tool in geometry definition; however, complex shapes require a large number of constraints, making their use prohibitive. It is therefore beneficial to represent them by using piecewise continuous Bèzier curves called B-spline curves . In fact B-spline curves are a widely utilized representation for geometrical entities in computer-aided geometric design (CAGD) systems. Their convex hull, local support, shape-preserving forms, affine invariance, and variation-diminishing properties are extremely attractive in engineering design applications .
A particular Bèzier curve is set up by its parametric representation; let be defined by
here, denotes the order or degree of the curve, are the Bernstein polynomials, defined as
and are the control points. Notice in Eq. (14) that Bernstein polynomials satisfy the barycentric property, meaning that they add up to 1, which explains why a given curve cannot be outside its control polygon that is the convex-hull property. The control points of a given BC can be calculated in several ways since the Bèzier curve evaluated in must provide the corresponding base point ; a linear system of equations can be formed for the unknown control points as
which is the well-known chord-length parametrization.
The last approach is a powerful tool in curve design, but it has a limitation: if the geometry that we model has a complex shape (i.e., a significant number of base points), then its Bèzier curve representation may be of a prohibitively high degree. Since the Bèzier curve is forced to satisfy several constraints according to Eq. (15), the resulting curve may experience inflection points and sudden changes in curvature (see Figure 4, where points in Eq. (15) are represented by circles). For practical purposes, degrees exceeding 10 are prohibitive [5, 6].
Such complex geometries can be modeled using piecewise polynomial curves named B-spline curves [5, 6] (see Figure 5). B-spline curves are a set of Bèzier curves of degree that must satisfy at least the continuity. A spline curve is the continuous mapping of a collection of global parameter values into , where each interval is mapped onto a polynomial curve segment as shown in Figure 5. We define as the computational space. A local coordinate for the interval can be defined by setting :
3.1. cubic curves
Let be a set of points defining the de Boor’s polygon that generates individual cubic curves as shown in Figure 6. The required Bèzier control points are calculated with the aid of and continuity criteria. conditions lead to
while conditions require that
where and . The end points are
This construction is due to W. Boehm .
Chord-length parametrization 
A parametrization proposed by the author 
This latest parametrization has the advantage that it yields to a symmetric curve if the control polygon is symmetric as well, which may be interesting for certain applications.
3.2. Inverse design and interpolation problems
A two-dimensional geometric description based on B-spline curves requires the definition of a control polygon (the de Boor’s polygon) that mimics the curve. Therefore two approaches are possible. The first method consists of providing the set of the de Boor’s points (i.e., define the de Boor’s polygon interactively from user’s input) that defines the composite curve, which is known as “inverse design,” and it has its application in “TrueType” font technology as shown in Figure 7. The second approach consists of defining the set of base points and then solves a linear system of equations for the de Boor’s points such that the resulting curve passes through them. This latter is known as the “interpolation” problem.
Figure 7 shows the difference between the above approaches; from left to right, it has the de Boor’s control polygon, inverse design, and interpolation problems both taking into account the same polygon as an argument with cubic curves.
3.3. Interpolation with cubic curves
In order to interpolate with cubic B-spline curves, we find unknown junction points such that they pass through a given set of data points and corresponding parameter values . A composite cubic curve , determined by its de Boor’s polygon [4, 5] such that , is required. The solution to this problem is obtained by finding the relationship between the data points and the control vertices . This leads to the following linear system of equations for the unknown de Boor’s points :
where (with )
The first and last vertices of the polygon are given by
The points and can be calculated from a given end condition. There are two possibilities for the choice of B-spline ending conditions. A natural spline requires that
Using the relations in Eq. (29), we obtain that
These two equations replace the first and last rows of the linear system of equations given in Eq. (26). Notice that in either case, natural ending conditions or prescribed tangent vectors, the linear system in Eq. (24) is a tridiagonal matrix. Since the coefficient matrix is real and scalar, and the left- and right-hand-side vectors are in fact hypervectors (i.e., an array of vectors), it is recommendable to use a type of Gaussian elimination method against multiple right-hand sides to achieve performance [6, 7]. If we recompute the interpolant in Figure 4 but this time with a B-spline cubic curve, we then obtain a smoother and steadier interpolant free of unwelcome inflection points.
3.4. NURBS curves
NURBS curves are useful when we require an exact geometrical representation of some entities, such as circles, parabolas, ellipses, spheres, cylinders, etc. This is precisely the case in various applications in aerospace and mechanical engineering where NURBS are quite popular [4, 5, 8, 9, 10]. For instance, a NURBS curve is defined by its rational representation in Eq. (31):
where the weights are positive real scalars. The usual B-spline definition is recovered if all those weights equal 1. Generally speaking, the weights play the role of attracting the curve toward its control polygon when we increase their values [5, 8]. It turns out that specific weights lead to exact representation of circles, for instance, as shown in Figure 8, where the real numbers depicted are the given weights. A circle can be exactly represented by three- or four-quadratic arcs (left and right side, respectively, in Figure 8); this latter alternative is more attractive, for instance, to generate a surface of revolution [4, 5, 8]. NURBS also provide exact representation for 3D surfaces such as spheres and cylinders as well as volumes [4, 10]. The reader may refer to [4, 5, 8, 9, 10, 11] for further details for this well-established area of computational geometry. We depict an example of practical interest, in the context of geomechanics, in Figure 9. Herein we precisely represent a near-borehole section. To describe the borehole geometry, we employ four line segments and a quadratic arc as shown.
In order to interpolate a set of points with a NURBS curve in two or three dimensions, one may follow the same procedure described with B-spline curves except that a mapping to must be carried out first. The new input points are given by
then the linear system in Eq. (24) with the right boundary conditions can be solved with replaced by , as defined by the Eq. (32). After solving this system, the solution control polygon is still in , which implies that a mapping back to is required:
The latter is a straightforward procedure to reuse the subroutines already developed for B-spline curves.
3.5. Numerical examples
We implemented the proposed approaches in a computer code named “LogProc” which is a graphical user interface application developed with the C++ programming language. LogProc is proprietary software, but a free community version will be available for download from
Figure 11 shows a zenith view of a pair of petroleum reservoir outlines. Few base points permit to approximate their complex shape. These shapes could be used as arguments to generate an unstructured mesh appropriate to simulate the flow in a porous media .
We interpolate a discrete blade geometry approximation in both Figures 12 and 13. A3K7 profiles are constructed in the same way as NACA 65, but they present rounded trailing edges . A3K7’s shape is represented by 46 base points and circle arcs both in leading and trailing edges. Therefore the curves that interpolate the data points are tangents to circle arcs to ensure continuity between these entities. The profiles in the left of Figure 12 were obtained applying Eq. (15) to construct single Bèzier curves (two curves, each of them approximating to suction and pressure sides, respectively). However, note that the resulting curves have unexpected inflection points on the suction side. The enforcement of continuity conditions implies a large number of constraints after Eq. (15), which explains this behavior. If we replace the former Bèzier curves by two B-spline curves, computed from Eqs. (24) and (29), we obtain a smooth A3K7 geometry description (see profiles in the right-hand side of Figure 12). Figure 14 zooms in to show that the resulting geometrical description is smooth and free of unwelcome features such as inflection points.
4. Interpolation surfaces
Let be a two-parameter mapping which represents a given surface. If structured data, i.e., tensor product data, needs to be interpolated, one may expect to come up with tensor product surfaces as well, where two parameters allow covering two different directions associated with the surface. In the computational space, i.e., in the plane , the domain, , is a rectangle, and its image is the surface in 3D as shown in Figure 15, where and are the number of curves in their respective directions.
B-spline tensor product surfaces allow interpolating structured data, and they are defined as the tensor product of two families of curves and , which is
4.1. Creating a surface of interpolation
The following steps describe creating a surface of interpolation:
1. An input control polygon, whose points are in , is provided. They correspond to data that is structured and ordered, which is usually a matrix-type array of points (see left side of Figure 16). For simplicity, points in the -direction are associated with the parameter while are associated with .
2. Create interpolation curves with constant values of , so-called -interpolants (see right side of Figure 16).
3. Proceed accordingly with previous step, interpolation curves with constant values of ; so-called -interpolants are created this time (see left side of Figure 17).
4. Compute the tensor product between - and -interpolants in order to get bicubic patches (see right side of Figure 17).
The right-hand side in Figure 17 shows typical bicubic patches as a chessboard surface emphasizing that we deal with a piecewise continuous entity. The computational cost associated with the above algorithm is reasonable because the most expensive part is computing the interpolants (see Section 3).
5. Translational surfaces
These surfaces are again a two-parameter mapping, , but their construction procedure is simpler than interpolation surfaces; see, for instance, [5, 8, 9]. The idea here is just literally translating a curve along another curve , which yields
This idea became very popular in CAGD systems long time ago. Those systems usually support a command which allows extruding a geometrical entity, for instance, a cylinder can be easily created by extruding a circle along a straight line. Figure 18 shows the above procedure applied to an aircraft wing where an NURBS airfoil profile was translated or extruded along a straight line accordingly. We interpolated an NACA 65 polyline with a NURBS curve as we mentioned in Section 3.4.
This procedure becomes very useful in the geometrical reconstruction of oil reservoirs (RS). Indeed, we reconstructed the geometry of RS in  by using B-spline surfaces. The technique exploits input mesh’s simplicity to build a robust piecewise continuous geometrical representation using Bèzier bicubic patches. We manage the reservoir’s topology with interpolation surfaces, while translational surfaces allow extrapolating it toward its side burdens. After that, transfinite interpolation can be applied to generate decent hexahedral meshes. Figure 19 shows a sample translational surface that we obtain by extruding a curve that interpolates the reservoir’s edge as shown. We render the surfaces in blue color with a white wireframe, while the RS is the color-contoured surface that represents the porosity, a scalar property. We tackle the RS itself after interpolating the control polygon that Figure 20 highlights in red color. The polygon is a array of points representing the RS topology. The procedure works well for a variety of so-called open-to-the-public RS data sets that we reconstructed in . It is also possible to utilize these NURBS curves and surfaces as interfaces for gluing nonmatching interfaces for the finite element method as we showed in .
6. Duchon splines
In the context of applications in statistical analysis involving very high dimensional data sets, response surfaces are growing popularity. By running the simulations at a set of points (e.g., experimental design) and fitting response surfaces, i.e., splines, for instance, to the resulting input-output data that is characterized by sparsity, we can obtain fast surrogates for the objective function for optimization purposes [14, 15]. The appeal of the latter approach goes beyond reducing runtime. Since the method begins with experimental design, statistical analyses can be done to identify which input variables are the most important, and thus we can create “main effect plots” to visualize input-output relationships . We must recognize interpolation methods in which the basis functions are fixed and those in which they have parameters that are tuned (e.g., kriging, which has a statistical interpretation that allows one to construct an estimate of the potential error in the interpolator). We refer the reader to [14, 15] for further reading.
There are different ways to approximate a function of several variables: multivariate piecewise polynomials, splines, and tensor product methods, among others. All these approaches have advantages and drawbacks, but if the rank of the linear system to solve may become large, a natural choice is radial basis functions, which are also useful in lower dimensional problems [14, 16, 17]. This may be particularly true if the input data is scattered, which excludes tensor product methods at first glance. Duchon splines are a class of positive definite and compactly supported radial functions, which consist of univariate polynomial within their support. It can be proven that they are of minimal degree and unique up to a constant factor, for given smoothness and space dimension . They are particularly suitable to compute interpolants for very large scatter datasets .
where is a linear polynomial in two or three dimensions:
Notice that and the polynomial coefficients are all scalar quantities. In order to guarantee existence and uniqueness for these splines, an orthogonality condition with respect to linear polynomials is enforced, for instance, in two dimensions this yields to
By considering this result, the interpolation problem becomes
which implies points plus orthogonality conditions; here, are the nodal values to be interpolated. The resultant linear system to solve for is of rank.
Duchon splines are certainly suitable to interpolate scattered data sets that we cannot tackle with the tensor product surfaces that we discussed before. Indeed, Figure 21 depicts such an application, in optimization, where an objective function that we wish to minimize was sampled randomly by Monte-Carlo (MC) realizations. To compute a minimum, we interpolate the black dots, and then we minimize the resulting spline with standard Newton stochastic techniques . It is true that Duchon splines are a valid choice for “surrogate” models for such applications.
7. Concluding remarks
We presented a concise introduction to scalar and parametric spline interpolants. We introduced cubic and tension splines for scalar functions, and then we generalized them for the parametric case via Bèzier, B-spline, and NURBS curves. These latter entities are of the particular interest for applications in CAGD. We thus elaborated on topics such as inverse design and interpolation. We extended the treatment also to cover interpolation and translational surfaces with examples in mechanical and petroleum engineering. We wrapped up with the topic of interpolating sparse very high dimensional data sets via Duchon splines which are a kind of response surfaces suitable for applications in statistical analysis and optimization.
We acknowledge the financial support of the project “Reduced-Order Parameter Estimation for Underbody Blasts” funded by Army Research Laboratory.