有理贝塞尔曲线切线计算出现NaN/无穷值问题的排查咨询

问题背景 #

正在练习贝塞尔曲线相关数学知识（基础较为薄弱），现有一段可正常运行的三次有理贝塞尔曲线计算代码（基于Pomax的《Bezier曲线入门教程》及Freya Holmér的视频教程实现），但编写切线计算函数时遇到问题：计算出的分母值常为0或极小值，导致返回NaN或Infinity。

原曲线计算代码：

float CalculateCubic(float t, float[] points, float[] weights)
{
    float t2 = t * t;
    float t3 = t2 * t;

    float[] f = {
        weights[0] * (-t3 + 3 * t2 - 3*t + 1),
        weights[1] * (3 * t3 - 6 * t2 + 3 * t),
        weights[2] * (-3 * t3 + 3 * t2),
        weights[3] * (t3)
    };
    
    float basis = f[0] + f[1] + f[2] + f[3];
    
    return (f[0] * points[0] + f[1] * points[1] + f[2] * points[2] + f[3] * points[3]) /basis;
}

错误的切线计算代码：

float CalculateCubicTangent(float t, float[] points, float[] weights)
{
    float t2 = t * t;

    float[] f = {
        weights[0] * (-3 * t2 + 6 * t - 3),
        weights[1] * (9 * t2 - 12 * t + 3),
        weights[2] * (-9 * t2 + 6 * t),
        weights[3] * (3 * t2)
    };

    float basis = f[0] + f[1] + f[2] + f[3];

    return (f[0] * points[0] + f[1] * points[1] + f[2] * points[2] + f[3] * points[3]) / basis;
}

错误分析 #

直接照搬原曲线的计算逻辑，对带权重的基函数求导后相除是错误的。有理贝塞尔曲线本质是带权点的贝塞尔曲线（分子N(t)）除以权重的贝塞尔曲线（分母D(t)），即P(t) = N(t)/D(t)。求导必须使用商的导数法则：
P’(t) = (N’(t)*D(t) - N(t)*D’(t)) / D(t)²

错误代码的核心问题：

1. 将sum(w_i*dB_i(t))（即D’(t)）作为分母，而D’(t)在很多t值下会为0（比如对称权重时t=0.5），直接除以它必然导致NaN/Infinity
2. 没有正确应用商的导数法则，忽略了N(t)和D(t)的交叉项

正确实现思路与代码 #

按照商的导数法则，步骤如下：

1. 计算原伯恩斯坦基函数B₀-B₃及其导数dB₀-dB₃
2. 计算分子N(t)、分母D(t)（与原曲线计算逻辑一致）
3. 计算分子的导数N’(t)、分母的导数D’(t)
4. 代入商的导数公式计算切线，同时加入极小值判断避免除以0

正确代码：

float CalculateCubicTangent(float t, float[] points, float[] weights)
{
    float t2 = t * t;
    float t3 = t2 * t;

    // 原三次伯恩斯坦基函数
    float B0 = -t3 + 3 * t2 - 3 * t + 1;
    float B1 = 3 * t3 - 6 * t2 + 3 * t;
    float B2 = -3 * t3 + 3 * t2;
    float B3 = t3;

    // 伯恩斯坦基函数的导数
    float dB0 = -3 * t2 + 6 * t - 3;
    float dB1 = 9 * t2 - 12 * t + 3;
    float dB2 = -9 * t2 + 6 * t;
    float dB3 = 3 * t2;

    // 计算N(t)和D(t)（对应原曲线的分子和分母）
    float N = weights[0] * B0 * points[0] + weights[1] * B1 * points[1] + weights[2] * B2 * points[2] + weights[3] * B3 * points[3];
    float D = weights[0] * B0 + weights[1] * B1 + weights[2] * B2 + weights[3] * B3;

    // 计算N'(t)和D'(t)
    float dN = weights[0] * dB0 * points[0] + weights[1] * dB1 * points[1] + weights[2] * dB2 * points[2] + weights[3] * dB3 * points[3];
    float dD = weights[0] * dB0 + weights[1] * dB1 + weights[2] * dB2 + weights[3] * dB3;

    // 应用商的导数法则，避免除以0
    float denominator = D * D;
    if (MathF.Abs(denominator) < 1e-8)
        return 0; // 可根据需求调整默认返回值

    return (dN * D - N * dD) / denominator;
}

学习资源推荐 #

• Pomax《Bezier曲线入门教程》：详细覆盖有理贝塞尔曲线的数学推导与导数计算
• Freya Holmér的贝塞尔曲线视频教程：通过可视化演示帮助理解几何意义
• 《计算机图形学原理及实践》：经典教材，系统讲解有理贝塞尔曲线的理论与应用

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