Я хочу преобразовать серию упорядоченных (довольно плотных) 2D точек, описывающих произвольные кривые, в представление NURBS, которое можно записать в файл IGES.
Я использую spprep scipy.interpolate, чтобы получить представление B-сплайна данной серии точек, а затем я предположил, что определение NURBS по сути будет этим плюсом, говоря, что все веса равны 1. Однако Я думаю, что я в корне неверно истолковываю вывод splprep, в частности, соотношение между «коэффициентами B-сплайна» и контрольными точками, необходимыми для ручного воссоздания сплайна в некотором пакете САПР (я использую 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 и, по крайней мере, еще один пакет не могут оценить его, ссылаясь на «недопустимые параметры усечения или параметры c значений для кривой B-сплайна».
Мне кажется, что есть как минимум три возможности:
- Требуется нетривиальное преобразование в go из нерациональных в рациональные B-сплайны
- Существует специфическая для приложения c интерпретация сплайнов IGES (т.е. моя интерпретация вывода splprep верно, но это упрощается / приближается NX при рисовании вручную / во время процедуры преобразования IGES). Кажется маловероятным.
- Коэффициенты из splprep не могут быть истолкованы как контрольные точки, как я описал.
Первую возможность я списал, сравнив уравнения для сципиона Б. -spline ( link ) и сплайн IGES NURBS со всеми весами = 1 ( link , стр. 14). Они выглядят одинаково, и именно это заставило меня поверить в то, что коэффициенты splprep = контрольные точки.
Любой , помощь в прояснении любого из вышеупомянутых пунктов была бы очень признательна!
NB, я хотел бы иметь возможность представления замкнутых кривых, поэтому, если возможно, хочу придерживаться 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]])
Я думаю, что это правильно, потому что значения y для P соответствуют коэффициентам из splrep, c. Интересно отметить, что значения х, по-видимому, представляют собой средние значения узлов (которые можно рассчитать отдельно, как показано ниже). Возможно, этот результат очевиден для кого-то, очень хорошо знакомого с математикой, это точно не для меня.
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, я подумал, что это был случай определения нормализованного вектора узла , установив все веса, равные единице, и включая контрольные точки P сверху. Я нормализовал это как ниже, и включил файл IGES ниже этого.
Однако, когда я пытаюсь импортировать файл в NX, он снова не в состоянии указать недопустимые параметры обрезки в определении. Может кто-нибудь сказать мне, если это правильное определение NURBS?
Или, может быть, это какое-то ограничение с NX? Например, я заметил, когда в интерактивном режиме рисования сплайнов вектор узла был вынужден быть (зажатым) равномерным (на что ссылается клык). Это ограничение (и weights all = 1) должно быть необходимо для однозначного определения кривой. Интересно, что если я заставлю splrep возвращать представление сплайна, используя равномерный вектор узла (то есть, зажатый, но в остальном равномерный), IGES будет считан. Я не должен думать, что это необходимо, хотя с точки зрения NXs - это противоречит цели иметь NURBS в первую очередь. Так что это маловероятно, и я oop задаюсь вопросом, верна ли моя интерпретация вывода 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