def GetTree(self):
		"""Returns a full vtkTree based on data loaded in LoadData()."""

		if self.data_loaded:

			vertex_id = vtk.vtkIdTypeArray()
			vertex_id.SetName('vertex_ids')
			for ii in range(len(self.cp)):
				vertex_id.InsertNextValue(ii)

			NINvtk = VN.numpy_to_vtk(self.NumberInNet, deep=True)
			NINvtk.SetName('num_in_vertex')
			SCALESvtk = VN.numpy_to_vtk(self.Scales, deep=True)
			SCALESvtk.SetName('scale')

			# This array will default to empty strings
			BLANKvtk = vtk.vtkStringArray()
			BLANKvtk.SetNumberOfComponents(1)
			BLANKvtk.SetNumberOfTuples(self.NumberInNet.shape[0])
			BLANKvtk.SetName('blank')

			# Build tree out of CP list of "is a child of"
			#	remembering that Matlab indices are 1-based and numpy/VTK 0-based
			print 'Building graph'
			dg = vtk.vtkMutableDirectedGraph()
			edge_id = vtk.vtkIdTypeArray()
			edge_id.SetName('edge_ids')
			for ii in range(self.cp.size):
				dg.AddVertex()
			for ii in range(self.cp.size):
				if self.cp[ii] > 0:		# CP already zero-based
					dg.AddGraphEdge(self.cp[ii],ii)		# Method for use with wrappers -- AddEdge() in C++
					edge_id.InsertNextValue(ii)

			dg.GetVertexData().AddArray(NINvtk)
			dg.GetVertexData().AddArray(SCALESvtk)
			dg.GetVertexData().AddArray(vertex_id)
			dg.GetVertexData().SetActiveScalars('scale')
			dg.GetVertexData().SetActivePedigreeIds('vertex_ids')
			dg.GetEdgeData().AddArray(edge_id)
			dg.GetEdgeData().SetActivePedigreeIds('edge_ids')

			tree = vtk.vtkTree()
			tree.CheckedShallowCopy(dg)

			return tree

		else:
			raise IOError, "Can't get tree until data is loaded successfully"

    def selectCell(self, cellId):
        if cellId in (None, -1):
            return
        ids = vtk.vtkIdTypeArray();
        ids.SetNumberOfComponents(1);
        ids.InsertNextValue(cellId);

        selectionNode = vtk.vtkSelectionNode();
        selectionNode.SetFieldType(vtk.vtkSelectionNode.CELL);
        selectionNode.SetContentType(vtk.vtkSelectionNode.INDICES);
        selectionNode.SetSelectionList(ids);

        selection = vtk.vtkSelection();
        selection.AddNode(selectionNode);

        extractSelection = vtk.vtkExtractSelection();

        extractSelection.SetInputData(0, self.actor.GetMapper().GetInput());
        extractSelection.SetInputData(1, selection);

        extractSelection.Update();

        selected = vtk.vtkUnstructuredGrid();
        selected.ShallowCopy(extractSelection.GetOutput());

        self.selectedMapper.SetInputData(selected);
        self.selectedMapper.Update()

def numpy_to_vtkIdTypeArray(num_array, deep=0):
    isize = vtk.vtkIdTypeArray().GetDataTypeSize()
    dtype = num_array.dtype
    if isize == 4:
        if dtype != numpy.int32:
            raise ValueError, \
             'Expecting a numpy.int32 array, got %s instead.' % (str(dtype))
    else:
        if dtype != numpy.int64:
            raise ValueError, \
             'Expecting a numpy.int64 array, got %s instead.' % (str(dtype))
        
    return numpy_to_vtk(num_array, deep, vtkConstants.VTK_ID_TYPE)

def numpy2vtkDataArrayInt(npa):
#    print npa[1][0]
    size0, size1 = npa.shape
    data = vtk.vtkIdTypeArray()
    data.SetNumberOfComponents(4)
