How to build a recursive spline function in C++

Question

At the moment I am working on a differential equation solving method called basis-spline collocation. What I am having trouble with is building a method to build an arbitrary order spline, with the relationship $$ B^{k+1}_{i}(x)= \frac{x-x_i}{x_{k+i}-x_i}B^k_i + \frac{x_{k+i+1}-x}{x_{k+i+1}-x_{i+1}}B^k_{i+1}(x) $$ with the initial condition $$B^1_i(x)=\begin{cases} 1 & \; \text{for } \; x_i \leq x < x_{i+1} \\ 0 & \; \text{otherwise} \end{cases}$$ and I am having trouble even starting with this problem, as it is recursive on could start from either the "top" or the "bottom", and I am running into a general writers block type of thing, where I can't get my mind around what I'm needing to do.

Answer 1

I can recommend consulting The NURBS book, which seems to be a classic text on this subject. The algorithm itself is given on page 72, it is available for online viewing.

Answer 2

I second the NURBS book and also would like to highlight that while this recursive formula is typically how one expresses the B-spline basis (as in a paper) it is not how you implement the basis functions. Since all higher order functions are based on lower order, you can compute all $p+1$ nonzero functions at once and reuse the lower order evaluations. That is how Piegl's algorithm does this.

Answer 3

I honestly don't know how efficient this is, but one way to do it is with c++ templates:

The order is k, t is the knot structure, and x is the value you want.

template <int k> 
real BSpline(real x, real *t)
{
    if (*t <= x && x < *(t+k))
    {
        real a = (x - *t) / (*(t+k-1) - *t);
        real b = (*(t+k) - x) / (*(t+k) - *(t+1));

        return a * BSpline<k-1>(x, t) + b * BSpline<k-1>(x, (t+1));
    }
    else
        return 0;
};

template <>
real BSpline<1>(real x, real *t)
{
    if (*t <= x && x < *(t+1))
        return 1.;
    else
        return 0.;
};