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 divn";
        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. There the spline is NOT clamped )