#    data.SetName("CELLS")
    for i in range(size0):
        n0 = int(npa[i, 0])
        n1 = int(npa[i, 1])
        n2 = int(npa[i, 2])
        n3 = int(npa[i, 3])
        data.InsertNextTuple4(n0, n1, n2, n3)
    
    return data

        def __init__(self, times, slider_repres):
            self._stimes = set(times)
            self._opacity = 1.0
            self._time_step = (max(self._stimes) - min(self._stimes)) \
                               / len(self._stimes)
            self._time = min(times)
            self._slider_repres = slider_repres
            self._current_id = vtk.vtkIdTypeArray()
            self._renderer = renderer
            self._renderer_window = renderer_window
            self._times = times
            self._image_counter = 0

            self._recording = False

 def WriteTecplotMeshFile(self):
     if (self.OutputFileName == ''):
         self.PrintError('Error: no OutputFileName.')
     self.PrintLog('Writing Tecplot file.')
     triangleFilter = vtk.vtkDataSetTriangleFilter()
     triangleFilter.SetInputData(self.Mesh)
     triangleFilter.Update()
     self.Mesh = triangleFilter.GetOutput()
     f=open(self.OutputFileName, 'w')
     line = "VARIABLES = X,Y,Z"
     arrayNames = []
     for i in range(self.Mesh.GetPointData().GetNumberOfArrays()):
         array = self.Mesh.GetPointData().GetArray(i)
         arrayName = array.GetName()
         if arrayName == None:
             continue
         if (arrayName[-1]=='_'):
             continue
         arrayNames.append(arrayName)
         if (array.GetNumberOfComponents() == 1):
             line = line + ',' + arrayName
         else:
             for j in range(array.GetNumberOfComponents()):
                 line = line + ',' + arrayName + str(j)
     line = line + '\n'
     f.write(line)
     tetraCellIdArray = vtk.vtkIdTypeArray()
     tetraCellType = 10
     self.Mesh.GetIdsOfCellsOfType(tetraCellType,tetraCellIdArray)
     numberOfTetras = tetraCellIdArray.GetNumberOfTuples()
     line = "ZONE " + "N=" + str(self.Mesh.GetNumberOfPoints()) + ',' + "E=" + str(numberOfTetras) + ',' + "F=FEPOINT" + ',' + "ET=TETRAHEDRON" + '\n'
     f.write(line)
     for i in range(self.Mesh.GetNumberOfPoints()):
         point = self.Mesh.GetPoint(i)
         line = str(point[0]) + ' ' + str(point[1]) + ' ' + str(point[2])
         for arrayName in arrayNames:
             array = self.Mesh.GetPointData().GetArray(arrayName)
             for j in range(array.GetNumberOfComponents()):
                 line = line + ' ' + str(array.GetComponent(i,j))
         line = line + '\n'
         f.write(line)
     for i in range(numberOfTetras):
         cellPointIds = self.Mesh.GetCell(tetraCellIdArray.GetValue(i)).GetPointIds()
         line = ''
         for j in range(cellPointIds.GetNumberOfIds()):
             if (j>0):
                 line = line + ' '
             line = line + str(cellPointIds.GetId(j)+1)
         line = line + '\n'
         f.write(line)

def GetSelection(ids, inverse=False):
    selNode = vtk.vtkSelectionNode()
    selNode.SetContentType(vtk.vtkSelectionNode.BLOCKS)

    idArray = vtk.vtkIdTypeArray()
    idArray.SetNumberOfTuples(len(ids))
    for i in range(len(ids)):
        idArray.SetValue(i, ids[i])
    selNode.SetSelectionList(idArray)

    if inverse:
        selNode.GetProperties().Set(vtk.vtkSelectionNode.INVERSE(), 1)

    sel = vtk.vtkSelection()
    sel.AddNode(selNode)

    return sel

	def picker_callback(self,obj,event):
		
		extract = vtk.vtkExtractSelectedFrustum()
	
		fPlanes=obj.GetFrustum() #collection of planes based on unscaled display
	
		#scale frustum to account for the zaspect
		scaledPlanes=vtk.vtkPlanes()
		scaledNormals=vtk.vtkDoubleArray()
		scaledNormals.SetNumberOfComponents(3)
		scaledNormals.SetNumberOfTuples(6)
		scaledOrigins=vtk.vtkPoints()
		for j in range(6):
			i=fPlanes.GetPlane(j)
			k=i.GetOrigin()
			q=i.GetNormal()
			scaledOrigins.InsertNextPoint(k[0],k[1],k[2]/float(self.Zaspect))
			scaledNormals.SetTuple(j,(q[0],q[1],q[2]*float(self.Zaspect)))
		scaledPlanes.SetNormals(scaledNormals)
		scaledPlanes.SetPoints(scaledOrigins)
			
		
		extract.SetFrustum(scaledPlanes)
		extract.SetInputData(self.vtkPntsPolyData)
		extract.Update()
		extracted = extract.GetOutput()
		
		ids = vtk.vtkIdTypeArray()
		ids = extracted.GetPointData().GetArray("vtkOriginalPointIds")

		
		if ids:
			#store them in an array for an undo operation
			self.lastSelectedIds=ids
			for i in range(ids.GetNumberOfTuples()):
				#turn them red
				self.colors.SetTuple(ids.GetValue(i),(255,0,0))
				self.bool_pnt[ids.GetValue(i)]=False
		
			self.vtkPntsPolyData.GetPointData().SetScalars(self.colors)
			self.vtkPntsPolyData.Modified()
		
		
		self.ui.vtkWidget.update()
		#set flag on ui to show that data has been modified
		self.unsaved_changes=True

    def setPoints( self, point_data ):
        self.vtkPoints.SetData( point_data )
        ncells = point_data.GetNumberOfTuples()               
        cells = vtk.vtkIdTypeArray()
        cell_data_array = np.empty( 2*ncells, dtype=np.int64 )  
        cell_data_array[1:2*ncells:2] = range( ncells )
        cell_data_array[0:2*ncells:2] = 1
