def as_vtkarray(a):
    if isinstance(a, npy.ndarray):
        return vnp.numpy_to_vtk(a)
    elif isinstance(a, rpy2.rinterface.SexpVector):
        return vnp.numpy_to_vtk(npy.asarray(a))
    else:
        return None 

    def render(self, **kwargs):
        """ Plots the volume and the control points. """
        # Calling parent function
        super(VisVolume, self).render(**kwargs)

        # Initialize a list to store VTK actors
        vtk_actors = []

        # Start plotting
        for plot in self._plots:
            # Plot control points
            if plot['type'] == 'ctrlpts' and self.vconf.display_ctrlpts:
                # Points as spheres
                pts = np.array(plot['ptsarr'], dtype=np.float)
                vtkpts = numpy_to_vtk(pts, deep=False, array_type=VTK_FLOAT)
                vtkpts.SetName(plot['name'])
                temp_actor = vtkh.create_actor_pts(pts=vtkpts, color=vtkh.create_color(plot['color']),
                                                   name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

            # Plot evaluated points
            if plot['type'] == 'evalpts' and self.vconf.display_evalpts:
                pts = np.array(plot['ptsarr'], dtype=np.float)
                vtkpts = numpy_to_vtk(pts, deep=False, array_type=VTK_FLOAT)
                vtkpts.SetName(plot['name'])
                temp_actor = vtkh.create_actor_pts(pts=vtkpts, color=vtkh.create_color(plot['color']),
                                                   name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

        # Render actors
        return vtkh.create_render_window(vtk_actors, dict(KeyPressEvent=(self.vconf.keypress_callback, 1.0)),
                                         figure_size=self.vconf.figure_size)

    def setPointHeights( self, ptheights ):
        try:
            if self.topo == PlotType.Planar:   
                self.np_points_data[2::3] =  ptheights
                vtk_points_data = numpy_support.numpy_to_vtk( self.np_points_data ) 
                vtk_points_data.SetNumberOfComponents( 3 )
                vtk_points_data.SetNumberOfTuples( len( self.np_points_data ) / 3 ) 
                self.vtk_planar_points.SetData( vtk_points_data )
                self.polydata.SetPoints( self.vtk_planar_points )  
                self.vtk_planar_points.Modified()
            elif self.topo == PlotType.Spherical:
                self.np_sp_grid_data[0::3] =  self.spherical_scaling * ptheights + self.earth_radius
                vtk_sp_grid_data = numpy_support.numpy_to_vtk( self.np_sp_grid_data ) 
                size = vtk_sp_grid_data.GetSize()                    
                vtk_sp_grid_data.SetNumberOfComponents( 3 )
                vtk_sp_grid_data.SetNumberOfTuples( size/3 )   
                vtk_sp_grid_points = vtk.vtkPoints()
                vtk_sp_grid_points.SetData( vtk_sp_grid_data )
                self.vtk_spherical_points = vtk.vtkPoints()
                self.shperical_to_xyz_trans.TransformPoints( vtk_sp_grid_points, self.vtk_spherical_points ) 
#                pt0 = self.vtk_spherical_points.GetPoint(0)
    #            print "VTK Set point Heights, samples: %s %s %s " % ( str( ptheights[0] ), str( self.np_sp_grid_data[0] ), str( pt0 ) )
                self.polydata.SetPoints( self.vtk_spherical_points ) 
                self.vtk_spherical_points.Modified()
            self.polydata.Modified()
        except Exception, err:
            self.printLogMessage( "Processing point heights: %s " % str( err ), error=True )

def Gen_uGrid(new_pt, new_fc):
    """ Generates a vtk unstructured grid given points and triangular faces"""    
    ints = np.ones(len(new_fc), 'int')*3
    cells = np.hstack((ints.reshape(-1, 1), np.vstack(new_fc)))
              
    # Generate vtk mesh
    
    # Convert points to vtkfloat object
    points = np.vstack(new_pt)
    vtkArray = VN.numpy_to_vtk(np.ascontiguousarray(points), deep=True)#, deep=True) 
    points = vtk.vtkPoints()
    points.SetData(vtkArray) 
                
    # Convert to vtk arrays
    tritype = vtk.vtkTriangle().GetCellType()*np.ones(len(new_fc), 'int')
    cell_type = np.ascontiguousarray(tritype).astype('uint8')
    cell_type = VN.numpy_to_vtk(cell_type, deep=True)
    
    offset = np.cumsum(np.hstack(ints + 1))
    offset = np.ascontiguousarray(np.delete(np.insert(offset, 0, 0), -1)).astype('int64')  # shift
    offset = VN.numpy_to_vtkIdTypeArray(offset, deep=True)
    
    cells = np.ascontiguousarray(np.hstack(cells).astype('int64'))
    vtkcells = vtk.vtkCellArray()
    vtkcells.SetCells(cells.shape[0], VN.numpy_to_vtkIdTypeArray(cells, deep=True))
    
    
    # Create unstructured grid
    uGrid = vtk.vtkUnstructuredGrid()
    uGrid.SetPoints(points)
    uGrid.SetCells(cell_type, offset, vtkcells)
    
    return uGrid

def putMaskOnVTKGrid(data,grid,actorColor=None,deep=True):
  #Ok now looking 
  msk = data.mask
  imsk =  VN.numpy_to_vtk(msk.astype(numpy.int).flat,deep=deep)
  mapper = None
  if msk is not numpy.ma.nomask:
      msk =  VN.numpy_to_vtk(numpy.logical_not(msk).astype(numpy.uint8).flat,deep=deep)
      if actorColor is not None:
          grid2 = vtk.vtkStructuredGrid()
          grid2.CopyStructure(grid)
          geoFilter = vtk.vtkDataSetSurfaceFilter()
          if grid.IsA("vtkStructuredGrid"):
              grid2.GetPointData().SetScalars(imsk)
              #grid2.SetCellVisibilityArray(imsk)
              p2c = vtk.vtkPointDataToCellData()
              p2c.SetInputData(grid2)
              geoFilter.SetInputConnection(p2c.GetOutputPort())
          else:
              grid2.GetCellData().SetScalars(imsk)
              geoFilter.SetInputData(grid)
          geoFilter.Update()
          mapper = vtk.vtkPolyDataMapper()
          mapper.SetInputData(geoFilter.GetOutput())
          lut = vtk.vtkLookupTable()
          lut.SetNumberOfTableValues(1)
          r,g,b = actorColor
          lut.SetNumberOfTableValues(2)
          lut.SetTableValue(0,r/100.,g/100.,b/100.)
          lut.SetTableValue(1,r/100.,g/100.,b/100.)
          mapper.SetLookupTable(lut)
          mapper.SetScalarRange(1,1)
      if grid.IsA("vtkStructuredGrid"):
          grid.SetPointVisibilityArray(msk)
          #grid.SetCellVisibilityArray(msk)
  return mapper

    def makeVTK(self, filename):
        nelements = len(self.data[filename]['elements'])
        n = np.array(self.data[filename]['nodes'])
        arr = numpy_support.numpy_to_vtk(
            n.ravel(), deep=True, array_type=vtk.VTK_DOUBLE)
        arr.SetNumberOfComponents(3)
        tetraPoints = vtk.vtkPoints()
        tetraPoints.SetData(arr)

        vtkMesh = vtk.vtkUnstructuredGrid()
        vtkMesh.Allocate(nelements, nelements)
        vtkMesh.SetPoints(tetraPoints)

        e = np.array(self.data[filename]['elements'], np.uint32) - 1
        e = np.hstack((np.ones((e.shape[0], 1), np.uint32) * 4, e))
        arr = numpy_support.numpy_to_vtk(e.ravel(), deep=True,
                                         array_type=vtk.VTK_ID_TYPE)
        tet = vtk.vtkCellArray()
        tet.SetCells(old_div(e.size, 5), arr)
        vtkMesh.SetCells(10, tet)

        centroids = (n[e[:, 1]] + n[e[:, 2]] + n[e[:, 3]] + n[e[:, 4]])
        centroids /= 4.0
        arr = numpy_support.numpy_to_vtk(
            centroids.ravel(), deep=True, array_type=vtk.VTK_FLOAT)
        arr.SetName("Centroids")
        arr.SetNumberOfComponents(3)
        vtkMesh.GetCellData().AddArray(arr)
        self.vtkMeshes[filename] = vtkMesh

	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 writeVTK(mesh, fileName, models=None):
        """
        Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.

        Input:
        :param str, path to the output vtk file
        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells

        """
        # Import
        from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION
        from vtk.util.numpy_support import numpy_to_vtk

        # Deal with dimensionalities
        if mesh.dim >= 1:
            vX = mesh.vectorNx
            xD = mesh.nNx
            yD,zD = 1,1
            vY, vZ = np.array([0,0])
        if mesh.dim >= 2:
            vY = mesh.vectorNy
            yD = mesh.nNy
        if mesh.dim == 3:
            vZ = mesh.vectorNz
            zD = mesh.nNz
        # Use rectilinear VTK grid.
        # Assign the spatial information.
        vtkObj = rectGrid()
        vtkObj.SetDimensions(xD,yD,zD)
        vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
        vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
        vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))

        # Assign the model('s) to the object
        if models is not None:
            for item in models.iteritems():
                # Convert numpy array
                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
                vtkDoubleArr.SetName(item[0])
                vtkObj.GetCellData().AddArray(vtkDoubleArr)
            # Set the active scalar
            vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
        # vtkObj.Update()

        # Check the extension of the fileName
        ext = os.path.splitext(fileName)[1]
        if ext is '':
            fileName = fileName + '.vtr'
        elif ext not in '.vtr':
            raise IOError('{:s} is an incorrect extension, has to be .vtr')
        # Write the file.
        vtrWriteFilter = rectWriter()
        if float(VTK_VERSION.split('.')[0]) >=6:
            vtrWriteFilter.SetInputData(vtkObj)
        else:
            vtuWriteFilter.SetInput(vtuObj)
        vtrWriteFilter.SetFileName(fileName)
        vtrWriteFilter.Update()

	def keyPressEvent(self, event):
		if event.key() == QtCore.Qt.Key_L:
			self.basis_idx += 1
			coeffs = VN.numpy_to_vtk(self.WavBases[self.basis_idx][::-1,0]*100, deep=True)
			coeffs.SetName('coefficient')
			c_sign = VN.numpy_to_vtk(N.sign(self.WavBases[self.basis_idx][::-1,0]), deep=True)
			c_sign.SetName('sign')
			
			self.table.RemoveColumn(2)
			self.table.RemoveColumn(1)
			self.table.AddColumn(coeffs)
			self.table.AddColumn(c_sign)
			self.WordleView.RemoveAllRepresentations()
			self.WordleView.AddRepresentationFromInput(self.table)

			self.table.Modified()
			self.WordleView.Update()
		
		if event.key() == QtCore.Qt.Key_Space:
			if event.modifiers() == QtCore.Qt.NoModifier:
				self.WordleView.Modified()
				self.WordleView.Update()
			elif event.modifiers() == QtCore.Qt.ShiftModifier:
				self.table.Modified()
				self.WordleView.Update()

		# Write PNG (n)
		# Trying to use a integer-based QImage
		if event.key() == QtCore.Qt.Key_N:
			scene = self.WordleView.GetScene()
			image = QtGui.QImage(256,256,QtGui.QImage.Format_ARGB32)
			painter = QtGui.QPainter(image)
			painter.setRenderHint(QtGui.QPainter.Antialiasing)
			scene.render(painter)
			painter.end()
			image.save("out.png")
		
		# Write PDF (p)
		if event.key() == QtCore.Qt.Key_P:
			scene = self.WordleView.GetScene()
			printer = QtGui.QPrinter()
			printer.setOutputFormat(QtGui.QPrinter.PdfFormat)
			printer.setOutputFileName("out.pdf")
			pdfPainter = QtGui.QPainter(printer)
			scene.render(pdfPainter)
			pdfPainter.end()
		
		# Write SVG (s)
		if event.key() == QtCore.Qt.Key_S:
			scene = self.WordleView.GetScene()
			svggen = QtSvg.QSvgGenerator()
			svggen.setFileName("out.svg")
			svggen.setSize(QtCore.QSize(600, 600))
			svggen.setViewBox(QtCore.QRect(0, 0, 600, 600))
			svggen.setTitle("SVG Generator Example Drawing")
			svggen.setDescription("An SVG drawing created by the SVG Generator")
			svgPainter = QtGui.QPainter(svggen)
			scene.render(svgPainter)
			svgPainter.end()

    def execute(self):
        res = self.input("in")
    
        #get primary variable..
        #need to ensure that results that go back to VisIt
        #are at least correct size..
        varr = self.context.mesh.GetCellData().GetArray(self.context.primary_var)
        if varr is None:
            varr = self.context.mesh.GetPointData().GetArray(self.context.primary_var)

        if res is None:
            return self.context.mesh

        if isinstance(res,vtk.vtkDataSet):
            return res

        if isinstance(res,npy.ndarray):
            res = npy.ascontiguousarray(res)
            res = vnp.numpy_to_vtk(res)

        if not isinstance(res, vtk.vtkDataArray):
            if isinstance(res,npy.ndarray) or isinstance(res, (list,tuple)):
                np_tmp = npy.ascontiguousarray(res)
            else:
                np_tmp = npy.ascontiguousarray([res])

            #ensure 1 dimension before putting it in vtk..
            np_tmp = npy.ravel(np_tmp)
            #pad with zeros if incorrect size..
            if varr is not None and varr.GetDataSize() > len(np_tmp):
                np_tmp = npy.pad(np_tmp,(0,len(np_tmp)-var.GetDataSize()),'constant')
            res = vnp.numpy_to_vtk(np_tmp)

        #if isinstance(res,npy.ndarray):
        #    # create
        #    vdata = vtk.vtkFloatArray()
        #    vdata.SetNumberOfComponents(1)
        #    vdata.SetNumberOfTuples(res.shape[0])
        #    npo = vnp.vtk_to_numpy(vdata)
        #    npo[:] = res
        #    res = vdata
        if isinstance(res,vtk.vtkDataArray):
            res.SetName(self.context.primary_var)
            rset = self.context.mesh.NewInstance()
            rset.ShallowCopy(self.context.mesh)

            #only handles scalar data right now. TODO: add more support
            vtk_data = rset.GetCellData().GetScalars(self.context.primary_var)
            if vtk_data :
                rset.GetCellData().RemoveArray(self.context.primary_var)
                rset.GetCellData().AddArray(res)
                rset.GetCellData().SetScalars(res)
            else:
                rset.GetPointData().RemoveArray(self.context.primary_var)
                rset.GetPointData().AddArray(res)
                rset.GetPointData().SetScalars(res)
        else: #should not get here..
            rset = res
        return rset

