B-spline basis function

Basic properties of B-spline basis function

Def. B-spline space

B-spline basis function is defined by Cox–de Boor recursion formula.

\[\begin{aligned} {B}_{(i,p,k)}(t) &= \frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t) +\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\ {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i}\le t< k_{i+1})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

If the denominator is zero, then the term is assumed to be zero.

The next figure shows the plot of B-spline basis functions. You can manipulate these plots on desmos graphing calculator!

Thm. Basis of B-spline space

The set of functions $\{B_{(i,p,k)}\}_i$ is a basis of B-spline space $\mathcal{P}[p,k]$.

julia> p = 22
julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])
julia> P = BSplineSpace{p}(k)BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))
julia> plot([t->bsplinebasis₊₀(P,i,t) for i in 1:dim(P)], 0, 10, ylims=(0,1), label=false)Plot{Plots.PlotlyBackend() n=5}

You can choose the first terms in different ways.

\[\begin{aligned} {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i} < t \le k_{i+1}) \\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

julia> p = 22
julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])
julia> P = BSplineSpace{p}(k)BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))
julia> plot([t->bsplinebasis₋₀(P,i,t) for i in 1:dim(P)], 0, 10, ylims=(0,1), label=false)Plot{Plots.PlotlyBackend() n=5}

In these cases, each B-spline basis function $B_{(i,2,k)}$ is coninuous, so bsplinebasis₊₀ and bsplinebasis₋₀ are equal.

Support of B-spline basis function

Thm. Support of B-spline basis function

If a B-spline space$\mathcal{P}[p,k]$ is non-degenerate, the support of its basis function is calculated as follows:

\[\operatorname{supp}(B_{(i,p,k)})=[k_{i},k_{i+p+1}]\]

[TODO: fig]

BasicBSpline.bsplinesupportFunction

Return the support of $i$-th B-spline basis function.

\[\operatorname{supp}(B_{(i,p,k)})=[k_{i},k_{i+p+1}]\]

Examples

julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])
KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])

julia> P = BSplineSpace{2}(k)
BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))

julia> bsplinesupport(P,1)
0.0..5.5

julia> bsplinesupport(P,2)
1.5..8.0
source

Partition of unity

Thm. Partition of unity

Let $B_{(i,p,k)}$ be a B-spline basis function, then the following equation is satisfied.

\[\begin{aligned} \sum_{i}B_{(i,p,k)}(t) &= 1 & (k_{p+1} \le t < k_{l-p}) \\ 0 \le B_{(i,p,k)}(t) &\le 1 \end{aligned}\]

julia> p = 22
julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])
julia> P = BSplineSpace{p}(k)BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))
julia> plot(t->sum(bsplinebasis₊₀(P,i,t) for i in 1:dim(P)), 0, 10, ylims=(0,1.1), label=false)Plot{Plots.PlotlyBackend() n=1}

To satisfy the partition of unity on whole interval $[1,8]$, sometimes more knots will be inserted to the endpoints of the interval.

julia> p = 22
julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]) + p * KnotVector([0,10])KnotVector([0.0, 0.0, 0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0, 10.0, 10.0])
julia> P = BSplineSpace{p}(k)BSplineSpace{2, Float64}(KnotVector([0.0, 0.0, 0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0, 10.0, 10.0]))
julia> plot(t->sum(bsplinebasis₊₀(P,i,t) for i in 1:dim(P)), 0, 10, ylims=(0,1.1), label=false)Plot{Plots.PlotlyBackend() n=1}

But, the sum $\sum_{i} B_{(i,p,k)}(t)$ is not equal to $1$ at $t=8$. Therefore, to satisfy partition of unity on closed interval $[k_{p+1}, k_{l-p}]$, the definition of first terms of B-spline basis functions are sometimes replaced:

\[\begin{aligned} {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i} \le t<k_{i+1})\\ &1\quad (k_{i} < t = k_{i+1}=k_{l})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

