code: def spline(xyz, s=3, k=2, nest=-1):
      """Generate B-splines as documented in
      https://scipy-cookbook.readthedocs.io/items/Interpolation.html
  
      The scipy.interpolate packages wraps the netlib FITPACK routines
      (Dierckx) for calculating smoothing splines for various kinds of
      data and geometries. Although the data is evenly spaced in this
      example, it need not be so to use this routine.
  
      Parameters
      ----------
      xyz : array, shape (N,3)
          array representing x,y,z of N points in 3d space
      s : float, optional
          A smoothing condition.  The amount of smoothness is determined by
          satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s
          where g(x) is the smoothed interpolation of (x,y).  The user can
          use s to control the tradeoff between closeness and smoothness of
          fit.  Larger satisfying the conditions: sum((w * (y -
          g))**2,axis=0) <= s where g(x) is the smoothed interpolation of
          (x,y).  The user can use s to control the tradeoff between
          closeness and smoothness of fit.  Larger s means more smoothing
          while smaller values of s indicate less smoothing. Recommended
          values of s depend on the weights, w.  If the weights represent
          the inverse of the standard-deviation of y, then a: good s value
          should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is
          the number of datapoints in x, y, and w.
      k : int, optional
          Degree of the spline.  Cubic splines are recommended.  Even
          values of k should be avoided especially with a small s-value.
          for the same set of data.  If task=-1 find the weighted least
          square spline for a given set of knots, t.
      nest : None or int, optional
          An over-estimate of the total number of knots of the spline to
          help in determining the storage space.  None results in value
          m+2*k. -1 results in m+k+1. Always large enough is nest=m+k+1.
          Default is -1.
  
  
      Returns
      -------
      xyzn : array, shape (M,3)
          array representing x,y,z of the M points inside the sphere
  
      Examples
      --------
      >>> import numpy as np
      >>> t=np.linspace(0,1.75*2*np.pi,100)# make ascending spiral in 3-space
      >>> x = np.sin(t)
      >>> y = np.cos(t)
      >>> z = t
      >>> x+= np.random.normal(scale=0.1, size=x.shape) # add noise
      >>> y+= np.random.normal(scale=0.1, size=y.shape)
      >>> z+= np.random.normal(scale=0.1, size=z.shape)
      >>> xyz=np.vstack((x,y,z)).T
      >>> xyzn=spline(xyz,3,2,-1)
      >>> len(xyzn) > len(xyz)
      True
  
      See Also
      --------
      scipy.interpolate.splprep
      scipy.interpolate.splev
      """
      # find the knot points
      tckp, u = splprep(
          [xyz[:, 0], xyz[:, 1], xyz[:, 2]], s=s, k=k, nest=nest
      )
      # evaluate spline, including interpolated points
      xnew, ynew, znew = splev(np.linspace(0, 1, 400), tckp)
      return np.vstack((xnew, ynew, znew)).T
  
imports: from dipy.utils.deprecator import deprecate_with_version
  from dipy.tracking.streamline import set_number_of_points
  import numpy
  from scipy.interpolate import splprep, splev
label_desc: Use image node_id 16 for calling a global function with example usage: spline(xyz, s, k, nest) and returns: np