Mathematica BSplineCurve in Julia

I am trying to achieve the same result as BSplineCurve—Wolfram Language Documentation in julia. I have used Dierckx.jl for more straightforward spline fitting but am stumped on how to place data points and knots to achieve this result.

image

Any help would be appreciated.

There is a more recent package, BasicBSpline.jl, with docs at https://hyrodium.github.io/BasicBSpline.jl/dev/, by @hyrodium. At the moment it only provides the metods to derive the objects of interest to plot a B-Spline manifold (i.e. a curve, surface or solid), and the user should define a function to get points on such a manifold and use his/her favourite graphics library to plot that manifold. I defined a function which returns a number of points on a B-spline curve, and plotted the curve as a PlotlyJS.jl figure:

using BasicBSpline, PlotlyJS

function bsplcurve(p::Int, knot:: KnotVector{Float64}, ctrl:: Vector{Vector{T}}, N::Int) where T
    """
    p: degree of B-spline basis functions
    knot: the sequence of m knots, m-p >p 
    ctrl: the vector of control points; length(ctrl)= n =  m-p-1
    N: the number of points to be computed on the B-spline curve
    returns N points on the B-spline curve defined by the knots and ctrl
    """
    m = length(knot)
    BS = BSplineSpace{p}(knot) # B-spline space
    n  = dim(BS)
    @assert length(ctrl)==n
    I = domain(BS)
    a, b =extrema(I)
    t = LinRange(a, b, N)
    bs_curve= Vector{Float64}[]
    for s in t
        basis= [bsplinebasis(BS, i, s) for i in 1:n]
        push!(bs_curve, sum(ctrl[i,]*basis[i] for i in 1:n))
    end
    return bs_curve
end;   
function figbspline(ctrl:: Vector{Vector{T}}, bs_pts::Vector{Vector{Float64}}; 
                    width=600, height=600) where T
    # function to display a B-spline curve as a PlotlyJS.jl plot
    fig = Plot(scatter(x=[v[1] for v in ctrl], y =[v[2] for v in ctrl], 
                       name="ctrl polygon", marker_size=8))
    addtraces!(fig, scatter(x= [v[1] for v in bs_pts], y=[v[2] for v in bs_pts],
                            name="b-spline", mode="lines", line_color="red"))
    relayout!(fig, width=width, height=height) 
    return fig
end;

p = 3 #degree of BS polynomials
#end knots must have multiplicity p+1 to ensure that the B-spline curve contains the end control points
knot = KnotVector(1, 1, 1, 1, 1.5, 2.3, 3., 3.5, 3.5, 3.5, 3.5) 
m = length(knot)
#Give m-p-1  control points
#println(m-p-1)
ctrl=[[0., 0.45], [-0.5, 0.75], [-1.,2.0], [0.5, 3.2], 
      [1.6, 1.0], [1.3, -0.3], [0.75, -1.0]];
bs_pts = bsplcurve(p, knot, ctrl, 200)   
fig=figbspline(ctrl, bs_pts; height=450)
display(fig)

cubic-B-spline

1 Like

I think this one supports them.

Using Dierckx’s weighted parametric smoothing splines the following result can be obtained (edited knots and weights):

CODE
using Dierckx, Plots; gr()

x, y = 0:5,  [0, 1, -1, 0, -2, 1]

