def BlenderToPolyData(me, uvlayer=None):

	pcoords = vtk.vtkFloatArray()
	pcoords.SetNumberOfComponents(3)
	pcoords.SetNumberOfTuples(len(me.verts))
	for i in range(len(me.verts)):
		p0 = me.verts[i].co[0]
		p1 = me.verts[i].co[1]
		p2 = me.verts[i].co[2]
		pcoords.SetTuple3(i, p0, p1, p2)

	points = vtk.vtkPoints()
	points.SetData(pcoords)

	polys = vtk.vtkCellArray()
	lines = vtk.vtkCellArray()
	for face in me.faces:
		if len(face.v) == 4:
			polys.InsertNextCell(4)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
			polys.InsertCellPoint(face.v[3].index)
		elif len(face.v) == 3:
			polys.InsertNextCell(3)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
		elif len(face.v) == 2:
			lines.InsertNextCell(2)
			lines.InsertCellPoint(face.v[0].index)
			lines.InsertCellPoint(face.v[1].index)

	for edge in me.edges:
		lines.InsertNextCell(2)
		lines.InsertCellPoint(edge.v1.index)
		lines.InsertCellPoint(edge.v2.index)

	pdata =vtk.vtkPolyData()
	pdata.SetPoints(points)
	pdata.SetPolys(polys)
	pdata.SetLines(lines)

	if me.faceUV:
		if uvlayer:
			uvnames = me.getUVLayerNames()
			if uvlayer in uvnames:
				me.activeUVLayer = uvlayer
		tcoords = vtk.vtkFloatArray()
		tcoords.SetNumberOfComponents(2)
		tcoords.SetNumberOfTuples(len(me.verts))
		for face in me.faces:
			for i in range(len(face.verts)):
				uv = face.uv[i]
				tcoords.SetTuple2(face.v[i].index, uv[0], uv[1])
		pdata.GetPointData().SetTCoords(tcoords);

	pdata.Update()

	return pdata

def generatePolyData(orientation,fillWith,factor):
   """
   Generate poly-data and point-scalars
   """
   poly = vtk.vtkPolyData()
   pts = vtk.vtkPoints()
   coords=[ (0,0,0),(1,0,0),(1,1,0),(0,1,0)]
   for coord in coords:
      pts.InsertNextPoint(coord[0],coord[1],coord[2])
   poly.SetPoints(pts)

   # Vertices at all corners
   # two 1-point vertices and 1 2-point poly-vertex
   vertices = [[0],[1],[2,3]]
   verts = vtk.vtkCellArray()
   for vertex in vertices:
      InsertCell(verts,vertex,orientation)
   poly.SetVerts(verts)

   # Lines at all sides of the quad
   # two 2-point lines and 1 3-point line
   edges = [ (0,1),(1,2),(2,3,0) ]
   lines = vtk.vtkCellArray()
   for edge in edges:
      InsertCell(lines,edge,orientation)
   poly.SetLines(lines)

   # Fill with one quad, two triangles or a triangle-strip
   if fillWith=='quad':
      quad = (0,1,2,3)
      polys = vtk.vtkCellArray()
      InsertCell(polys,quad,orientation)
      poly.SetPolys(polys)
   elif fillWith=='triangles':
      triangles=[(0,1,3),(3,1,2)]
      strips = vtk.vtkCellArray()
      for triangle in triangles:
         InsertCell(strips,triangle,orientation)
      poly.SetStrips(strips)
   elif fillWith=='strip':
      strip=(0,1,3,2)
      strips = vtk.vtkCellArray()
      InsertCell(strips,strip,orientation)
      poly.SetStrips(strips)

   # Scalars for contouring
   values = [ 0.0, 0.5, 1.5, 1.0 ]
   array=vtk.vtkDoubleArray()
   for v in values:
      array.InsertNextValue(factor*v)
   poly.GetPointData().SetScalars(array)

   return poly

