How to use the splipy.surface_factory function in Splipy

elif kinks[2]-kinks[1] > kinks[3]-kinks[2]:
            corners = [0, kinks[1], (kinks[1]+kinks[2])/2], kinks[2]
        else:
            corners = [0] + kinks[1:3] + [(kinks[2]+kinks[3])/2]

    else:
        while len(kinks) > 5:
            max_span   = 0
            max_span_i = 0
            for i in range(1,len(kinks)-1):
                max_span   = max(max_span, kinks[i+1]-kinks[i-1])
                max_span_i = i
            del kinks[max_span_i]
        corners = kinks[0:4]

    return surface_factory.edge_curves(crv.split(corners))

dim      = int(     self.read_next_non_whitespace().strip())
        r        = float(   next(self.fstream).strip())
        center   = np.array(next(self.fstream).split(), dtype=float)
        z_axis   = np.array(next(self.fstream).split(), dtype=float)
        x_axis   = np.array(next(self.fstream).split(), dtype=float)
        finite   =          next(self.fstream).strip() != '0'
        param_u  = np.array(next(self.fstream).split(), dtype=float)
        if finite:
            param_v=np.array(next(self.fstream).split(), dtype=float)
        else:
            param_v=[-state.unlimited, state.unlimited]
        swap     =          next(self.fstream).strip() != '0'

        center = center + z_axis*param_v[0]
        h      = param_v[1] - param_v[0]
        result = SurfaceFactory.cylinder(r=r, center=center, xaxis=x_axis, axis=z_axis, h=h)
        result.reparam(param_u, param_v)
        if swap:
            result.swap()
        return result

# normalize all values to be in range [0, 1]
    u     = [float(i)/u[-1]     for i in u]
    v     = [float(i)/v[-1]     for i in v]
    knot1 = [float(i)/knot1[-1] for i in knot1]
    knot2 = [float(i)/knot2[-1] for i in knot2]

    # flip and reverse image so coordinate (0,0) is at lower-left corner
    imGrey = imGrey.T  / 255.0
    imGrey = np.flip(imGrey, axis=1)
    x,y    = np.meshgrid(u,v, indexing='ij')
    pts   = np.stack([x,y,imGrey], axis=2)

    basis1 = BSplineBasis(p[0], knot1)
    basis2 = BSplineBasis(p[1], knot2)

    return surface_factory.least_square_fit(pts,[basis1, basis2], [u,v])

def surface_of_linear_extrusion(self):
        dim      = int(      self.read_next_non_whitespace().strip())
        crv      = self.splines(1)
        normal   = np.array(self.read_next_non_whitespace().split(), dtype=float)
        finite   =          next(self.fstream).strip() != '0'
        param_u  = np.array(next(self.fstream).split(), dtype=float)
        if finite:
            param_v=np.array(next(self.fstream).split(), dtype=float)
        else:
            param_v=[-state.unlimited, +state.unlimited]
        swap     =          next(self.fstream).strip() != '0'

        result = SurfaceFactory.extrude(crv + normal*param_v[0], normal*(param_v[1]-param_v[0]))
        result.reparam(param_u, param_v)

        if swap:
            result.swap()
        return result

def trefoil(t):
    t = np.array(t)
    return np.array([sin(t) + 2*sin(2*t), cos(t) - 2*cos(2*t), -sin(3*t)]).T

### do an adaptive best cubic spline-fit of this function
path  = curves.fit(trefoil, 0, 2*pi)

### since we know it is a closed curve, enforce this on the path
path  = path.make_periodic(0,0)

### create a sweeping curve (either a circle or square)
shape = curves.circle(r=0.2)
# shape = 16*curves.n_gon(4)

### sweep *shape along *path
srf = surfaces.sweep(path, shape)

### write results to file. Use meshlab (www.meshlab.net) to view stl-files
with STL('trefoil.stl') as f:
    f.write(srf, n=(150,30))

c = c1.append(c2).append(c3)
c -= c.center()

# plot the reuleaux triangle
t = np.linspace(c.start(0), c.end(0), 151) # 151 parametric evaluation points
x = c(t)                                   # evaluate (x,y)-coordinates
plt.plot(x[:,0], x[:,1])
plt.axis('equal')
plt.show()

# split the triangle in two, and align this with the y-axis
two_parts = c.split((c.start(0) + c.end(0)) / 2.0)
half_curve = two_parts[0].rotate(2*pi/3)

# create a surface by revolving around its symmetry axis (y-axis)
surf = surfaces.revolve(half_curve, axis=[0,1,0])

# plot the resulting surface on a uniform grid of 71x71 evaluation points
u = np.linspace(surf.start(0), surf.end(0), 71)
v = np.linspace(surf.start(1), surf.end(1), 71)
x = surf(u,v)
ax = plt.figure().gca(projection='3d')
ax.plot_wireframe(x[:,:,0], x[:,:,1], x[:,:,2])
plt.show()

