def Create_Topo(Slope,Plane):
	ugrid = vtk.vtkUnstructuredGrid()
	gridreader=vtk.vtkUnstructuredGridReader()
	gridreader.SetFileName(Slope)
	gridreader.Update()
	ugrid = gridreader.GetOutput()
	GeomFilt1 = vtk.vtkGeometryFilter()
	GeomFilt1.SetInput(ugrid)
	GeomFilt1.Update()
	x = GeomFilt1.GetOutput()

	u = vtk.vtkUnstructuredGrid()
	bgridreader=vtk.vtkXMLUnstructuredGridReader()
	bgridreader.SetFileName(Plane)
	bgridreader.Update()
	u = bgridreader.GetOutput() 
	GeomFilt2 = vtk.vtkGeometryFilter()
	GeomFilt2.SetInput(u)
	GeomFilt2.Update()
	y = GeomFilt2.GetOutput()

	append = vtk.vtkAppendPolyData()
	append.AddInput(x)
	append.AddInput(y)
	data = append.GetOutput()

	d = vtk.vtkDelaunay3D()
	d.SetInput(data)
	d.Update
	d.BoundingTriangulationOff()  
	d.SetTolerance(0.00001)
	d.SetAlpha(0)
	d.Update()
	z = d.GetOutput()
	return z

    def removeOldGeometry(self, fileName):
        if fileName is None:
            self.grid = vtk.vtkUnstructuredGrid()
            self.gridResult = vtk.vtkFloatArray()
            #self.emptyResult = vtk.vtkFloatArray()
            #self.vectorResult = vtk.vtkFloatArray()
            self.grid2 = vtk.vtkUnstructuredGrid()
            self.scalarBar.VisibilityOff()
            skipReading = True
        else:
            self.TurnTextOff()
            self.grid.Reset()
            self.grid2.Reset()
            self.gridResult = vtk.vtkFloatArray()
            self.gridResult.Reset()
            self.gridResult.Modified()
            self.eidMap = {}
            self.nidMap = {}

            self.resultCases = {}
            self.nCases = 0
            try:
                del self.caseKeys
                del self.iCase
                del self.iSubcaseNameMap
            except:
                print "cant delete geo"
                pass
            ###
            #print dir(self)
            skipReading = False
        self.scalarBar.VisibilityOff()
        self.scalarBar.Modified()
        return skipReading

    def test_contours(self):
        cell = vtk.vtkUnstructuredGrid()
        cell.ShallowCopy(self.Cell)

        np = self.Cell.GetNumberOfPoints()
        ncomb = pow(2, np)

        scalar = vtk.vtkDoubleArray()
        scalar.SetName("scalar")
        scalar.SetNumberOfTuples(np)
        cell.GetPointData().SetScalars(scalar)

        incorrectCases = []
        for i in range(1,ncomb-1):
            c = Combination(np, i)
            for p in range(np):
                scalar.SetTuple1(p, c[p])

            gradientFilter = vtk.vtkGradientFilter()
            gradientFilter.SetInputData(cell)
            gradientFilter.SetInputArrayToProcess(0,0,0,0,'scalar')
            gradientFilter.SetResultArrayName('grad')
            gradientFilter.Update()

            contourFilter = vtk.vtkContourFilter()
            contourFilter.SetInputConnection(gradientFilter.GetOutputPort())
            contourFilter.SetNumberOfContours(1)
            contourFilter.SetValue(0, 0.5)
            contourFilter.Update()

            normalsFilter = vtk.vtkPolyDataNormals()
            normalsFilter.SetInputConnection(contourFilter.GetOutputPort())
            normalsFilter.SetConsistency(0)
            normalsFilter.SetFlipNormals(0)
            normalsFilter.SetSplitting(0)

            calcFilter = vtk.vtkArrayCalculator()
            calcFilter.SetInputConnection(normalsFilter.GetOutputPort())
            calcFilter.SetAttributeTypeToPointData()
            calcFilter.AddVectorArrayName('grad')
            calcFilter.AddVectorArrayName('Normals')
            calcFilter.SetResultArrayName('dir')
            calcFilter.SetFunction('grad.Normals')
            calcFilter.Update()

            out = vtk.vtkUnstructuredGrid()
            out.ShallowCopy(calcFilter.GetOutput())

            numPts = out.GetNumberOfPoints()
            if numPts > 0:
                dirArray = out.GetPointData().GetArray('dir')
                for p in range(numPts):
                    if(dirArray.GetTuple1(p) > 0.0): # all normals are reversed
                        incorrectCases.append(i)
                        break

        self.assertEquals(','.join([str(i) for i in incorrectCases]), '')

    def __init__(self, *args, **kwargs):

        self.grid = vtk.vtkUnstructuredGrid()
        #gridResult = vtk.vtkFloatArray()
        #self.emptyResult = vtk.vtkFloatArray()
        #self.vectorResult = vtk.vtkFloatArray()
        self.grid2 = vtk.vtkUnstructuredGrid()

        # edges
        self.edgeActor = vtk.vtkActor()
        self.edgeMapper = vtk.vtkPolyDataMapper()

        self.is_edges = kwargs['is_edges']
        self.is_nodal = kwargs['is_nodal']
        self.is_centroidal = kwargs['is_centroidal']
        self.magnify = kwargs['magnify']
        self.debug = True

        gui_parent = kwargs['gui_parent']
        if 'log' in kwargs:
            self.log = kwargs['log']
            del kwargs['log']
        del kwargs['gui_parent']
        del kwargs['is_edges']
        del kwargs['is_nodal']
        del kwargs['is_centroidal']
        del kwargs['magnify']
        #del kwargs['rotation']
        wx.Panel.__init__(self, *args, **kwargs)
        NastranIO.__init__(self)
        Cart3dIO.__init__(self)
        LaWGS_IO.__init__(self)
        PanairIO.__init__(self)
        STL_IO.__init__(self)
        TetgenIO.__init__(self)
        Usm3dIO.__init__(self)


        self.nCases = 0

        #print("is_edges = %r" % self.is_edges)
        self.widget = pyWidget(self, -1)

        window = self.widget.GetRenderWindow()
        self.iren = vtk.vtkRenderWindowInteractor()

        if platform.system == 'Windows':
            self.iren.SetRenderWindow(window)

        self.iText = 0
        self.textActors = {}