def draw_lines(nodes, color):
    colors = vtk.vtkUnsignedCharArray()
    colors.SetNumberOfComponents(3)
    colors.SetName("Colors")
    cnt = 0
    noderange = 100
    mod = 1
    edges = nodes[0].getedges()
    points = vtk.vtkPoints()
    lines = vtk.vtkCellArray()
    nodecnt = 0

    while cnt < len(nodes):
        node = nodes[cnt]
        cnt += 1
        edges = node.getedges()
        for edge in edges:
            x0,y0,z0 = edge[0]
            x1,y1,z1 = edge[1]

            points.InsertNextPoint(edge[0])
            points.InsertNextPoint(edge[1])

            line = vtk.vtkLine()
            line.GetPointIds().SetId(0,nodecnt)
            line.GetPointIds().SetId(1,nodecnt+1)
            lines.InsertNextCell(line)
            nodecnt += 2
            colors.InsertNextTupleValue(color)

        if cnt % mod == 0:
            print "noderange", noderange, "cnt", cnt, "mod",mod
            # Create a polydata to store everything in
            linesPolyData = vtk.vtkPolyData()
            # Add the points to the dataset
            linesPolyData.SetPoints(points)
            # Add the lines to the dataset
            linesPolyData.SetLines(lines)
            linesPolyData.GetCellData().SetScalars(colors)

            mapper = vtk.vtkPolyDataMapper()
            mapper.SetInput(linesPolyData)
            actor = vtk.vtkActor()
            actor.SetMapper(mapper)
            renderer.AddActor(actor)
            points = vtk.vtkPoints()
            lines = vtk.vtkCellArray()
            nodecnt = 0    
            renderWindow.Render()
            camera = renderer.GetActiveCamera()
            camera.Azimuth(0.1)

    print "done!"

def BlenderToPolyData(me):

	pcoords = vtk.vtkFloatArray()
	pcoords.SetNumberOfComponents(3)
	pcoords.SetNumberOfTuples(len(me.verts))
	for i in range(len(me.verts)):
		p0 = me.verts[i][0]
		p1 = me.verts[i][1]
		p2 = me.verts[i][2]
		pcoords.SetTuple3(i, p0, p1, p2)

	points = vtk.vtkPoints()
	points.SetData(pcoords)

	polys = vtk.vtkCellArray()
	lines = vtk.vtkCellArray()
	for face in me.faces:
		if len(face.v) == 4:
			polys.InsertNextCell(4)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
			polys.InsertCellPoint(face.v[3].index)
		elif len(face.v) == 3:
			polys.InsertNextCell(3)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
		elif len(face.v) == 2:
			lines.InsertNextCell(2)
			lines.InsertCellPoint(face.v[0].index)
			lines.InsertCellPoint(face.v[1].index)

	pdata =vtk.vtkPolyData()
	pdata.SetPoints(points)
	pdata.SetPolys(polys)
	pdata.SetLines(lines)

#------------------------------------------------------------------
# CODIGO INTRODUCIDO PARA PERMITIR LA INCLUSION DE LA             #
# INFORMACION DE LA TEXTURA EN EL POLYDATA GENERADO               #

	if me.hasFaceUV():
		pdata.GetPointData().SetTCoords(calc_coords(me))	# Se insertan las coordenadas en el polydata
#------------------------------------------------------------------

	pdata.Update()

	return pdata