def write_vtu(filename, atoms, data=None):
    from vtk import VTK_MAJOR_VERSION, vtkUnstructuredGrid, vtkPoints, vtkXMLUnstructuredGridWriter
    from vtk.util.numpy_support import numpy_to_vtk

    if isinstance(atoms, list):
        if len(atoms) > 1:
            raise ValueError('Can only write one configuration to a VTI file!')
        atoms = atoms[0]

    # Create a VTK grid of structured points
    ugd = vtkUnstructuredGrid()

    # add atoms as vtk Points
    p = vtkPoints()
    p.SetNumberOfPoints(len(atoms))
    p.SetDataTypeToDouble()
    for i,pos in enumerate(atoms.get_positions()):
        p.InsertPoint(i,pos[0],pos[1],pos[2])
    ugd.SetPoints(p)

    # add atomic numbers
    numbers = numpy_to_vtk(atoms.get_atomic_numbers(), deep=1)
    ugd.GetPointData().AddArray(numbers)
    numbers.SetName("atomic numbers")

    # add tags
    tags = numpy_to_vtk(atoms.get_tags(), deep=1)
    ugd.GetPointData().AddArray(tags)
    tags.SetName("tags")

    # add covalent radii
    from ase.data import covalent_radii
    radii = numpy_to_vtk(np.array([covalent_radii[i] for i in atoms.get_atomic_numbers()]), deep=1)
    ugd.GetPointData().AddArray(radii)
    radii.SetName("radii")

    # Save the UnstructuredGrid dataset to a VTK XML file.
    w = vtkXMLUnstructuredGridWriter()

    if fast:
        w.SetDataModeToAppend()
        w.EncodeAppendedDataOff()
    else:
        w.GetCompressor().SetCompressionLevel(0)
        w.SetDataModeToAscii()

    if isinstance(filename, str):
        w.SetFileName(filename)
    else:
        w.SetFileName(filename.name)
    if VTK_MAJOR_VERSION <= 5:
        w.SetInput(ugd)
    else:
        w.SetInputData(ugd)
    w.Write()

        def method(self):

            # multiblock += contact points
            output_a = self._contact_source_a.GetPolyDataOutput()
            output_b = self._contact_source_b.GetPolyDataOutput()

            id_f = numpy.where(
                abs(self._data[:, 0] - self._time) < 1e-15)[0]

            self.cpa_export = self._data[
                id_f, 2:5].copy()

            self.cpb_export = self._data[
                id_f, 5:8].copy()

            self.cn_export = self._data[
                id_f, 8:11].copy()

            self.cf_export = self._data[
                id_f, 11:14].copy()

            self.cpa_ = numpy_support.numpy_to_vtk(
                self.cpa_export)
            self.cpa_.SetName('contact_positions_A')

            self.cpb_ = numpy_support.numpy_to_vtk(
                self.cpb_export)
            self.cpb_.SetName('contact_positions_B')

            self.cn_ = numpy_support.numpy_to_vtk(
                self.cn_export)
            self.cn_.SetName('contact_normals')

            self.cf_ = numpy_support.numpy_to_vtk(
                self.cf_export)
            self.cf_.SetName('contact_forces')

            output_a.Allocate(len(self.cpa_export), 1)

            cpa_points = vtk.vtkPoints()
            cpa_points.SetNumberOfPoints(len(self.cpa_export))
            cpa_points.SetData(self.cpa_)
            output_a.SetPoints(cpa_points)

            # normal and forces are attached to A points
            output_a.GetPointData().AddArray(self.cn_)
            output_a.GetPointData().AddArray(self.cf_)

            output_b.Allocate(len(self.cpb_export), 1)

            cpb_points = vtk.vtkPoints()
            cpb_points.SetNumberOfPoints(len(self.cpb_export))
            cpb_points.SetData(self.cpb_)

            output_b.SetPoints(cpb_points)

    def attach_pos(self):

        pos = self._pos
        rad = self._rad

        if pos is not None and rad is not None:
            positions_vtk = numpy_support.numpy_to_vtk(num_array=pos.ravel(), deep=True, array_type=vtk.VTK_FLOAT)
            positions_vtk.SetName("positions")

            radius_vtk = numpy_support.numpy_to_vtk(num_array=rad.ravel(), deep=True, array_type=vtk.VTK_FLOAT)
            radius_vtk.SetName("radius")

            sphere = vtk.vtkSphereSource()
            sphere.SetRadius(1.0)

            ballGlyph = vtk.vtkGlyph3D()

            if vtk.VTK_MAJOR_VERSION <= 5:
                ballGlyph.SetSource(sphere.GetOutput())
            else:
                ballGlyph.SetSourceConnection(sphere.GetOutputPort())

            polydata = vtk.vtkPolyData()
            polydata.SetPoints(self._points)
            polydata.GetPointData().AddArray(radius_vtk)
            polydata.GetPointData().SetActiveScalars("radius") # this scales the source (sphere) radius (1.0)
            
            ballGlyph.SetInputData(polydata)

            #ballGlyph.SetScaleModeToDataScalingOn() 
            mapper = vtk.vtkPolyDataMapper()

            if vtk.VTK_MAJOR_VERSION <= 5:
               mapper.SetInput(ballGlyph.GetOutput())
            else:
               mapper.SetInputConnection(ballGlyph.GetOutputPort())

            # Set colors depending on the color transfer functions
            # mapper.SetLookupTable(self.colorTransferFunction)

            # actor
            ballActor = vtk.vtkActor()
            ballActor.GetProperty().SetColor(0,0,1)
            ballActor.SetMapper(mapper)

            #self._ren.AddActor(ballActor)
            self._marker2 = vtk.vtkOrientationMarkerWidget()
            self._marker2.SetInteractor( self.widget._Iren )
            self._marker2.SetOrientationMarker( ballActor )
            self._marker2.SetViewport(0.75,0,1,0.25)
            self._marker2.SetEnabled(1)

        else:
            print("No particles found. Make sure the particles loaded have positions and radii.")

