title: Converting output of scipy.interpolate.splprep into NURBS format for IGES display
body: <p>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.  </p>
  
  <p>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). </p>
  
  <p>I've tried a simple example of approximating the function y = x^3 from a sparse set of points:</p>
  
  <pre class="lang-py prettyprint-override"><code>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()
  
  </code></pre>
  
  <p>Which gives the following parameters:</p>
  
  <pre><code>&gt;&gt;&gt; t
  array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.39084883,
          0.5       ,  0.60915117,  1.        ,  1.        ,  1.        ,  1.        ])
  &gt;&gt;&gt; c_x
  array([ -1.00000000e+00,  -9.17992269e-01,  -6.42403598e-01,
          -2.57934892e-16,   6.42403598e-01,   9.17992269e-01,
           1.00000000e+00])
  &gt;&gt;&gt; c_y
  array([ -1.00000000e+00,  -7.12577481e-01,  -6.82922469e-03,
          -1.00363771e-18,   6.82922469e-03,   7.12577481e-01,
           1.00000000e+00])
  &gt;&gt;&gt; k
  3
  &gt;&gt;&gt; u
  array([ 0.        ,  0.25341516,  0.39084883,  0.5       ,  0.60915117,
          0.74658484,  1.        ])
  &gt;&gt;&gt; 
  </code></pre>
  
  <p>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).</p>
  
  <pre><code>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])
  </code></pre>
  
  <p>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'.</p>
  
  <p>There seem to me to be at least three possibilities:</p>
  
  <ol>
  <li>A non-trivial conversion is necessary to go from non-rational to rational B-splines </li>
  <li>There is an application-specific interpretation of IGES splines (i.e. my interpretation of splprep output is correct, but this is simplified/approximated by NX when manually drawn/during the IGES conversion routine).  Seems unlikely.</li>
  <li>The coefficients from splprep cannot be interpreted as control points in the manner I've described</li>
  </ol>
  
  <p>I had written off the first possibility by comparing the equations for a scipy B-spline (<a href="https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.interpolate.BSpline.html#scipy.interpolate.BSpline" rel="nofollow noreferrer">link</a>) and an IGES NURBS spline with all weights = 1 (<a href="https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19950005054.pdf" rel="nofollow noreferrer">link</a>, page 14).  They look identical, and it was this that led me to believe splprep coefficients = control points.</p>
  
  <p><strong>Any</strong> help clarifying any of the above points would be very much appreciated!</p>
  
  <p>NB, I would like the possibility of representing closed curves, so want to stick to splprep if possible.</p>
  
  <hr>
  
  <p><strong>EDIT:</strong>
  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 <a href="https://math.stackexchange.com/questions/3219717/b-spline-curve-how-to-extract-the-control-points-if-knot-vector-is-given">this</a> 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 <a href="https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.interpolate.BSpline.html#scipy.interpolate.BSpline" rel="nofollow noreferrer">here</a>.  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:</p>
  
  <pre class="lang-py prettyprint-override"><code>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] &lt;= x &lt; 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 &gt;= k+1) and (len(c) &gt;= 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)
  </code></pre>
  
  <p>Resulting in:</p>
  
  <pre><code>&gt;&gt;&gt; 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]])
  </code></pre>
  
  <p>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.</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>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.  </p>
  
  <p>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? </p>
  
  <p>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?</p>
  
  <pre><code># Original knot vector
  &gt;&gt;&gt; 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:
  &gt;&gt;&gt; norm_t
  array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.33333333,
          0.5       ,  0.66666667,  1.        ,  1.        ,  1.        ,  1.        ])
  </code></pre>
  
  <p>IGES definition:</p>
  
  <pre><code>                                                                        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
  </code></pre>