My tool of choice to investigate, learn about and toy with splines has always been GeoGebra, because it is easy to drag control points and sliders around to see what happens. Splines like B'ezier curves are easily created there as long as you remembered the basis functions by heart. Luckily these basis functions are easy to come by:

\[B_i^k(t) = {k\choose i} (1-t)^{k-i} t^i,\]

here \(k\) is the degree and \(i\) the index so that the i-th basis function is multiplied with the $i$-th control point.

When considering B-splines or even non-uniform rational B-splines (NURBS) this becomes a little harder to do. Using these in GeoGebra requires evaluating the B-spline through De Boor’s algorithm:

\[N_{i, 0}(t) = \begin{cases} 1 & t_i \leq t < t_{i+1} \\ 0 & \text{otherwise} \end{cases}\] \[N_{i, k}(t) = \frac{t - t_i}{t_{i+k} - t_i} N_{i, k-1}(t) + \frac{t_{i+k+1} - t}{t_{i+k+1} - t_{i+1}} N_{i+1, k - 1}(t)\]

As anyone who has used any of the newer (browser-based) Geogebra versions knows, inputting formulas is an excruciatingly painful ordeal. Moreover, with higher degree splines come increasing levels of recursion, and therefore unbearable slowness in GeoGebra. For this I usually resorted to not use the recursion based formulation of the basis functions but their closed form expressions. Again, for higher degrees this becomes a very error prone and tedious process.

For this purpose I made this simple Python script that can be used expand and write out the B-spline basis functions given a degree \(k\), a knot vector \(\mathbf{t} = [t_0, \ldots, t_{2(k + 1)}]\), the index \(i\) (with respect to the knot vector) and a knot span \(j\) (meaning \([t_j, t_{j+1}]\)):

from sympy import symbols, simplify
import numpy as np

T = symbols('t')

def N(i, k, knots, span):
   if(k == 0):
      if(i == span):   
         return 1
      else:
         return 0
   
   a = 0
   D_1 = knots[i + k] - knots[i]
   if(D_1 != 0):
      a = (T - knots[i]) / D_1

   b = 0
   D_2 = knots[i + k + 1] - knots[i+1]
   if(D_2 != 0):
      b = (knots[i + k + 1] - T) / D_2

   return a * N(i, k - 1, knots, span) + b * N(i + 1, k - 1, knots, span)

Then outputting the uniform cubic B-spline basis functions on the interval \([0, 1]\) is as easy as:

knots = np.array([-3, -2, -1, 0, 1, 2, 3, 4])
span = 3
print(simplify(N(0, 3, knots, span)))
print(simplify(N(1, 3, knots, span)))
print(simplify(N(2, 3, knots, span)))
print(simplify(N(3, 3, knots, span)))

Which returns:

-(t - 1)**3/6
t**3/2 - t**2 + 2/3
-t**3/2 + t**2/2 + t/2 + 1/6
t**3/6