def CreatePolyData( pts, faces ):
    """
    Creates vtkPolyData from vertices and faces
    
    pts numpy.array: Nx3 array of vertices
    faces numpy.array: Mx3 array of faces

    Return vtkPolyData
    """
    (nv,mv) = pts.shape
    (nf,mf) = faces.shape
    cells = vtk.vtkCellArray()
    for j in range(nf):
        cell = vtk.vtkTriangle()
        cell.GetPointIds().SetNumberOfIds(3)
        cell.GetPointIds().SetId( 0, faces[j,0] )
        cell.GetPointIds().SetId( 1, faces[j,1] )
        cell.GetPointIds().SetId( 2, faces[j,2] )
        cells.InsertNextCell( cell )
    
    
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(nv)
    for j in range(nv):
        points.SetPoint( j, pts[j,0], pts[j,1], pts[j,2] )
        
    new_mesh = vtk.vtkPolyData()
    new_mesh.SetPoints( points )
    new_mesh.SetPolys( cells )
    new_mesh.BuildCells()	
    
    return new_mesh

    def __init__(self, pointlist=[]):
        points = vtk.vtkPoints()
        cellArr = vtk.vtkCellArray()
        Colors = vtk.vtkUnsignedCharArray()
        Colors.SetNumberOfComponents(3)
        Colors.SetName("Colors")
        
        n=0
        for p in pointlist:
            vert = vtk.vtkVertex()
            points.InsertNextPoint(p.x, p.y, p.z)
            vert.GetPointIds().SetId(0,n)
            cellArr.InsertNextCell( vert )
            col = clColor(p.cc())
            Colors.InsertNextTuple3( float(255)*col[0], float(255)*col[1], float(255)*col[2] )
            n=n+1
            
        polydata= vtk.vtkPolyData()
        polydata.SetPoints(points)
        polydata.SetVerts( cellArr )
        polydata.GetPointData().SetScalars(Colors)

        polydata.Modified()
        polydata.Update()
        self.src=polydata
        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInput(self.src)
        self.SetMapper(self.mapper)

 def get_cell(self, name):
     try:
         return getattr(self, 'cell_'+name)
     except AttributeError:
         cellarray = vtk.vtkCellArray()
         setattr(self, 'cell_'+name, cellarray)
         return cellarray

    def _update_vtk_objects(self):
        """When n is changed the thus the number of coordinates this function is needed
        to update the vtk objects with the new number of points."""
        # self._vtk_points.SetNumberOfPoints(len(self._points))
        # for i, c in enumerate(self._points):
        #     self._vtk_points.InsertPoint(i, c[0], c[1], c[2])
        self._vtk_points = _vtk.vtkPoints()
        for coordinates in self._points:
            self._vtk_points.InsertNextPoint(coordinates[0], coordinates[1], coordinates[2])

        self._vtk_polygons = _vtk.vtkCellArray()
        for polygon in self._polygons:
            vtk_polygon = _vtk.vtkPolygon()
            vtk_polygon.GetPointIds().SetNumberOfIds(3)
            for local_index, global_index in enumerate(polygon):
                vtk_polygon.GetPointIds().SetId(local_index, global_index)
            self._vtk_polygons.InsertNextCell(vtk_polygon)

        self._vtk_poly_data.SetPoints(self._vtk_points)
        self._vtk_poly_data.SetPolys(self._vtk_polygons)

        self._vtk_scalars = _vtk.vtkFloatArray()
        self._vtk_scalars.SetNumberOfValues(self._vtk_poly_data.GetPoints().GetNumberOfPoints())
        for i in range(self._vtk_scalars.GetNumberOfTuples()):
            self._vtk_scalars.SetValue(i, 0.)

        self._vtk_poly_data.GetPointData().SetScalars(self._vtk_scalars)
        self._vtk_poly_data.Modified()