def Creating_z(Plane,Slope,Sediment,Start_time,End_time,Time_step):

	vtuObject = vtktools.vtu(Plane)
	vtuObject.GetFieldNames()
	gradient = vtuObject.GetScalarField('u')
	ugrid = vtk.vtkUnstructuredGrid()
	gridreader=vtk.vtkXMLUnstructuredGridReader()
	gridreader.SetFileName(Plane)
	gridreader.Update()
	ugrid = gridreader.GetOutput()
	points = ugrid.GetPoints()

	nPoints = ugrid.GetNumberOfPoints()
	for p in range(0,nPoints):
		x = (points.GetPoint(p)[:2] + (gradient[p],))
		points.SetPoint(p,x)
    
	ugrid.Update()
###################################################################################################################
	t = Start_time
	dt = Time_step
	et = End_time
	while t <= et:

		Import = Sediment + str(t) +'000000.vtu'
		NewSave = Sediment + str(t) + '_sed_slope.pvd'
		vtuObjectSed = vtktools.vtu(Import)
		vtuObjectSed.GetFieldNames()
		gradientSed = vtuObjectSed.GetScalarField('u')
		sedgrid = vtk.vtkUnstructuredGrid()
		sedgridreader=vtk.vtkXMLUnstructuredGridReader()
		sedgridreader.SetFileName(Import)
		sedgridreader.Update()
		sedgrid = sedgridreader.GetOutput()
		s = sedgrid.GetPoints()
	
		for p in range(0,nPoints):
			x = ((s.GetPoint(p)[0],) + (s.GetPoint(p)[1],) + ((gradientSed[p]+gradient[p]),))
			s.SetPoint(p,x)

		writer = vtk.vtkUnstructuredGridWriter()
		writer.SetFileName(NewSave)
		writer.SetInput(sedgrid)
		writer.Update()
		writer.Write()
		t += dt
	writer = vtk.vtkUnstructuredGridWriter()
	writer.SetFileName(Slope)
	writer.SetInput(ugrid)
	writer.Update()
	writer.Write()

    def convert_to_vtk(self, pvar_list):
        """
        Convert particle data into vtk format and return the vtk objects.

        call signature:

        convert_to_vtk(self, pvar_list)

        Keyword arguments:

        *pvar_list*:
          List of particle objects.
        """

        import numpy as np
        import vtk
        
        for pvar in pvar_list:
            points = np.vstack([pvar.xp, pvar.yp, pvar.zp])
            
            vtk_grid_data = vtk.vtkUnstructuredGrid()
            vtk_points = vtk.vtkPoints()
            
            # Add the data to the vtk points.
            for point_idx in range(points.shape[1]):
                vtk_points.InsertNextPoint(points[:, point_idx])
            # Add the vtk points to the vtk grid.
            vtk_grid_data.SetPoints(vtk_points)
            
            self.vtk_grid_data.append(vtk_grid_data)

 def pca(self):
   """Performs Principle Component Analysis on the alignedGrids. Also calculates the mean shape.
   """
   
   logging.info("running pca")
   
   size = len(self.alignedGrids)
   
   self.filterPCA.SetNumberOfInputs(size)
   
   for id, grid in enumerate(self.alignedGrids):    
     self.filterPCA.SetInput(id, grid)
   
   self.filterPCA.Update()
      
   #Get the eigenvalues
   evals = self.filterPCA.GetEvals()
   
   #Now let's get mean ^^
   b = vtk.vtkFloatArray()
   b.SetNumberOfComponents(0)
   b.SetNumberOfTuples(0)
   mean = vtk.vtkUnstructuredGrid()
   mean.DeepCopy(self.alignedGrids[0])
   #Get the mean shape:
   self.filterPCA.GetParameterisedShape(b, mean)
   self.meanShape.append(mean)   
   
   #get the meanpositions
   for pos in range(self.meanShape[0].GetNumberOfCells()):
     bounds = self.meanShape[0].GetCell(pos).GetBounds()
     self.meanPositions.append((bounds[0],bounds[2]))
     
   logging.info("done")

    def doNonLinear(self, ghosts, ncells):
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(10)
        pts.InsertPoint(0, (0, 0, 0))
        pts.InsertPoint(1, (1, 0, 0))
        pts.InsertPoint(2, (0.5, 1, 0))
        pts.InsertPoint(3, (0.5, 0.5, 1))
        pts.InsertPoint(4, (0.5, 0, 0))
        pts.InsertPoint(5, (1.25, 0.5, 0))
        pts.InsertPoint(6, (0.25, 0.5, 0))
        pts.InsertPoint(7, (0.25, 0.25, 0.5))
        pts.InsertPoint(8, (0.75, 0.25, 0.5))
        pts.InsertPoint(9, (0.5, 0.75, 0.5))

        te = vtk.vtkQuadraticTetra()
        ptIds = te.GetPointIds()
        for i in range(10):
            ptIds.SetId(i, i)

        grid = vtk.vtkUnstructuredGrid()
        grid.Allocate(1, 1)
        grid.InsertNextCell(te.GetCellType(), te.GetPointIds())
        grid.SetPoints(pts)
        grid.GetPointData().AddArray(ghosts)

        ugg = vtk.vtkUnstructuredGridGeometryFilter()
        ugg.SetInputData(grid)

        dss = vtk.vtkDataSetSurfaceFilter()
        dss.SetNonlinearSubdivisionLevel(2)
        dss.SetInputConnection(ugg.GetOutputPort())
        dss.Update()
        self.assertEqual(dss.GetOutput().GetNumberOfCells(), ncells)

    def asVtkUnstructuredGrid(self):
        '''
        Return object as a vtk.vtkUnstructuredMesh instance.

        .. note:: This method uses the compiled vtk module (which is a wrapper atop the c++ VTK library) -- in contrast to :obj:`UnstructuredMesh.getVTKRepresentation`, which uses the pyvtk module (python-only implementation of VTK i/o supporting only VTK File Format version 2).
        '''
        import vtk
        # vertices
        pts=vtk.vtkPoints()
        for ip in range(self.getNumberOfVertices()): pts.InsertNextPoint(self.getVertex(ip).getCoordinates())
        # cells
        cells,cellTypes=vtk.vtkCellArray(),[]
        for ic in range(self.getNumberOfCells()):
            c=self.getCell(ic)
            cgt=c.getGeometryType()
            cellGeomTypeMap={
                CellGeometryType.CGT_TRIANGLE_1: (vtk.vtkTriangle,vtk.VTK_TRIANGLE),
                CellGeometryType.CGT_QUAD:       (vtk.vtkQuad,vtk.VTK_QUAD),
                CellGeometryType.CGT_TETRA:      (vtk.vtkTetra,vtk.VTK_TETRA),
                CellGeometryType.CGT_HEXAHEDRON: (vtk.vtkHexahedron,vtk.VTK_HEXAHEDRON),
                CellGeometryType.CGT_TRIANGLE_2: (vtk.vtkQuadraticTriangle,vtk.VTK_QUADRATIC_TRIANGLE)
            }
            c2klass,c2type=cellGeomTypeMap[cgt] # instantiate the VTK cell with the correct type
            c2=c2klass()
            verts=c.getVertices() # those should be all instances of Vertex...? Hopefully so.
            for i,v in enumerate(verts): c2.GetPointIds().SetId(i,v.getNumber())
            cells.InsertNextCell(c2)
            cellTypes.append(c2type)
        ret=vtk.vtkUnstructuredGrid()
        ret.SetPoints(pts)
        ret.SetCells(cellTypes,cells)
        return ret

    def __init__(self, pos, cell):
        """Construct basic VTK-representation of a set of atomic positions.

        pos: NumPy array of dtype float and shape ``(n,3)``
            Cartesian positions of the atoms.
        cell: Instance of vtkUnitCellModule of subclass thereof
            Holds information equivalent to that of atoms.get_cell().

        """
        # Make sure position argument is a valid array
        if not isinstance(pos, np.ndarray):
            pos = np.array(pos)

        assert pos.dtype == float and pos.shape[1:] == (3,)

        vtkBaseGrid.__init__(self, len(pos), cell)

        # Convert positions to VTK array
        npy2da = vtkDoubleArrayFromNumPyArray(pos)
        vtk_pda = npy2da.get_output()
        del npy2da

        # Transfer atomic positions to VTK points
        self.vtk_pts = vtkPoints()
        self.vtk_pts.SetData(vtk_pda)

        # Create a VTK unstructured grid of these points
        self.vtk_ugd = vtkUnstructuredGrid()
        self.vtk_ugd.SetWholeBoundingBox(self.cell.get_bounding_box())
        self.vtk_ugd.SetPoints(self.vtk_pts)

        # Extract the VTK point data set
        self.set_point_data(self.vtk_ugd.GetPointData())

