I'm looking to convert a series of ordered (pretty dense) 2D points describing arbitrary curves into a NURBS representation, which can be written into an IGES file.
I'm using scipy.interpolate's splprep to get a B-spline representation of the given series of points, and then I had presumed the NURBS definition would essentially be this plus saying all weights are equal to 1. However I think I am fundamentally misinterpreting the output of splprep, specifically the relation between 'B-spline coefficients' and the control points needed to manually recreate the spline in some CAD package (I am using Siemens NX11).
I've tried a simple example of approximating the function y = x^3 from a sparse set of points:
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()
Which gives the following parameters:
>>> 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. ])
>>>
I've assumed that the two sets of coefficients (c_x, c_y) describe the (x,y) coordinates of poles needed to construct the spline. Trying this manually in NX gives a similar spline, though not quite the same, with other points in the interval being evaluated differently than in Python. When I export this manual spline to IGES format, NX changes the knots to the below (while obviously keeping the same control points/poles and setting all weights = 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])
Going the other way and writing the splprep knots (t) into the IGES definition (with said 'control points' and weights = 1) does not seem to give a valid spline. NX and at least one other package cannot evaluate it, citing 'invalid trim or parametric values for B-spline curve'.
There seem to me to be at least three possibilities:
I had written off the first possibility by comparing the equations for a scipy B-spline (link) and an IGES NURBS spline with all weights = 1 (link, page 14). They look identical, and it was this that led me to believe splprep coefficients = control points.
Any help clarifying any of the above points would be very much appreciated!
NB, I would like the possibility of representing closed curves, so want to stick to splprep if possible.
EDIT: I thought it would be simpler to try this process first using splrep, as the outputs seemed more intuitive to me. I assumed the coefficients returned were the y-values of the control points, but didn't know to what x position they corresponded. I therefore tried to calculate them from the spline definition and input data using 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. Then all that remains is to invert N, and pre-multiply C by it to get the control points. The code and result is below:
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)
Resulting in:
>>> 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]])
I think it's correct because the y-values of P match the coefficients from splrep, c. Interestingly the x-values seem to be the knot averages (which could be separately calculated as below). Perhaps this result is obvious to someone very familiar with the maths, it certainly wasn't to me.
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)
Now, to convert these to IGES NURBS I thought it was a case of defining the normalised knot vector, setting the weights all equal to one, and including the P control points from above. I normalised it as below, and have included the IGES file below that.
However when I try to import the file into NX, it again fails stating invalid trim parameters in the definition. Can anyone tell me if this is a valid NURBS definition?
Or perhaps this is some limitation with NX? For instance, I noticed when interactively drawing studio splines the knot vector was forced to be (clamped) uniform (as alluded to by fang). This constraint (and weights all = 1) must be required to uniquely define the curve. Interestingly if I force splrep to return a spline representation using a uniform knot vector (that is, clamped but otherwise uniform), the IGES is read in. I shouldn't think this is necessary though from NXs point of view - it defeats the purpose of having a NURBS in the first place. So it doesn't seem likely and I loop round wondering if my interpretation of the output of splrep is correct...can someone please point out where I've gone wrong?
# 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 definition:
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
On the off chance this niche query helps anyone else- it turns out the problem was incorrect formatting of the parameter data section in the IGES. The data describing the spline can't take up > 64 characters per line. The interpretation of splprep output was correct, the (c_x, c_y) arrays describe the (x,y) coordinates of successive poles. The equivalent NURBS definition just requires specification of all weights = 1.