#        cell_data_array = np.array( range( ncells ), dtype=np.int64 )  
        cells.SetVoidArray( cell_data_array, 2*ncells, 1 )
        self.vtkCells.SetCells ( ncells, cells )
#        self.vtkCells.InsertNextCell( cell )
#        for iCell in range( ncells ):
#            self.vtkCells.InsertNextCell( 1 )
#            self.vtkCells.InsertCellPoint( iCell )
        cellData = self.vtkCells.GetData ()
        print " Cell Data:\n %s " % str( [ cellData.GetValue(iCell) for iCell in range(cellData.GetNumberOfTuples())]  )
        self.vtkCells.Modified()
        self.vtkPoints.Modified()

def trimesh_to_vtk(trimesh):
    r"""Return a `vtkPolyData` representation of a :map:`TriMesh` instance

    Parameters
    ----------
    trimesh : :map:`TriMesh`
        The menpo :map:`TriMesh` object that needs to be converted to a
        `vtkPolyData`

    Returns
    -------
    `vtk_mesh` : `vtkPolyData`
        A VTK mesh representation of the Menpo :map:`TriMesh` data

    Raises
    ------
    ValueError:
        If the input trimesh is not 3D.
    """
    import vtk
    from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray
    if trimesh.n_dims != 3:
        raise ValueError('trimesh_to_vtk() only works on 3D TriMesh instances')

    mesh = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    points.SetData(numpy_to_vtk(trimesh.points, deep=1))
    mesh.SetPoints(points)

    cells = vtk.vtkCellArray()

    # Seemingly, VTK may be compiled as 32 bit or 64 bit.
    # We need to make sure that we convert the trilist to the correct dtype
    # based on this. See numpy_to_vtkIdTypeArray() for details.
    isize = vtk.vtkIdTypeArray().GetDataTypeSize()
    req_dtype = np.int32 if isize == 4 else np.int64
    cells.SetCells(trimesh.n_tris,
                   numpy_to_vtkIdTypeArray(
                       np.hstack((np.ones(trimesh.n_tris)[:, None] * 3,
                                  trimesh.trilist)).astype(req_dtype).ravel(),
                       deep=1))
    mesh.SetPolys(cells)
    return mesh

                neighb = find_celledge_neighbors(tri_id, tri)

                #print('Neighbors ', neighb)

                for j in range( len(neighb) ):
                    if trilabel[neighb[j]]==0:
                        #see if the triangles tri_id and neighb[j] are on the different sides of the line
                        #i.e. if they share any pair of points of the line
                        if not triangles_on_any_line(tri_id, neighb[j], tri, lines):
                            tri_stack.append(neighb[j])


# In[11]:


array = vtk.vtkIdTypeArray()
array.SetName(array_name)
array.SetNumberOfComponents(1)
array.SetNumberOfTuples(trilabel.shape[0])

for i in range(trilabel.shape[0]):
    array.SetTuple1(i, trilabel[i])
    
mesh.GetCellData().AddArray(array)