def SaveParentArtery(centerlines):
   numberOfCells = centerlines.GetNumberOfCells()

   cell0 = centerlines.GetCell(0)
   numberOfArteryPoints = centerlines.GetNumberOfPoints()-cell0.GetNumberOfPoints()
  
   artery = vtk.vtkPolyData()
   arteryPoints = vtk.vtkPoints()
   arteryCellArray = vtk.vtkCellArray()

   radiusArray = vtk.vtkDoubleArray()
   radiusArray.SetName(radiusArrayName)
   radiusArray.SetNumberOfComponents(1)
   radiusArray.SetNumberOfTuples(numberOfArteryPoints)
   radiusArray.FillComponent(0,0.0)

   count = 0
   for i in range(1,numberOfCells): # cell0 is the one that goes to the aneurysm dome
      cell = vtk.vtkGenericCell()
      centerlines.GetCell(i,cell)
   
      arteryCellArray.InsertNextCell(cell.GetNumberOfPoints())
   
      for j in range(cell.GetNumberOfPoints()):
         arteryPoints.InsertNextPoint(cell.GetPoints().GetPoint(j))
         radiusArray.SetTuple1(count,centerlines.GetPointData().GetArray(radiusArrayName).GetTuple1(cell.GetPointId(j)))
         arteryCellArray.InsertCellPoint(count)
         count+=1
         
   artery.SetPoints(arteryPoints)
   artery.SetLines(arteryCellArray)
   artery.GetPointData().AddArray(radiusArray)
   return artery 

 def setEdgesPolydata(self, vd):
     self.edges = []
     self.edges = vd.getEdgesGenerators()
     self.epts = vtk.vtkPoints()
     nid = 0
     lines=vtk.vtkCellArray()
     for e in self.edges:
         p1 = self.scale*e[0]
         p2 = self.scale*e[1] 
         self.epts.InsertNextPoint( p1.x, p1.y, p1.z)
         self.epts.InsertNextPoint( p2.x, p2.y, p2.z)
         line = vtk.vtkLine()
         line.GetPointIds().SetId(0,nid)
         line.GetPointIds().SetId(1,nid+1)
         nid = nid+2
         lines.InsertNextCell(line)
     
     linePolyData = vtk.vtkPolyData()
     linePolyData.SetPoints(self.epts)
     linePolyData.SetLines(lines)
     
     mapper = vtk.vtkPolyDataMapper()
     mapper.SetInput(linePolyData)
     
     self.edge_actor = vtk.vtkActor()
     self.edge_actor.SetMapper(mapper)
     self.edge_actor.GetProperty().SetColor( camvtk.cyan )
     myscreen.addActor( self.edge_actor )
     myscreen.render() 

def draw_lines(nodes, color):
    points = vtk.vtkPoints()
    lines = vtk.vtkCellArray()
    nodecnt = 0
    colors = vtk.vtkUnsignedCharArray()
    colors.SetNumberOfComponents(3)
    colors.SetName("Colors")
    
    for node in nodes: 
        edges = node.getedges()
       
        for edge in edges:
            x0,y0,z0 = edge[0]
            x1,y1,z1 = edge[1]

            points.InsertNextPoint(edge[0])
            points.InsertNextPoint(edge[1])

            line = vtk.vtkLine()
            line.GetPointIds().SetId(0,nodecnt)
            line.GetPointIds().SetId(1,nodecnt+1)
            lines.InsertNextCell(line)
            nodecnt += 2
            colors.InsertNextTupleValue(color)
            
    # Create a polydata to store everything in
    linesPolyData = vtk.vtkPolyData()
    # Add the points to the dataset
    linesPolyData.SetPoints(points)
    # Add the lines to the dataset
    linesPolyData.SetLines(lines)
    linesPolyData.GetCellData().SetScalars(colors)
    return linesPolyData

    def __init__(self, dataShape, interactor):
        self.dataShape = dataShape
        self.planes = []
        self.coordinate = [0,0,0]
        self.lastChangedAxis = -1
        for i in range(3):
            p = vtkImplicitPlaneRepresentation()
            p.SetPlaceFactor(1.0)
            p.OutsideBoundsOn()
            p.ScaleEnabledOff()
            p.SetOrigin(0.25,0.25,0.25)
            p.PlaceWidget([0.1,dataShape[0],0.1,dataShape[1],0.1,dataShape[2]])

            args = [0, 0, 0]
            args[i] = 1
            p.SetNormal(*args)
            p.GetSelectedPlaneProperty().SetColor(*args)
            p.GetEdgesProperty().SetColor(*args) #bug in VTK

            p.GetPlaneProperty().SetOpacity(0.001)
            #do not draw outline
            p.GetOutlineProperty().SetColor(0,0,0)
            p.GetOutlineProperty().SetOpacity(0.0)
            #do not draw normal
            p.GetSelectedNormalProperty().SetOpacity(0.0)
            p.GetNormalProperty().SetOpacity(0.0)
            p.OutlineTranslationOff()
            p.TubingOff()
            
            self.cross = vtkPolyData()
            points = vtkPoints()
            polys = vtkCellArray()
            points.SetNumberOfPoints(6)
            for i in range(3):
                polys.InsertNextCell(2)
                polys.InsertCellPoint(2*i); polys.InsertCellPoint(2*i+1)
            self.cross.SetPoints(points)
            self.cross.SetLines(polys)
            
            pw = vtkImplicitPlaneWidget2()
            pw.SetRepresentation(p)
            pw.SetInteractor(interactor)
            pw.AddObserver("InteractionEvent", self.__PlanePositionCallback)
            
            self.planes.append(pw)
            
        tubes = vtkTubeFilter()
        tubes.SetNumberOfSides(16)
        tubes.SetInput(self.cross)
        tubes.SetRadius(1.0)
        
        crossMapper = vtkPolyDataMapper()
        crossMapper.SetInput(self.cross)
        crossActor = vtkActor()
        crossActor.SetMapper(crossMapper)
        crossActor.GetProperty().SetColor(0,0,0)
        self.AddPart(crossActor)
        
        #initially invoke the event!
        self.InvokeEvent("CoordinatesEvent")