def polydata_from_numpy(coords):
    """
    Converts Numpy Array to vtkPolyData

    Parameters
    ----------
    coords: array, shape= [number of points, point features]
        array containing the point cloud in which each point has three coordinares [x, y, z]
        and none, one or three values corresponding to colors.

    Returns:
    --------
    PolyData: vtkPolyData
        concrete dataset represents vertices, lines, polygons, and triangle strips

    """

    Npts, Ndim = np.shape(coords)

    Points = vtkPoints()
    ntype = get_numpy_array_type(Points.GetDataType())
    coords_vtk = numpy_to_vtk(np.asarray(coords[:,0:3], order='C',dtype=ntype), deep=1)
    Points.SetNumberOfPoints(Npts)
    Points.SetData(coords_vtk)

    Cells = vtkCellArray()
    ids = np.arange(0,Npts, dtype=np.int64).reshape(-1,1)
    IDS = np.concatenate([np.ones(Npts, dtype=np.int64).reshape(-1,1), ids],axis=1)
    ids_vtk = numpy_to_vtkIdTypeArray(IDS, deep=True)

    Cells.SetNumberOfCells(Npts)
    Cells.SetCells(Npts,ids_vtk)

    if Ndim == 4:
        color = [128]*len(coords)
        color = np.c_[color, color, color]
    elif Ndim == 6:
        color = coords[:,3:]
    else:
        color = [[128, 128, 128]]*len(coords)

    color_vtk = numpy_to_vtk(
            np.ascontiguousarray(color, dtype=get_vtk_to_numpy_typemap()[VTK_UNSIGNED_CHAR]),
            deep=True)

    color_vtk.SetName("colors")

    PolyData = vtkPolyData()
    PolyData.SetPoints(Points)
    PolyData.SetVerts(Cells)
    PolyData.GetPointData().SetScalars(color_vtk)
    return PolyData

  def extract_particles(self,ids):
    data=vtk.vtkPolyData()
    points=vtk_to_numpy(self._in_vtk.GetPoints().GetData())

    s_points=vtk.vtkPoints()
    cell_arr=vtk.vtkCellArray()
    
    #s_points.SetData(numpy_to_vtk(points[ids,:]))
    #s_points.SetNumberOfPoints(s_points.GetData().GetNumberOfTuples())
    s_p = points[ids,:]
    s_points.SetNumberOfPoints(s_p.shape[0])
    cell_arr.SetNumberOfCells(s_p.shape[0])
    for kk in xrange(s_p.shape[0]):
      s_points.SetPoint(kk,s_p[kk,0],s_p[kk,1],s_p[kk,2])
      cell_arr.InsertNextCell(1)
      cell_arr.InsertCellPoint(kk)
    
    data.SetPoints(s_points)
    data.SetVerts(cell_arr)
    
    #Transfer point data and field data
    for pd,out_pd in zip([self._in_vtk.GetPointData(),self._in_vtk.GetFieldData()],[data.GetPointData(),data.GetFieldData()]):
      for k in xrange(pd.GetNumberOfArrays()):
        arr=vtk_to_numpy(pd.GetArray(pd.GetArrayName(k)))
        if len(arr.shape) == 1:
          s_vtk_arr=numpy_to_vtk(arr[ids],1)
        else:
          s_vtk_arr=numpy_to_vtk(arr[ids,:],1)
        
        s_vtk_arr.SetName(pd.GetArrayName(k))
        out_pd.AddArray(s_vtk_arr)
    
    return data
                      
    #Method to do the extraction using a vtk pipeline (experimental with seg fault)
    def extract_using_vtk(self,ids):
      node=vtk.vtkSelectionNode()
      sel = vtk.vtkSelection()
      node.GetProperties().Set(vtk.vtkSelectionNode.CONTENT_TYPE(),\
                               vtk.vtkSelectionNode.INDICES)
      node.GetProperties().Set(vtk.vtkSelectionNode.FIELD_TYPE(),\
                               vtk.vtkSelectionNode.POINT)
      
      #Create Id Array with point Ids for each cluster
      vtk_ids=numpy_to_vtkIdTypeArray(ids)
      node.SetSelectionList(vtk_ids)
      #sel_filter = vtk.vtkExtractSelectedPolyDataIds()
      sel_filter = vtk.vtkExtractSelection()
      sel_filter.SetInput(0,self._in_vtk)
      sel_filter.SetInput(1,sel)
      sel_filter.Update()
      return sel_filter.GetOutput()

	def GetNodeOneScaleCoeffTable(self, node_id):
		"""Returns a table of the wavelet coefficients at a single node at a single
		scale for plotting on a scatter plot. Relying on icicle_view already having
		called GetWaveletCoeffImages() with correct positions of leaf nodes in view,
		otherwise just using original Matlab-saved LeafNodes ordering.
		This version supports category labels."""

		if self.data_loaded:

			# For a given node_id, concatenate wavelet coeffs in proper order
			# (according to leaf node positions in icicle view if it was set already)
			# Columns of table will be rows of the wavelet coeffs image

			scale = self.Scales[node_id]

			leaf_offspring = N.array(list(self.get_leaf_children(self.cp, node_id)))
			offspring_idxs = self.LeafNodesImap[leaf_offspring]
			offspring_pos = self.mapped_leaf_pos[offspring_idxs]
			sorted_offspring_pos_idxs = N.argsort(offspring_pos)
			sorted_offspring_idxs = offspring_idxs[sorted_offspring_pos_idxs]
			img_tuple = tuple(pp for pp in self.CelCoeffs[sorted_offspring_idxs, scale])
			# The image comes out with shape (npts, ndims)
			# May need to reorder (reverse) this...?
			img = N.concatenate(img_tuple, axis=0)

			table = vtk.vtkTable()
			for ii in range(img.shape[1]):
				column = VN.numpy_to_vtk(img[:,ii].copy(), deep=True)
				column.SetName(str(scale) + '.' + str(ii))
				table.AddColumn(column)

			IDtuple = tuple(self.PointsInNet[xx] for xx in self.LeafNodes[sorted_offspring_idxs])
			IDarray = N.concatenate(IDtuple)
			# Trying to set PedigreeIds to that parallel coords selections have correct IDs
			IDvtk = VN.numpy_to_vtk(IDarray, deep=True)
			IDvtk.SetName('pedigree_ids')
			table.AddColumn(IDvtk)
			table.GetRowData().SetActivePedigreeIds('pedigree_ids')

			# Adding in category labels
			# Name needs to end in _ids so plots will ignore it
			if self.hasLabels:
				for ii in range(self.cat_labels.shape[0]):
					CATvtk = VN.numpy_to_vtk(self.cat_labels[ii,IDarray], deep=True)
					CATvtk.SetName(self.label_names[ii])
					table.AddColumn(CATvtk)

			return table


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

