Consider the following inverse triangular formula
$$\left( \begin{array}{ccccc} & & & & N_{i-p,p}\left(u_0\right) \\ & & N_{i-2,2}\left(u_0\right) & & \\ & N_{i-1,1}\left(u_0\right) & & & \\ N_{i,0}\left(u_0\right) & & N_{i-1,2}\left(u_0\right) & \cdots & \vdots \\ & N_{i,1}\left(u_0\right) & & & \\ & & N_{i,2}\left(u_0\right) & & \\ & & & & N_{i,p}\left(u_0\right) \end{array} \right)$$
where, $N_{i,0}=1$, and ![enter image description here][1]
In addition, $u_0 \in [u_i,u_{i+1})$ knots = $\{u_0,u_1, \cdots, u_m\}$
Here is a procedural implementaion
Search the index $i$ by the auxiliary function
searchSpan![enter image description here][2]
In the code, I use the following local array to store the values of $N_{m,n}$
$$ \left( \begin{array}{ccccc} & & & & N_{i-p,p}\left(u_0\right) \\ & & & & \\ & & N_{i-2,2}\left(u_0\right) & & \vdots \\ & N_{i-1,1}\left(u_0\right) & N_{i-1,2}\left(u_0\right) & \cdots & \\ N_{i,0}\left(u_0\right) & N_{i,1}\left(u_0\right) & N_{i,2}\left(u_0\right) & \cdots & N_{i,p}\left(u_0\right) \end{array} \right)_{(p+1)\times (p+1)} $$
where $N_{m,n}$ was stored in the position $(p+1-i+m,n+1)$
Search the index of span $[u_i,u_{i+1})$
searchSpan[{deg_, knots_}, u0_] :=
Module[{biSearch},
biSearch =
Function[{low, high},
With[{mid = Floor[(low + high)/2]},
If[u0 < knots[[mid]], {low, mid}, {mid, high}]]
];(*Do bisection search*)
First@
NestWhile[
biSearch[Sequence @@ #, u0] &,
{deg + 1, Length@knots - deg}, Subtract @@ # != -1 &] - 1
]
NonzeroBasis[{deg_, knots_}, u0_] :=
Module[{coeff, basis, i},
coeff =
(u0 - knots[[#1 + 1]])/(knots[[#1 + #2 + 1]] - knots[[#1 + 1]]) &;
basis = ConstantArray[1, {deg + 1, deg + 1}];
i = searchSpan[{deg, knots}, u0];
Do[
basis[[deg + 1 - k, k + 1]] =
(1 - coeff[i - k + 1, k]) basis[[deg + 2 - k, k]];
With[{m = deg + 1 - i},
Do[
basis[[m + j, k + 1]] =
{coeff[j, k], 1 - coeff[j + 1, k]}.{basis[[m + j, k]], basis[[m + j + 1, k]]},
{j, i - k + 1, i - 1}]
];
basis[[deg + 1, k + 1]] =
coeff[i, k] basis[[deg + 1, k]],
{k, deg}];
basis
]
###Test
knots = {0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 5, 6, 7, 8, 9, 10, 10, 10, 10, 10};
deg = 4;
NonzeroBasis[{deg, knots}, 5/2] // MatrixForm
$\left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & \frac{1}{288} \\ 1 & 1 & 1 & \frac{1}{48} & \frac{227}{1152} \\ 1 & 1 & \frac{1}{8} & \frac{23}{48} & \frac{205}{384} \\ 1 & \frac{1}{2} & \frac{3}{4} & \frac{15}{32} & \frac{25}{96} \\ 1 & \frac{1}{2} & \frac{1}{8} & \frac{1}{32} & \frac{1}{192} \end{array} \right)$
###Question
- How to implement this
triangular formulain a non-procedural(like functional or rule-based) method?
Related question