def findPointsInCell(points,
                     cell,
                     verbose=1):

    ugrid_cell = vtk.vtkUnstructuredGrid()
    ugrid_cell.SetPoints(cell.GetPoints())
    cell = vtk.vtkHexahedron()
    for k_point in xrange(8): cell.GetPointIds().SetId(k_point, k_point)
    cell_array_cell = vtk.vtkCellArray()
    cell_array_cell.InsertNextCell(cell)
    ugrid_cell.SetCells(vtk.VTK_HEXAHEDRON, cell_array_cell)

    geometry_filter = vtk.vtkGeometryFilter()
    geometry_filter.SetInputData(ugrid_cell)
    geometry_filter.Update()
    cell_boundary = geometry_filter.GetOutput()

    pdata_points = vtk.vtkPolyData()
    pdata_points.SetPoints(points)

    enclosed_points_filter = vtk.vtkSelectEnclosedPoints()
    enclosed_points_filter.SetSurfaceData(cell_boundary)
    enclosed_points_filter.SetInputData(pdata_points)
    enclosed_points_filter.Update()

    points_in_cell = [k_point for k_point in xrange(points.GetNumberOfPoints()) if enclosed_points_filter.GetOutput().GetPointData().GetArray('SelectedPoints').GetTuple(k_point)[0]]
    return points_in_cell

    def __init__(self, parent = None):
        super(VTKFrame, self).__init__(parent)

        self.vtkWidget = QVTKRenderWindowInteractor(self)
        vl = QtGui.QVBoxLayout(self)
        vl.addWidget(self.vtkWidget)
        vl.setContentsMargins(0, 0, 0, 0)
 
        self.ren = vtk.vtkRenderer()
        self.ren.SetBackground(0.1, 0.2, 0.4)
        self.vtkWidget.GetRenderWindow().AddRenderer(self.ren)
        self.iren = self.vtkWidget.GetRenderWindow().GetInteractor()
 
        # Create source
        # Create five points. 
        origin = [0.0, 0.0, 0.0]
        p0 = [1.0, 0.0, 0.0]
        p1 = [0.0, 1.0, 0.0]
        p2 = [0.0, 1.0, 2.0]
        p3 = [1.0, 2.0, 3.0]
        p4 = [1.0, 2.0, 8.0]
         
        # Create a vtkPoints object and store the points in it
        points = vtk.vtkPoints()
        points.InsertNextPoint(origin)
        points.InsertNextPoint(p0)
        points.InsertNextPoint(p1)
        points.InsertNextPoint(p2)
        points.InsertNextPoint(p3)
        points.InsertNextPoint(p4)
         
        # Create a cell array to store the lines in and add the lines to it
        lines = vtk.vtkCellArray()
         
        for i in range(4):
          line = vtk.vtkLine()
          line.GetPointIds().SetId(0,i)
          line.GetPointIds().SetId(1,i+1)
          lines.InsertNextCell(line)
         
        # Create a polydata to store everything in
        linesPolyData = vtk.vtkPolyData()
         
        # Add the points to the dataset
        linesPolyData.SetPoints(points)
         
        # Add the lines to the dataset
        linesPolyData.SetLines(lines)
         
        # Setup actor and mapper
        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInput(linesPolyData)
         
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
 
        self.ren.AddActor(actor)
        self.ren.ResetCamera()

        self._initialized = False

