sostock/bsplines.jl

A Julia package for working with B-splines

BSplines

This package provides data types and functions for working with B-splines as a means to approximate real functions. For its usage, see the documentation.

This package is not related to the package at https://github.com/gusl/BSplines.jl. While gusl/BSplines.jl is only compatible with Julia ≤ 0.4, this package works with Julia ≥ 1.0.

Installation

This package is compatible with Julia ≥ 1.0. It can be installed by typing

] add BSplines

in the Julia REPL.

Acknowledgments

The algorithms used for evaluating B-splines and their derivatives are adapted from the Fortran code found in Carl de Boor’s book A practical Guide to Splines [1], in particular from the subroutines BSPLVB, BSPLVD and BVALUE. This package is published with his permission.

[1] Carl de Boor, A Practical Guide to Splines, New York, N.Y.: Springer-Verlag, 1978.

bsplines.jl's Issues

Number of spline basis functions

Hi

Thank you for this package. I consider using it in another package, but there is a thing I don't understand -- from simple examples I expect to have fewer functions.

Consider the knot vector

using BSplines
knots = [0; 0; 0; 1; 1; 1]
basis = BSplineBasis(1, knots)

Here, length(basis) returns 5 as expected and plotting the function look right. But when increasing the order I run into trouble

basis = BSplineBasis(2, knots)

Here, length(basis) returns 6, but I expect 4. Consider the output of plot(basis)

I would index the functions from 0 through 3. I expect the 0th and 3rd functions to be constantly 0. The 1st and 2nd should be the two straight lines.
What are the rest?

This becomes a problem when constructing a Spline because it insists on having a larger coeffs vector than the number of basis functions I anticipate.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Slightly confused about the notation used in this package

The order k in this package appears to differ from the corresponding usage in other packages such as Dierckx or in scipy.interpolate (which also uses Dierckx I believe). It appears that k here is somewhat equivalent to k+1 in these packages. For example, in scipy:

In [1]: from scipy import interpolate

In [2]: knots = [0,0,1,2,3,4,4]; k = 2; c = [1,0,0,0,0];

In [3]: interpolate.BSpline(knots, c, k)(1)
Out[3]: array(0.5)

In [4]: c = [0,1,0,0,0];

In [5]: interpolate.BSpline(knots, c, k)(1)
Out[5]: array(0.5)

In Dierckx,

julia> k = 2; breakpoints = 0:4;

julia> t = vcat([breakpoints[1] for i in 1:k], breakpoints[2:end-1], [breakpoints[end] for i = 1:k])
7-element Array{Int64,1}:
 0
 0
 1
 2
 3
 4
 4

julia> c = [1,0,0,0,0];

julia> Dierckx.Spline1D(t, c, k, 0, 0)(1) # matches with scipy
0.5

# comparing with this package
julia> basis = BSplineBasis(k, breakpoints); knots(basis) == t
true

whereas in this package

julia> basis = BSplineBasis(2, 0:4)
5-element BSplineBasis{UnitRange{Int64}}:
 order: 2
 breakpoints: 0:4

julia> knots(basis)
7-element BSplines.KnotVector{Int64,UnitRange{Int64}}:
 0
 0
 1
 2
 3
 4
 4

julia> bsplines(basis, 1)
2-element OffsetArray(::Array{Float64,1}, 2:3) with eltype Float64 with indices 2:3:
 1.0
 0.0

julia> basis = BSplineBasis(3, 0:4)
6-element BSplineBasis{UnitRange{Int64}}:
 order: 3
 breakpoints: 0:4

julia> knots(basis)
9-element BSplines.KnotVector{Int64,UnitRange{Int64}}:
 0
 0
 0
 1
 2
 3
 4
 4
 4

julia> bsplines(basis, 1)
3-element OffsetArray(::Array{Float64,1}, 2:4) with eltype Float64 with indices 2:4:
 0.5
 0.5
 0.0

Using order=3 in this package leads to values that seem to match the results for k=2 in scipy, although the knots differ as expected at the edges. I wonder if these differences could be explained in the documentation?

Roadmap

• Simplify tests. Currently, running the tests takes a lot of time. We probably don’t need all of them. (#3)
• Allow knot vectors where the first and last breakpoints are repeated less than k times. (#5)
  • Add a new constructor with mandatory keyword argument:
    • BSplineBasis(k, knots=...) for supplying the knot sequence directly.
    • BSplineBasis(k, breakpoints=...) which creates the knot vector from the breakpoints by repeating the first and last breakpoint k times, equivalent to the old constructor without keyword argument.
  • Indexing a BSplineBasis with an AbstractUnitRange creates a new BSplineBasis with an appropriate knot vector. view(::BSplineBasis, ::AbstractUnitRange) creates an equivalent BSplineBasis whose knot vector is a view. The indices keyword argument can then be removed.
  • Since the BSplineBasis is changed to contain the knot vector instead of the breakpoints, breakpoints(basis) is changed to return unique(knots(basis)). Alternatively, the breakpoints function can be removed completely, since the user can just call unique(knots(basis)) themself.
• Add a StaticBSplineBasis{K, T} type where the order K is a type parameter for potentially better performance.
  • Specialized implementations for small K are possible.
  • Add AbstractBSplineBasis{T} as supertype for BSplineBasis{T} and StaticBSplineBasis{K, T}.