julia> p = 22
julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]) + p * KnotVector([0,10])KnotVector([0.0, 0.0, 0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0, 10.0, 10.0])
julia> P = BSplineSpace{p}(k)BSplineSpace{2, Float64}(KnotVector([0.0, 0.0, 0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0, 10.0, 10.0]))
julia> plot(t->sum(bsplinebasis(P,i,t) for i in 1:dim(P)), 0, 10, ylims=(0,1.1), label=false)Plot{Plots.PlotlyBackend() n=1}
BasicBSpline.bsplinebasis₊₀Function

$i$-th B-spline basis function. Right-sided limit version.

\[\begin{aligned} {B}_{(i,p,k)}(t) &= \frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t) +\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\ {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i}\le t< k_{i+1})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

Examples

julia> P = BSplineSpace{0}(KnotVector(1:6))
BSplineSpace{0, Int64}(KnotVector([1, 2, 3, 4, 5, 6]))

julia> bsplinebasis₊₀.(P,1:5,(1:6)')
5×6 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0
source
BasicBSpline.bsplinebasis₋₀Function

$i$-th B-spline basis function. Left-sided limit version.

\[\begin{aligned} {B}_{(i,p,k)}(t) &= \frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t) +\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\ {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i}< t\le k_{i+1})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

Examples

julia> P = BSplineSpace{0}(KnotVector(1:6))
BSplineSpace{0, Int64}(KnotVector([1, 2, 3, 4, 5, 6]))

julia> bsplinebasis₋₀.(P,1:5,(1:6)')
5×6 Matrix{Float64}:
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0
source
BasicBSpline.bsplinebasisFunction

$i$-th B-spline basis function. Modified version.

\[\begin{aligned} {B}_{(i,p,k)}(t) &= \frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t) +\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\ {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i} \le t<k_{i+1})\\ &1\quad (k_{i} < t = k_{i+1}=k_{l})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}\]

Examples

julia> P = BSplineSpace{0}(KnotVector(1:6))
BSplineSpace{0, Int64}(KnotVector([1, 2, 3, 4, 5, 6]))

julia> bsplinebasis.(P,1:5,(1:6)')
5×6 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0
source

B-spline basis functions at specific point

Sometimes, you may need the non-zero values of B-spline basis functions at specific point. The bsplinebasisall function is much more efficient than evaluating B-spline functions one by one with bsplinebasis function.

julia> using BenchmarkTools, BasicBSpline
julia> P = BSplineSpace{2}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))
julia> t = 6.36.3
julia> (bsplinebasis(P, 2, t), bsplinebasis(P, 3, t), bsplinebasis(P, 4, t))(0.2101818181818182, 0.7166753246753247, 0.0731428571428571)
julia> bsplinebasisall(P, 2, t)3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3): 0.21018181818181822 0.7166753246753247 0.0731428571428571
julia> @benchmark (bsplinebasis($P, 2, $t), bsplinebasis($P, 3, $t), bsplinebasis($P, 4, $t))BenchmarkTools.Trial: 10000 samples with 993 evaluations. Range (minmax): 34.038 ns 2.280 μs GC (min … max): 0.00% … 0.00% Time (median): 38.973 ns GC (median): 0.00% Time (mean ± σ): 42.714 ns ± 29.357 ns GC (mean ± σ): 0.00% ± 0.00% ▁▆▇██ ▇▇ ▂▂▃▃▂▂▁ ▂ ███████▄▅▆▄▆▅▆▇████████▇██▆▆▅▂▅▇▇▇█▇▇▇▆▆▆▆▆▅▅▃▅▅▅▄▃▅▃▅▅▂▄ █ 34 ns Histogram: log(frequency) by time 96.4 ns < Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark bsplinebasisall($P, 2, $t)BenchmarkTools.Trial: 10000 samples with 999 evaluations. Range (minmax): 8.809 ns241.143 ns GC (min … max): 0.00% … 0.00% Time (median): 9.910 ns GC (median): 0.00% Time (mean ± σ): 10.905 ns ± 6.521 ns GC (mean ± σ): 0.00% ± 0.00% ▅██▇▆ ▁ ▂ █████▃▅▄▃▄▃▃▃▃▁▁▁▁▄▃▃▁▁▁▁▁▃▁▃▃▁▁▄▁▁▄▆▁▁▁▄▇██████▇▇▆▆▃▆▅▅▅▅ █ 8.81 ns Histogram: log(frequency) by time 31.5 ns < Memory estimate: 0 bytes, allocs estimate: 0.
BasicBSpline.intervalindexFunction

