B-splines

BSplines

Module defining B-spline bases and B-spline related functions.

AbstractBSplineBasis{k,T}

Abstract type defining a B-spline basis, or more generally, a functional basis defined from B-splines.

The basis is represented by a B-spline order k and a knot element type T.

(B::AbstractBSplineBasis)(
    x::Real, [op = Derivative(0)], [T = float(typeof(x))];
    [ileft = nothing],
) -> (i, bs)

Evaluates all basis functions which are non-zero at x.

This is a convenience alias for evaluate_all. See evaluate_all for details on optional arguments and on the returned values.

Examples

julia> B = BSplineBasis(BSplineOrder(4), -1:0.1:1)
23-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]

julia> i, bs = B(0.42)
(18, (0.0013333333333333268, 0.28266666666666657, 0.6306666666666667, 0.08533333333333339))

julia> sum(bs)
1.0

julia> bs[1] - B[i](0.42)
0.0

julia> bs[2] - B[i - 1](0.42)
-5.551115123125783e-17

julia> B(0.44; ileft = i)
(18, (0.01066666666666666, 0.4146666666666667, 0.5386666666666665, 0.03599999999999999))

julia> B(0.42, Float32)
(18, (0.0013333336f0, 0.28266668f0, 0.6306667f0, 0.085333325f0))

julia> B(0.42, Derivative(1))
(18, (0.19999999999999937, 6.4, -3.3999999999999977, -3.200000000000001))

BSplineBasis{k, T}

B-spline basis for splines of order k and knot element type T <: Real.

The basis is defined by a set of knots and by the B-spline order.

BSplineBasis(order::BSplineOrder{k}, ts::AbstractVector; augment = Val(true))
BSplineBasis(order::BSplineOrder{k}, ts::NTuple; augment = Val(true))

Create B-spline basis of order k with breakpoints ts.

If augment = Val(true), breakpoints will be "augmented" so that boundary knots have multiplicity k. Note that, if they are passed as a regular Vector, the input may be modified. See augment_knots! for details.

Examples

julia> breaks = range(-1, 1; length = 21)
-1.0:0.1:1.0

julia> B = BSplineBasis(BSplineOrder(5), breaks)
24-element BSplineBasis of order 5, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5  …  0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0]

Note that first and last knots are repeated $k = 5$ times.

If augment = Val(false), input breakpoints are taken without modification as the knots $t_i$ of the B-spline basis. Note that the valid domain is reduced to $[-0.6, 0.6]$. The domain is always defined as the range $[t_k, t_{N + 1}]$, where $N$ is the length of the basis (below, $N = 16$).

julia> Bn = BSplineBasis(5, breaks, augment = Val(false))
16-element BSplineBasis of order 5, domain [-0.6, 0.6]
 knots: -1.0:0.1:1.0

Statically-sized bases

To define a basis with static size (i.e. size known at compile time), the breakpoints ts should be passed as a tuple or as an SVector (from the StaticArrays package):

julia> breaks = (0.0, 0.1, 0.2, 0.6, 1.0);

julia> B = BSplineBasis(BSplineOrder(3), breaks)
6-element BSplineBasis of order 3, domain [0.0, 1.0]
 knots: [0.0, 0.0, 0.0, 0.1, 0.2, 0.6, 1.0, 1.0, 1.0]

julia> knots(B)
9-element StaticArraysCore.SVector{9, Float64} with indices SOneTo(9):
 0.0
 0.0
 0.0
 0.1
 0.2
 0.6
 1.0
 1.0
 1.0

boundaries(B::AbstractBSplineBasis)

Returns (xmin, xmax) tuple with the boundaries of the domain supported by the basis.

order(::Type{AbstractBSplineBasis}) -> Int
order(::AbstractBSplineBasis) -> Int
order(::Type{Spline}) -> Int
order(::BSplineOrder) -> Int

Returns order of B-splines as an integer.

knots(B::AbstractBSplineBasis)
knots(g::Spline)

Returns the knots of the B-spline basis.

getindex(B::AbstractBSplineBasis, i, [T = Float64])

Get $i$-th basis function.

This is an alias for BasisFunction(B, i, T) (see BasisFunction for details).

The returned object can be evaluated at any point within the boundaries defined by the basis.

Examples

julia> B = BSplineBasis(BSplineOrder(4), -1:0.1:1)
23-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]

julia> B[6]
Basis function i = 6
  from 23-element BSplineBasis of order 4, domain [-1.0, 1.0]
  support: [-0.8, -0.4) (knots 6:10)

