# 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 = 22julia> 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{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 = 22julia> 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{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{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 = 22julia> 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{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 = 22julia> 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{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 julia> p = 22julia> 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{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{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{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{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
BasicBSpline.bsplinebasis₋₀IFunction

$i$-th B-spline basis function. Modified version (2).

\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_{1+p} \le k_{i} < t \le k_{i+1}) \\ &1\quad (t = k_{1+p} = k_{i} < k_{i+1}) \\ &1\quad (k_{i} \le t < k_{i+1} \le k_{1+p}) \\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}

Examples

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

julia> BasicBSpline.bsplinebasis₋₀I.(P,1:5,(1:6)')
5×6 Matrix{Float64}:
1.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

## 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, BasicBSplinejulia> P = BSplineSpace{2}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))BSplineSpace{2, Float64, KnotVector{Float64}}(KnotVector([0.0, 1.5, 2.5, 5.5, 8.0, 9.0, 9.5, 10.0]))julia> t = 6.36.3julia> (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.0731428571428571julia> @benchmark (bsplinebasis($P, 2,$t), bsplinebasis($P, 3,$t), bsplinebasis($P, 4,$t))BenchmarkTools.Trial: 10000 samples with 996 evaluations.
Range (min … max):  25.389 ns … 50.938 ns  ┊ GC (min … max): 0.00% … 0.00%
Time  (median):     25.701 ns              ┊ GC (median):    0.00%
Time  (mean ± σ):   25.771 ns ±  0.851 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

█   ▆▁
▂▂▂▁▁▂▂▂▂▃▃███▄███▅▄▄▃▂▂▂▂▂▂▂▂▂▂▂▂▃▂▂▂▂▂▂▂▁▁▁▂▁▁▁▂▁▂▂▂▂▂▂▂▂ ▃
25.4 ns         Histogram: frequency by time        26.6 ns <

Memory estimate: 0 bytes, allocs estimate: 0.julia> @benchmark bsplinebasisall($P, 2,$t)BenchmarkTools.Trial: 10000 samples with 999 evaluations.
Range (min … max):  8.333 ns … 25.092 ns  ┊ GC (min … max): 0.00% … 0.00%
Time  (median):     8.354 ns              ┊ GC (median):    0.00%
Time  (mean ± σ):   8.390 ns ±  0.473 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

▁▆  ▁█  ▅
▃▁▁██▁▁██▁▁█▁▁▁▂▁▁▁▂▁▁▁▄▁▁▃▄▁▁▃▄▁▁▂▁▁▁▂▁▁▁▃▁▁▁▄▁▁▂▃▁▁▂▂▁▁▂ ▃
8.33 ns        Histogram: frequency by time        8.48 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{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)