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.107311, 0.159983, 0.213316, 0.268371, 0.321046, 0.36973 … 0.625816, 0.681882, 0.732199, 0.783289, 0.837281, 0.886254, 1.0, 1.0, 1.0, 1.0]
coefficients: [0.0183804, 0.15391, 0.315331, 0.496405, 0.598167, 0.712468, 0.930592, 0.875045, 0.990095, 1.00157, 0.981121, 1.0346, 0.933992, 0.850968, 0.735805, 0.666223, 0.459099, 0.407906, -0.123718, -0.0129226]
interpolation points: [0.0, 0.0518416, 0.109298, 0.160794, 0.209856, 0.269298, 0.325959, 0.367882, 0.415349, 0.465107, 0.524201, 0.570411, 0.635975, 0.671061, 0.738611, 0.786926, 0.824332, 0.900587, 0.933845, 1.0]
Let's plot the result:
lines!(0..1, S; label = "k = 4", color = Cycled(4 - 3))
We can also choose other interpolation orders for comparison:
for k ∈ (5, 6, 8)
local S = interpolate(xs, ys, BSplineOrder(k))
lines!(0..1, S; 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 BCsSplineInterpolation 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.0518416, 0.109298 … 0.900587, 0.933845, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
coefficients: [0.0137583, 0.142842, 0.381383, 0.482671, 0.653978, 0.576594, 1.10858, 0.718314, 1.09114, 0.975321, 0.940988, 1.12046, 0.863254, 0.931811, 0.613254, 0.83853, 0.206394, 0.679701, 0.046989, -0.0095833]
interpolation points: [0.0, 0.0518416, 0.109298, 0.160794, 0.209856, 0.269298, 0.325959, 0.367882, 0.415349, 0.465107, 0.524201, 0.570411, 0.635975, 0.671061, 0.738611, 0.786926, 0.824332, 0.900587, 0.933845, 1.0]
Let's look at the result:
scatter(xs, ys; label = "Data", color = :black)
lines!(0..1, S; label = "k = $k (original)", linewidth = 2)
lines!(0..1, Snat; label = "k = $k (natural)", linestyle = :dash, linewidth = 4)
axislegend()
Clearly, the spurious oscillations are strongly suppressed near the boundaries.
Smoothing cubic splines
One can use smoothing splines to fit noisy data. A smoothing spline is a curve which passes close to the input data, while avoiding strong fluctuations due to possible noise. The smoothing strength is controlled by a regularisation parameter $λ$. Setting $λ = 0$ corresponds to a regular interpolation (the obtained spline passes through all the points), while increasing $λ$ leads to a smoother curve which roughly approximates the data.
Given a set of data points $(x_i, y_i)$, the idea is to construct a spline $S(x)$ that minimises:
\[∑_{i = 1}^N w_i |y_i - S(x_i)|^2 + λ ∫_{x_1}^{x_N} \left[ S''(x) \right]^2 \, \mathrm{d}x\]
Here $w_i$ are optional weights that may be used to give "priority" to certain data points.
Note that only cubic splines (order $k = 4$) are currently supported.
rng = MersenneTwister(42)
Ndata = 20
xs = sort!(rand(rng, Ndata))
xs[begin] = 0; xs[end] = 1; # not strictly necessary; just to set the data limits
ys = cospi.(2 .* xs) .+ 0.04 .* randn(rng, Ndata);Create smoothing spline from data:
λ = 1e-3
S_fit = fit(xs, ys, λ)20-element Spline{Float64}:
basis: 20-element RecombinedBSplineBasis of order 4, domain [0.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
order: 4
knots: [0.0, 0.0, 0.0, 0.0, 0.0644853, 0.107498, 0.143591, 0.177709, 0.194236, 0.208728 … 0.561824, 0.69722, 0.710824, 0.734579, 0.802455, 0.806648, 1.0, 1.0, 1.0, 1.0]
coefficients: [0.95914, 0.66248, 0.76206, 0.583945, 0.432862, 0.312774, 0.0755456, -0.187885, -0.481406, -0.705676, -0.841713, -0.887651, -0.91701, -0.828688, -0.49396, -0.210413, -0.0337114, 0.130158, 0.415708, 0.805077]If we want the spline to pass very near a single data point, we can assign a larger weight to that point:
weights = fill!(similar(xs), 1)
weights[12] = 100 # larger weight to point i = 12
S_fit_weight = fit(xs, ys, λ; weights)20-element Spline{Float64}:
basis: 20-element RecombinedBSplineBasis of order 4, domain [0.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
order: 4
knots: [0.0, 0.0, 0.0, 0.0, 0.0644853, 0.107498, 0.143591, 0.177709, 0.194236, 0.208728 … 0.561824, 0.69722, 0.710824, 0.734579, 0.802455, 0.806648, 1.0, 1.0, 1.0, 1.0]
coefficients: [0.961671, 0.664922, 0.76546, 0.586503, 0.433923, 0.312021, 0.0703509, -0.201631, -0.514269, -0.758976, -0.914636, -0.965504, -0.995151, -0.881144, -0.516941, -0.219777, -0.0360762, 0.131485, 0.421835, 0.80903]Plot results and compare with natural cubic spline interpolation:
S_interp = interpolate(xs, ys, BSplineOrder(4), Natural())
scatter(xs, ys; label = "Data", color = :black)
lines!(0..1, S_interp; label = "Interpolation", linewidth = 2)
lines!(0..1, S_fit; label = "Fit (λ = $λ)", linewidth = 2)
lines!(0..1, S_fit_weight; label = "Fit (λ = $λ) with weight", linestyle = :dash, linewidth = 2)
axislegend(position = (0.5, 1))
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.2:0.2:1.2
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.2, 1.2]
order: 4
knots: [0.2, 0.2, 0.2, 0.2, 0.6, 0.8, 1.2, 1.2, 1.2, 1.2]
coefficients: [1.61803, 1.13682, -0.301508, -2.08779, -2.14709, -1.61803]
interpolation points: 0.2:0.2:1.2
One can directly evaluate these interpolations outside of the domain $[0, 1]$, but the result will always be zero:
S(-0.32)0.0To 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_linear = extrapolate(S, Linear())
E_smooth = extrapolate(S, Smooth())SplineExtrapolation containing the 6-element Spline{Float64}:
basis: 6-element BSplineBasis of order 4, domain [0.2, 1.2]
order: 4
knots: [0.2, 0.2, 0.2, 0.2, 0.6, 0.8, 1.2, 1.2, 1.2, 1.2]
coefficients: [1.61803, 1.13682, -0.301508, -2.08779, -2.14709, -1.61803]
extrapolation method: Smooth()fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1])
scatter!(ax, xs, ys; label = "Data", color = :black)
lines!(ax, -0.5..1.5, S; label = "No extrapolation", linewidth = 2)
lines!(ax, -0.5..1.5, E_smooth; label = "Smooth", linestyle = :dash, linewidth = 2)
lines!(ax, -0.5..1.5, E_linear; label = "Linear", linestyle = :dashdot, linewidth = 2)
lines!(ax, -0.5..1.5, E_flat; label = "Flat", linestyle = :dot, linewidth = 2)
axislegend(ax)
fig
Periodic data
It is also possible to interpolate or fit data which is expected to result from a periodic function, such that $f(x + L) = f(x)$ for some period $L$. For this, one can pass Periodic(L) as a boundary condition to interpolate or fit.
The following example starts from data at points $x_j$ within the $[-1, 1]$ interval, and assumes the resulting spline can be extended periodically outside of this interval.
We start by generating some data:
N = 40
f_slow(x) = cospi(x)
f_fast(x) = 0.2 * sinpi(40x)
xs = [-cospi(n / N) for n = 0:(N - 1)] # in [-1, 1) // NOTE: the endpoint (x = 1) must be excluded!!
ys = @. f_slow(xs) + f_fast(xs);Interpolate the data:
L = 2
S_interp = interpolate(xs, copy(ys), BSplineOrder(4), Periodic(L))SplineInterpolation containing the 40-element Spline{Float64}:
basis: 40-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
order: 4
knots: [..., -1.01231, -1.00308, -1.0, -0.996917, -0.987688, -0.97237, -0.951057, -0.92388, -0.891007, -0.85264 … 0.809017, 0.85264, 0.891007, 0.92388, 0.951057, 0.97237, 0.987688, 0.996917, 1.0, 1.00308, ...]
coefficients: [..., -1.00002, -0.871461, -0.698456, -1.24032, -0.914953, -1.10651, -0.590139, -1.0634, -1.00603, -0.929873 … -0.428381, -0.529331, -0.644522, -0.726095, -1.2938, -0.836605, -1.06145, -0.752128, -1.30002, -1.12841, ...]
interpolation points: [-1.0, -0.996917, -0.987688, -0.97237, -0.951057, -0.92388, -0.891007, -0.85264, -0.809017, -0.760406 … 0.707107, 0.760406, 0.809017, 0.85264, 0.891007, 0.92388, 0.951057, 0.97237, 0.987688, 0.996917]
Create a periodic cubic smoothing spline. Note that BSplineOrder(4) is assumed (it's currently the only supported choice). We also compare with a smoothing spline which doesn't assume periodic boundary conditions.
λ = 0.001 # smoothing parameter
S_fit_natural = fit([xs; xs[begin] + L], [ys; ys[begin]], λ) # for comparison, compute a natural spline (no implied periodicity)
S_fit_periodic = fit(xs, ys, λ, Periodic(L))40-element Spline{Float64}:
basis: 40-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
order: 4
knots: [..., -1.01231, -1.00308, -1.0, -0.996917, -0.987688, -0.97237, -0.951057, -0.92388, -0.891007, -0.85264 … 0.809017, 0.85264, 0.891007, 0.92388, 0.951057, 0.97237, 0.987688, 0.996917, 1.0, 1.00308, ...]
coefficients: [..., -0.998215, -0.999056, -1.00003, -1.00021, -0.997255, -0.988569, -0.972036, -0.947159, -0.901331, -0.81504 … -0.529776, -0.637428, -0.74265, -0.836333, -0.906633, -0.949966, -0.975088, -0.988484, -0.994806, -0.997216, ...]Plot the results:
fig = Figure()
ax = Axis(fig[1, 1])
scatter!(ax, xs, ys; label = "Data")
lines!(ax, -1..1, S_interp; label = "Interpolation", color = (:grey, 0.5))
lines!(ax, -1..1, S_fit_periodic; linewidth = 3, label = "Smoothing (periodic)")
lines!(ax, -1..1, S_fit_natural; linewidth = 3, linestyle = :dash, label = "Smoothing (natural)")
axislegend(ax)
fig
As can be expected, the smoothing spline with periodic boundary conditions mostly differs from the "natural" smoothing spline near the boundaries. In particular, the natural smoothing spline simply does not satisfy periodic boundary conditions (even at the level of the zero-th derivative!), so it's not very adapted for constructing periodic functions.
This page was generated using Literate.jl.