The BSplineBasis type and related functions

To work with B-spline bases, this package defines the parametric BSplineBasis{T} type which represents B-spline bases with breakpoint vector of type T.

A B-spline basis is completely characterized by its order $k$ and knot vector $t$. In the case of the BSplineBasis type, the knot vector of a basis is generated from its breakpoint vector by repeating the first and last breakpoints so that they appear $k$ times.

Properties of B-spline bases – order, breakpoints, and knots

A BSplineBasis is constructed from its order and breakpoint vector. The order and breakpoints functions return these properties. The length function returns the number of B-splines in the basis:

julia> basis = BSplineBasis(4, 0:5)
8-element BSplineBasis{UnitRange{Int64}}:
 order: 4
 breakpoints: 0:5

julia> order(basis)
4

julia> breakpoints(basis)
0:5

julia> length(basis)
8

The knots function returns the knot vector that is generated from the breakpoints. In order to not allocate memory for a new array, a wrapper type around the original breakpoint vector is used:

julia> knots(basis)
12-element BSplines.KnotVector{Int64, UnitRange{Int64}}:
 0
 0
 0
 0
 1
 2
 3
 4
 5
 5
 5
 5

The interval over which a BSplineBasis is defined (i.e., the interval between the first and the last breakpoint) can be obtained (as a Tuple) with the support function:

julia> support(basis)
(0, 5)

Evaluating B-Splines and their derivatives

The bsplines and bsplines! functions can be used to obtain the values (or derivatives) of all B-splines that are non-zero at a given point.

Evaluating B-splines

If x is within the support of the B-spline basis, bsplines(basis, x) returns an OffsetArray which contains the value of the i-th B-spline at the index i. The array always has the length order(basis):

julia> basis = BSplineBasis(4, 0:5);

julia> bsplines(basis, 3.2)
4-element OffsetArray(::Vector{Float64}, 4:7) with eltype Float64 with indices 4:7:
 0.08533333333333327
 0.6306666666666667
 0.28200000000000014
 0.0020000000000000052

