B-Splines

# Constructing Explicit B-Spline

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

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

 x=v13+b1(v11−v13)+b2(v12−v13), y=v23+b1(v21−v23)+b2(v22−v23), (1.1)

or, equivalently,

 x=b1v11+b2v12+(1−b1−b2)v13, y=b1v21+b2v22+(1−b1−b2)v23. (1.2)

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

 0≤b1+b2≤1. (1.3)

We define

 b3≡1−b1−b2. (1.4)

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

 b1+b2+b3=1. (1.5)

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

 x=b1v11+b2v12+b3v13, y=b1v21+b2v22+b3v23, (1.6)

or, equivalently,

 (xy)=[v1v2v3]⎛⎜⎝b1b2b3⎞⎟⎠. (1.7)

Say, we have two barycentric coordinate systems and :

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

 (xy)=[v1v2v3]⎛⎜⎝b1b2b3⎞⎟⎠,(x,y)∈T1, (1.8)

and

 (xy)=[~v1~v2~v3]⎛⎜ ⎜⎝~b1~b2~b3⎞⎟ ⎟⎠,(x,y)∈T2, (1.9)

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

 (b1b2)=[v11−v13v12−v13v21−v23v22−v23]−1(x−v13y−v23), (1.10)

and, by construction, (1.4):

 b3=1−b1−b2.

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

 p(x,y)=∑0≤i+j≤dαijxiyj (1.11)

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):

 p(b1,b2,b3)=∑0≤i+j≤dγijbi1bj2bd−i−j3. (1.12)

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

 p(~b1,~b2,~b3)=∑0≤i+j≤d~γij~bi1~bj2~bd−i−j3. (1.13)

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

 u=a1v1+a2v2+a3v3. (1.14)

Then the barycentric coordinates of in are

 a=(a1,a2,a3)T. (1.15)

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:

 D(0)up(x,y)=p(x,y)=p(b1,b2,b3)=D(0)ap(b1,b2,b3). (1.16)

The first-order directional derivative is defined as:

 D(1)up(x,y) =uT∇p(x,y) =aT∇p(b1,b2,b3) (1.17) =D(1)ap(b1,b2,b3).

where is defined, respectively, as

 ∇=(∂∂x,∂∂y),and∇=(∂∂b1,∂∂b2,∂∂b3),

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):

 D(0)ap(b1,b2,b3) =∑0≤i+j≤2γijbi1bj2bd−i−j3 =γ00b23+γ10b1b3+γ01b2b3+γ11b1b2+γ20b21+γ02b22

and

 D(1)ap(b1,b2,b3) =aT∇∑0≤i+j≤2γijbi1bj2bd−i−j3 =aT⎛⎜⎝γ10b3+γ11b2+2γ20b1γ01b3+γ11b1+2γ02b22γ00b3+γ10b1+γ01b2⎞⎟⎠ (1.20) =γ002a3b3+γ10(a1b3+a3b1)+γ01(a2b3+a3b2) +γ11(a1b2+a2b1)+γ202a1b1+γ022a2b2.

and

 D(2)ap(b1,b2,b3) =aT∇[D(1)ap(b1,b2,b3)] =aT⎛⎜⎝γ10a3+γ11a2+γ202a1γ01a3+γ11a1+γ022a2γ002a3+γ10a1+γ01a2⎞⎟⎠ (1.21) =γ002a23+γ102a1a3+γ012a2a3 +γ112a1a2+γ202a21+γ022a22.

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:

 p(x,y)=α1+α2x+α3y+α4xy+α5x2+α6y2. (2.1)

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 :

 p1(b1,b2,b3)=γ00b23+γ10b1b3+γ01b2b3+γ11b1b2+γ20b21+γ02b22. (2.2)

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 :

 p2(~b1,~b2,~b3)=~γ00~b23+~γ10~b1~b3+~γ01~b2~b3+~γ11~b1~b2+~γ20~b21+~γ02~b22. (2.3)

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:

 (2.4)

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

 q1=b1=~b1,q2=b3=~b2,b2=~b3=0. (2.5)

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

 D(0)ap1(q1,0,q2)=γ00q22+γ10q1q2γ20q21+γ02. (2.6)

And for the piecewise polynomial defined on , we find:

 D(0)~ap(q1,q2,0)=γ11q1q2+γ20q21+γ02q22. (2.7)

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

 Q=[q22q1q200q210000000000000000q1q2q21q22]. (2.8)

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:

 γ20=~γ20,γ00=~γ02,γ10=~γ11. (2.9)

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

 p1(b1,b2,b3)=γ00b23+γ10b1b3+γ20b21+γ01b2b3+γ11b1b2+γ02b22. (2.10)

