Generating a clamped B-Spline

Question

I am trying to display a clamped B-Spline curve. ( That is, one where the curve starts on the first control point and ends on last control point )

Here is the input:

Control points

    xcp = {
        5, 5, 16, 31, 22, 33, 44, 42, 51, 50, 59};
    ycp = {
        27, 12, 29, 18, 9, 9, 20, 29, 28, 10, 10};

Knots

    vknot = {
        0.0,        0.0,        0.0,        0.0,        
        1.0,        2.0,        3.0,        4.0,        5.0,        6.0,        7.0,
        8.0,        8.0,        8.0,        8.0};

Note the repeated knot values at beginning and end - that is what should clamp the curve. See https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html

I have tried implementing the described algorithm here https://mathworld.wolfram.com/B-Spline.html

Here is my code

#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
#include <vector>
#include <algorithm>


class CSpline
{
public:
    int m_ControlPointCount;

    std::vector<double> xcp, ycp, vknot;

    double current;

    void generateInput();

    bool getNextPoint(int &xp, int &yp);

private:
    double BSN(int i, int j, double t);
};

CSpline theSpline;

void CSpline::generateInput()
{
    vknot = {
        0.0,
        0.0,
        0.0,
        0.0,
        1.0,
        2.0,
        3.0,
        4.0,
        5.0,
        6.0,
        7.0,
        8.0,
        8.0,
        8.0,
        8.0};
    if (fabs(1 - vknot.back()) > 0.01)
    {
        // normalize knot values
        for (double &v : vknot)
            v /= vknot.back();
    }
    xcp = {
        5, 5, 16, 31, 22, 33, 44, 42, 51, 50, 59};
    ycp = {
        27, 12, 29, 18, 9, 9, 20, 29, 28, 10, 10};

    current = 0;
}

bool CSpline::getNextPoint(int &xp, int &yp)
{
    // number of points drawn along curve
    const int Ndiv = 100;
    if (current == Ndiv)
        return false;
    double t = current / Ndiv;
    int degree = vknot.size() - xcp.size() - 1;

    double x,y;
    x = y = 0;
    for (int i = 0; i < xcp.size(); i++)
    {
        double N = BSN(i, degree, t);
        x += xcp[i] * N;
        y += ycp[i] * N;
    }

    std::cout << t <<" "<< x <<" "<< y << "\n";

    current++;

    return true;
}

double CSpline::BSN(int i, int j, double t)
{
    if (j == 0)
    {
        if (vknot[i] <= t && t < vknot[i + 1] &&
            vknot[i] < vknot[i + 1])
        {
            return 1;
        }
        return 0;
    }
    double adiv = vknot[i + j] - vknot[i];
    double bdiv = vknot[i + j + 1] - vknot[i + 1];
    if (fabs(adiv) < 0.00001 || fabs(bdiv) < 0.00001) {
        std::cout << "warning zero div\n";
        return 0;
    }
    double a = (t - vknot[i]) / adiv;
    double b = (vknot[i + j + 1] - t) / bdiv;
    return a * BSN(i, j - 1, t) + b * BSN(i + 1, j - 1, t);
}

main()
{
    theSpline.generateInput();

    int x, y;
    while( theSpline.getNextPoint(x,y))
        ;

    cGUI theGUI;
    return 0;
}

The code causes numerous divide by zero problems

This happens because

double adiv = vknots[i+j] - vknots[i];

is zero when i = 0 and j = 3

because

I have tried increasing the degree ( more knots ) but the same thing happens deeper in the recursive calls to BSN as J becomes small enough.

I have tried returning 0 or 1 instead of throwing an exception, but the resulting points on the curve are nonsense rather than clamping.

( Side note 1: the code seems to work OK if the spline is not clamped )

( Side note 2: I am aware of this question ( B-Splines in C++ ). There the spline is NOT clamped )

Answer 1

The best algorithm for calculation a good BSpline curve does seem to be deBoor. I have looked at several descriptions of this algorithm on the web. IMHO they are all inadequate in that they do not describe their parameters and usage. Some have errors in the sample code.

For me, the best is on Wikipedia ( it still does not properly explain the input variables but gives enough clues that it can be fixed with some work )

python code from https://en.wikipedia.org/wiki/De_Boor%27s_algorithm#Example_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]

Note 1: input x "the position" is declared to be an integer. This makes no sense.
IMHO this should be a double in the range [0,1] where 0 is the first point and 1 is the last. [] indicates end points included in range

Note 2: input c: array of control points is an array of values for one dimension of the control points for a 2D curve this must be called twice, once for the x values and once for the y values

Here is my port to C++

    /// @brief deBoor algorithm for B-Spline generation
    /// @param k index of greatest knot less than x
    /// @param x position along the curve in the range [0,1] where 0 is the first point and 1 is the last.
    /// @param t knots in the range [0,1] where 0 is the first point and 1 is the last
    /// @param c values for one dimension of the control points ( e.g x or y for 2D curve)
    /// @param p degree of spline 
    /// @return value of spline curve in dimenson provided by c
    ///
    /// call this once for each dimension of the curve

    double deBoor(
        int k,
        double x,
        const std::vector<double> &t,
        const std::vector<double> &c,
        int p)
    {
        std::vector<double> d;
        for (int j1 = 0; j1 < p + 1; j1++)
            d.push_back(c[j1 + k - p]);

        for (int r = 1; r < p + 2; r++)
        {
            for (int j = p; j > r - 1; j--)
            {
                double alpha = (x - t[j + k - p]) / (t[j + 1 + k - r] - t[j + k - p]);
                d[j] = (1 - alpha) * d[j - 1] + alpha * d[j];
            }
        }
        return d[p];
    }

This gives a decent result with the input data given in my question