Constructing Explicit B-Spline

R.O. Linger H.R.N. van Erp  and  P.H.A.J.M. van Gelder

We introduce here a direct method to construct multivariate explicit B-spline bases. B-splines are piecewise polynomials, which are defined on adjacent tetrahedra and which are continuous throughout. The continuity is enforced by making sure that all directional derivatives of order , and lower, on the boundaries of adjacent tetrahedra give the same values for both tetrahedra. The method presented here is explicit, in that we will provide an algorithm with which one can analytically construct the B-spline base that enforces continuity for a given geometry.

1. Some Preliminaries

In order to give a direct method for constructing B-splines, we must first give some preliminaries of these B-splines. In what follows, we will introduce the concepts of barycentric coordinates, barycentric polynomials, and directional derivatives on barycentric polynomials. This is done for two-dimensional functions.

Say, we have three points , , , in a two-dimensional Cartesian coordinate system with coordinates. Then the difference vectors and are the axes of a coordinate system with origin and coordinates, say, and . We then have that


or, equivalently,


If a given coordinate lies in the triangle, say, , with origin and spanned by the vectors and , then


We define


By way of (1.3) and (1.4), we have that the coordinates , for a given point in triangle , are all greater than zero. Furthermore, these coordinates sum to


The coordinates are called barycentric coordinates. Using these coordinates, we may rewrite (1) as


or, equivalently,


Say, we have two barycentric coordinate systems and :

Figure 1. Simple Triangulation

Then a given Cartesian point in either or may be written as, respectively,




We call the collection of triangles and a triangulation of the Cartesian plane.

Note that for a given point in, say, the corresponding barycentric coordinates may be found by way of


and, by construction, (1.4):

If we have some th-order polynomial defined on the whole Cartesian plane , say:


Then we may make a Jacobian transformation from to , by way of (1.10) and (1.4), in order to obtain a polynomial in barycentric coordinates, equivalent to (1.11):


Likewise, we may make a Jacobian transformation from to , in order to obtain the equivalent polynomial:


Note that in [1], a more complicated proof of the equivalence of, say, (1.11) and (1.12) is given. But such a proof is unnecessary, in that a simple change of variable argument will suffice.

In what follows, we will have to make use of the fact that th order continuity of two different piecewise polynomials defined, respectively, on two connected triangles and , implies that all th-order directional derivatives of these polynomials be equal on their shared boundary. So, we introduce here the concept of directional derivatives of polynomials in barycentric coordinates.

Let be the directional vector


Then the barycentric coordinates of in are


Let be the first-order directional derivative operator in the Cartesian plane . Then we may define the zeroth-order directional derivative, in both the Cartesian and barycentric coordinate systems, as:


The first-order directional derivative is defined as:


where is defined, respectively, as

depending on the coordinate system on which the polynomial is defined. By way of (1.16) and (1), we may define the th-order directional derivative recursively as:


For example, the zeroth-, first-, and second-order directional derivatives for arbitrary and a second-order piecewise polynomial may be written down as, (1.12) and (1.16):





Having introduced th-order directional derivatives for barycentric polynomials (1.12) and (1.13), we may now proceed to the construction of B-spline bases that are th-order continuous throughout.

2. Enforcing Continuity

We now will construct B-spline bases that enforce th-order continuity throughout. We will show how to enforce zeroth- and first-order continuity for second-order polynomial functions defined on the plane. It is left to the reader to generalize to higher variate functions and higher order polynomials.

Say, we have the second-order polynomial:


Then we may make a change of variable from (2.1) to a polynomial which takes as its arguments the barycentric coordinates relative to the sides of the triangle :


So, if , then we have , for . Likewise, may make a change of variable from (2.1) to a polynomial which takes as its arguments the barycentric coordinates relative to the sides of the triangle :


So, if , then , for .

The transformed polynomials (2.2) and (2.3) are valid on the whole -plane. However, we will constrain the polynomials and to the triangles and , respectively. The basis for the unconnected piecewise polynomials of the triangulation , then may be given as:


where the columns of the basis correspond, respectively, with the coefficients

If , then the corresponding barycentric coordinates, found by way of (1.4) and (1.10), can be fed into the first row. But if , then the corresponding barycentric coordinates are fed into the second row.

Now, if we look at the triangulation in Figure 1 and equations (1.8) and (1.9), then we may see that all points on the boundary of are equivalent to the points on the boundary of , if and . Stated differently, in order for the polynomials (2.2) and (2.3) to be connected at their shared boundary, we must have that the zeroth-order directional derivative of these polynomials generate the same values for and , where


By substituting the values (2.5) into (1), we find for the piecewise polynomial defined on :


And for the piecewise polynomial defined on , we find:


Substituting (2.6) and (2.7) into (2.4), we obtain the constraint matrix, say, :