def edges(indices):
    """
    Maps a numpy ndarray to an vtkCellArray of vtkLines

    Args:
        indices (numpy.ndarray<int>): A numpy.ndarray of shape (n,2) of indices that define n edges

    Returns:
        vtk_lines (vtk.vtkCellArray): VTK representation of the edges
    """
    if not isinstance(indices, numpy.ndarray):
        raise Numpy2VtkFormatException(
            'lines needs numpy array as input'
        )
    if len(indices.shape) != 2 or indices.shape[1] != 2:
        raise Numpy2VtkFormatException(
            'lines needs a nx2 ndarray as input'
        )
    if indices.dtype != numpy.int:
        raise Numpy2VtkFormatException(
            'lines needs to be numpy array of type numpy.int'
        )
    vtk_lines = vtk.vtkCellArray()
    for e in indices:
        line = vtk.vtkLine()
        line.GetPointIds().SetId(0, e[0])
        line.GetPointIds().SetId(1, e[1])
        vtk_lines.InsertNextCell(line)
    return vtk_lines

def displayClickPoints(clickPoints):
	
	points = vtk.vtkPoints()
	lines = vtk.vtkCellArray()
	polygon = vtk.vtkPolyData()
	polygonMapper = vtk.vtkPolyDataMapper()
	polygonActor = vtk.vtkActor()
		
	points.SetNumberOfPoints(4)
	points.SetPoint(0, 0.0, -1.0, 0.0)
	points.SetPoint(1, -0.7, -0.5, 0.0)
	points.SetPoint(2,   0.7,  0.5, 0.0)
	points.SetPoint(3, 0.0,  -1.0, 0.0)

		
	lines.InsertNextCell(4)
	lines.InsertCellPoint(0)
	lines.InsertCellPoint(1)
	lines.InsertCellPoint(2)
	lines.InsertCellPoint(3)
		
	polygon.SetPoints(points)
	polygon.SetLines(lines)
		
	polygonMapper.SetInputConnection(polygon.GetProducerPort())
	polygonActor.SetMapper(polygonMapper)	
	
	mpv_renderer[1].ResetCamera()
	mpv_renderer[1].AddActor(polygonActor)
	rw.Render()

	def makeEdgeVTKObject(mesh,model):
		"""
		Make and return a edge based VTK object for a simpeg mesh and model.

		Input:
		:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
		:param model, dictionary of numpy.array - Name('s) and array('s).
			Property array must be order hstack(Ex,Ey,Ez)

		Output:
        :rtype: vtkUnstructuredGrid object
        :return: vtkObj
		"""

		## Convert simpeg mesh to VTK properties
		# Convert mesh nodes to vtkPoints
		vtkPts = vtk.vtkPoints()
		vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1))

		# Define the face "cells"
		# Using VTK_QUAD cell for faces (see VTK file format)
		nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M')
		def edgeR(mat,length):
			return mat.T.reshape((length,1))
		# First direction
		nTEx = np.prod(mesh.nEx)
		ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)])
		# Second direction
		if mesh.dim >= 2:
			nTEy = np.prod(mesh.nEy)
			EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)])
		# Third direction
		if mesh.dim == 3:
			nTEz = np.prod(mesh.nEz)
			EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)])
		# Cells -cell array
		ECellArr = vtk.vtkCellArray()
		ECellArr.SetNumberOfCells(mesh.nE)
		ECellArr.SetCells(mesh.nE,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1))
		# Cell type
		ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1)
		# Cell location
		ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nE*3,3,dtype='int64'),deep=1)

		## Make the object
		vtkObj = vtk.vtkUnstructuredGrid()
		# Set the objects properties
		vtkObj.SetPoints(vtkPts)
		vtkObj.SetCells(ECellType,ECellLoc,ECellArr)

		# Assign the model('s) to the object
		for item in model.iteritems():
			# Convert numpy array
			vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
			vtkDoubleArr.SetName(item[0])
			vtkObj.GetCellData().AddArray(vtkDoubleArr)

		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
		vtkObj.Update()
		return vtkObj

