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.
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)
https://github.com/JuliaMath/Interpolations.jl
I think this one supports them.
Using Dierckx’s weighted parametric smoothing splines the following result can be obtained (edited knots and weights):
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)
@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))
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!
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
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]
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)]
@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)]
Thanks! I need this code today