创建满足曲率条件的Bezier曲线

Question

在一次作业中，我正在尝试创建贝塞尔曲线（有4个点），并希望它的最大曲率不超过K_x。曲率是用一个公式计算的，公式是k = |y''|/(1+y'^2)^1.5。

我有4个点，其中p0=(0,0)和p3=(0.3,0)，而p1和p2是可以自由选择的位置。

我尝试的方法是先在[(0,0), (0,0.15), (0.3,0.15), (0,-0.15)]这个区域内进行扫描。对于我选择的每一个p1点（由于对称性，我只选择右上角的点），我会扫描所有可能的p2点，然后检查每种组合的曲线是否符合曲率的要求，如果符合，就保存这些点。

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, :]
                        
                    #fig=plt.figure(1,figsize=(6,6))
                    #ax=plt.subplot(111);
                    #plt.plot(x_vec, y_vec)
                        
                    # 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。我希望能得到一些建议，看看应该添加什么条件或限制来解决这个问题。也许我应该限制选择p2的区域？

我希望能确保最大曲率小于K_max。我使用了数值导数，但当我用解析方法计算曲线时，结果没有变化。

编辑： 我在我添加的图片中可以看到这个问题——标记为'x'的控制点导致曲线在左侧的曲率很大。 查看图片