wr = vtk.vtkPolyDataWriter()
wr.SetFileName(output_mesh)
wr.SetInputData(mesh)
wr.Write()

       xM = data2[:,1].max()
       ym = data2[:,0].min()
       yM = data2[:,0].max()
       N = data2.shape[0]
       m2 = numpy.ascontiguousarray(numpy.transpose(data2,(0,2,1)))
       nVertices = m2.shape[-2]
       m2.resize((m2.shape[0]*m2.shape[1],m2.shape[2]))
       m2=m2[...,::-1]
       # here we add dummy levels, might want to reconsider converting "trimData" to "reOrderData" and use actual levels?
       m3=numpy.concatenate((m2,numpy.zeros((m2.shape[0],1))),axis=1)
 except Exception,err: # Ok no mesh on file, will do with lat/lon
   pass
 if m3 is not None:
   #Create unstructured grid points
   vg = vtk.vtkUnstructuredGrid()
   lst = vtk.vtkIdTypeArray()
   cells = vtk.vtkCellArray()
   numberOfCells = N
   lst.SetNumberOfComponents(nVertices + 1)
   lst.SetNumberOfTuples(numberOfCells)
   for i in range(N):
     tuple = [None] * (nVertices + 1)
     tuple[0] = nVertices
     for j in range(nVertices):
       tuple[j + 1] = i * nVertices + j
     lst.SetTuple(i, tuple)
     ## ??? TODO ??? when 3D use CUBE?
   cells.SetCells(numberOfCells, lst)
   vg.SetCells(vtk.VTK_POLYGON, cells)
 else:
   #Ok a simple structured grid is enough

    def WriteLifeVMeshFile(self):
        if (self.OutputFileName == ''):
            self.PrintError('Error: no OutputFileName.')
        self.PrintLog('Writing LifeV file.')

        self.Mesh.BuildLinks()

        cellEntityIdsArray = vtk.vtkIntArray()
        cellEntityIdsArray.DeepCopy(self.Mesh.GetCellData().GetArray(self.CellEntityIdsArrayName))

        tetraCellType = 10
        triangleCellType = 5

        f=open(self.OutputFileName, 'w')
        line = "MeshVersionFormatted 1\n\n"
        line += "Dimension\n"
        line += "3\n\n"
        line += "Vertices\n"
        line += "%d\n" % self.Mesh.GetNumberOfPoints()
        f.write(line)
        for i in range(self.Mesh.GetNumberOfPoints()):
            point = self.Mesh.GetPoint(i)
            pointCells = vtk.vtkIdList()
            self.Mesh.GetPointCells(i,pointCells)
            minTriangleCellEntityId = -1
            tetraCellEntityId = -1
            for j in range(pointCells.GetNumberOfIds()):
                cellId = pointCells.GetId(j)
                if self.Mesh.GetCellType(cellId) == triangleCellType:
                    cellEntityId = cellEntityIdsArray.GetValue(cellId)
                    if cellEntityId < minTriangleCellEntityId or minTriangleCellEntityId == -1:
                        minTriangleCellEntityId = cellEntityId
                else:
                    tetraCellEntityId = cellEntityIdsArray.GetValue(cellId)
            cellEntityId = tetraCellEntityId
            if minTriangleCellEntityId != -1:
                cellEntityId = minTriangleCellEntityId
            line = "%f  %f  %f  %d\n" % (point[0], point[1], point[2], cellEntityId)
            f.write(line)
        line = "\n"

        tetraCellIdArray = vtk.vtkIdTypeArray()
        self.Mesh.GetIdsOfCellsOfType(tetraCellType,tetraCellIdArray)
        numberOfTetras = tetraCellIdArray.GetNumberOfTuples()
        line += "Tetrahedra\n"
        line += "%d\n" % numberOfTetras
        f.write(line)
        for i in range(numberOfTetras):
            tetraCellId = tetraCellIdArray.GetValue(i) 
            cellPointIds = self.Mesh.GetCell(tetraCellId).GetPointIds()
            line = ''
            for j in range(cellPointIds.GetNumberOfIds()):
                if j>0:
                    line += '  '
                line += "%d" % (cellPointIds.GetId(j)+1)
            line += '  %d\n' % 1
            f.write(line)
        line = "\n"

        triangleCellIdArray = vtk.vtkIdTypeArray()
        self.Mesh.GetIdsOfCellsOfType(triangleCellType,triangleCellIdArray)
        numberOfTriangles = triangleCellIdArray.GetNumberOfTuples()
        line += "Triangles\n"
        line += "%d\n" % numberOfTriangles
        f.write(line)
        for i in range(numberOfTriangles):
            triangleCellId = triangleCellIdArray.GetValue(i)
            cellPointIds = self.Mesh.GetCell(triangleCellId).GetPointIds()
            line = ''
            for j in range(cellPointIds.GetNumberOfIds()):
                if j>0:
                    line += '  '
                line += "%d" % (cellPointIds.GetId(j)+1)
            cellEntityId = cellEntityIdsArray.GetValue(triangleCellId)
            line += '  %d\n' % cellEntityId
            f.write(line)

