从点集推导三次贝塞尔曲线控制点与控制柄的问题求助

我尝试从一系列点中拟合三次贝塞尔曲线并获取对应的控制点与控制柄，目前用scipy.interpolate.CubicSpline实现了预期的曲线拟合效果，但手动设置的橙色控制点（控制柄）明显不正确。当前代码如下（致谢Python Discord的Zero Zero）：

import numpy as np
import scipy as sp

def fit_curve(points):
    # Fit a cubic bezier curve to the points
    curve = sp.interpolate.CubicSpline(points[:, 0], points[:, 1], bc_type=((1, 0.0), (1, 0.0)))

    # Get 4 control points for the curve
    p = np.zeros((4, 2))
    p[0, :] = points[0, :]
    p[3, :] = points[-1, :]
    p[1, :] = points[0, :] + 0.3 * (points[-1, :] - points[0, :])
    p[2, :] = points[-1, :] - 0.3 * (points[-1, :] - points[0, :])

    return p, curve

ypoints = [0.0, 0.03771681353260319, 0.20421680080883106, 0.49896111463402026, 0.7183501026981503, 0.8481517096346528, 0.9256128196832564, 0.9705404287079152, 0.9933297674379904, 1.0]
xpoints = [x for x in range(len(ypoints))]

points = np.array([xpoints, ypoints]).T

from scipy.interpolate import splprep, splev
tck, u = splprep([xpoints, ypoints], s=0)
#print(tck, u)
xnew, ynew = splev(np.linspace(0, 1, 100), tck)

# Plot the original points and the Bézier curve
import matplotlib.pyplot as plt
#plt.plot(xpoints, ypoints, 'x', xnew, ynew, xpoints, ypoints, 'b')
plt.axis([0, 10, -0.05, 1.05])
plt.legend(['Points', 'Bézier curve', 'True curve'])
plt.title('Bézier curve fitting')

# Get the curve
p, curve = fit_curve(points)

# Plot the points and the curve
plt.plot(points[:, 0], points[:, 1], 'o')
plt.plot(p[:, 0], p[:, 1], 'o')
plt.plot(np.linspace(0, 9, 100), curve(np.linspace(0, 9, 100)))

plt.show()

现在需要解决的是：如何获取与拟合的三次样条曲线匹配的正确三次贝塞尔控制点（包括两个控制柄）？

解决方案：从样条端点导数推导贝塞尔控制点 #

你的问题出在手动设置p[1]和p[2]的方式——这只是简单的线性比例，和实际拟合的样条曲线没有关联。要得到匹配三次样条的贝塞尔控制点，需要利用三次贝塞尔曲线的端点导数与控制点的关系：

对于三次贝塞尔曲线 ( B(t) = (1-t)^3P_0 + 3(1-t)^2tP_1 + 3(1-t)t^2P_2 + t^3P_3 )：

• 起点 ( t=0 ) 的一阶导数：( B'(0) = 3(P_1 - P_0) )
• 终点 ( t=1 ) 的一阶导数：( B'(1) = 3(P_3 - P_2) )

而CubicSpline可以直接计算端点的导数，我们可以用这个导数来反推控制柄 ( P_1 ) 和 ( P_2 )：

修改后的代码 #

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def fit_curve(points):
    # 拟合三次样条曲线
    curve = sp.interpolate.CubicSpline(points[:, 0], points[:, 1], bc_type=((1, 0.0), (1, 0.0)))
    
    # 获取端点导数
    start_deriv = curve(points[0, 0], 1)  # 起点一阶导数
    end_deriv = curve(points[-1, 0], 1)    # 终点一阶导数
    
    # 计算贝塞尔控制点
    p = np.zeros((4, 2))
    p[0, :] = points[0, :]
    p[3, :] = points[-1, :]
    # 从导数反推控制柄：P1 = P0 + B'(0)/3，P2 = P3 - B'(1)/3
    p[1, :] = p[0, :] + start_deriv / 3
    p[2, :] = p[3, :] - end_deriv / 3
    
    return p, curve

# 原始数据点
ypoints = [0.0, 0.03771681353260319, 0.20421680080883106, 0.49896111463402026, 0.7183501026981503, 0.8481517096346528, 0.9256128196832564, 0.9705404287079152, 0.9933297674379904, 1.0]
xpoints = [x for x in range(len(ypoints))]
points = np.array([xpoints, ypoints]).T

# 获取贝塞尔控制点和样条曲线
p, curve = fit_curve(points)

# 绘制结果
plt.axis([0, 10, -0.05, 1.05])
plt.title('Bézier curve fitting with correct control points')

# 原始数据点
plt.plot(points[:, 0], points[:, 1], 'o', label='Original points')
# 样条拟合曲线
plt.plot(np.linspace(0, 9, 100), curve(np.linspace(0, 9, 100)), label='Cubic Spline curve')
# 贝塞尔控制点及控制柄连线
plt.plot(p[:, 0], p[:, 1], 'o', label='Bezier control points')
plt.plot([p[0,0], p[1,0]], [p[0,1], p[1,1]], '--', color='orange')
plt.plot([p[2,0], p[3,0]], [p[2,1], p[3,1]], '--', color='orange')

plt.legend()
plt.show()

关键说明 #

1. 导数与控制点的关系：利用三次贝塞尔曲线的导数公式，直接从样条的端点导数计算控制柄位置，保证贝塞尔曲线的端点切线与样条曲线完全一致。
2. 拟合一致性：如果你需要的是单个三次贝塞尔曲线拟合所有点，上述方法能保证端点匹配，但中间点可能有偏差。如果要让贝塞尔曲线更贴合所有点，需要用最小二乘法拟合四个控制点，可以参考scipy.optimize.curve_fit实现：

最小二乘法拟合完整贝塞尔曲线的示例

from scipy.optimize import curve_fit

def bezier_curve(t, p0x, p0y, p1x, p1y, p2x, p2y, p3x, p3y):
    """计算三次贝塞尔曲线上的点"""
    return (
        (1-t)**3 * p0x + 3*(1-t)**2*t * p1x + 3*(1-t)*t**2 * p2x + t**3 * p3x,
        (1-t)**3 * p0y + 3*(1-t)**2*t * p1y + 3*(1-t)*t**2 * p2y + t**3 * p3y
    )

# 固定起点和终点，只拟合中间两个控制点
def constrained_bezier(t, p1x, p1y, p2x, p2y):
    p0x, p0y = points[0]
    p3x, p3y = points[-1]
    return bezier_curve(t, p0x, p0y, p1x, p1y, p2x, p2y, p3x, p3y)

# 生成参数t（对应x坐标归一化到0-1）
t_data = (xpoints - xpoints[0]) / (xpoints[-1] - xpoints[0])
# 拟合参数
initial_guess = [p[1,0], p[1,1], p[2,0], p[2,1]]  # 用之前的结果做初始值
params, _ = curve_fit(constrained_bezier, t_data, (xpoints, ypoints), p0=initial_guess)

# 提取拟合后的控制点
p_fit = np.zeros((4,2))
p_fit[0] = points[0]
p_fit[1] = params[:2]
p_fit[2] = params[2:]
p_fit[3] = points[-1]

# 绘制拟合的贝塞尔曲线
t_new = np.linspace(0,1,100)
x_bez, y_bez = bezier_curve(t_new, *p_fit.flatten())
plt.plot(x_bez, y_bez, label='Fitted Bezier curve', linestyle='--')

这个方法会让贝塞尔曲线尽可能贴合所有数据点，适合需要单个贝塞尔曲线的场景。

内容的提问来源于stack exchange，提问作者Dr. Pontchartrain