def ExtractSingleLine(centerlines,id):
   cell = vtk.vtkGenericCell()
   centerlines.GetCell(id,cell) 

   line = vtk.vtkPolyData()
   points = vtk.vtkPoints()
   cellArray = vtk.vtkCellArray()
   cellArray.InsertNextCell(cell.GetNumberOfPoints())

   radiusArray = vtk.vtkDoubleArray()
   radiusArray.SetName(radiusArrayName)
   radiusArray.SetNumberOfComponents(1)
   radiusArray.SetNumberOfTuples(cell.GetNumberOfPoints())
   radiusArray.FillComponent(0,0.0)

   for i in range(cell.GetNumberOfPoints()):
      point = [0.0,0.0,0.0]
      point = cell.GetPoints().GetPoint(i)

      points.InsertNextPoint(point)
      cellArray.InsertCellPoint(i)
      radius = centerlines.GetPointData().GetArray(radiusArrayName).GetTuple1(cell.GetPointId(i))
      radiusArray.SetTuple1(i,radius)

   line.SetPoints(points)   
   line.SetLines(cellArray)
   line.GetPointData().AddArray(radiusArray)

   return line 

def VTKPoints2PolyData( pts ):
    """
    Transforms numpy point array into  vtkPolyData
    
    pts numpy.array: Nx3 array of vertices

    Return vtkPolyData
    """
    (nv,mv) = pts.shape
    
    
    vertices = vtk.vtkCellArray()
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(nv)
    for j in range(nv):
        points.SetPoint( j, pts[j,0], pts[j,1], pts[j,2] )
        vertices.InsertNextCell(1)
        vertices.InsertCellPoint(j)    
 
    # Create a polydata object
    mesh = vtk.vtkPolyData()
 
    # Set the points and vertices we created as the geometry and topology of the polydata
    mesh.SetPoints(points)
    mesh.SetVerts(vertices)

    mesh.BuildCells()	
    
    return mesh

def skippoints(polydata,nskippoints):
    """Generate a single cell line from points in idlist."""

    # derive number of nodes
    numberofnodes = polydata.GetNumberOfPoints() - nskippoints

    # define points and line
    points = vtk.vtkPoints()
    polyline = vtk.vtkPolyLine()
    polyline.GetPointIds().SetNumberOfIds(numberofnodes)

    # assign id and x,y,z coordinates
    for i in range(nskippoints,polydata.GetNumberOfPoints()):
        pointid = i - nskippoints
        polyline.GetPointIds().SetId(pointid,pointid)
        point = polydata.GetPoint(i)
        points.InsertNextPoint(point)


    # define cell
    cells = vtk.vtkCellArray()
    cells.InsertNextCell(polyline)

    # add to polydata
    polyout = vtk.vtkPolyData()
    polyout.SetPoints(points)
    polyout.SetLines(cells)

    return polyout