julia> bsplines(basis, 7//3)
4-element OffsetArray(::Vector{Rational{Int64}}, 3:6) with eltype Rational{Int64} with indices 3:6:
  4//81
 31//54
 10//27
  1//162

The type of the elements of the array depends on the type of the breakpoints and the type of x. If the point x is outside of the support of basis, nothing is returned instead:

julia> bsplines(basis, 6) # returns `nothing`, which is not printed in the REPL

Evaluating derivatives

The bsplines function can also calculate derivatives of the B-splines. The N-th derivative is specified via an optional third argument of the singleton type Derivative{N}. Instead of Derivative{N}(), the shorter constructor Derivative(N) can be used:

julia> basis = BSplineBasis(4, 0:5);

julia> bsplines(basis, 7//3, Derivative(1)) # calculate first derivatives of non-zero B-splines
4-element OffsetArray(::Vector{Rational{Int64}}, 3:6) with eltype Rational{Int64} with indices 3:6:
 -2//9
 -1//2
  2//3
  1//18

To calculate the value of the non-zero B-splines as well as all of their derivatives up to the N-1-th, the AllDerivatives{N} singleton type can be used instead. AllDerivatives(N) is equivalent to AllDerivatives{N}(). If x is within the support of the B-spline basis, bsplines(basis, x, AllDerivatives(N)) returns an OffsetArray which contains the j-th derivative of the i-th B-spline at the index i, j (the zeroth derivative is the value itself):

julia> bsplines(basis, 4, AllDerivatives(3)) # calculate zeroth, first and second derivatives
4×3 OffsetArray(::Matrix{Float64}, 5:8, 0:2) with eltype Float64 with indices 5:8×0:2:
 0.166667  -0.5    1.0
 0.583333  -0.25  -2.5
 0.25       0.75   1.5
 0.0        0.0    0.0

Pre-allocating output arrays

The bsplines function allocates a new array to return the values (except when it returns nothing). In order to use a pre-allocated array instead, the bsplines! function can be used: bsplines!(dest, args...) behaves like bsplines(args...), but the calculations are done in dest and the returned OffsetArray is a wrapper around dest.

julia> basis = BSplineBasis(4, 0:5);

julia> destvec = zeros(4);

julia> bsplines!(destvec, basis, 3.2)
4-element OffsetArray(::Vector{Float64}, 4:7) with eltype Float64 with indices 4:7:
 0.08533333333333327
 0.6306666666666667
 0.28200000000000014
 0.0020000000000000052

julia> parent(ans) === destvec
true

julia> bsplines!(destvec, basis, 7//3, Derivative(2))
4-element OffsetArray(::Vector{Float64}, 3:6) with eltype Float64 with indices 3:6:
  0.6666666666666665
 -0.9999999999999996
 -4.440892098500626e-16
  0.3333333333333335

julia> parent(ans) === destvec
true

julia> destmat = zeros(Rational{Int}, 4, 2);

julia> bsplines!(destmat, basis, 4, AllDerivatives(2))
4×2 OffsetArray(::Matrix{Rational{Int64}}, 5:8, 0:1) with eltype Rational{Int64} with indices 5:8×0:1:
 1//6   -1//2
 7//12  -1//4
 1//4    3//4
 0//1    0//1

julia> parent(ans) === destmat
true

When calculating values of B-splines or their derivatives via the Derivative{N} argument, dest must be a vector of length order(basis). In the case of the AllDerivatives{N} argument, it must be a matrix of size (order(basis), N). In any case, dest must not have offset axes itself.

Improving performance

The bsplines and bsplines! functions accept two optional keyword arguments (one of them only when evaluating derivatives) can be used to speed up the evaluation of splines:

• derivspace: When evaluating derivatives, coefficients are stored in a matrix of size (order(basis), order(basis)). In order to avoid allocating a new matrix every time, a pre-allocated matrix can be supplied with the derivspace argument. It can only be used when calculating derivatives, i.e., with Derivative(N) where N ≥ 1 or AllDerivatives(N) where N ≥ 2.
• leftknot: Evaluating B-splines at a point x requires finding the largest index i such that t[i] ≤ x and t[i] < t[end] where t is the knot vector. The bsplines and bsplines! functions use the intervalindex function to find this index. If the index is already known, it can be specified with the leftknot keyword argument.

julia> using BenchmarkTools

julia> basis = BSplineBasis(4, 0:5);

julia> dest = zeros(order(basis));

julia> space = zeros(order(basis), order(basis));

julia> left = intervalindex(basis, 2.5);

julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1));
  136.674 ns (2 allocations: 240 bytes)

julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1), derivspace=$space);
  105.235 ns (1 allocation: 32 bytes)

julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1), derivspace=$space, leftknot=$left);
  84.636 ns (0 allocations: 0 bytes)

Constructing knot vectors

As stated above, the knots of a BSplineBasis are generated from its breakpoints by repeating the first and last breakpoints so that they appear order(basis) times. However, in many situations it may be desirable to

• repeat other knots than the first and last or
• have the first and/or last knot repeated less than $k$ times (e.g. to implement certain boundary conditions).

In order to have repeated interior knots, it is sufficient to include it in the breakpoint vector multiple times when constructing the basis:

julia> basis = BSplineBasis(3, [1, 2, 3, 4, 5, 5, 6, 7, 8]) # 5 appears twice
10-element BSplineBasis{Vector{Int64}}:
 order: 3
 breakpoints: [1, 2, 3, 4, 5, 5, 6, 7, 8]

julia> knots(basis)
13-element BSplines.KnotVector{Int64, Vector{Int64}}:
 1
 1
 1
 2
 3
 4
 5
 5
 6
 7
 8
 8
 8

There is currently no way to create a BSplineBasis where the first and last knots appear less than $k$ times. However, some functions provided by this package (see Higher-level functions) provide a indices keyword argument to select only a certain range of B-splines from a basis, thus achieving the same result as if the knot vector had been shortened.

Because of repeated knots, not every pair of knots (t[i], t[i+1]) describes an actual interval. The intervalindices function helps with finding all indices which describe relevant intervals. See its docstrings for more information.