将 scipy.interpolate.splprep 的输出转换为 NURBS 格式以供 IGES 显示
Converting output of scipy.interpolate.splprep into NURBS format for IGES display
我正在寻找将描述任意曲线的一系列有序（非常密集）2D 点转换为 NURBS 表示，可以将其写入 IGES 文件。
我正在使用 scipy.interpolate 的 splprep 来获得给定点系列的 B 样条表示，然后我假设 NURBS 定义基本上是这样加上说所有权重都等于1. 但是我认为我从根本上误解了 splprep 的输出，特别是 'B-spline coefficients' 和在某些 CAD 包中手动重新创建样条曲线所需的控制点之间的关系（我使用的是西门子 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 将结更改为以下（同时显然保持相同的控制 points/poles 并设置所有权重 = 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 定义（'control points' 和权重 = 1）似乎无法给出有效的样条曲线。 NX 和至少一个其他软件包无法对其进行评估，引用 'invalid trim or parametric values for B-spline curve'.
在我看来至少有三种可能性：
- 从非有理 B 样条到有理 B 样条需要一个不平凡的转换
- 有特定于应用程序的 IGES 样条解释（即我对 splprep 输出的解释是正确的，但这是 simplified/approximated 由 NX 手动 drawn/during IGES 转换例程）。似乎不太可能。
- splprep 的系数不能按照我描述的方式解释为控制点
我通过比较 scipy B 样条方程（link) and an IGES NURBS spline with all weights = 1 (link，第 14 页）排除了第一种可能性。它们看起来一模一样，正是这一点让我相信 splprep 系数 = 控制点。
任何帮助澄清以上任何一点将不胜感激！
注意，我想要表示闭合曲线的可能性，所以如果可能的话我想坚持使用 splprep。
编辑：
我认为首先使用 splrep 尝试这个过程会更简单，因为输出对我来说似乎更直观。我假设系数 returned 是控制点的 y 值，但不知道它们对应的 x 位置。因此，我尝试使用 this matrix approach. The C matrix is just the input data. The N matrix is the evaluation of each basis function for each x-value, I did this using the (slightly modified) recursive functions shown here 从样条定义和输入数据中计算它们。然后剩下的就是反转 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 时，它再次失败，说明定义中的 trim 参数无效。谁能告诉我这是否是有效的 NURBS 定义？
或者这可能是 NX 的一些限制？例如，我注意到在交互式绘制工作室样条曲线时，节点矢量被迫（钳制）均匀（如 fang 所暗示的）。必须要求此约束（并且权重全部 = 1）才能唯一定义曲线。有趣的是，如果我强制 splrep 到 return 使用统一结向量（即夹紧但其他方面统一）的样条表示，则读入 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
这个小众查询对其他人有帮助的机会很小——事实证明问题是 IGES 中参数数据部分的格式不正确。描述样条曲线的数据每行不能超过 64 个字符。 splprep 输出的解释是正确的，(c_x, c_y) 数组描述了连续极点的 (x,y) 坐标。等效的 NURBS 定义只需要指定所有权重 = 1。
我正在寻找将描述任意曲线的一系列有序（非常密集）2D 点转换为 NURBS 表示，可以将其写入 IGES 文件。
我正在使用 scipy.interpolate 的 splprep 来获得给定点系列的 B 样条表示，然后我假设 NURBS 定义基本上是这样加上说所有权重都等于1. 但是我认为我从根本上误解了 splprep 的输出，特别是 'B-spline coefficients' 和在某些 CAD 包中手动重新创建样条曲线所需的控制点之间的关系（我使用的是西门子 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 将结更改为以下（同时显然保持相同的控制 points/poles 并设置所有权重 = 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 定义（'control points' 和权重 = 1）似乎无法给出有效的样条曲线。 NX 和至少一个其他软件包无法对其进行评估，引用 'invalid trim or parametric values for B-spline curve'.
在我看来至少有三种可能性：
- 从非有理 B 样条到有理 B 样条需要一个不平凡的转换
- 有特定于应用程序的 IGES 样条解释（即我对 splprep 输出的解释是正确的，但这是 simplified/approximated 由 NX 手动 drawn/during IGES 转换例程）。似乎不太可能。
- splprep 的系数不能按照我描述的方式解释为控制点
我通过比较 scipy B 样条方程（link) and an IGES NURBS spline with all weights = 1 (link，第 14 页）排除了第一种可能性。它们看起来一模一样，正是这一点让我相信 splprep 系数 = 控制点。
任何帮助澄清以上任何一点将不胜感激！
注意，我想要表示闭合曲线的可能性，所以如果可能的话我想坚持使用 splprep。
编辑： 我认为首先使用 splrep 尝试这个过程会更简单，因为输出对我来说似乎更直观。我假设系数 returned 是控制点的 y 值，但不知道它们对应的 x 位置。因此，我尝试使用 this matrix approach. The C matrix is just the input data. The N matrix is the evaluation of each basis function for each x-value, I did this using the (slightly modified) recursive functions shown here 从样条定义和输入数据中计算它们。然后剩下的就是反转 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 时，它再次失败，说明定义中的 trim 参数无效。谁能告诉我这是否是有效的 NURBS 定义？
或者这可能是 NX 的一些限制？例如，我注意到在交互式绘制工作室样条曲线时，节点矢量被迫（钳制）均匀（如 fang 所暗示的）。必须要求此约束（并且权重全部 = 1）才能唯一定义曲线。有趣的是，如果我强制 splrep 到 return 使用统一结向量（即夹紧但其他方面统一）的样条表示，则读入 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
这个小众查询对其他人有帮助的机会很小——事实证明问题是 IGES 中参数数据部分的格式不正确。描述样条曲线的数据每行不能超过 64 个字符。 splprep 输出的解释是正确的，(c_x, c_y) 数组描述了连续极点的 (x,y) 坐标。等效的 NURBS 定义只需要指定所有权重 = 1。