def ExtractSelection(fileList):
    # Report our CWD just for testing purposes.
    print "CWD:", os.getcwd()

    ringOffset = numQuadsPerRing * numSMCsPerCol * 4
    branchOffset = numRingsPerBranch * ringOffset

    # Prepare selection id array.    
    selectionIds = vtk.vtkIdTypeArray()
    
    for branch in branches:
        thisBranchOffset = branch * branchOffset
        # print thisBranchOffset,
    
        for ring in range(numRingsPerBranch):
            thisRingOffset = ring * ringOffset
            # print thisRingOffset, 
            
            for cell in range(numSMCsPerCol):
                thisOffset = start + thisBranchOffset + thisRingOffset + cell
                # print thisOffset,
                selectionIds.InsertNextValue(thisOffset)
    
    # Create selection node.
    selectionNode = vtk.vtkSelectionNode()
    selectionNode.SetFieldType(selectionNode.CELL)
    selectionNode.SetContentType(selectionNode.INDICES)
    selectionNode.SetSelectionList(selectionIds)
    
    # Create selection.
    selection = vtk.vtkSelection()
    selection.AddNode(selectionNode)
    
    sortNicely(fileList)
    
    # Process every file by extracting selection.
    for inFile in fileList:
        print 'Reading file', inFile
        reader = vtk.vtkXMLUnstructuredGridReader()
        reader.SetFileName(inFile)
        
        print '\tExtracting ', selectionIds.GetNumberOfTuples(), 'cells...'
        selectionExtractor = vtk.vtkExtractSelection()
        selectionExtractor.SetInputConnection(reader.GetOutputPort())
        if vtk.vtkVersion().GetVTKMajorVersion() > 5:
            selectionExtractor.SetInputData(1, selection)
        else:
            selectionExtractor.SetInput(1, selection)
        
        baseName = os.path.basename(inFile)
        number = [int(s) for s in re.split('([0-9]+)', baseName) if s.isdigit()]
        outFile = outputPattern + str(number[0]) + '.vtu'
        
        print '\t\tSaving file', outFile
        writer = vtk.vtkXMLUnstructuredGridWriter()
        writer.SetInputConnection(selectionExtractor.GetOutputPort())
        writer.SetFileName(outFile)
        writer.Update()
        
        print '\n'
    
    print 'All done...'

 - You need to make sure you hold a reference to a Numpy array you want
   to import into VTK.  If not you'll get a segfault (in the best case).
   The same holds in reverse when you convert a VTK array to a numpy
   array -- don't delete the VTK array.