# Examples on how to write geometries to file
#
# Author:    Kjetil Andre Johannessen
# Institute: Norwegian University of Science and Technology (NTNU)
# Date:      March 2016
#

from sys import path
path.append('../')
from splipy.io import *
import splipy.surface_factory as surfaces

# create a NURBS torus
torus = surfaces.torus(minor_r=1, major_r=4)


# STL files are tesselated linear triangles. View with i.e. meshlab
with STL('torus.stl') as my_file:
    my_file.write(torus, n=(50, 150)) # specify resolution of 50x150 evaluation pts


# G2 files are native GoTools (http://www.sintef.no/projectweb/geometry-toolkits/gotools/)
with G2('torus.g2') as my_file:
    my_file.write(torus)

def loft(*surfaces):
    if len(surfaces) == 1:
        surfaces = surfaces[0]

    # clone input, so we don't change those references
    # make sure everything has the same dimension since we need to compute length
    surfaces = [s.clone().set_dimension(3) for s in surfaces]
    if len(surfaces)==2:
        return SurfaceFactory.edge_curves(surfaces)
    elif len(surfaces)==3:
        # can't do cubic spline interpolation, so we'll do quadratic
        basis3 = BSplineBasis(3)
        dist  = basis3.greville()
    else:
        x = [s.center() for s in surfaces]

        # create knot vector from the euclidian length between the surfaces
        dist = [0]
        for (x1,x0) in zip(x[1:],x[:-1]):
            dist.append(dist[-1] + np.linalg.norm(x1-x0))

        # using "free" boundary condition by setting N'''(u) continuous at second to last and second knot
        knot = [dist[0]]*4 + dist[2:-2] + [dist[-1]]*4
        basis3 = BSplineBasis(4, knot)

# ngeom    = np.floor(self.n / (p-1))
        # ntexture = np.floor(self.n * n)
        # ngeom    = ngeom.astype(np.int32)
        # ntexture = ntexture.astype(np.int32)
        ngeom    = ensure_listlike(ngeom, 3)
        ntexture = ensure_listlike(ntexture, 3)
        p        = ensure_listlike(p, 3)


        # Create the geometry
        ngx, ngy, ngz = ngeom
        b1 = BSplineBasis(p[0], [0]*(p[0]-1) + [i/ngx for i in range(ngx+1)] + [1]*(p[0]-1))
        b2 = BSplineBasis(p[1], [0]*(p[1]-1) + [i/ngy for i in range(ngy+1)] + [1]*(p[1]-1))
        b3 = BSplineBasis(p[2], [0]*(p[2]-1) + [i/ngz for i in range(ngz+1)] + [1]*(p[2]-1))

        l2_fit = surface_factory.least_square_fit
        vol = self.get_c0_mesh()

        i = slice(irange[0], irange[1], None)
        j = slice(jrange[0], jrange[1], None)
        # special case number of evaluation points for full domain
        if irange[1] == None: irange[1] = vol.shape[0]
        if jrange[1] == None: jrange[1] = vol.shape[1]
        if irange[0] == None: irange[0] = 0
        if jrange[0] == None: jrange[0] = 0

        nu = np.diff(irange)
        nv = np.diff(jrange)
        nw = vol.shape[2]

        u = np.linspace(0, 1, nu)
        v = np.linspace(0, 1, nv)

def sphere(r=1, center=(0,0,0), type='radial'):
    """  Create a solid sphere

    :param float r: Radius
    :param array-like center: Local origin of the sphere
    :param string type: The type of parametrization ('radial' or 'square')
    :return: A solid ball
    :rtype: Volume
    """
    if type == 'radial':
        shell    = SurfaceFactory.sphere(r, center)
        midpoint = shell*0 + center
        return edge_surfaces(shell, midpoint)
    elif type == 'square':
        # based on the work of James E.Cobb: "Tiling the Sphere with Rational Bezier Patches"
        # University of Utah, July 11, 1988. UUCS-88-009
        b = BSplineBasis(order=5)
        sr2 = sqrt(2)
        sr3 = sqrt(3)
        sr6 = sqrt(6)
        cp = [[      -4*(sr3-1),       4*(1-sr3),      4*(1-sr3),   4*(3-sr3)  ], # row 0
              [           -sr2 ,     sr2*(sr3-4),    sr2*(sr3-4), sr2*(3*sr3-2)],
              [              0 ,  4./3*(1-2*sr3), 4./3*(1-2*sr3),4./3*(5-sr3)  ],
              [            sr2 ,     sr2*(sr3-4),    sr2*(sr3-4), sr2*(3*sr3-2)],
              [       4*(sr3-1),       4*(1-sr3),      4*(1-sr3),   4*(3-sr3)  ],
              [    -sr2*(4-sr3),            -sr2,    sr2*(sr3-4), sr2*(3*sr3-2)], # row 1
              [    -(3*sr3-2)/2,     (2-3*sr3)/2,     -(sr3+6)/2,     (sr3+6)/2],