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.110915, 0.275742, 0.70627, 0.894336, 0.422503, 0.452695, 0.173176, 6.12789e-5]

julia> basis(spl) === b
true

julia> coeffs(spl) === c
true

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.41251124266445405

julia> splinevalue(spl, 2.5)
0.41251124266445405

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

julia> spl(2.5, Derivative(1))
-0.0013034002533023892

julia> splinevalue(spl, 2.5, Derivative(2))
-0.209661397174141

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);
  93.046 ns (1 allocation: 112 bytes)

julia> @btime $spl(2.5, workspace=$work);
  72.501 ns (0 allocations: 0 bytes)

julia> @btime $spl(2.5, workspace=$work, leftknot=$left);
  54.493 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 * spl
Spline{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 / 4
Spline{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 == ans
true

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> spl
Spline{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 + spl2
Spline{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 - spl2
Spline{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 == bspl
true

julia> support(spl) # support of the basis
(0, 5)

julia> support(bspl) # actual support of the B-spline
(0, 3)