julia> B[6](-0.5)
0.16666666666666666

julia> B[6, Float32](-0.5)
0.16666667f0

julia> B[6](-0.5, Derivative(1))
-5.000000000000001

length(g::BSplineBasis) -> Int

Returns the number of B-splines composing a spline.

PeriodicBSplineBasis{k, T}

B-spline basis for splines of order k and knot element type T <: Real.

The basis is defined by a set of knots and by the B-spline order.

PeriodicBSplineBasis(order::BSplineOrder{k}, ts::AbstractVector)

Create periodic B-spline basis of order k with knots ts.

The knot period is taken to be L = ts[end] - ts[begin].

Examples

Create B-spline basis on periodic domain with period $L = 2$.

julia> ts = range(-1, 1; length = 11)
-1.0:0.2:1.0

julia> B = PeriodicBSplineBasis(BSplineOrder(4), ts)
10-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
 knots: [..., -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, ...]

julia> period(B)
2.0

julia> length(B)
10

julia> boundaries(B)
(-1.0, 1.0)

julia> B(-0.42)
(5, (0.12150000000000002, 0.6571666666666667, 0.22116666666666668, 0.00016666666666666563))

julia> B(-0.42 + 2)
(15, (0.12150000000000015, 0.6571666666666667, 0.22116666666666657, 0.00016666666666666674))

julia> knots(B)
14-element PeriodicKnots{Float64, 4, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}:
 -1.4
 -1.2
 -1.0
 -0.8
 -0.6
 -0.4
 -0.2
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2

PeriodicKnots{T} <: AbstractVector{T}

Describes an infinite vector of knots with periodicity L.

PeriodicKnots(ξs::AbstractVector{T}, ::BSplineOrder{k})

Construct a periodic knot sequence from breakpoints ξs.

The knot period is taken to be L = ξs[end] - ξs[begin]. In other words, the breakpoints ξs are expected to include the endpoint ξs[begin] + L.

The breakpoints should be given in non-decreasing order.

Note that the indices of the returned knots ts are offset with respect to the input ξs according to ts[i] = ξs[i + offset] where offset = k ÷ 2.

period(bc::Periodic) -> Real
period(B::PeriodicBSplineBasis) -> Real
period(ts::PeriodicKnots) -> Real

Returns the period L associated to a periodic boundary condition or B-spline basis.

BasisFunction{B <: AbstractBSplineBasis, T}

Describes a single basis function.

