B-spline Curves: Computing the Coefficients

Although de Boor's algorithm is a standard way for computing the point on a B-spline curve that corresponds to a given u, we really need these coefficients in many cases (e.g., curve interpolation and approximation). We shall illustrate a simple way to do this.

Given a clamped B-spline curve of degree p defined by n+1 control points P₀, P₁, ..., P_n, and m+1 knots u₀=u₁=...=u_p=0, u_p+1, u_p+2, ..., u_m-p-1, u_m-p=u_m-p+1=...= u_m=1, we want to compute the coefficients N_0,p(u), N_1,p(u), ..., N_n,p(u) for any given u in [0,1]. A naive way is the use of the following recurrence relation:

However, this is a very time consuming process. To compute N_i,p(u), we need to compute N_i,p-1(u) and N_i+1,p-1(u). To compute N_i,p-1(u), we need to compute N_i,p-2(u) and N_i+1,p-2(u). To compute N_i+1,p-1(u), we need N_i+1,p-2(u) and N_i+2,p-2(u). As you can see, N_i+1,p-2(u) appears twice, and, as a result, its recursive computations will also be repeated. As the recursion goes deeper, more duplicated computations will occur. This is very similar to the recursive version of de Casteljau's algorithm discussed on a previous page. Consequently, the computation speed is very slow.

There is an easy way. Suppose u is in knot span [u_k,u_k+1). An important property of B-spline basis functions states that at most p+1 basis functions of degree p are non-zero on [u_k,u_k+1), namely: N_k-p,p(u), N_k-p+1,p(u), N_k-p+2,p(u), ..., N_k-1,p(u), N_k,p(u). By definition, the only non-zero basis function of degree 0 on [u_k,u_k+1) is N_k,0(u). As a result, the coefficients can be computed from N_k,0(u) in a "fan-out" triangular form as shown below:

Since N_k,0(u) = 1 on knot span [u_k,u_k+1) and other B-spline basis functions of degree 0 are zero on [u_k,u_k+1), we can start from N_k,0(u) and compute the basis functions of degree 1 N_k-1,1(u) and N_k,1(u). From these two values, we can compute the basis functions of degree 2 N_k-2,2(u), N_k-1,2(u) and N_k,2(u). This process repeats until all p+1 non-zero coefficients are computed.

In this computation, "internal" values such as N_k-1,2(u) has a north-west predecessor (i.e., N_k-1,1(u)) and a south-west predecessor (i.e., N_k,1(u)); values on the north-east direction boundary of the above triangle such as N_k-1,1(u) has only a south-west predecessor (i.e., N_k,0(u)); and values on the south-east direction boundary of this triangle such as N_k,2(u) has only a north-west predecessor (i.e., N_k,1(u)). Thus, values that are on the north-east (resp., south-east) direction boundary use the second (resp., first) term of the recurrence relation in the definition. Only the internal values use both terms. Based on this observation, we have the following algorithm:

Input

₀

Output

_0,p

_1,p

_n,p

Algorithm

₀

then

return

else

then

return

end if

// now u is between u₀ and u_m

Let u be in knot span [u_k,u_k+1);
N[k] := 1.0 // degree 0 coefficient
for d :=1 to p do // degree d goes from 1 to p

begin

k-d

for

k-d

end

The above is not a very efficient algorithm. Its purpose is to illustrate the idea in an intuitive and easy to understand way. Array N[ ] holds all intermediate and the final results. For a degree d, N[i] stores the value of N_i,d(u), and, at the end, N[k-d], N[k-d+1], ..., N[k] contain the non-zero coefficients. The computation starts with d=1 because we know that the only non-zero basis function is N_k,0(u) if u is in knot span [u_k,u_k+1). The outer loop lets degree d go from 1 to p. The first assignment following begin computes N_k-d,d(u) using only one term (i.e., its south-west term in the triangle, N_k-d+1,d-1(u)), the inner for loop computes the "internal" terms, and the last statement in the outer loop computes N_k,d(u) using only one term (i.e., its north-west term in the triangle, N_k,d-1(u)).

Can you make this algorithm more efficient?