Created by Prabhu Ramachandran in Feb. 2008.
"""

import vtk
import vtkConstants
import numpy

# Useful constants for VTK arrays.
VTK_ID_TYPE_SIZE = vtk.vtkIdTypeArray().GetDataTypeSize()
if VTK_ID_TYPE_SIZE == 4:
    ID_TYPE_CODE = numpy.int32
elif VTK_ID_TYPE_SIZE == 8:
    ID_TYPE_CODE = numpy.int64

VTK_LONG_TYPE_SIZE = vtk.vtkLongArray().GetDataTypeSize()
if VTK_LONG_TYPE_SIZE == 4:
    LONG_TYPE_CODE = numpy.int32
    ULONG_TYPE_CODE = numpy.uint32
elif VTK_LONG_TYPE_SIZE == 8:
    LONG_TYPE_CODE = numpy.int64
    ULONG_TYPE_CODE = numpy.uint64


def get_vtk_array_type(numpy_array_type):

    def GenerateTetrInFile(self):

        self.PrintLog('Generating Tetr .in file.')
        
        f=open(self.OutputFileName, 'w')

        line = '$title' + '\n'
        f.write(line)
        line = self.OutputFileName + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$compress' + '\n'
        f.write(line)
        line = 'gzip -f' + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        self.NormalizationTransform.Identity()

        if (self.NormalizationEntity != ''):
            sectionBoundaryPointIds = self.GetSectionBoundaryPointIds(self.NormalizationEntity)
            sectionProperties = SectionProperties()
            sectionProperties.Mesh = self.Mesh
            sectionProperties.NormalizationTransform = self.NormalizationTransform
            sectionProperties.SectionBoundaryPointIds = sectionBoundaryPointIds
            sectionProperties.Execute()
            self.NormalizationRadius = sectionProperties.Radius

        if (self.NormalizationRadius != 0.0) & (self.NormalizationRadius != 1.0):
            self.NormalizationTransform.Scale(1.0/self.NormalizationRadius,1.0/self.NormalizationRadius,1.0/self.NormalizationRadius)

        line = '$radius' + '\n'
        f.write(line)
        line = str(self.NormalizationRadius) + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$viscosity' + '\n'
        f.write(line)
        line = str(0.035) + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$density' + '\n'
        f.write(line)
        line = str(1.06) + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$red' + '\n'
        f.write(line)
        line = str(0) + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$alpha' + '\n'
        f.write(line)
        line = str(0) + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$nsolve' + '\n'
        f.write(line)
        line = '22' + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$Uzawa_PC' + '\n'
        f.write(line)
        line = '1' + '\n'
        f.write(line)
        line = '\n'
        f.write(line)

        line = '$node' + '\n'
        f.write(line)
        line = str(self.Mesh.GetNumberOfPoints()) + '\n'
        f.write(line)
        for i in range(self.Mesh.GetNumberOfPoints()):
            point = [0.0,0.0,0.0]
            self.NormalizationTransform.TransformPoint(self.Mesh.GetPoint(i),point)
            line = str(i+1) + ' ' + str(point[0]) + ' ' + str(point[1]) + ' ' + str(point[2]) + '\n'
            f.write(line)

        line = '\n'
        f.write(line)
            
        line = '$elem' + '\n'
        f.write(line)
        quadratidTetraCellType = 24
        quadraticTetraCellIds = vtk.vtkIdTypeArray()
#.........这里部分代码省略.........

def main():
    pointSource = vtk.vtkPointSource()
    pointSource.SetNumberOfPoints(50)
    pointSource.Update()

    print("There are %s input points\n" % pointSource.GetOutput().GetNumberOfPoints())

    ids = vtk.vtkIdTypeArray()
    ids.SetNumberOfComponents(1)

    # Set values
    i = 10
    while i < 20:
        ids.InsertNextValue(i)
        i += 1

    selectionNode = vtk.vtkSelectionNode()
    selectionNode.SetFieldType(1) # POINT
    #  CELL_DATA = 0
    #  POINT_DATA = 1
    #  FIELD_DATA = 2
    #  VERTEX_DATA = 3
    #  EDGE_DATA = 4

    selectionNode.SetContentType(4) # INDICES
    #SELECTIONS = 0
    #GLOBALIDS = 1
    #PEDIGREEIDS = 2
    #VALUES = 3
    #INDICES = 4
    #FRUSTUM = 5
    #LOCATIONS = 6
    #THRESHOLDS = 7
    #BLOCKS = 8
    selectionNode.SetSelectionList(ids)

    selection = vtk.vtkSelection()
    selection.AddNode(selectionNode)

    extractSelection = vtk.vtkExtractSelection()

    extractSelection.SetInputConnection(0, pointSource.GetOutputPort())
    if vtk.VTK_MAJOR_VERSION <= 5:
        extractSelection.SetInput(1, selection)
    else:
        extractSelection.SetInputData(1, selection)
    extractSelection.Update()

    # In selection
    selected = vtk.vtkUnstructuredGrid()
    selected.ShallowCopy(extractSelection.GetOutput())

    print("There are %s points in the selection" % selected.GetNumberOfPoints())
    print("There are %s cells in the selection" %selected.GetNumberOfCells())

    # Get points that are NOT in the selection
    # invert the selection
    selectionNode.GetProperties().Set(vtk.vtkSelectionNode.INVERSE(), 1)
    extractSelection.Update()

    notSelected = vtk.vtkUnstructuredGrid()
    notSelected.ShallowCopy(extractSelection.GetOutput())

    print("There are %s points NOT in the selection" % notSelected.GetNumberOfPoints())
    print("There are %s cells NOT in the selection" % notSelected.GetNumberOfCells())

    inputMapper = vtk.vtkDataSetMapper()
    inputMapper.SetInputConnection(pointSource.GetOutputPort())
    inputActor = vtk.vtkActor()
    inputActor.SetMapper(inputMapper)

    selectedMapper = vtk.vtkDataSetMapper()
    if vtk.VTK_MAJOR_VERSION <= 5:
        selectedMapper.SetInputConnection(selected.GetProducerPort())
    else:
        selectedMapper.SetInputData(selected)
    selectedActor = vtk.vtkActor()
    selectedActor.SetMapper(selectedMapper)

    notSelectedMapper = vtk.vtkDataSetMapper()
    if vtk.VTK_MAJOR_VERSION <= 5:
        notSelectedMapper.SetInputConnection(notSelected.GetProducerPort())
    else:
        notSelectedMapper.SetInputData(notSelected)
    notSelectedActor = vtk.vtkActor()
    notSelectedActor.SetMapper(notSelectedMapper)


    # There will be one render window
    renderWindow = vtk.vtkRenderWindow()
    renderWindow.SetSize(900, 300)

    # And one interactor
    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(renderWindow)

    # Define viewport ranges
    # (xmin, ymin, xmax, ymax)
    leftViewport = [0.0, 0.0, 0.33, 1.0]
    centerViewport = [0.33, 0.0, .66, 1.0]
#.........这里部分代码省略.........

def convert2VPD(M,clean=False):
    """Convert pyFormex data to vtkPolyData.
    
    Convert a pyFormex Mesh or Coords object into vtkPolyData.
    This is limited to vertices, lines, and polygons.
    Lines should already be ordered (with connectedLineElems for instance).

    Parameters:

    - `M`: a Mesh or Coords type. If M is a Coords type it will be saved as
      VERTS. Else...
    - `clean`: if True, the resulting vtkdata will be cleaned by calling
      cleanVPD.

    Returns a vtkPolyData.
    """
    from vtk import vtkPolyData,vtkPoints,vtkIdTypeArray,vtkCellArray
    
    print('STARTING CONVERSION FOR DATA OF TYPE %s '%type(M))
    
    if type(M) == Coords:
        M = Mesh(M,arange(M.ncoords()))
    
    Nelems = M.nelems() # Number of elements
    Ncxel = M.nplex() # # Number of nodes per element
    
    # create a vtkPolyData variable
    vpd=vtkPolyData()
    
    # creating  vtk coords
    pts = vtkPoints()
    ntype=gnat(pts.GetDataType())
    coordsv = n2v(asarray(M.coords,order='C',dtype=ntype),deep=1) #.copy() # deepcopy array conversion for C like array of vtk, it is necessary to avoid memry data loss
    pts.SetNumberOfPoints(M.ncoords())
    pts.SetData(coordsv)
    vpd.SetPoints(pts)
    
    
    # create vtk connectivity
    elms = vtkIdTypeArray()
    ntype=gnat(vtkIdTypeArray().GetDataType())
    elmsv = concatenate([Ncxel*ones(Nelems).reshape(-1,1),M.elems],axis=1)
    elmsv = n2v(asarray(elmsv,order='C',dtype=ntype),deep=1) #.copy() # deepcopy array conversion for C like array of vtk, it is necessary to avoid memry data loss
    elms.DeepCopy(elmsv)

    # set vtk Cell data
    datav = vtkCellArray()
    datav.SetCells(Nelems,elms)
    if Ncxel == 1:
        try:
            print("setting VERTS for data with %s maximum number of point for cell "%Ncxel)
            vpd.SetVerts(datav)
        except:
            raise ValueError,"Error in saving  VERTS"

    elif Ncxel == 2:
        try:
            print ("setting LINES for data with %s maximum number of point for cell "%Ncxel)
            vpd.SetLines(datav)
        except:
            raise  ValueError,"Error in saving  LINES"
            
    else:
        try:
            print ("setting POLYS for data with %s maximum number of point for cell "%Ncxel)
            vpd.SetPolys(datav)
        except:
            raise ValueError,"Error in saving  POLYS"
            
    vpd.Update()
    if clean:
        vpd=cleanVPD(vpd)
    return vpd

def writeHdf5():
    
    if timeStep < 0.01:
        exit("Timestep is too small, choose 0.01 or larger")
    
    # This is where the data is for testing purposes.
    print "Current working directory:", os.getcwd()
    print taskMeshIn
    
    numQuadsPerRing = circQuads

    # Working with the task mesh junt to figure out the quads and rows numbers.
    taskMeshReader = vtk.vtkXMLPolyDataReader()
    taskMeshReader.SetFileName(taskMeshIn)
    taskMeshReader.Update()

    taskMesh = taskMeshReader.GetOutput()
    print taskMesh.GetNumberOfPoints()

    
    # Get the range of branch labels.
    labelRange = [0, 0]
    taskMesh.GetCellData().GetScalars().GetRange(labelRange, 0)

    # Convert label range to a list of labels.
    labelRange = range(int(labelRange[0]), int(labelRange[1]) + 1)
    print "Labels found in task mesh:", labelRange

    # Store the number of rings for each label. 
    numRingsPerLabel = {}   

    # For every label in the range of labels we want to extract all cells/quads.
    for label in labelRange:
        
        # Use this filter to extract the cells for a given label value.
        branchSelector = vtk.vtkThreshold()
        branchSelector.SetInputData(taskMesh)
        branchSelector.ThresholdBetween(label,label);
        branchSelector.Update()

        taskMeshBranch = branchSelector.GetOutput()

        numQuadRowsPerBranch = taskMeshBranch.GetNumberOfCells() / numQuadsPerRing;
        numRingsPerLabel[label] = numQuadRowsPerBranch
        
        
    atpFiles = glob.glob(atpMeshPattern)
    
    parentFile = h5py.File(atpHdf5Files[0], 'w')
    leftBranchFile = h5py.File(atpHdf5Files[1], 'w')
    rightBranchFile = h5py.File(atpHdf5Files[2], 'w')
    
    for atpIndex in range(0, len(atpFiles), int(timeStep / originalTimeStep)):
        print "Time step " + str(atpIndex * timeStep )
        

        atpMeshReader = vtk.vtkXMLPolyDataReader()
        atpMeshReader.SetFileName(atpFiles[atpIndex])
        atpMeshReader.Update()
    
        atpMesh = atpMeshReader.GetOutput()

        # For every label in the range of labels we want to extract all ECs.
        for label in labelRange:
            
    
            # Keep track of how many branches we need to skip.
            numECsPerLabel = numQuadsPerRing * numRingsPerLabel[label] * numECsPerQuad
            atpCellOffset = label * numECsPerLabel
    
    
            # Collect cell ids to select.
            selectionIds = vtk.vtkIdTypeArray()
            for sId in range(0, numECsPerLabel):
                selectionIds.InsertNextValue(atpCellOffset + sId)
    
            # Create selecion node.
            selectionNode = vtk.vtkSelectionNode()
            selectionNode.SetFieldType(selectionNode.CELL)
            selectionNode.SetContentType(selectionNode.INDICES)
            selectionNode.SetSelectionList(selectionIds)
    
            # Create selection.
            selection = vtk.vtkSelection()
            selection.AddNode(selectionNode)
    
            # Use vtkSelection filter.
            selectionExtractor = vtk.vtkExtractSelection()
            selectionExtractor.SetInputData(0, atpMesh)
            selectionExtractor.SetInputData(1, selection)
            selectionExtractor.Update()
    
            extractedCells = selectionExtractor.GetOutput()
    
            # Ring ids list for traversal.
            ringIds = range(0, numRingsPerLabel[label])
            ringIds.reverse()
    
            # Number of ECs rows is the number of ECs per quad.
            rowIds = range(0, numECsPerCol)
#.........这里部分代码省略.........

def _generate_vtk_mesh(points, cells):
    import vtk
    from vtk.util import numpy_support

    mesh = vtk.vtkUnstructuredGrid()

    # set points
    vtk_points = vtk.vtkPoints()
    # Not using a deep copy here results in a segfault.
    vtk_array = numpy_support.numpy_to_vtk(points, deep=True)
    vtk_points.SetData(vtk_array)
    mesh.SetPoints(vtk_points)

    # Set cells.
    meshio_to_vtk_type = {y: x for x, y in vtk_to_meshio_type.items()}

    # create cell_array. It's a one-dimensional vector with
    # (num_points2, p0, p1, ... ,pk, numpoints1, p10, p11, ..., p1k, ...
    cell_types = []
    cell_offsets = []
    cell_connectivity = []
    len_array = 0
    for meshio_type, data in cells.items():
        numcells, num_local_nodes = data.shape
        if meshio_type[:7] == "polygon":
            vtk_type = meshio_to_vtk_type[meshio_type[:7]]
        else:
            vtk_type = meshio_to_vtk_type[meshio_type]
        # add cell types
        cell_types.append(numpy.empty(numcells, dtype=numpy.ubyte))
        cell_types[-1].fill(vtk_type)
        # add cell offsets
        cell_offsets.append(
            numpy.arange(
                len_array,
                len_array + numcells * (num_local_nodes + 1),
                num_local_nodes + 1,
                dtype=numpy.int64,
            )
        )
        cell_connectivity.append(
            numpy.c_[
                num_local_nodes * numpy.ones(numcells, dtype=data.dtype), data
            ].flatten()
        )
        len_array += len(cell_connectivity[-1])

    cell_types = numpy.concatenate(cell_types)
    cell_offsets = numpy.concatenate(cell_offsets)
    cell_connectivity = numpy.concatenate(cell_connectivity)

    connectivity = vtk.util.numpy_support.numpy_to_vtkIdTypeArray(
        cell_connectivity.astype(numpy.int64), deep=1
    )

    # wrap the data into a vtkCellArray
    cell_array = vtk.vtkCellArray()
    cell_array.SetCells(len(cell_types), connectivity)

    # Add cell data to the mesh
    mesh.SetCells(
        numpy_support.numpy_to_vtk(
            cell_types, deep=1, array_type=vtk.vtkUnsignedCharArray().GetDataType()
        ),
        numpy_support.numpy_to_vtk(
            cell_offsets, deep=1, array_type=vtk.vtkIdTypeArray().GetDataType()
        ),
        cell_array,
    )

    return mesh