def ply_exporter(mesh, file_path, binary=False, **kwargs):
    r"""
    Given a file path to write in to write out the mesh data.
    No value is returned. Only file paths are supported and if a file handle
    is passed it will be ignored and a warning will be raised.

    Note that this does not save out textures of textured images, and so should
    not be used in isolation.

    Parameters
    ----------
    file_path : `str`
        The full path where the obj will be saved out.
    mesh : :map:`TriMesh`
        Any subclass of :map:`TriMesh`. If :map:`TexturedTriMesh` texture
        coordinates will be saved out. Note that :map:`ColouredTriMesh`
        will only have shape data saved out, as .PLY doesn't robustly support
        per-vertex colour information.
    binary: `bool`, optional
        Specify whether to format output in binary or ascii, defaults to False
    """
    import vtk
    from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray

    file_path = _enforce_only_paths_supported(file_path, 'PLY')

    polydata = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    points.SetData(numpy_to_vtk(mesh.points))
    polydata.SetPoints(points)

    cells = vtk.vtkCellArray()
    counts = np.empty((mesh.trilist.shape[0], 1), dtype=np.int)
    counts.fill(3)
    tris = np.concatenate((counts, mesh.trilist), axis=1)
    cells.SetCells(mesh.trilist.shape[0], numpy_to_vtkIdTypeArray(tris))
    polydata.SetPolys(cells)

    if isinstance(mesh, TexturedTriMesh):
        pointdata = polydata.GetPointData()
        pointdata.SetTCoords(numpy_to_vtk(mesh.tcoords.points))

    ply_writer = vtk.vtkPLYWriter()
    ply_writer.SetFileName(str(file_path))
    ply_writer.SetInputData(polydata)
    if not binary:
        ply_writer.SetFileTypeToASCII()
    else:
        ply_writer.SetFileTypeToBinary()
    ply_writer.Update()
    ply_writer.Write()

