Derivatives of B-Spline w.r.t their degrees of freedom

[touste]

Hi all,

I’m trying to approximate a shape described by a set of 3D points with a smooth interpolant, I’ve settled on a parametric B-Spline.

I’ve got a set of 3D points P_i:

using GeometryBasics
P₁ = Point3(0.0, 0.0, 0.0)
P₂ = Point3(1.0, 1.0, 0.1)
P₃ = Point3(1.0, 1.2, 0.5)
P₄ = Point3(2.0, 0.9, 3.0)

and I’m using Interpolations.jl to build the interpolant as:

using Interpolations
A = [P₁, P₂, P₃, P₄]
Q = interpolate(A, BSpline(Cubic(Free(OnGrid()))))

Then I’m able to use the interpolation function Q(t), t \in [1,4]

Q(1.5)

However, I’m interested in using this interpolating approximation Q(t) within an optimization problem, by changing the degrees of freedom of the B-Spline until I’ve minimized a functional. As a consequence, I’m looking for the jacobian of a given set of points on the curve w.r.t the degrees of freedom of the B-Spline \partial Q(t)/\partial P_i.
I can get this using FiniteDiff.jl as:

using FiniteDiff

function interp(P₁, P₂, P₃, P₄, t)
	A = [P₁, P₂, P₃, P₄]
	itp = interpolate(A, BSpline(Cubic(Free(OnGrid()))))
	return itp(t)
end

# dQ(1.5)/dP₁
FiniteDiff.finite_difference_jacobian(x->interp(x, P₂, P₃, P₄, 1.5), P₁)

However I’m struggling a bit to simplify this process and get the analytical derivatives.

Usually, B-Splines are defined using a set of control points P'_i that are different from the interpolating points (knots) P_i:

Q(t) = \sum P'_i N_{i,k}(t).

It may be wise to use these control points as degrees of freedom instead, as the derivatives are simply the basis functions N_{i,k}. Can I retrieve these control points and basis functions using the Interpolations.jl package? I’ve looked into the documentation but I’m a bit lost by the way interpolations are computed.

I may be on the wrong track here though, I’m not quite sure how to get these derivatives in a easier way.

Thanks!

[Alec_Loudenback]

This isn’t really answering your question directly, but have you seen the package BSplineKit.jl?

I have used that package with ForwardDIff and not had issues getting the derivative from the curve (EconomicScnearioGenerators.jl in places uses the derivative of yield curves from Yields.jl, and those curves use BSplineKit interpolated curves underneath).

[jipolanco]

As @Alec_Loudenback suggests, you should be able to do this using BSplineKit.jl.

Here is how you would get a parametric spline interpolation (the example is in 2D because it’s easier to visualise, but the same works in 3D):

using BSplineKit
using GeometryBasics

P₁ = Point2(0.0, 0.0)
P₂ = Point2(1.0, 1.0)
P₃ = Point2(1.0, 1.2)
P₄ = Point2(0.3, 0.9)

t = 0:3
A = [P₁, P₂, P₃, P₄]

Q = interpolate(t, A, BSplineOrder(4))

using CairoMakie
tfine = 0:0.1:3
fig = Figure()
ax = Axis(fig[1, 1])
lines!(ax, Q.(tfine))
scatter!(ax, A)
fig

Note that I’m interpolating using B-splines of order k = 4, which means cubic splines in the BSplineKit convention. To get the very same result as in Interpolations.jl, you could also try using Natural boundary conditions as explained here, which I’m not doing here to keep things simple.

As you mentioned, and keeping your notation, the interpolating spline can be written as Q(t) = \sum_i P_i' N_i(t), where the N_i(t) are the B-spline basis functions (i \in [1, 4] in the example).

In my example, one can easily evaluate each basis function as:

julia> B = basis(Q)  # get underlying B-spline basis
4-element BSplineBasis of order 4, domain [0.0, 3.0]
 knots: [0.0, 0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 3.0]

julia> B[2](2.13)  # either evaluate B-spline N₂ at point t = 2.13...
0.17913300000000004

julia> B(2.13)   # ...or evaluate all B-splines which are non-zero at t = 2.13
(4, (0.357911, 0.43856700000000004, 0.17913300000000004, 0.024389000000000008))

The last method is more efficient, especially if you want to evaluate all non-zero B-splines at a point. As explained here, what this returns is the index i of the last B-spline, and then the evaluated B-splines in reverse order. In other words, the result above gives you (N_4(t), N_3(t), N_2(t), N_1(t)). I hope this helps!

And in case you need them, the interpolation coefficients P'_i themselves can be obtained as:

julia> coefficients(Q)
4-element Vector{Point2{Float64}}:
 [0.0, 0.0]
 [1.5999999999999996, 1.4999999999999998]
 [1.2500000000000002, 1.3499999999999999]
 [0.3, 0.9]

[touste]

Thank you very much for the detailed answer! BSplineKit.jl seems indeed much better suited in this case, and I was able to retrieve the basis functions as you describe and use the coefficients as degrees of freedom.
Really appreciate it!