使用法向量计算三次贝塞尔曲线偏移结果异常求助

问题概述 #

我尝试通过以下步骤实现三次贝塞尔曲线的偏移：

• 求解x和y方向的导数得到切向量(dx, dy)
• 将切向量旋转90度得到法向量(-dy, dx)
• 将法向量归一化为单位向量
• 单位法向量乘以偏移量后与原曲线上的点相加

对应的实现代码：

def cubic_offset(P0, P1, P2, P3, t, dist): 
    initx, inity = cubic(P0, P1, P2, P3, t) 
    dx, dy = cubic_dt(P0, P1, P2, P3, t) # 切向量
    
    normx, normy = -dy, dx
    mag = sqrt(normx**2 + normy**2)
    normx, normy = normx/mag, normy/mag
    
    return initx + dist*normx, inity + dist*normy

但该实现无法正常工作，尤其在极值点附近问题明显（附图中红色点为极值点）。我用于求解极值和包围盒的导数功能是正常的，且参考了相关方法验证可行，但代码与示例几乎一致却结果不符。

可能的问题点 #

• 法向量方向翻转：旋转90度得到的法向量有两个方向，若曲线参数t的走向导致极值点附近切向量方向突变，会使法向量突然反转，造成偏移轨迹拧结。
• 极值点数值不稳定：极值点处切向量模长趋近于0，直接归一化会触发除以0风险，导致单位法向量异常，偏移点位置错误。
• 曲线与导数实现不一致：若cubic（曲线求值）和cubic_dt（导数求值）的参数化逻辑不符合三次贝塞尔标准公式，会导致点与切向量不匹配，偏移结果自然出错。

解决方案 #

1. 固定法向量方向 #

通过二阶导数判断曲率方向，确保法向量始终指向曲线的同一侧，避免极值点附近的方向翻转：

def cubic_offset(P0, P1, P2, P3, t, dist): 
    initx, inity = cubic(P0, P1, P2, P3, t) 
    dx, dy = cubic_dt(P0, P1, P2, P3, t)
    d2x, d2y = cubic_d2t(P0, P1, P2, P3, t) # 实现二阶导数计算
    
    normx, normy = -dy, dx
    mag = sqrt(normx**2 + normy**2) + 1e-8 # 加epsilon避免除以0
    normx, normy = normx/mag, normy/mag
    
    # 用二阶导数校验法向量方向
    dot = normx*d2x + normy*d2y
    if dot < 0:
        normx, normy = -normx, -normy
    
    return initx + dist*normx, inity + dist*normy

2. 处理极值点数值问题 #

在计算模长时加入极小值epsilon，避免除以0的情况；同时可在模长过小时，用相邻点的法向量插值替代当前点计算，减少突变。

3. 校验曲线与导数的标准实现 #

确保cubic和cubic_dt完全符合三次贝塞尔的标准公式：

# 标准三次贝塞尔曲线求值
def cubic(P0, P1, P2, P3, t):
    mt = 1 - t
    mt2 = mt * mt
    t2 = t * t
    x = mt2 * mt * P0[0] + 3 * mt2 * t * P1[0] + 3 * mt * t2 * P2[0] + t2 * t * P3[0]
    y = mt2 * mt * P0[1] + 3 * mt2 * t * P1[1] + 3 * mt * t2 * P2[1] + t2 * t * P3[1]
    return x, y

# 标准三次贝塞尔曲线一阶导数求值
def cubic_dt(P0, P1, P2, P3, t):
    mt = 1 - t
    mt2 = mt * mt
    t2 = t * t
    dx = 3 * mt2 * (P1[0] - P0[0]) + 6 * mt * t * (P2[0] - P1[0]) + 3 * t2 * (P3[0] - P2[0])
    dy = 3 * mt2 * (P1[1] - P0[1]) + 6 * mt * t * (P2[1] - P1[1]) + 3 * t2 * (P3[1] - P2[1])
    return dx, dy

# 标准三次贝塞尔曲线二阶导数求值
def cubic_d2t(P0, P1, P2, P3, t):
    mt = 1 - t
    d2x = 6 * mt * (P2[0] - 2*P1[0] + P0[0]) + 6 * t * (P3[0] - 2*P2[0] + P1[0])
    d2y = 6 * mt * (P2[1] - 2*P1[1] + P0[1]) + 6 * t * (P3[1] - 2*P2[1] + P1[1])
    return d2x, d2y

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