n = length(x)
δ = 0.05
x2 = [x[1]; (1-δ)*x[1] + δ*x[2]; x[2:end-1]; δ*x[end-1] + (1-δ)*x[end]; x[end]]
y2 = [y[1]; (1-δ)*y[1] + δ*y[2]; y[2:end-1]; δ*y[end-1] + (1-δ)*y[end]; y[end]]
l = [1; 1+δ; 2:n-1; n-δ; n]
w = ones(n+2)
w[1] = w[2]  = w[end-1] = w[end] = 100
spl = ParametricSpline(l, [x2 y2]', w=w, s=1.5) 
xy  = evaluate(spl, range(1, n, 100))
plot(x, y, marker=:circle, ms=3, c=:lime, mc=:red, frame=:none)
plot!(xy[1,:], xy[2,:], c=:black)
2 Likes

@rafael.guerra Your example seems to be the graph of a function. Is it possible to define a parametric B-spline curve (not necessarily a graph), with Dierkx, similar as with scipy?

@empet, the code should be general and can be extended to 3D parametric curves too. The reason why the example looks like a “graph” of a function is because the input data used honors OP’s Mathematica link provided. See an example here.

In the example at the link you posted, the given points are interpolated. I asked for free form B-spline curves, not interpolatory ones. A free form B-spline approximates, i e imitates the form of the control polygon (the line that connects consecutive control points). I came across BasicBSpline.jl, when searched for free form B- splines.

This is very close to what I was hoping to achieve. Is it possible to ensure the resulting spline is tangent to the input line segments at the end points as in empet’s solution?

Thank you for your solution! This is what I was hoping to achieve. Can you explain how you chose your knot values? I may need to read a bit more and explore the math underlying the bspline in order to make this work for me in general.

@arnief3, see edited post above with revised knots and weights.

With BasicBSpline.jl and BasicBSplineExpoerter.jl, you can draw the B-spline curve like this:

using StaticArrays
using BasicBSpline
using BasicBSplineExporter

pts = [SVector(0, 0), SVector(1, 1), SVector(2, -1), SVector(3, 0), SVector(4, -2), SVector(5, 1)]
p = 3
k = KnotVector(0:3)+p*KnotVector(0,3)
P = BSplineSpace{p}(k)
M = BSplineManifold(pts, (P,))
save_png("a.png", M, xlims=(-1,6), ylims=(-3,2))

1 Like

Can you explain how you chose your knot values?

@arnief3 The minimal number of control points of B-spline curve is 4. (assuming the polynomial degree is 3.) The B-spline curve with 4 control points is almost the same as a Bézier curve, and its domain is closed interval [0,1]. We have 6 control points here, so we need to add two knots and replace (0,1) with (0,1,2,3).

We would like to have the endpoints of the curve are located in the same place as the ends of the control points, so we need to duplicated knots in the endpoints. This will replace (0,1,2,3) to (0,0,0,0,1,2,3,3,3,3). This type of knot vector is called “open knot vector”. (B-splines)

WIth BasicBSpline.jl, the addition+ and multiplication* is supported for modifying knot vectors:

julia> using BasicBSpline

julia> p = 3
3

julia> KnotVector(0:3)
KnotVector([0.0, 1.0, 2.0, 3.0])

julia> KnotVector(0:3) + p*KnotVector(0,3)
KnotVector([0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0])

The following books are recommended for learning B-splines.

Thank you very much for introducing my package!

1 Like

Yuto,
Thank you for the outstanding package and the resources. My ultimate goal is to parameterize the resulting bspline (similar to what I can achieve with the ParametricSpline function in the Dierckx.jl package). Is there a recommended way to convert the BSplineManifold object to points along the resulting curve?

Thanks for your help,
Arnold

1 Like

The BSplineManifold instance M is a function-like object, so you can input parameter like M(1.3):

julia> using StaticArrays

julia> using BasicBSpline

julia> using BasicBSplineExporter

julia> pts = [SVector(0, 0), SVector(1, 1), SVector(2, -1), SVector(3, 0), SVector(4, -2), SVector(5, 1)]
6-element Vector{SVector{2, Int64}}:
 [0, 0]
 [1, 1]
 [2, -1]
 [3, 0]
 [4, -2]
 [5, 1]

julia> p = 3
3

julia> k = KnotVector(0:3)+p*KnotVector(0,3)
KnotVector([0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0])

julia> P = BSplineSpace{p}(k)
BSplineSpace{3, Float64}(KnotVector([0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0]))

julia> M = BSplineManifold(pts, (P,))
BSplineManifold{1, (3,), SVector{2, Int64}, Tuple{BSplineSpace{3, Float64}}}((BSplineSpace{3, Float64}(KnotVector([0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0])),), SVector{2, Int64}[[0, 0], [1, 1], [2, -1], [3, 0], [4, -2], [5, 1]])

julia> M(0)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0

julia> M(3)
2-element SVector{2, Float64} with indices SOneTo(2):
 5.0
 1.0

julia> M(-1)
ERROR: DomainError with -1:
The input -1 is out of range.
Stacktrace:
 [1] (::BSplineManifold{1, (3,), SVector{2, Int64}, Tuple{BSplineSpace{3, Float64}}})(t1::Int64)
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_BSplineManifold.jl:86
 [2] top-level scope
   @ REPL[13]:1

julia> M.(0:0.1:3)
31-element Vector{SVector{2, Float64}}:
 [0.0, 0.0]
 [0.2854166666666667, 0.2426666666666667]
 [0.5433333333333336, 0.38133333333333347]
 [0.7762499999999998, 0.43199999999999994]
 [0.9866666666666668, 0.41066666666666674]
 [1.1770833333333335, 0.3333333333333333]
 [1.35, 0.21599999999999997]
 [1.5079166666666666, 0.07466666666666671]
 [1.6533333333333333, -0.07466666666666677]
 [1.7887500000000003, -0.21599999999999997]
 [1.9166666666666665, -0.33333333333333326]
 [2.0393333333333334, -0.4146666666666668]
 [2.158, -0.46399999999999997]
 [2.2736666666666667, -0.48933333333333334]
 [2.3873333333333333, -0.49866666666666665]
 [2.5, -0.5]
 [2.6126666666666667, -0.5013333333333332]
 [2.7263333333333333, -0.5106666666666666]
 [2.8419999999999996, -0.536]
 [2.9606666666666666, -0.5853333333333333]
 [3.083333333333333, -0.6666666666666666]
 [3.21125, -0.7820000000000001]
 [3.3466666666666667, -0.9093333333333335]
 [3.4920833333333334, -1.0206666666666666]
 [3.65, -1.0879999999999999]
 [3.8229166666666665, -1.0833333333333333]
 [4.013333333333334, -0.9786666666666665]
 [4.223750000000001, -0.7459999999999996]
 [4.456666666666666, -0.35733333333333417]
 [4.714583333333333, 0.2153333333333327]
 [5.0, 1.0]
1 Like

Using all of the input from this thread, here is the solution I was hoping to find. It gives me the curve I want as a parametric function with input range [0,1]. It seems to work for all p and number of points assuming length(pts) - p > 0. It also works for 2d and 3d points. Thanks again for everyone’s help.

using StaticArrays
using BasicBSpline

function bsplcurve(p, pts)
    max_k = length(pts)-p
    @assert max_k > 0
    k = KnotVector(0:max_k)+p*KnotVector(0,max_k)
    P = BSplineSpace{p}(k)
    M = BSplineManifold(pts, (P,))
    return t->M(t*max_k)
end

pts = [SVector(0, 0), SVector(1, 1), SVector(2, -1), SVector(3, 0), SVector(4, -2), SVector(5, 1)]
p = 2

bspl = bsplcurve(p, pts)
bspl_pts = [bspl(i) for i = range(0, 1, 100)]
1 Like

@arnief3
If the knot vector k is defined like this, then we don’t need the scaling t->M(t*max_k).

using StaticArrays
using BasicBSpline

function bsplcurve(p, pts)
    @assert length(pts) > p
    k = KnotVector(range(0,1,length(pts)-p+1))+p*KnotVector(0,1)
    P = BSplineSpace{p}(k)
    M = BSplineManifold(pts, (P,))
    return M
end

pts = [SVector(0, 0), SVector(1, 1), SVector(2, -1), SVector(3, 0), SVector(4, -2), SVector(5, 1)]
p = 2

bspl = bsplcurve(p, pts)
bspl_pts = [bspl(i) for i = range(0, 1, 100)]
1 Like