def MakeHexagonalPrism():
    '''
      3D: hexagonal prism: a wedge with an hexagonal base.
      Be careful, the base face ordering is different from wedge.
    '''
 
    numberOfVertices = 12
 
    points = vtk.vtkPoints()
 
    points.InsertNextPoint(0.0, 0.0, 1.0)
    points.InsertNextPoint(1.0, 0.0, 1.0)
    points.InsertNextPoint(1.5, 0.5, 1.0)
    points.InsertNextPoint(1.0, 1.0, 1.0)
    points.InsertNextPoint(0.0, 1.0, 1.0)
    points.InsertNextPoint(-0.5, 0.5, 1.0)
 
    points.InsertNextPoint(0.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 0.0, 0.0)
    points.InsertNextPoint(1.5, 0.5, 0.0)
    points.InsertNextPoint(1.0, 1.0, 0.0)
    points.InsertNextPoint(0.0, 1.0, 0.0)
    points.InsertNextPoint(-0.5, 0.5, 0.0)
 
    hexagonalPrism = vtk.vtkHexagonalPrism()
    for i in range(0, numberOfVertices):
        hexagonalPrism.GetPointIds().SetId(i, i)
 
    ug = vtk.vtkUnstructuredGrid()
    ug.InsertNextCell(hexagonalPrism.GetCellType(),
                       hexagonalPrism.GetPointIds())
    ug.SetPoints(points)
 
    return ug

 def __init___(self, scene, V, opacity=1, color='gray'):
     self.scene = scene
     self.frame = scene.frame
     self.V = V
     n_voxels = len(V)
     size = V.bin_size
     pts = vtk.vtkPoints()
     pts.SetNumberOfPoints(8 * n_voxels)
     grid = vtk.vtkUnstructuredGrid()
     grid.Allocate(n_voxels, 1)
     vx = vtk.vtkVoxel()
     for i, q in enumerate(V):
         pos = q * size + V.low_range
         pts.InsertPoint(i * 8 + 0, *pos)
         pts.InsertPoint(i * 8 + 1, *(pos + (size,0,0)))
         pts.InsertPoint(i * 8 + 2, *(pos + (0,size,0)))
         pts.InsertPoint(i * 8 + 3, *(pos + (size,size,0)))
         pts.InsertPoint(i * 8 + 4, *(pos + (0,0,size)))
         pts.InsertPoint(i * 8 + 5, *(pos + (size,0,size)))
         pts.InsertPoint(i * 8 + 6, *(pos + (0,size,size)))
         pts.InsertPoint(i * 8 + 7, *(pos + (size,size,size)))
         for j in range(8):
             vx.GetPointIds().SetId(j, i * 8 + j)
         grid.InsertNextCell(vx.GetCellType(), vx.GetPointIds())
     grid.SetPoints(pts)
     mapper = vtk.vtkDataSetMapper()
     mapper.SetInput(grid)
     self.actor = vtk.vtkActor()
     self.actor.SetMapper(mapper)
     #self.actor.GetProperty().SetDiffuseColor(*name_to_rgb_float(color))
     self.actor.GetProperty().SetColor(*name_to_rgb_float(color))
     self.actor.GetProperty().SetOpacity(opacity)
     self.frame.ren.AddActor(self.actor)

