3D空间三点B样条曲线构建失败原因及替代方案咨询

在3D空间中有三个点，希望构建一条经过所有点的平滑曲线（如样条曲线）。当点数量为4个及以上时，可使用SciPy的B样条实现，代码如下：

import numpy as np
from scipy.interpolate import splev, splprep
import matplotlib.pyplot as plt


dim = 3
# doesn't work for 3 points
num_points = 4
X = np.random.rand(dim, num_points)

tck, u = splprep(X, s=0)
new_points = splev(np.linspace(0, 1, 100), tck)

# plot
ax = plt.figure().add_subplot(projection="3d")
ax.scatter(*X)
ax.plot(*new_points)
plt.show()

但该方法无法处理三个点，以下是原因分析和替代方案：

为什么splprep无法处理3个点？ #

SciPy的splprep默认使用三次B样条（阶数k=3），而B样条的阶数k对控制点数量n有硬性要求：n > k。三次B样条需要至少4个控制点才能生成有效的节点向量，3个点不满足这个条件，因此会直接报错。即使强制尝试，也无法构造出严格经过所有三点的三次B样条曲线。

替代方案 #

1. 降低B样条阶数 #

将B样条的阶数设为2（二次B样条），二次B样条仅要求控制点数量≥3，刚好适配三点场景，且能生成严格经过所有点的平滑曲线。修改代码如下：

import numpy as np
from scipy.interpolate import splev, splprep
import matplotlib.pyplot as plt

dim = 3
num_points = 3
X = np.random.rand(dim, num_points)

# 指定k=2（二次B样条），满足3个控制点的要求
tck, u = splprep(X, s=0, k=2)
new_points = splev(np.linspace(0, 1, 100), tck)

# 绘图
ax = plt.figure().add_subplot(projection="3d")
ax.scatter(*X)
ax.plot(*new_points)
plt.show()

2. 使用三次Hermite样条 #

如果希望保持三次样条的高阶平滑性，可以使用三次Hermite样条。该方法需要端点的导数信息，若没有实际导数数据，可通过三点的差分近似端点导数：

import numpy as np
from scipy.interpolate import CubicHermiteSpline
import matplotlib.pyplot as plt

dim = 3
num_points = 3
X = np.random.rand(dim, num_points)
u = np.linspace(0, 1, num_points)  # 均匀分布的参数化变量

# 用相邻点差分近似端点导数
dx_start = (X[:,1] - X[:,0]) / (u[1] - u[0])
dx_end = (X[:,2] - X[:,1]) / (u[2] - u[1])

# 对每个维度分别构建Hermite样条
splines = []
for i in range(dim):
    spl = CubicHermiteSpline(u, X[i], [dx_start[i], dx_end[i]])
    splines.append(spl)

# 生成曲线上的采样点
u_new = np.linspace(0, 1, 100)
new_points = np.array([spl(u_new) for spl in splines])

# 绘图
ax = plt.figure().add_subplot(projection="3d")
ax.scatter(*X)
ax.plot(*new_points)
plt.show()

3. 使用二次Bezier曲线 #

二次Bezier曲线天然适配3个控制点，曲线平滑且严格经过首尾两点，中间点作为控制点影响曲线的弯曲形状。实现代码如下：

import numpy as np
import matplotlib.pyplot as plt

dim = 3
num_points = 3
X = np.random.rand(dim, num_points)

# 生成参数t的采样点
t = np.linspace(0, 1, 100)
# 二次Bezier曲线公式计算
new_points = (1-t)**2 * X[:,0:1] + 2*(1-t)*t * X[:,1:2] + t**2 * X[:,2:3]

# 绘图
ax = plt.figure().add_subplot(projection="3d")
ax.scatter(*X)
ax.plot(*new_points)
plt.show()

内容的提问来源于stack exchange，提问作者Nico Schlömer