Drawing a B-spline

Question

From the Wikipedia article:

In the mathematical subfield of numerical analysis, de Boor's algorithm is a polynomial-time and numerically stable algorithm for evaluating spline curves in B-spline form. It is a generalization of de Casteljau's algorithm for Bézier curves. The algorithm was devised by German-American mathematician Carl R. de Boor. Simplified, potentially faster variants of the de Boor algorithm have been created but they suffer from comparatively lower stability.

That article offers what it calls a "naive" implementation

def deBoor(k: int, x: int, t, c, p: int):
    """Evaluates S(x).

    Arguments
    ---------
    k: Index of knot interval that contains x.
    x: Position.
    t: Array of knot positions, needs to be padded as described above.
    c: Array of control points.
    p: Degree of B-spline.
    """
    d = [c[j + k - p] for j in range(0, p + 1)] 

    for r in range(1, p + 1):
        for j in range(p, r - 1, -1):
            alpha = (x - t[j + k - p]) / (t[j + 1 + k - r] - t[j + k - p]) 
            d[j] = (1.0 - alpha) * d[j - 1] + alpha * d[j]

    return d[p]

I have a couple of questions:

1. It looks like the c array goes out of bounds in the first line of code

   d = [c[j + k - p] for j in range(0, p + 1)] 
   

   when generating the first point the first loop through range(0, p + 1) for a 3rd degree spline

   j = 0;
   k = 0;
   p = 3;
   

   so d[0] will be set to c[ -3 ]

   Question: How do I prevent attempts to access negative indices of c?

   Note 1: I tried checking for negative indices and accessing c[0] instead. The results look OK, but can this be correct?

2. In this line

   alpha = (x - t[j + k - p]) / (t[j + 1 + k - r] - t[j + k - p]) 
   

   should there not be a check to prevent divide by zero?

Answer 1

The knot values must be in the range 0 to 1, with 0 being the start of the spline curve and 1 the end of the curve.

If the input knot values are not in the expected range, they must be normalized.

void CSpline::BSplineInit()
{
    // check that there are some knots
    if (!vknots.size())
        return;

    // normalize the knots to the range 0 to 1 inclusive
    double min = vknots[0];
    double max = vknots.back();
    if (max <= min)
        throw std::runtime_error("BSpline has invalid knot values");
    double scale = max - min;
    double off = min;
    for (auto &t : vknots)
        t = (t - off) / scale;
}

When this is done, it is sufficient to check for valid values, and throw an exception on failure.

This C++ code works well for me ( no exceptions are thrown and the spline curves show as expected ):

double CSpline::deBoor(
    int k,
    double x,
    const std::vector<double> &c)
{

    std::vector<double> d;
    for (int j1 = 0; j1 < myDegree + 1; j1++)
    {
        int index = j1 + k - myDegree;
        if (index < 0)
            throw std::runtime_error("CSpline::deBoor bad knot index");

        d.push_back(c[index]);
    }
    for (int r = 1; r < myDegree + 2; r++)
    {
        for (int j = myDegree; j > r - 1; j--)
        {
            double div = vknots[j + 1 + k - r] - vknots[j + k - myDegree];
            if (div < 0.0000001)
                throw std::runtime_error("CSpline::deBoor divide by zero");
            double alpha = (x - vknots[j + k - myDegree]) / div;
            d[j] = (1 - alpha) * d[j - 1] + alpha * d[j];
        }
    }

    if (std::isnan(d[myDegree]))
        throw std::runtime_error("CSpline::deBoor nan");

    return d[myDegree];
}