The basis function may belong to a BSplineBasis (in which case it's effectively a B-spline), or to a basis derived from a B-spline basis (such as a RecombinedBSplineBasis).

BasisFunction(basis::AbstractBSplineBasis, i::Int, [T = Float64])

Construct i-th basis function of the given basis.

The constructed function can be evaluated as b(x), returning a value of type T.

(b::BasisFunction)(x, [op::AbstractDifferentialOp])

Evaluate basis function at coordinate x.

To evaluate a derivative, pass Derivative(n) as the op argument, with n the derivative order.

To evaluate multiple derivatives, pass a derivative range Derivative(m:n). In particular, Derivative(m:n) evaluates the basis function itself and its first n derivatives.

More general differential operators, such as Derivative(n) + λ Derivative(m), are also supported.

support(b::BasisFunction) -> UnitRange{Int}

Get range of knots supported by the basis function.

Returns the knot range i:j such that the basis function support is $t ∈ [t_i, t_j)$.

support(B::AbstractBSplineBasis, i::Integer) -> UnitRange{Int}

Get range of knots supported by the $i$-th basis function.

common_support(b1::BasisFunction, b2::BasisFunction, ...) -> UnitRange{Int}

Get range of knots commonly supported by different basis functions.

If the supports don't intersect, an empty range is returned (e.g. 6:5), following the behaviour of intersect. The lack of intersection can be checked using isempty, which returns true for such a range.

find_knot_interval(ts::AbstractVector, x::Real, [ileft = nothing]) -> (i, zone)

Finds the index $i$ corresponding to the knot interval $[t_i, t_{i + 1}]$ that should be used to evaluate B-splines at location $x$.

The knot vector is assumed to be sorted in non-decreasing order.

It also returns a zone integer, which is:

• 0 if x is within the knot domain (ts[begin] ≤ x ≤ ts[end]),
• -1 if x < ts[begin],
• 1 if x > ts[end].

This function is functionally equivalent to de Boor's INTERV routine (de Boor 2001, p. 74).

If one already knows the location i associated to the knot interval, then one can pass it as the optional ileft argument, in which case only the zone needs to be computed.

evaluate_all(
    B::AbstractBSplineBasis, x::Real,
    [op = Derivative(0)], [T = float(typeof(x))];
    [ileft = nothing],
) -> i, bs

Evaluate all B-splines which are non-zero at coordinate x.

Returns a tuple (i, bs), where i is an index identifying the basis functions that were computed, and bs is a tuple with the actual values.

More precisely:

• i is the index of the first B-spline knot $t_{i}$ when going from $x$ towards the left. In other words, it is such that $t_{i} ≤ x < t_{i + 1}$.

  It can be effectively computed as i = searchsortedlast(knots(B), x). If the correct value of i is already known, one can avoid this computation by manually passing this index via the optional ileft keyword argument.

• bs is a tuple of B-splines evaluated at $x$:

  \[(b_i(x), b_{i - 1}(x), …, b_{i - k + 1}(x)).\]

  It contains $k$ values, where $k$ is the order of the B-spline basis. Note that values are returned in backwards order starting from the $i$-th B-spline.

Computing derivatives

One can pass the optional op argument to compute B-spline derivatives instead of the actual B-spline values.

Examples

See AbstractBSplineBasis for some examples using the alternative evaluation syntax B(x, [op], [T]; [ileft]), which calls this function.

evaluate(B::AbstractBSplineBasis, i::Integer, x,
         [op::AbstractDifferentialOp], [T=Float64])

Evaluate $i$-th basis function in the given basis at x (can be a coordinate or a vector of coordinates).

To evaluate a derivative, pass Derivative(n) as the op argument, with n the derivative order.

More general differential operators, such as Derivative(n) + λ Derivative(m), are also supported.

See also evaluate!.

evaluate!(b::AbstractVector, B::BSplineBasis, i::Integer,
          x::AbstractVector, args...)

Evaluate i-th basis function at positions x and write result to b.

See also evaluate.

nonzero_in_segment(B::AbstractBSplineBasis, n::Int) -> UnitRange{Int}

Returns the range of basis functions that are non-zero in a given knot segment.

The $n$-th knot segment is defined by $Ω_n = [t_n, t_{n + 1}]$.

For BSplineBasis and RecombinedBSplineBasis, the number of non-zero functions in any given segment is generally equal to the B-spline order $k$. This number decreases near the borders, but this is not significant when B-spline knots have multiplicity $k$ there (as is the default).

For B-spline bases, and excepting the borders, the non-zero B-splines are $\left\{ b_i \right\}_{i = n - k + 1}^{n}$. This function thus returns (n - k + 1):N when B is a BSplineBasis.

See also support for the inverse operation.

BSplineOrder(k::Integer)

Specifies the B-spline order k.

AugmentedKnots{T,k} <: AbstractVector{T}

Pads from both sides a vector of B-spline breakpoints, making sure that the first and last values are repeated k times.

augment_knots!(breaks::AbstractVector, k::Union{Integer,BSplineOrder})

Modifies the input breakpoints to make sure that the first and last knot have multiplicity k for splines of order k.

To prevent allocations, this function will modify the input when this is a standard Vector. Otherwise, the input will be wrapped inside an AugmentedKnots object.

It is assumed that the input breakpoints have multiplicity 1 at the borders. That is, border coordinates should not be repeated in the input.

basis_to_array_index(B::AbstractBSplineBasis, axs, i::Int) -> Int

Converts from a basis index i (for basis function $b_i$) to a valid index in an array A with indices axs, where axs is typically the result of axes(A, dim) for a given dimension dim.

This is mainly relevant for periodic bases, in which case i may be outside of the "main" period of the basis.

It is assumed that length(B) == length(axs) (this is not checked by this function).

has_parent_basis(::Type{<:AbstractBSplineBasis}) -> Bool
has_parent_basis(::AbstractBSplineBasis) -> Bool

Trait determining whether a basis has a parent B-spline basis.

This is notably the case for RecombinedBSplineBasis, which are derived from regular B-spline bases.

static_length(::Type{<:AbstractBSplineBasis}) -> Union{Int, Nothing}
static_length(::AbstractBSplineBasis) -> Union{Int, Nothing}

Return the basis' length if it is statically known (i.e. at compile time); return nothing otherwise.

Typically, bases with statically-known length are those constructed using an SVector (from the StaticArrays package) to describe the basis breakpoints.

multiplicity(knots, i)

Determine multiplicity of knot knots[i].