and

 p2(~b1,~b2,~b3)=γ00~b22+γ10~b1~b2+γ20~b21+~γ00~b23+~γ10~b1~b3+~γ01~b2~b3. (2.11)

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

 B=[b23b1b3b21b1b2b2b3b22000~b22~b1~b2~b21000~b23~b1~b3~b2~b3], (2.12)

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

 (γ00γ10γ20γ01γ11γ02~γ00~γ10~γ01).

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:

 u=(1,0)T=v2,or, % equivalently,a=(0,1,0)T, (2.13)

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 :

 ~a=(1,1,−1)T. (2.14)

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

 ~v1+~v2−~v3=(1,0)T=u.

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

 D(1)ap1(b1,b2,b3)=γ01b3+γ11b1+γ022b2. (2.15)

The directional derivative of (2.11) is:

 D(1)~ap2(~b1,~b2,~b3) =~aT⎛⎜ ⎜⎝γ10~b2+γ202~b1+~γ10~b3γ002~b2+γ10~b1+~γ01~b3~γ002~b3+~γ10~b1+~γ01~b2⎞⎟ ⎟⎠ =γ002~a2~b2+γ10(~a1~b2+~a2~b1)+γ202~a1~b1 +~γ002~a3~b3+~γ10(~a1~b3+~a3~b1)+~γ01(~a2~b3+~a3~b2).

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

 D(1)~ap2(~b1,~b2,~b3) =γ002~b2+γ10(~b2+~b1)+γ202~b1 −~γ002~b3+~γ10(~b3−~b1)+~γ01(~b3−~b2). (2.17)

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

 D(1)ap1(q1,0,q2)=γ01q2+γ11q1 (2.18)

and

 D(1)~ap2(q1,q2,0)=γ002q2+γ10(q1+q2)+γ202q1−~γ10q1−~γ01q2. (2.19)

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

 Q=[000q1q200002q2q1+q12q10000−q1−q2]. (2.20)

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

 qk11qk22⋯qknn, (2.21)

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):

 {q1,q2},

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:

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣02q20q1+q202q1q10q200−q10−q2⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⟶⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0002001100201000010000−10000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (2.22)

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

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣000210000000011010000000200010000100000010000100000010000−100000010000−10000001⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (2.23)

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

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1000000100001000000100001000000−100001000000−1000010000020000010001100000010020⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (2.24)

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

 PT=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0001000000010000000−10000000−1100000201000110010020⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (2.25)

The non-zero columns of (2.20) are

 Q1=[000q1q2002q2q1+q12q100−q1−q2]. (2.26)

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

 Q1P=[q1q20000000q1q2000]. (2.27)

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):

 B1=[b23b1b3b21b1b2b2b300~b22~b1~b2~b2100~b1~b3~b2~b3] (2.28)

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

 B2 =B1P =[b1b2b2b300b23b1b3b2100−~b1~b3−~b2~b3~b22+2~b2~b3~b1~b2+~b1~b3+~b2~b3~b21+2~b1~b3].

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

 B3=[b1b2b2b3b23b1b3b21−~b1~b3−~b2~b3~b22+2~b2~b3~b1~b2+~b1~b3+~b2~b3~b21+2~b1~b3]. (2.30)

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

 B4=[b1b2b2b3b23b1b3b21b220−~b1~b3−~b2~b3~b22+2~b2~b3~b1~b2+~b1~b3+~b2~b3~b21+2~b1~b30~b23]. (2.31)

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):

 B=[b1b2b2b3b23b1b3b21b220−~b1~b3−~b2~b3~b22+2~b2~b3~b1~b2+~b1~b3+~b2~b3~b21+2~b1~b30~b23]. (3.1)

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

If we have a small data set of observations :

 (x1,y1,z1) =(0.2,0.1,1.0), (x2,y2,z2) =(0.2,0.7,3.0), (x3,y3,z3) =(0.1,0.3,2.0), (3.2) (x4,y4,z4) =(0.5,0.1,1.0), (x5,y5,z5) =(0.7,0.8,4.0).

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

 (x1,y1) =(0.2,0.1)∈T1, (b11,b12,b13)=(0.8,0.1,0.1), (x2,y2) =(0.2,0.7)∈T2, (~b21,~b22,~b23)=(0.3,0.2,0.5), (x3,y3) =(0.1,0.3)∈T2, (~b31,~b32,~b33)=(0.7,0.1,0.2), (3.3) (x4,y4) =(0.5,0.1)∈T1, (b41,b42,b43)=(0.5,0.4,0.1), (x5,y5) =(0.7,0.8)∈T2, (~b51,~b52,~b53)=(0.2,