实现Clamped B-Spline曲线时遭遇除零错误的技术求助

需求与输入参数 #

我需要绘制一条Clamped B-Spline曲线（曲线起始于第一个控制点，终止于最后一个控制点），输入参数如下：

控制点 #

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

节点向量 #

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};

注：首尾重复的节点值用于实现曲线的clamp效果。

问题现象 #

我按照B-Spline算法实现了代码，但运行时出现大量除零错误。具体原因是当i=0且j=3时，double adiv = vknot[i+j] - vknot[i];的计算结果为零。

我尝试过增加阶数（添加更多节点），但递归调用BSN函数时，当j足够小时仍会出现同样问题；也尝试返回0或1替代异常，但生成的曲线结果完全不符合预期，无法实现clamp效果。

注1：若不使用clamped样条，代码运行正常。
注2：我了解相关的B-Spline实现问题，但未找到针对clamped样条的解决方案。

实现代码 #

#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)
    {
        // 归一化节点值
        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)
{
    // 曲线上生成的点数
    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);
}

int main()
{
    theSpline.generateInput();

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

    // cGUI theGUI;
    return 0;
}

内容的提问来源于stack exchange，提问作者ravenspoint