The Spline type and related functions

A spline consists of a B-spline basis and a coefficient vector of the same length. The Spline{B, C} type represents a spline with a basis of type B and a coefficient vector of type C. The constructor Spline(basis, coeffs) returns a spline with the given B-spline basis and coefficients. The basis on which the spline is defined can be obtained with the basis function. The coeffs function returns the coefficient vector of a spline.

julia> b = BSplineBasis(4, 0:5);
julia> c = rand(length(b));
julia> spl = Spline(b, c)Spline{BSplineBasis{UnitRange{Int64}}, Vector{Float64}}:
 basis: 8-element BSplineBasis{UnitRange{Int64}}:
  order: 4
  breakpoints: 0:5
 coeffs: [0.736361, 0.975439, 0.658518, 0.97874, 0.571551, 0.0275178, 0.849618, 0.107628]
julia> basis(spl) === btrue
julia> coeffs(spl) === ctrue

Evaluating Splines and their derivatives

To evaluate a Spline, it can be called as a function. Alternatively, the splinevalue function can be used:

julia> spl = Spline(BSplineBasis(4, 0:5), rand(8));
julia> spl(2.5)0.6595129083193956
julia> splinevalue(spl, 2.5)0.6595129083193956

To evaluate derivatives of a spline, the optional Derivative(N) argument is used:

julia> spl(2.5, Derivative(1))-0.10093211466441107
julia> splinevalue(spl, 2.5, Derivative(2))-0.24898941775511785

Improving performance

Two optional keyword arguments can be used to speed up the evaluation of splines:

• workspace: To evaluate a spline, a vector of length order(spline) is used to store intermediate values. By default, a new vector is allocated for each evaluation. To avoid this, a pre-allocated vector can be specified with the workspace keyword argument. Note that in this case, the calculations are performed with numbers of type eltype(workspace).
• leftknot: Evaluating a spline 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 of the basis. By default, the intervalindex function is used to find this index. If the index is already known, it can be specified with the leftknot keyword argument.

julia> using BenchmarkTools
julia> spl = Spline(BSplineBasis(4, 0:5), rand(8));
julia> work = zeros(order(spl));
julia> left = intervalindex(basis(spl), 2.5);
julia> @btime $spl(2.5);  91.145 ns (1 allocation: 96 bytes)
julia> @btime $spl(2.5, workspace=$work);  68.540 ns (0 allocations: 0 bytes)
julia> @btime $spl(2.5, workspace=$work, leftknot=$left);  52.728 ns (0 allocations: 0 bytes)

Arithmetic with Splines

A Spline can be multiplied or divided by a real number using *, /, or \. These operations return a new spline:

julia> spl = Spline(BSplineBasis(3, 0:5), [1:7;])Spline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [1, 2, 3, 4, 5, 6, 7]
julia> 3.0 * splSpline{BSplineBasis{UnitRange{Int64}}, Vector{Float64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0]
julia> spl / 4Spline{BSplineBasis{UnitRange{Int64}}, Vector{Float64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75]
julia> 4 \ spl == anstrue

Using the lmul!/rmul! and ldiv!/rdiv! functions from the LinearAlgebra stdlib, a Spline can be multiplied/divided by a real number in-place, which modifies the spline (if its coefficient vector is mutable):

julia> using LinearAlgebra: rmul!
julia> spl = Spline(BSplineBasis(3, 0:5), [1:7;])Spline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [1, 2, 3, 4, 5, 6, 7]
julia> rmul!(spl, 3);
julia> splSpline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [3, 6, 9, 12, 15, 18, 21]

If two splines are defined on the same basis, they can be added and subtracted from each other using + and -:

julia> spl1 = Spline(BSplineBasis(3, 0:5), [1:7;]);
julia> spl2 = Spline(BSplineBasis(3, 0:5), [2:8;]);
julia> spl1 + spl2Spline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [3, 5, 7, 9, 11, 13, 15]
julia> spl1 - spl2Spline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 7-element BSplineBasis{UnitRange{Int64}}:
  order: 3
  breakpoints: 0:5
 coeffs: [-1, -1, -1, -1, -1, -1, -1]

The BSpline type

Indexing into or iterating over a B-spline basis yields BSplines. The parametric type BSpline{T} represents a single B-spline from a B-spline basis of type T:

julia> b = BSplineBasis(4, 0:5);
julia> b[3]BSpline{BSplineBasis{UnitRange{Int64}}}:
 basis: 8-element BSplineBasis{UnitRange{Int64}}:
  order: 4
  breakpoints: 0:5
 index: 3 (knots: [0, 0, 1, 2, 3])

BSplines can be used like any other Spline. However, they behave differently in two ways:

• BSplines are printed differently: instead of the coefficient vector, its index within the basis is shown.
• Calling support on a BSpline yields its actual support (the interval on which it is non-zero) instead of the support of the basis.

The second point is illustrated by the following example: even though spl and bspl compare equal, their support differs:

julia> spl = Spline(b, [0,0,1,0,0,0,0,0])Spline{BSplineBasis{UnitRange{Int64}}, Vector{Int64}}:
 basis: 8-element BSplineBasis{UnitRange{Int64}}:
  order: 4
  breakpoints: 0:5
 coeffs: [0, 0, 1, 0, 0, 0, 0, 0]
julia> bspl = b[3]BSpline{BSplineBasis{UnitRange{Int64}}}:
 basis: 8-element BSplineBasis{UnitRange{Int64}}:
  order: 4
  breakpoints: 0:5
 index: 3 (knots: [0, 0, 1, 2, 3])
julia> spl == bspltrue
julia> support(spl) # support of the basis(0, 5)
julia> support(bspl) # actual support of the B-spline(0, 3)