def gocad2simpegMeshIndex(gcFile,mesh,extractBoundaryCells=True,extractInside=True):
    """"
    Function to read gocad polystructure file and output indexes of mesh with in the structure.

    """

    # Make the polydata
    polyData = gocad2vtp(gcFile)

    # Make implicit func
    ImpDistFunc = vtk.vtkImplicitPolyDataDistance()
    ImpDistFunc.SetInput(polyData)

    # Convert the mesh
    vtkMesh = vtk.vtkRectilinearGrid()
    vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz)
    vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx,deep=1))
    vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy,deep=1))
    vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz,deep=1))
    # Add indexes cell data to the object
    vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC),deep=1)
    vtkInd.SetName('Index')
    vtkMesh.GetCellData().AddArray(vtkInd)

    # Define the extractGeometry
    extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor
    extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) #
    extractImpDistRectGridFilt.SetInputData(vtkMesh)

    # Set extraction type
    if extractBoundaryCells is True:
        extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
    else:
        extractImpDistRectGridFilt.ExtractBoundaryCellsOff()

    if extractInside is True:
        extractImpDistRectGridFilt.ExtractInsideOn()
    else:
        extractImpDistRectGridFilt.ExtractInsideOff()

    print "Extracting indices from grid..."
    # Executing the pipe
    extractImpDistRectGridFilt.Update()

    # Get index inside
    insideGrid = extractImpDistRectGridFilt.GetOutput()
    insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))


    # Return the indexes inside
    return insideGrid

	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 VertFacetoPoly(new_pt, new_fc):
    """ Creates a vtk polydata object given points and triangular faces """
    
    # Convert points to vtkfloat object
    points = np.vstack(new_pt)
    vtkArray = VN.numpy_to_vtk(np.ascontiguousarray(points), deep=True)#, deep=True) 
    points = vtk.vtkPoints()
    points.SetData(vtkArray) 
    
    # Convert faces to vtk cells object
    ints = np.ones(len(new_fc), 'int')*3
    cells = np.hstack((ints.reshape(-1, 1), np.vstack(new_fc)))
    cells = np.ascontiguousarray(np.hstack(cells).astype('int64'))
    vtkcells = vtk.vtkCellArray()
    vtkcells.SetCells(cells.shape[0], VN.numpy_to_vtkIdTypeArray(cells, deep=True))

    # Create polydata object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    pdata.SetPolys(vtkcells)
    
    # Remove duplicate verticies
    clean = vtk.vtkCleanPolyData()
    clean.ConvertPolysToLinesOff()
    clean.ConvertLinesToPointsOff()
    clean.ConvertStripsToPolysOff()
    if vtk.vtkVersion().GetVTKMajorVersion() > 5:
        clean.SetInputData(pdata)
    else:
        clean.SetInput(pdata) 
    clean.Update()

    return clean.GetOutput()