def create_vessel_actor(ref_data):
    vessel_ref_data = ref_data
    points = vtk.vtkPoints()

    # insert each properties of points into obj.
    for i in range(vessel_ref_data.get_len_of_vassel_value()):
        x = vessel_ref_data.get_abscissa_value_at(i)
        y = vessel_ref_data.get_ordinate_value_at(i)
        z = vessel_ref_data.get_iso_value_at(i)
        points.InsertPoint(i, x, y, z)

    poly = vtk.vtkPolyVertex()
    poly.GetPointIds().SetNumberOfIds(vessel_ref_data.get_len_of_vassel_value())

    cont = 0
    while cont < vessel_ref_data.get_len_of_vassel_value():
        poly.GetPointIds().SetId(cont, cont)
        cont += 1

    grid = vtk.vtkUnstructuredGrid()
    grid.SetPoints(points)
    grid.InsertNextCell(poly.GetCellType(), poly.GetPointIds())

    mapper = vtk.vtkDataSetMapper()
    if sys.hexversion == 34015984:
        mapper.SetInputData(grid)
    elif sys.hexversion == 34015728:
        mapper.SetInput(grid)

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    return actor

def makeTetraGrid(nCubes, flip=0):
    max_x=nCubes+1
    max_y=nCubes+1
    max_z=nCubes+1
    scale=20./nCubes #*(19./20)
    meshPoints = vtk.vtkPoints()
    meshPoints.SetNumberOfPoints(max_x*max_y*max_z)
    #i*(max_y)*(max_z)+j*(max_z)+k
    for i in xrange(max_x):
        for j in xrange(max_y):
            for k in xrange(max_z):
                meshPoints.InsertPoint(i*(max_y)*(max_z)+j*(max_z)+k,scale*i,scale*j,scale*k)

    nelements = 5*(max_x-1)*(max_y-1)*(max_z-1)
    meshGrid = vtk.vtkUnstructuredGrid()
    meshGrid.Allocate(nelements, nelements)
    meshGrid.SetPoints(meshPoints)

    for i in range(max_x-1):                  
        for j in range(max_y-1):              
            for k in range(max_z-1):
                ulf = (i+1)*(max_y)*(max_z)+j*(max_z)+k        # upper left front
                urf = (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k    # upper right front
                lrf = i*(max_y)*(max_z)+(j+1)*(max_z)+k        # lower right front
                llf = i*(max_y)*(max_z)+j*(max_z)+k            # lower left front 
                ulb = (i+1)*(max_y)*(max_z)+j*(max_z)+k+1      # upper left back
                urb = (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k+1  # upper right back 
                lrb = i*(max_y)*(max_z)+(j+1)*(max_z)+k+1      # lower right back
                llb = i*(max_y)*(max_z)+j*(max_z)+k+1          # lower left back
                
                point_order = [  # not flip
                    [[llf,urf,lrf,lrb],
                     [llf,ulf,urf,ulb],
                     [lrb,urf,urb,ulb],
                     [llf,lrb,llb,ulb],
                     [llf,ulb,urf,lrb]],
                    # flip
                    [[llf,ulf,lrf,llb],
                     [ulf,urf,lrf,urb],
                     [ulf,ulb,urb,llb],
                     [lrf,urb,lrb,llb],
                     [ulf,urb,lrf,llb]
                     ]]
                
                for o in point_order[flip]:
                    cell = vtk.vtkTetra()
                    id=0
                    for p in o:
                        cell.GetPointIds().SetId(id,p)
                        id+=1
                    meshGrid.InsertNextCell(cell.GetCellType(),cell.GetPointIds())

                flip = not flip

            if (max_z-1)%2==0:
                flip = not flip
        if (max_y-1)%2==0:
            flip = not flip

    return meshGrid

	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 makeVoxelGrid(nCubes):
    max_x=nCubes+1
    max_y=nCubes+1
    max_z=nCubes+1
    scale=20./nCubes #*(19./20)
    meshPoints = vtk.vtkPoints()
    meshPoints.SetNumberOfPoints(max_x*max_y*max_z)
    #i*(max_y)*(max_z)+j*(max_z)+k
    for i in xrange(max_x):
        for j in xrange(max_y):
            for k in xrange(max_z):
                meshPoints.InsertPoint(i*(max_y)*(max_z)+j*(max_z)+k,scale*i,scale*j,scale*k)

    nelements = (max_x-1)*(max_y-1)*(max_z-1)
    meshGrid = vtk.vtkUnstructuredGrid()
    meshGrid.Allocate(nelements, nelements)
    meshGrid.SetPoints(meshPoints)

    for i in range(max_x-1):                  
        for j in range(max_y-1):              
            for k in range(max_z-1):          
                cell = vtk.vtkHexahedron()
                # This is the order we need such that the normals point out.
                cell.GetPointIds().SetId(0, (i+1)*(max_y)*(max_z)+j*(max_z)+k)       # upper left front
                cell.GetPointIds().SetId(1, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k)   # upper right front
                cell.GetPointIds().SetId(2, i*(max_y)*(max_z)+(j+1)*(max_z)+k)       # lower right front
                cell.GetPointIds().SetId(3, i*(max_y)*(max_z)+j*(max_z)+k)           # lower left front node index
                cell.GetPointIds().SetId(4, (i+1)*(max_y)*(max_z)+j*(max_z)+k+1)     # upper left back
                cell.GetPointIds().SetId(5, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k+1) # upper right back
                cell.GetPointIds().SetId(6, i*(max_y)*(max_z)+(j+1)*(max_z)+k+1)     # lower right back
                cell.GetPointIds().SetId(7, i*(max_y)*(max_z)+j*(max_z)+k+1)         # lower left back
                meshGrid.InsertNextCell(cell.GetCellType(), cell.GetPointIds())

    return meshGrid

def create_vessel_actor(ref_data):
    # vessel_ref_data = ref_data

    points = vtk.vtkPoints()

    # insert each properties of points into obj.
    for i in xrange(len(ref_data)):
        x = ref_data[i][0]
        y = ref_data[i][1]
        z = ref_data[i][2]
        points.InsertPoint(i, x, y, z)

    poly = vtk.vtkPolyVertex()
    poly.GetPointIds().SetNumberOfIds(len(ref_data))

    cont = 0
    while cont < len(ref_data):
        poly.GetPointIds().SetId(cont, cont)
        cont += 1

    grid = vtk.vtkUnstructuredGrid()
    grid.SetPoints(points)
    grid.InsertNextCell(poly.GetCellType(), poly.GetPointIds())

    mapper = vtk.vtkDataSetMapper()
    if sys.hexversion == 34015984:
        mapper.SetInputData(grid)
    if sys.hexversion == 34015728:
        mapper.SetInput(grid)

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    return actor

def show_beta_sheets(renderer):
    for sheet in sheet_defs:
#    helix = helix_defs[1]
        aPolyLineGrid = vtk.vtkUnstructuredGrid()
        aPolyLineGrid.Allocate(5, 1)

        (polyLinePoints, aPolyLine) = make_polyline(sheet)

#         # Create a tube filter to represent the lines as tubes.  Set up the
#         # associated mapper and actor.
#         tuber = vtk.vtkTubeFilter()
#         tuber.SetInputConnection(appendF.GetOutputPort())
#         tuber.SetRadius(0.1)
#         lineMapper = vtk.vtkPolyDataMapper()
#         lineMapper.SetInputConnection(tuber.GetOutputPort())
#         lineActor = vtk.vtkActor()
#         lineActor.SetMapper(lineMapper)

        aPolyLineGrid.InsertNextCell(aPolyLine.GetCellType(),
                                     aPolyLine.GetPointIds())
        aPolyLineGrid.SetPoints(polyLinePoints)

        aPolyLineMapper = vtk.vtkDataSetMapper()
        aPolyLineMapper.SetInput(aPolyLineGrid)
        aPolyLineActor = vtk.vtkActor()
        aPolyLineActor.SetMapper(aPolyLineMapper)
        aPolyLineActor.GetProperty().SetDiffuseColor(1, 1, 1)

        renderer.AddActor(aPolyLineActor)

  def addPointSet(self, data):
    """Add a set of points to the UnstructuredGrid vertexGrid (the source of our calculations)
    """
    
    logging.info("adding point set")
    
    size = len(data)

    # Create some landmarks, put them in UnstructuredGrid
    # Start off with some landmarks
    vertexPoints = vtk.vtkPoints()
    vertexPoints.SetNumberOfPoints(size)
    
    vertexGrid = vtk.vtkUnstructuredGrid()
    vertexGrid.Allocate(size, size)
    
    for id, (x, y, z) in enumerate(data):
      vertexPoints.InsertPoint(id, x, y, z)
      
      # Create vertices from them
      vertex = vtk.vtkVertex()
      vertex.GetPointIds().SetId(0, id)
      
      #Create an unstructured grid with the landmarks added
      vertexGrid.InsertNextCell(vertex.GetCellType(), vertex.GetPointIds())
    
    vertexGrid.SetPoints(vertexPoints)
    self.vertexGrids.append(vertexGrid)
        
    logging.info("done")

def MakeVoxel():
    '''
      A voxel is a representation of a regular grid in 3-D space.
    '''
    numberOfVertices = 8
 
    points = vtk.vtkPoints()
    points.InsertNextPoint(0, 0, 0)
    points.InsertNextPoint(1, 0, 0)
    points.InsertNextPoint(0, 1, 0)
    points.InsertNextPoint(1, 1, 0)
    points.InsertNextPoint(0, 0, 1)
    points.InsertNextPoint(1, 0, 1)
    points.InsertNextPoint(0, 1, 1)
    points.InsertNextPoint(1, 1, 1)
 
    voxel = vtk.vtkVoxel()
    for i in range(0, numberOfVertices):
        voxel.GetPointIds().SetId(i, i)
 
    ug = vtk.vtkUnstructuredGrid()
    ug.SetPoints(points)
    ug.InsertNextCell(voxel.GetCellType(), voxel.GetPointIds())
 
    return ug