Spline interpolations

Interpolating data

BSplineKit can interpolate evenly and unevenly-distributed data.

For example, we can try interpolating random data on randomly-distributed data points:

using Random
rng = MersenneTwister(42)

Ndata = 20
xs = range(0, 1; length = Ndata) .+ 0.01 .* randn(rng, Ndata)
sort!(xs)  # make sure coordinates are sorted
xs[begin] = 0; xs[end] = 1;   # not strictly necessary; just to set the data limits
ys = sinpi.(xs) .+ 0.02 .* randn(rng, Ndata);

Let's start by plotting the generated data:

using CairoMakie
scatter(xs, ys; label = "Data", color = :black)

To interpolate the data using splines, we may choose any arbitrary B-spline order $k$ (in particular, $k = 4$ corresponds to cubic splines). The main interpolation function is interpolate.

using BSplineKit
S = interpolate(xs, ys, BSplineOrder(4))

SplineInterpolation containing the 20-element Spline{Float64}:
 basis: 20-element BSplineBasis of order 4, domain [0.0, 1.0]
 order: 4
 knots: [0.0, 0.0, 0.0, 0.0, 0.102874, 0.162913, 0.211638, 0.263706, 0.310932, 0.358573  …  0.638901, 0.692134, 0.740085, 0.790027, 0.836131, 0.893086, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-0.0210198, 0.147291, 0.24597, 0.508261, 0.416719, 0.841692, 0.801509, 0.972899, 0.870355, 0.999667, 1.0157, 0.959612, 0.931631, 0.825159, 0.742533, 0.641917, 0.431575, 0.325183, 0.127764, 0.0140354]
 interpolation points: [0.0, 0.0481877, 0.105535, 0.1549, 0.228305, 0.251709, 0.311103, 0.369982, 0.394633, 0.483717, 0.53714, 0.580818, 0.63676, 0.699124, 0.740518, 0.780612, 0.848951, 0.878831, 0.951475, 1.0]

Let's plot the result:

lines!(0..1, x -> S(x); label = "k = 4", color = Cycled(4 - 3))

We can also choose other interpolation orders for comparison:

for k ∈ (5, 6, 8)
    S = interpolate(xs, ys, BSplineOrder(k))
    lines!(0..1, x -> S(x); label = "k = $k", color = Cycled(k - 3))
end
axislegend()

We see larger and larger oscillations, especially near the boundaries, as the spline order increases. We can try to fix this using natural splines, which impose some derivatives to be zero at the boundaries.

Natural splines

Natural splines usually refer to cubic splines (order $k = 4$) with the additional constraint $S''(a) = S''(b) = 0$ at the boundaries ($x ∈ \{a, b\}$).

In BSplineKit, this concept is generalised for all even spline orders $k$, by setting all derivatives of order $2, 3, …, k / 2$ to be zero at the boundaries. For instance, for $k = 6$ (quintic splines), the constraint is $S'' = S''' = 0$. We achieve this by using basis recombination to implicitly impose the wanted boundary conditions.

The natural boundary condition not only allows to suppress oscillations near the boundaries, but is also quite convenient for interpolations, as it reduces the number of degrees of freedom such that the number of unique spline knots is equal to the number of data points. This simply means that one can set the knots to be equal to the data points. All of this is done internally when using this boundary condition.

To impose natural boundary conditions, one just needs to pass Natural to interpolate, as illustrated below.

k = 8
S = interpolate(xs, ys, BSplineOrder(k))  # without BCs
Snat = interpolate(xs, ys, BSplineOrder(k), Natural())  # with natural BCs

SplineInterpolation containing the 20-element Spline{Float64}:
 basis: 20-element RecombinedBSplineBasis of order 8, domain [0.0, 1.0], BCs {left => (D{2}, D{3}, D{4}), right => (D{2}, D{3}, D{4})}
 order: 8
 knots: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0481877, 0.105535  …  0.878831, 0.951475, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-0.0160231, 0.139829, 0.164684, 0.876983, -0.080688, 1.21806, 0.476603, 1.28948, 0.578772, 1.14528, 0.980904, 0.963455, 0.950701, 0.825614, 0.730136, 0.698789, 0.338693, 0.409011, 0.132431, 0.0105342]
 interpolation points: [0.0, 0.0481877, 0.105535, 0.1549, 0.228305, 0.251709, 0.311103, 0.369982, 0.394633, 0.483717, 0.53714, 0.580818, 0.63676, 0.699124, 0.740518, 0.780612, 0.848951, 0.878831, 0.951475, 1.0]

Let's look at the result:

scatter(xs, ys; label = "Data", color = :black)
lines!(0..1, x -> S(x); label = "k = $k (original)", linewidth = 2)
lines!(0..1, x -> Snat(x); label = "k = $k (natural)", linestyle = :dash, linewidth = 4)
axislegend()

Clearly, the spurious oscillations are strongly suppressed near the boundaries.

Extrapolations

One can use extrapolation to evaluate splines outside of their domain of definition. A few different extrapolation strategies are implemented in BSplineKit. See Extrapolation methods for details.

Below we compare a few possible extrapolation strategies included in BSplineKit.

First, we generate and interpolate new data:

xs = 0:0.2:1
ys = 2 * cospi.(xs)
S = interpolate(xs, ys, BSplineOrder(4))

SplineInterpolation containing the 6-element Spline{Float64}:
 basis: 6-element BSplineBasis of order 4, domain [0.0, 1.0]
 order: 4
 knots: [0.0, 0.0, 0.0, 0.0, 0.4, 0.6, 1.0, 1.0, 1.0, 1.0]
 coefficients: [2.0, 2.02957, 1.10398, -1.10398, -2.02957, -2.0]
 interpolation points: 0.0:0.2:1.0

One can directly evaluate these interpolations outside of the domain $[0, 1]$, but the result will always be zero:

S(-0.32)

0.0

To enable extrapolations, one must call extrapolate with the desired extrapolation strategy (see Extrapolation methods for a list). Here we compare both Flat and Smooth methods:

E_flat   = extrapolate(S, Flat())
E_smooth = extrapolate(S, Smooth())

SplineExtrapolation containing the 6-element Spline{Float64}:
 basis: 6-element BSplineBasis of order 4, domain [0.0, 1.0]
 order: 4
 knots: [0.0, 0.0, 0.0, 0.0, 0.4, 0.6, 1.0, 1.0, 1.0, 1.0]
 coefficients: [2.0, 2.02957, 1.10398, -1.10398, -2.02957, -2.0]
 extrapolation method: Smooth()

fig = Figure(resolution = (600, 400))
ax = Axis(fig[1, 1])
scatter!(ax, xs, ys; label = "Data", color = :black)
lines!(ax, -0.5..1.5, x -> S(x); label = "No extrapolation", linewidth = 2)
lines!(ax, -0.5..1.5, x -> E_smooth(x); label = "Smooth", linestyle = :dash, linewidth = 2)
lines!(ax, -0.5..1.5, x -> E_flat(x); label = "Flat", linestyle = :dot, linewidth = 2)
axislegend(ax)
fig

This page was generated using Literate.jl.