满足曲率约束的三次Bezier曲线控制点筛选问题

任务描述 #

需构造最大曲率不超过$K_{max}$的三次Bezier曲线（含4个控制点），曲率计算公式为：
$$k = \frac{|y''|}{(1+y'²⁾{1.5}}$$
已知固定控制点$p_0=(0,0)$、$p_3=(0.3,0)$，$p_1$、$p_2$为可调整控制点。当前采用网格扫描方案：利用对称性仅在右上象限选取$p_1$，遍历所有可能的$p_2$，检查曲线是否满足曲率约束，保存符合条件的控制点组合。

现有实现代码 #

import numpy as np
import scipy as sc
from matplotlib import pyplot as plt

def bezzier_curve(p0, p1, p2, p3, t):
    b0 = b_in_calc(3, 0, t) * np.transpose(p0)
    b1 = b_in_calc(3, 1, t) * np.transpose(p1)
    b2 = b_in_calc(3, 2, t) * np.transpose(p2)
    b3 = b_in_calc(3, 3, t) * np.transpose(p3)
    b = b0 + b1 + b2 + b3
    return b


def b_in_calc(n, i, t):
    b = sc.special.binom(n, i)
    t1 = np.power(1-t, n-i)
    t2 = np.power(t, i)
    return b*t1*t2


def numerical_deriv(y, x):
    # finding numerical derivative of the finction
    dx = np.gradient(x)
    dy = np.gradient(y)
    d = dy/dx
    return d


k_max_allowed = 5
d = 0.3
p0 = np.zeros([1, 2])
p3 = np.array([[0.3, 0]])
p1 = np.zeros([1, 2])
p2 = np.zeros([1, 2])
t = np.linspace(0, 1, 10000)
P1_space = np.linspace(0, d/2, num=40)
P2_xspace = np.linspace(0, d, endpoint=False, num=40)
P2_yspace = np.linspace(-d/2, d/2, 40)

output_name = 'desired_coordinates_K_max_5.txt'

condition = 0 # if we started writing P1 or not
condition1 = 0 # if we already wrote at least 1 P2 point

with open(output_name, 'w',encoding='utf-8') as f:
    for i in range(1, len(P1_space)):
        p1[0, 0] = P1_space[i]
        for j in range(1, len(P1_space)):
            p1[0, 1] = P1_space[j]            
            
            for n in range(1, len(P2_xspace)):
                p2[0,0] = P2_xspace[n]
                for m in range(len(P2_yspace)):
                    p2[0,1] = P2_yspace[m]
                        
                    # create bezier curve
                    bezz = bezzier_curve(p0, p1, p2, p3, t)
                    x_vec = bezz[0, :]
                    y_vec = bezz[1, :]
                        
                    # interpolate points to spread them
                    x_interp = np.linspace(0, d, len(t), endpoint=True)
                    y_interp = np.interp(x_interp, x_vec, y_vec)  
                        
                    # derive and find curvature vector and max k
                    y_deriv = numerical_deriv(y_interp, x_interp)
                    y_deriv2 = numerical_deriv(y_deriv, x_interp)
                    k_vec = y_deriv2/(1 + y_deriv**2)**1.5
                    k_max = np.amax(k_vec)
                        
                    # check if meets condition, if so write it to .txt file
                    if k_max <= k_max_allowed:
                        if condition == 0:
                            f.write('P1:' + str(P1_space[int(i)]) + ',' + str(P1_space[int(j)]) + '\n')
                            condition = 1
                        if condition1 ==0:
                            f.write('P2:')
                            condition1 = 1
                        f.write(str(P2_xspace[n]) + ',' + str(P2_yspace[m]) + '|')
                        
            condition = 0
            condition1 = 0
        f.write('\n')
f.close()

plt.text(5,5 , 'Complete', fontsize = 22)
plt.xlim(0, 15)
plt.ylim(0, 10)
plt.show()

遇到的问题 #

读取保存的控制点组合并验证时，部分组合生成的曲线曲率超过设定的$K_{max}$。即使改用解析求导，问题仍存在，部分控制点会导致曲线左侧出现大曲率区域。

改进建议 #

1. 修正曲率计算的绝对值处理 #

曲率是标量，需取绝对值。当前代码中k_vec未计算绝对值，导致仅考虑正的二阶导数对应的曲率，漏掉了负二阶导数的大曲率情况（比如左侧区域）。修改为：

k_vec = np.abs(y_deriv2)/(1 + y_deriv**2)**1.5
k_max = np.amax(k_vec)

2. 替换数值求导为解析导数（更稳定） #

三次Bezier曲线的导数有解析表达式，避免插值和数值求导的误差：

• 一阶导数（对参数t）：
  $$B'(t) = 3(1-t)^2(p_1-p_0) + 6(1-t)t(p_2-p_1) + 3t^2(p_3-p_2)$$
• 二阶导数（对参数t）：
  $$B''(t) = 6(1-t)(p_2-2p_1+p_0) + 6t(p_3-2p_2+p_1)$$
• 转换为对x的导数：
  $$\frac{dy}{dx} = \frac{dy/dt}{dx/dt}, \quad \frac{d^2y}{dx2} = \frac{d/dt(dy/dx)}{dx/dt}$$
  直接基于t计算这些值，避免x插值带来的误差，尤其是曲线x非单调时。

3. 限制p2的选取范围 #

• 排除x非单调的组合：当$dx/dt=0$时，曲线出现垂直切线或回折，曲率会趋于无穷大。提前检查$dx/dt$在$t∈[0,1]$内是否有零点，若有则直接跳过该控制点组合。
• 利用对称性约束：若目标是对称曲线，可强制$p_2=(0.3-p_1.x, -p_1.y)$，缩小搜索范围；非对称场景也可限制$p_2$的y值范围，使其与$p_1$的y绝对值相当，避免曲线起伏过大。

4. 精准捕捉曲率极值 #

• 仅靠采样点的最大值可能漏掉真实极值，可使用数值优化方法（如scipy.optimize.minimize_scalar）寻找$|k(t)|$在$t∈[0,1]$内的最大值，确保准确检测到最大曲率。
• 增加关键区域采样密度：在$dx/dt$变化大的区域（即曲线斜率变化快的地方）增加t的采样点，提高极值检测概率。

5. 添加后验证步骤 #

在保存控制点组合前，用更高精度的方法（解析导数+极值优化）重新验证一次，确认最大曲率严格小于$K_{max}$后再写入文件，避免误判。

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