It follows that zeroth-order continuity may be enforced by merging Columns 1 and 12, Columns 2 and 10, and Columns 5 and 11 of the bases (2.4). In other words, if we constrain the coefficients in (2.2) and (2.3) to adhere to:


Substituting the constraints (2.9) into, say, (2.3), and rearranging the terms in both (2.2) and (2.3) so that equal coefficients are placed beneath each other, we obtain the barycentric polynomials




So, by way of (2.10) and (2.11), we find the basis that enforces zeroth-order continuity, that is, , to be:


where the columns of the basis correspond, respectively, with the coefficients

We now proceed to find the basis that enforces first-order continuity. For the triangles and , any directional vector which is non-parallel to the boundary shared by and will suffice:


where we have used (1.14) to express in terms of barycentric coordinates. By substituting both and the vertices of triangle into (1.10) and (1.4), we may obtain an equivalent directional vector in the barycentric coordinates of :


It may be checked, by way of (1.14), that

If we take the directional vector from (2.13) and substitute it in (1), we obtain:


The directional derivative of (2.11) is:

If we substitute the directional vector from (2.14) in (2), we obtain:


Further, substituting (2.5) into (2.15) and (2), we obtain, respectively,




Equations (2.18) and (2.19) correspond with the constraint matrix:


Now, any permutation of the columns of will leave intact any constraints already in place. For example, the columns of (2.12) enforce , and any permutation of these columns will still enforce this continuity constraint. So, for some polynomial function defined on a -dimensional hyperplane, we are free to permutate the columns of any constraint matrix , in order to obtain rows in which we are left with terms of the type


where , and where is the order of the polynomial function and is the order of the directional derivative or, equivalently, the order of the continuity constraint we wish to enforce.

If we thus reduce the rows of our constraint matrix , then we can instantly see how to enforce our continuity constraint. Just merge the permutated columns having identical elements. For example, in (2.8) we have the situation that no permutation matrix is needed, seeing that the elements of the constraint matrix are already in the form (2.21). This holds generally for zeroth-order constraint matrices .

In what follows, we give the steps needed to derive a permutation matrix that will deliver us (2.20) in the desired form. First, we observe in (2.20) that the columns 6 and 7 of basis (2.12) play no role in the enforcement of the first-order continuity. So, for the moment we will set these columns aside. This leaves us with seven non-zero columns. Now, the dimensionality of our original polynomial (2.1) is , the order of this polynomial is , and the continuity constraint we wish to enforce is . So, we wish to find that permutation matrix that only leaves us with the terms, (2.21):

in every row of (2.20).

If we transpose (2.20), then we may designate for each row and each distinct element in the set a separate column:


Then we append an identity matrix to the right-hand side of (2.22):


If we row reduce (2.23), by way of Gaussian elimination [2], then we get:


By dropping the first four columns of (2.24), we are left with the transposed permutation matrix:


The non-zero columns of (2.20) are


It may be checked that multiplying with the permutation matrix , brings the non-zero columns of (2.20) in the desired form:


Now that we have found the needed permutation matrix , we may proceed to the construction of a first-order continuous B-spline.

We select the columns in (2.12) which are non-zero (2.20):


The selected columns of (2.12) are then multiplied with the permutation matrix :

We then, because of (2.27), merge Columns 1 and 3, and Columns 2 and 4:


Finally, we add the columns in (2.12) which were dropped in (2.28):


Seeing that the columns of this basis are a linear combination of the columns of the basis (2.12), which was constructed to be zeroth-order continuous throughout, that is, , we have that the basis (2.31) is first-order continuous throughout, that is, .

In order to construct a B-spline basis that is second-order continuous throughout, that is, , we simply repeat the steps that led us to (2.31).

First, we take the second-order directional derivative of (2.31). By doing so, we may identify and put aside the columns in (2.31) which do not participate in the second-order continuity constraint, those being the zero-columns of the second-order directional derivative. This then gives us a reduced first-order explicit base. We then proceed to construct a permutation matrix for the remaining, that is, non-zero, columns of the second-order directional derivative. If we multiply the reduced first-order explicit base with this permutation matrix , and adding again those columns of (2.31) that were put aside, then we obtain the second order B-spline basis that is second-order continuous throughout, that is, .

By an iteration of these steps, we may arrive at the explicit B-spline basis that is , for . If we expand the triangulation in Figure 1, in that we add more triangles, we must, for a given , repeat these steps for every pair of triangles which share a side.

Note that for -variate B-splines, which have -variate triangulations, all tetrahedra which share -dimensional sides are considered to be connected.

3. A Simple B-Spline Analysis

Say, we have a triangulation as in Figure 1, a polynomial of degree and continuity throughout, then the corresponding B-spline base is (2.31):


The first and second row of correspond, respectively, with and of Figure 1.

If we have a small data set of observations :


Then we have that the input values, both in Cartesian and barycentric coordinates, (1.10) and (1.4), are given as: