我希望将描述任意曲线的一系列有序（非常密集）2D 点转换为 NURBS 表示，可以将其写入 IGES 文件。
我正在使用 scipy.interpolate 的 splprep 来获得给定一系列点的 B 样条表示，然后我假设 NURBS 定义基本上就是这样，加上所有权重都等于 1。但是我认为我从根本上误解了splprep 的输出，特别是“B 样条系数”与在某些 CAD 包中手动重新创建样条所需的控制点之间的关系（我使用的是 Siemens NX11）。
我尝试了一个从稀疏点集逼近函数 y = x^3 的简单示例：
import scipy.interpolate as si
import numpy as np
import matplotlib.pyplot as plt
# Sparse points defining cubic
x = np.linspace(-1,1,7)
y = x**3
# Get B-spline representation
tck, u = si.splprep([x,y],s=0.0)
# Get (x,y) coordinates of control points
c_x = tck[1][0]
c_y = tck[1][1]
# Plotting
u_fine = np.linspace(0,1,1000)
x_fine, y_fine = si.splev(u_fine, tck)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y, 'o', x_fine, y_fine)
ax.axis('equal')
plt.show()
它给出了以下参数：
>>> t
array([ 0. , 0. , 0. , 0. , 0.39084883,
0.5 , 0.60915117, 1. , 1. , 1. , 1. ])
>>> c_x
array([ -1.00000000e+00, -9.17992269e-01, -6.42403598e-01,
-2.57934892e-16, 6.42403598e-01, 9.17992269e-01,
1.00000000e+00])
>>> c_y
array([ -1.00000000e+00, -7.12577481e-01, -6.82922469e-03,
-1.00363771e-18, 6.82922469e-03, 7.12577481e-01,
1.00000000e+00])
>>> k
3
>>> u
array([ 0. , 0.25341516, 0.39084883, 0.5 , 0.60915117,
0.74658484, 1. ])
>>>
我假设两组系数 (c_x, c_y) 描述了构造样条曲线所需的极点 (x,y) 坐标。在 NX 中手动尝试会得到类似的样条曲线，但并不完全相同，区间中的其他点的评估方式与 Python 中的不同。当我将此手动样条曲线导出为 IGES 格式时，NX 将节点更改为以下（同时显然保持相同的控制点/极点并设置所有权重 = 1）。
t_nx = np.array([0.0, 0.0, 0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0, 1.0, 1.0])
反过来，将 splprep 结 (t) 写入 IGES 定义（具有所述“控制点”和权重 = 1）似乎没有给出有效的样条曲线。NX 和至少一个其他软件包无法评估它，理由是“B 样条曲线的修剪或参数值无效”。
在我看来，至少有三种可能性：
- 从非有理 B 样条到有理 B 样条需要进行非平凡的转换
- IGES 样条有一个特定于应用程序的解释（即我对 splprep 输出的解释是正确的，但是在手动绘制/在 IGES 转换例程期间，NX 对此进行了简化/近似）。似乎不太可能。
- splprep 的系数不能以我描述的方式解释为控制点
通过比较所有权重 = 1 的 scipy B 样条（链接）和 IGES NURBS 样条的方程，我已经取消了第一种可能性（链接，第 14 页）。它们看起来相同，正是这一点让我相信 splprep 系数 = 控制点。
任何澄清上述任何点的帮助将不胜感激！
注意，我希望可以表示闭合曲线，所以如果可能的话，我想坚持使用 splprep。
编辑： 我认为首先使用 splrep 尝试这个过程会更简单，因为输出对我来说似乎更直观。我假设返回的系数是控制点的 y 值，但不知道它们对应的 x 位置。因此，我尝试使用这种矩阵方法从样条定义和输入数据中计算它们。C 矩阵只是输入数据。N 矩阵是对每个 x 值的每个基函数的评估，我使用此处显示的（略微修改的）递归函数进行了此操作。然后剩下的就是将 N 取反，并将 C 与它预乘以得到控制点。代码和结果如下：
import numpy as np
import scipy.interpolate as si
# Functions to evaluate B-spline basis functions
def B(x, k, i, t):
if k == 0:
return 1.0 if t[i] <= x < t[i+1] else 0.0
if t[i+k] == t[i]:
c1 = 0.0
else:
c1 = (x - t[i])/(t[i+k] - t[i]) * B(x, k-1, i, t)
if t[i+k+1] == t[i+1]:
c2 = 0.0
else:
c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t)
return c1 + c2
def bspline(x, t, c, k):
n = len(t) - k - 1
assert (n >= k+1) and (len(c) >= n)
cont = []
for i in range(n):
res = B(x, k, i, t)
cont.append(res)
return cont
# Input data
x = np.linspace(-1,1,7)
y = x**3
# B-spline definition
t, c, k = si.splrep(x,y)
# Number of knots = m + 1 = n + k + 2
m = len(t) - 1
# Number of kth degree basis fcns
n = m - k - 1
# Define C and initialise N matrix
C_mat = np.column_stack((x,y))
N_mat = np.zeros(((n+1),(n+1)))
# Calculate basis functions for each x, store in matrix
for i, xs in enumerate(x):
row = bspline(xs, t, c, k)
N_mat[i,:] = row
# Last value must be one...
N_mat[-1,-1] = 1.0
# Invert the matrix
N_inv = np.linalg.inv(N_mat)
# Now calculate control points
P = np.dot(N_inv, C_mat)
导致：
>>> P
array([[ -1.00000000e+00, -1.00000000e+00],
[ -7.77777778e-01, -3.33333333e-01],
[ -4.44444444e-01, -3.29597460e-17],
[ -3.12250226e-17, 8.67361738e-18],
[ 4.44444444e-01, -2.77555756e-17],
[ 7.77777778e-01, 3.33333333e-01],
[ 1.00000000e+00, 1.00000000e+00]])
我认为这是正确的，因为 P 的 y 值与 splrep、c 的系数相匹配。有趣的是，x 值似乎是节点平均值（可以单独计算如下）。也许这个结果对非常熟悉数学的人来说是显而易见的，但对我来说肯定不是。
def knot_average(knots, degree):
"""
Determines knot average vector from knot vector.
:knots: A 1D numpy array describing knots of B-spline.
(NB expected from scipy.interpolate.splrep)
:degree: Integer describing degree of B-spline basis fcns
"""
# Chop first and last vals off
knots_to_average = knots[1:-1]
num_averaged_knots = len(knots_to_average) - degree + 1
knot_averages = np.zeros((num_averaged_knots,))
for i in range(num_averaged_knots):
avg = np.average(knots_to_average[i: i + degree])
knot_averages[i] = avg
return(knot_averages)
现在，要将这些转换为 IGES NURBS，我认为这是定义归一化节点向量、将权重设置为 1 并包括上方的 P 个控制点的情况。我将其标准化如下，并在其下方包含了 IGES 文件。
但是，当我尝试将文件导入 NX 时，它再次无法在定义中说明无效的修剪参数。谁能告诉我这是否是一个有效的 NURBS 定义？
或者这可能是 NX 的一些限制？例如，我注意到当交互式绘制工作室样条曲线时，节点向量被迫（夹紧）均匀（正如 fang 所暗示的那样）。必须要求此约束（并且权重全部 = 1）才能唯一地定义曲线。有趣的是，如果我强制 splrep 使用统一的节点向量（即，夹紧但否则统一）返回样条表示，则会读入 IGES。尽管从 NX 的角度来看，我不认为这是必要的 - 它违背了目的首先要有一个 NURBS。所以这似乎不太可能，我循环想知道我对 splrep 输出的解释是否正确......有人可以指出我哪里出错了吗？
# Original knot vector
>>> t
array([-1. , -1. , -1. , -1. , -0.33333333,
0. , 0.33333333, 1. , 1. , 1. , 1. ])
mini = min(t)
maxi = max(t)
r = maxi - mini
norm_t = (t-mini)/r
# Giving:
>>> norm_t
array([ 0. , 0. , 0. , 0. , 0.33333333,
0.5 , 0.66666667, 1. , 1. , 1. , 1. ])
IGES 定义：
S 1
,,11Hspline_test,13Hsome_path.igs,19HSpline to iges v1.0,4H 0.1,,,,,,, G 1
1.0, 2,2HMM,,,8H 8:58:19,,,,; G 2
126 1 1 1 0 0 0D 1
126 27 4 0 Spline1 1D 2
126,6,3,0,0,1,0,0.0,0.0,0.0,0.0,0.33333,0.5,0.6666666,1.0,1.0,1.0,1.0, 1P 1
1.0,1.0,1.0,1.0,1.0,1.0,1.0,-1.0,-1.0,0.0,-0.7777,-0.33333,0.0, 1P 2
-0.444444,0.0,0.0,0.0,0.0,0.0,0.4444444,0.0,0.0,0.777777777,0.33333, 1P 3
0.0,1.0,1.0,0.0,0.0,1.0,0.0,0.0,0.0,0; 1P 4
S 1G 2D 2P 4 T 1