Return an index of a interval in the domain of B-spline space

Examples

julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]);

julia> P = BSplineSpace{2}(k);

julia> domain(P)
2.5..9.0

julia> intervalindex(P,2.6)
1

julia> intervalindex(P,5.6)
2

julia> intervalindex(P,8.5)
3

julia> intervalindex(P,9.5)
3
source
BasicBSpline.bsplinebasisallFunction

B-spline basis functions at point t on i-th interval.

Examples

julia> k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])
KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])

julia> p = 2
2

julia> P = BSplineSpace{p}(k)
BSplineSpace{2, Float64}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))

julia> t = 5.7
5.7

julia> i = intervalindex(P,t)
2

julia> bsplinebasisall(P,i,t)
3-element SVector{3, Float64} with indices SOneTo(3):
 0.3847272727272727
 0.6107012987012989
 0.00457142857142858

julia> bsplinebasis.(P,i:i+p,t)
3-element Vector{Float64}:
 0.38472727272727264
 0.6107012987012989
 0.00457142857142858
source

The next figures illustlates the relation between domain(P), intervalindex(P,t) and bsplinebasisall(P,i,t).

k = KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0])

for p in 1:3
    P = BSplineSpace{p}(k)
    plot(P, legend=:topleft, label="B-spline basis (p=1)")
    plot!(t->intervalindex(P,t),0,10, label="Interval index")
    plot!(t->sum(bsplinebasis(P,i,t) for i in 1:dim(P)),0,10, label="Sum of B-spline basis")
    scatter!(k.vector,zero(k.vector), label="knot vector")
    plot!([t->bsplinebasisall(P,1,t)[i] for i in 1:p+1],0,10, color=:black, label="bsplinebasisall (i=1)", ylim=(-1,8-2p))
end

Uniform B-spline basis and uniform distribution

Let $X_1, \dots, X_n$ be i.i.d. random variables with $X_i \sim U(0,1)$, then the probability density function of $X_1+\cdots+X_n$ can be obtained via BasicBSpline.uniform_bsplinebasis_kernel(Val(n-1),t).

N = 100000
# polynomial degree 0
plot1 = histogram([rand() for _ in 1:N], normalize=true, label=false)
plot!(t->BasicBSpline.uniform_bsplinebasis_kernel(Val(0),t), label=false)
# polynomial degree 1
plot2 = histogram([rand()+rand() for _ in 1:N], normalize=true, label=false)
plot!(t->BasicBSpline.uniform_bsplinebasis_kernel(Val(1),t), label=false)
# polynomial degree 2
plot3 = histogram([rand()+rand()+rand() for _ in 1:N], normalize=true, label=false)
plot!(t->BasicBSpline.uniform_bsplinebasis_kernel(Val(2),t), label=false)
# polynomial degree 3
plot4 = histogram([rand()+rand()+rand()+rand() for _ in 1:N], normalize=true, label=false)
plot!(t->BasicBSpline.uniform_bsplinebasis_kernel(Val(3),t), label=false)
# plot all
plot(plot1,plot2,plot3,plot4)
"/home/runner/work/BasicBSpline.jl/BasicBSpline.jl/docs/build/histogram-uniform.html"