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

This might not an answer your question but it should lead to your answer.

• Common inputs for creating B-Spline are control point list, knot list and degree. Your degree got by calculate this seem isn't correct.
• To create clamped B-Spline we need adjust input by repeat last knots and add first control points to the end of list (should be same number with degree). So your control point list seem isn't for create clamped B-Spline.
• We can use De Boor's algorithm for create common B-Spline de-Boor.

I was success to implement B-Spline here mine.

 struct Spline {
      std::vector<Vec3> control_points;
      std::vector<double> knots;
      std::vector<double> weights;
      int16_t degree = 2;

  };
 struct SplineGen {
        using point_type = Vec3;
        using iterator = math::SegmentIter<SplineGen, point_type, point_type>;

        SplineGen(const Spline& _sp)
            :sp{ _sp }
        {
            gen_component(10);
        }

        //--start deboor algorithm
        double basis(double u, int i, int p) const
        {
            if (p == 0) {
                return (sp.knots[i] <= u && u < sp.knots[i + 1]) ? 1. : 0.;
            }
            double a = 0.0;
            double b = 0.0;
            if (sp.knots[i + p] != sp.knots[i]) {
                a = (u - sp.knots[i]) / (sp.knots[i + p] - sp.knots[i]) * basis(u, i, p - 1);
            }
            if (sp.knots[i + p + 1] != sp.knots[i + 1]) {
                b = (sp.knots[i + p + 1] - u) / (sp.knots[i + p + 1] - sp.knots[i + 1]) * basis(u, i + 1, p - 1);
            }
            return a + b;
        }

        point_type at(double u) const
        {
            point_type result{};
            for (int i = 0;i <= sp.control_points.size(); ++i) {
                auto w = sp.weights.empty() ? 1. : sp.weights[i];
                result += basis(u, i, sp.degree) * w * sp.control_points[i];
            }
            return result;
        }
        //--end deboor algorithm

        point_type operator()(int t) const
        {
            return at(sp.knots.front() + (static_cast<double>(t) / seg) * knot_len);
        }


        void gen_component(int knot_seg)
        {
            seg = sp.knots.size() * knot_seg;
            knot_len = sp.knots.back() - sp.knots.front();
        }


        iterator begin() const { return iterator{ 0, *this }; }
        iterator end() const
        {
            return iterator{ seg, *this }; 
        }
        int size() const
        {
            return seg;
        }
    private:
        const Spline& sp;
        int seg;
        double knot_len;
    };

My Vec3

struct Vec3
{
    Vec3() = default;
    double x;
    double y;
    double z;
};
Vec3& operator+=(Vec3& a, const Vec3& b)
{
   a.x += b.x;
   a.y += b.y;
   a.z += b.z;
   return a;
}
Vec3 operator+(const Vec3& a, const Vec3& b)
{
   return { a.x + b.x, a.y + b.y, a.z + b.z };
}
Vec3 operator*(const Vec3& a, double s)
{
   return { a.x * s, a.y * s, a.z * s };
}
Vec3 operator*(double s, const Vec3& v)
{
   return v * s;
}

It's look similar to your. My spline require vector3D and weights you can ignore them. My Spline generator is implement c++ iterator but it isn't necessary.