10
$\begingroup$

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.

$\endgroup$
7
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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.;
};
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.