Python vtk.vtkCellArray() Examples

def create_pointcloud_polydata(points, colors=None):
    """https://github.com/lmb-freiburg/demon
    Creates a vtkPolyData object with the point cloud from numpy arrays
    
    points: numpy.ndarray
        pointcloud with shape (n,3)
    
    colors: numpy.ndarray
        uint8 array with colors for each point. shape is (n,3)

    Returns vtkPolyData object
    """
    vpoints = vtk.vtkPoints()
    vpoints.SetNumberOfPoints(points.shape[0])
    for i in range(points.shape[0]):
        vpoints.SetPoint(i, points[i])
    vpoly = vtk.vtkPolyData()
    vpoly.SetPoints(vpoints)
    
    if not colors is None:
        vcolors = vtk.vtkUnsignedCharArray()
        vcolors.SetNumberOfComponents(3)
        vcolors.SetName("Colors")
        vcolors.SetNumberOfTuples(points.shape[0])
        for i in range(points.shape[0]):
            vcolors.SetTuple3(i ,colors[i,0],colors[i,1], colors[i,2])
        vpoly.GetPointData().SetScalars(vcolors)

    vcells = vtk.vtkCellArray()
    
    for i in range(points.shape[0]):
        vcells.InsertNextCell(1)
        vcells.InsertCellPoint(i)
        
    vpoly.SetVerts(vcells)
    
    return vpoly 

def build_globe(self):
        """Generates the globe as ``vtkPolyData``"""
        # NOTE: https://gitlab.kitware.com/paraview/paraview/issues/19417
        from scipy.spatial import Delaunay
        pos, tex = self.create_sphere()
        pts = self.spherical_to_cartesian(pos[:,0], pos[:,1])
        points = interface.points_to_poly_data(pts).GetPoints()
        texcoords = interface.convert_array(tex, name='Texture Coordinates')
        # Now generate triangles
        cell_connectivity = Delaunay(pos).simplices.astype(int)
        cells = vtk.vtkCellArray()
        cells.SetNumberOfCells(cell_connectivity.shape[0])
        cells.SetCells(cell_connectivity.shape[0], interface.convert_cell_conn(cell_connectivity))
        # Generate output
        output = vtk.vtkPolyData()
        output.SetPoints(points)
        output.GetPointData().SetTCoords(texcoords)
        output.SetPolys(cells)
        return output 

def createDEM_v2():

    ds = xr.open_dataset('../data/output/Peru_20160601-20180530_comp4.nc')
    
    points = vtk.vtkPoints()
    quad = vtk.vtkQuad()
    cells = vtk.vtkCellArray()
    
    numPoints = ds.south_north.size*ds.west_east.size
   
    print('Write points \n')
    for i,j in product(ds.south_north.values,ds.west_east.values):
            points.InsertNextPoint(ds.lat.isel(south_north=i,west_east=j), ds.lon.isel(south_north=i,west_east=j), ds.HGT.sel(south_north=i,west_east=j).values/6370000.0)
    
    print('Write cells \n')
    for idx in range(points.GetNumberOfPoints()-ds.west_east.size):
        if (idx%ds.west_east.size != 0):
            quad.GetPointIds().SetId(0,idx)
            quad.GetPointIds().SetId(1,idx+1)
            quad.GetPointIds().SetId(2,idx+ds.west_east.size+1)
            quad.GetPointIds().SetId(3,idx+ds.west_east.size)
            cells.InsertNextCell(quad)

    print('Create unstructured grid \n') 
    grid = vtk.vtkUnstructuredGrid()
    grid.SetPoints(points)
    grid.SetCells(vtk.VTK_QUAD, cells)
    
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName('cosipy.vtu')
    writer.SetInputData(grid)
    writer.Write() 

def CreatePolygon(self, points):
        "Create a 3D plane from Nx3 numpy array"
        self.verts = vtk.vtkPoints()
        polygon = vtk.vtkPolygon()

        polygon_pid = polygon.GetPointIds()
        for i, p in enumerate(points):
            self.verts.InsertNextPoint(*tuple(p))
            polygon_pid.InsertNextId(i)

        polygon_cell = vtk.vtkCellArray()
        polygon_cell.InsertNextCell(polygon)

        self.pd = vtk.vtkPolyData()
        self.pd.SetPoints(self.verts)
        self.pd.SetPolys(polygon_cell)

        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInputData(self.pd)
        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper) 

def poly_data_builder(state, zf, register):
    instance = vtk.vtkPolyData()
    register.update({state['id']: instance})
    # geometry
    if 'points' in state['properties']:
        points = state['properties']['points']
        vtkpoints = vtk.vtkPoints()
        data_arr = ARRAY_TYPES[points['dataType']]()
        fill_array(data_arr, points, zf)
        vtkpoints.SetData(data_arr)
        instance.SetPoints(vtkpoints)
    for cell_type in ['verts', 'lines', 'polys', 'strips']:
        if cell_type in state['properties']:
            cell_arr = vtk.vtkCellArray()
            fill_array(cell_arr.GetData(), state['properties'][cell_type], zf)
            getattr(instance, 'Set' + capitalize(cell_type))(cell_arr)
    
    # datasets
    fields = state['properties']['fields']
    for dataset in fields:
        data_arr = ARRAY_TYPES[dataset['dataType']]()
        fill_array(data_arr, dataset, zf)
        location = getattr(instance, 'Get' + capitalize(dataset['location']))()
        getattr(location, capitalize(dataset['registration']))(data_arr) 

def get_polydata_from(points, tr_re):

    numberPoints = len(points)
    Points = vtkPoints()
    ntype = get_numpy_array_type(Points.GetDataType())
    points_vtk = numpy_to_vtk(np.asarray(points, order='C',dtype=ntype), deep=1)
    Points.SetNumberOfPoints(numberPoints)
    Points.SetData(points_vtk)

    Triangles = vtkCellArray()
    for item in tr_re:
        Triangle = vtkTriangle()
        Triangle.GetPointIds().SetId(0,item[0])
        Triangle.GetPointIds().SetId(1,item[1])
        Triangle.GetPointIds().SetId(2,item[2])
        Triangles.InsertNextCell(Triangle)

    polydata = vtkPolyData()
    polydata.SetPoints(Points)
    polydata.SetPolys(Triangles)

    polydata.Modified()
    polydata.Update()

    return polydata 

def create_surface_triangles(simplices):
        """
        Method to create the Triangles that form the surfaces
        Args:
            simplices (numpy.array): 2D array with the value of the vertices that form every single triangle

        Returns:
            vtk.vtkTriangle
        """

        Triangles = vtk.vtkCellArray()
        Triangle = vtk.vtkTriangle()

        for s in simplices:
            Triangle.GetPointIds().SetId(0, s[0])
            Triangle.GetPointIds().SetId(1, s[1])
            Triangle.GetPointIds().SetId(2, s[2])

            Triangles.InsertNextCell(Triangle)
        return Triangles 

def face_points2actor(fps_df):

    cell_array = vtk.vtkCellArray()
    points = vtk.vtkPoints()
    point_id = 0
    for i in range(fps_df.shape[0]):
        polygon = vtk.vtkPolygon()
        polygon.GetPointIds().SetNumberOfIds(3)
        for n in range(3):
            points.InsertNextPoint(fps_df.ix[i, 0+3*n], fps_df.ix[i, 1+3*n], fps_df.ix[i, 2+3*n])
            polygon.GetPointIds().SetId(n, point_id)
            point_id += 1
        cell_array.InsertNextCell(polygon)
    polydata = vtk.vtkPolyData()
    polydata.SetPoints(points)
    polydata.SetPolys(cell_array)

    # mapper
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInput(polydata)
    # actor
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    return actor, polydata 

def CreatePolyLine(self, points):
        "Create a 3D line from Nx3 numpy array"
        self.verts = vtk.vtkPoints()
        polyline = vtk.vtkPolyLine()

        polyline_pid = polyline.GetPointIds()
        for i, p in enumerate(points):
            self.verts.InsertNextPoint(*tuple(p))
            polyline_pid.InsertNextId(i)

        polyline_cell = vtk.vtkCellArray()
        polyline_cell.InsertNextCell(polyline)

        self.pd = vtk.vtkPolyData()
        self.pd.SetPoints(self.verts)
        self.pd.SetLines(polyline_cell)

        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInputData(self.pd)
        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper)
        self.actor.GetProperty().SetLineStipplePattern(0xF0F0)
        self.actor.GetProperty().SetLineStippleRepeatFactor(1)
        self.actor.GetProperty().SetLineWidth(1.5) 

def CreateTriangles(self, pc, triangles):
        "Create a mesh from triangles"

        self.verts = vtk.vtkPoints()
        self.points_npy = pc[:, :3].copy()
        self.verts.SetData(numpy_support.numpy_to_vtk(self.points_npy))

        nTri = len(triangles)
        self.cells = vtk.vtkCellArray()
        self.pd = vtk.vtkPolyData()
        # - Note that the cell array looks like this: [3 vtx0 vtx1 vtx2 3 vtx3 ... ]
        self.cells_npy = np.column_stack(
            [np.full(nTri, 3, dtype=np.int64), triangles.astype(np.int64)]
        ).ravel()
        self.cells.SetCells(nTri, numpy_support.numpy_to_vtkIdTypeArray(self.cells_npy))

        self.pd.SetPoints(self.verts)
        self.pd.SetPolys(self.cells)

        self.SetupPipelineMesh() 

def points2actor(xyz, apoint_size):
    import vtk
    points = vtk.vtkPoints()
    # Create the topology of the point (a vertex)
    vertices = vtk.vtkCellArray()
    # Add points
    for i in range(0, len(xyz)):
        p = xyz.loc[i].values.tolist()
        point_id = points.InsertNextPoint(p)
        vertices.InsertNextCell(1)
        vertices.InsertCellPoint(point_id)
    # Create a poly data object
    polydata = vtk.vtkPolyData()
    # Set the points and vertices we created as the geometry and topology of the polydata
    polydata.SetPoints(points)
    polydata.SetVerts(vertices)
    polydata.Modified()
    # Mapper for points
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInput(polydata)
    # ACTOR for points
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetPointSize(apoint_size)
    return actor, polydata 

def draw_box(x):
    cube = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    polys = vtk.vtkCellArray()
    scalars = vtk.vtkFloatArray()

    for i in range(8):
        points.InsertPoint(i, x[i])
    for i in range(6):
        polys.InsertNextCell(mkVtkIdList(pts[i]))
    for i in range(8):
        scalars.InsertTuple1(i, i)

    cube.SetPoints(points)
    del points
    cube.SetPolys(polys)
    del polys
    cube.GetPointData().SetScalars(scalars)
    del scalars

    cubeMapper = vtk.vtkPolyDataMapper()
    if vtk.VTK_MAJOR_VERSION <= 5:
        cubeMapper.SetInput(cube)
    else:
        cubeMapper.SetInputData(cube)
    cubeMapper.SetScalarRange(0, 7)
    # cubeMapper.SetScalarVisibility(2)
    cubeActor = vtk.vtkActor()
    cubeActor.SetMapper(cubeMapper)
    cubeActor.GetProperty().SetOpacity(0.4)
    return cubeActor 

def CreateAxesSimplified(self, length):
        "Create a simplified coordinate axes system with 3 lines"
        points = vtk.vtkPoints()
        lines = vtk.vtkCellArray()
        colors = vtk.vtkUnsignedCharArray()
        colors.SetNumberOfComponents(3)
        colors.SetName("Colors")

        ptIds = [None, None]
        end_points = length * np.eye(3)
        line_colors = 255 * np.eye(3, dtype="i")
        for i in range(3):
            ptIds[0] = points.InsertNextPoint([0, 0, 0])
            ptIds[1] = points.InsertNextPoint(end_points[i])
            lines.InsertNextCell(2, ptIds)

            colors.InsertNextTuple3(*line_colors[i])
            colors.InsertNextTuple3(*line_colors[i])

        # Add the lines to the polydata container
        self.pd = vtk.vtkPolyData()
        self.pd.SetPoints(points)
        self.pd.GetPointData().SetScalars(colors)
        self.pd.SetLines(lines)

        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInputData(self.pd)
        self.mapper.Update()
        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper) 

def clearPoints(self):
        """
        Clear all points from the point cloud
        :return: None
        """
        self.vtkPoints = vtk.vtkPoints()
        self.vtkCells = vtk.vtkCellArray()
        self.vtkDepth = vtk.vtkDoubleArray()
        self.vtkDepth.SetName('DepthArray')
        self.vtkPolyData.SetPoints(self.vtkPoints)
        self.vtkPolyData.SetVerts(self.vtkCells)
        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray') 

def create_empty_mesh(self):
        self.pcoords = vtk.vtkFloatArray()
        self.pcoords.SetNumberOfComponents(3)
        self.points = vtk.vtkPoints()
        self.polys = vtk.vtkCellArray()
        self.poly_data = vtk.vtkPolyData() 

def clear_points(self):
        self.vtkPoints = vtk.vtkPoints()
        self.vtkCells = vtk.vtkCellArray()
        self.vtkDepth = vtk.vtkDoubleArray()
        self.vtkDepth.SetName('DepthArray')
        self.vtkPolyData.SetPoints(self.vtkPoints)
        self.vtkPolyData.SetVerts(self.vtkCells)
        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray') 

def clear_points(self):
        self.vtkPoints = vtk.vtkPoints()
        self.vtkCells = vtk.vtkCellArray()
        self.vtkDepth = vtk.vtkDoubleArray()
        self.vtkDepth.SetName('DepthArray')
        self.vtkPolyData.SetPoints(self.vtkPoints)
        self.vtkPolyData.SetVerts(self.vtkCells)
        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray') 

def __init__(self, nparray):
        super(VTKActorWrapper, self).__init__()

        self.nparray = nparray

        nCoords = nparray.shape[0]
        nElem = nparray.shape[1]

        self.verts = vtk.vtkPoints()
        self.cells = vtk.vtkCellArray()
        self.scalars = None

        self.pd = vtk.vtkPolyData()
        self.verts.SetData(vtk_np.numpy_to_vtk(nparray))
        self.cells_npy = np.vstack([np.ones(nCoords,dtype=np.int64),
                               np.arange(nCoords,dtype=np.int64)]).T.flatten()
        self.cells.SetCells(nCoords,vtk_np.numpy_to_vtkIdTypeArray(self.cells_npy))
        self.pd.SetPoints(self.verts)
        self.pd.SetVerts(self.cells)

        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInputDataObject(self.pd)

        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper)
        self.actor.GetProperty().SetRepresentationToPoints()
        self.actor.GetProperty().SetColor(0.0,1.0,0.0) 

def create_pointcloud_polydata(points, colors=None):
    """Creates a vtkPolyData object with the point cloud from numpy arrays

    points: numpy.ndarray
        pointcloud with shape (n,3)

    colors: numpy.ndarray
        uint8 array with colors for each point. shape is (n,3)

    Returns vtkPolyData object
    """
    import vtk
    vpoints = vtk.vtkPoints()
    vpoints.SetNumberOfPoints(points.shape[0])
    for i in range(points.shape[0]):
        vpoints.SetPoint(i, points[i])
    vpoly = vtk.vtkPolyData()
    vpoly.SetPoints(vpoints)

    if not colors is None:
        vcolors = vtk.vtkUnsignedCharArray()
        vcolors.SetNumberOfComponents(3)
        vcolors.SetName("Colors")
        vcolors.SetNumberOfTuples(points.shape[0])
        for i in range(points.shape[0]):
            vcolors.SetTuple3(i ,colors[i,0],colors[i,1], colors[i,2])
        vpoly.GetPointData().SetScalars(vcolors)

    vcells = vtk.vtkCellArray()

    for i in range(points.shape[0]):
        vcells.InsertNextCell(1)
        vcells.InsertCellPoint(i)

    vpoly.SetVerts(vcells)

    return vpoly 

def _create_vtp_mesh(self):
        '''Internal function to creat a VTP mesh from the imported mesh data
        '''
        if self.VTK_installed is True:

            self.vtp_mesh    = vtk.vtkPolyData()
            points  = vtk.vtkPoints()
            polys   = vtk.vtkCellArray()

            for i in range(self.points.shape[0]):

                points.InsertPoint(i, self.points[i])


            for i in range(self.faces.shape[0]):

                polys.InsertNextCell(_mk_vtk_id_list(self.faces[i]))

            self.vtp_mesh.SetPoints(points)
            self.vtp_mesh.SetPolys(polys)

    # def _collapse(self,plane=2,value=0.0,direction=1):
    #This function is not yet working 100%

    #     '''Collapse points
    #     '''
    #     for face,face_n in xrange(self.faces.shape[0]):
    #
    #         for j in xrange(self.faces[i].size):
    #
    #             p = int(self.faces[i][j])
    #
    #             if self.points[p][plane] > value*direction:
    #
    #                 self.points[p][plane] = value 

def _create_vtp_mesh(self):
        '''Internal function to creat a VTP mesh from the imported mesh data
        '''
        if self.VTK_installed is True:

            self.vtp_mesh    = vtk.vtkPolyData()
            points  = vtk.vtkPoints()
            polys   = vtk.vtkCellArray()

            for i in range(self.points.shape[0]):

                points.InsertPoint(i, self.points[i])


            for i in range(self.faces.shape[0]):

                polys.InsertNextCell(_mk_vtk_id_list(self.faces[i]))

            self.vtp_mesh.SetPoints(points)
            self.vtp_mesh.SetPolys(polys)

    # def _collapse(self,plane=2,value=0.0,direction=1):
    #This function is not yet working 100%

    #     '''Collapse points
    #     '''
    #     for face,face_n in xrange(self.faces.shape[0]):
    #
    #         for j in xrange(self.faces[i].size):
    #
    #             p = int(self.faces[i][j])
    #
    #             if self.points[p][plane] > value*direction:
    #
    #                 self.points[p][plane] = value 

def clearPoints(self):
        self.vtkPoints = vtk.vtkPoints()
        self.vtkCells = vtk.vtkCellArray()
        self.vtkDepth = vtk.vtkDoubleArray()
        self.vtkDepth.SetName('DepthArray')
        self.vtkPolyData.SetPoints(self.vtkPoints)
        self.vtkPolyData.SetVerts(self.vtkCells)
        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray')

        self.Colors = vtk.vtkUnsignedCharArray()
        self.Colors.SetNumberOfComponents(3) 
        self.Colors.SetName("Colors") 

def surface_geom_to_vtk(surfgeom):
    """Convert the triangulated surface to a :class:`pyvista.UnstructuredGrid`
    object

    Args:
        surfgeom (:class:`omf.surface.SurfaceGeometry`): the surface geomotry to
            convert
    """

    output = vtk.vtkUnstructuredGrid()
    pts = vtk.vtkPoints()
    cells = vtk.vtkCellArray()

    # Generate the points
    pts.SetNumberOfPoints(surfgeom.num_nodes)
    pts.SetData(nps.numpy_to_vtk(surfgeom.vertices))

    # Generate the triangle cells
    cellConn = surfgeom.triangles.array
    cellsMat = np.concatenate(
        (np.ones((cellConn.shape[0], 1), dtype=np.int64)*cellConn.shape[1], cellConn),
        axis=1).ravel()
    cells = vtk.vtkCellArray()
    cells.SetNumberOfCells(cellConn.shape[0])
    cells.SetCells(cellConn.shape[0],
            nps.numpy_to_vtk(cellsMat, deep=True, array_type=vtk.VTK_ID_TYPE))

    # Add to output
    output.SetPoints(pts)
    output.SetCells(vtk.VTK_TRIANGLE, cells)
    return pyvista.wrap(output) 

def point_set_to_vtk(pse):
    """Convert the point set to a :class:`pyvista.PolyData` data object.

    Args:
        pse (:class:`omf.pointset.PointSetElement`): The point set to convert

    Return:
        :class:`pyvista.PolyData`
    """

    points = pse.geometry.vertices
    npoints = pse.geometry.num_nodes

    # Make VTK cells array
    cells = np.hstack((np.ones((npoints, 1)),
                       np.arange(npoints).reshape(-1, 1)))
    cells = np.ascontiguousarray(cells, dtype=np.int64)
    cells = np.reshape(cells, (2*npoints))
    vtkcells = vtk.vtkCellArray()
    vtkcells.SetCells(npoints, nps.numpy_to_vtk(cells, deep=True, array_type=vtk.VTK_ID_TYPE))

    # Convert points to vtk object
    pts = vtk.vtkPoints()
    pts.SetNumberOfPoints(pse.geometry.num_nodes)
    pts.SetData(nps.numpy_to_vtk(points))

    # Create polydata
    output = vtk.vtkPolyData()
    output.SetPoints(pts)
    output.SetVerts(vtkcells)

    # Now add point data:
    add_data(output, pse.data)

    add_textures(output, pse.textures, pse.name)

    return pyvista.wrap(output) 

def _create_vtp_mesh(self):
        '''Internal function to creat a VTP mesh from the imported mesh data
        '''
        if self.VTK_installed is True:

            self.vtp_mesh    = vtk.vtkPolyData()
            points  = vtk.vtkPoints()
            polys   = vtk.vtkCellArray()

            for i in range(self.points.shape[0]):

                points.InsertPoint(i, self.points[i])


            for i in range(self.faces.shape[0]):

                polys.InsertNextCell(_mk_vtk_id_list(self.faces[i]))

            self.vtp_mesh.SetPoints(points)
            self.vtp_mesh.SetPolys(polys)

    # def _collapse(self,plane=2,value=0.0,direction=1):
    #This function is not yet working 100%

    #     '''Collapse points
    #     '''
    #     for face,face_n in xrange(self.faces.shape[0]):
    #
    #         for j in xrange(self.faces[i].size):
    #
    #             p = int(self.faces[i][j])
    #
    #             if self.points[p][plane] > value*direction:
    #
    #                 self.points[p][plane] = value 

def clear_points(self):
        self.vtkPoints = vtk.vtkPoints()
        self.vtkCells = vtk.vtkCellArray()
        self.vtkDepth = vtk.vtkDoubleArray()
        self.vtkDepth.SetName('DepthArray')
        self.vtkPolyData.SetPoints(self.vtkPoints)
        self.vtkPolyData.SetVerts(self.vtkCells)
        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray') 

def save_vtk_labels(filename, tracts, scalars, lines_indices=None):

	lengths = [len(p) for p in tracts]
	line_starts = ns.numpy.r_[0, ns.numpy.cumsum(lengths)]
	if lines_indices is None:
		lines_indices = [ns.numpy.arange(length) + line_start for length, line_start in izip(lengths, line_starts)]

	ids = ns.numpy.hstack([ns.numpy.r_[c[0], c[1]] for c in izip(lengths, lines_indices)])
	vtk_ids = ns.numpy_to_vtkIdTypeArray(ids, deep=True)

	cell_array = vtk.vtkCellArray()
	cell_array.SetCells(len(tracts), vtk_ids)
	points = ns.numpy.vstack(tracts).astype(ns.get_vtk_to_numpy_typemap()[vtk.VTK_DOUBLE])
	points_array = ns.numpy_to_vtk(points, deep=True)

	poly_data = vtk.vtkPolyData()
	vtk_points = vtk.vtkPoints()
	vtk_points.SetData(points_array)
	poly_data.SetPoints(vtk_points)
	poly_data.SetLines(cell_array)
	poly_data.GetPointData().SetScalars(ns.numpy_to_vtk(scalars))
	poly_data.BuildCells()
#    poly_data.SetScalars(scalars)

	if filename.endswith('.xml') or filename.endswith('.vtp'):
		writer = vtk.vtkXMLPolyDataWriter()
		writer.SetDataModeToBinary()
	else:
		writer = vtk.vtkPolyDataWriter()
		writer.SetFileTypeToBinary()


	writer.SetFileName(filename)
	if hasattr(vtk, 'VTK_MAJOR_VERSION') and vtk.VTK_MAJOR_VERSION > 5:
		writer.SetInputData(poly_data)
	else:
		writer.SetInput(poly_data)
	writer.Write() 

def save_vtk(filename, tracts, lines_indices=None, scalars = None):

	lengths = [len(p) for p in tracts]
	line_starts = ns.numpy.r_[0, ns.numpy.cumsum(lengths)]
	if lines_indices is None:
		lines_indices = [ns.numpy.arange(length) + line_start for length, line_start in izip(lengths, line_starts)]

	ids = ns.numpy.hstack([ns.numpy.r_[c[0], c[1]] for c in izip(lengths, lines_indices)])
	vtk_ids = ns.numpy_to_vtkIdTypeArray(ids, deep=True)

	cell_array = vtk.vtkCellArray()
	cell_array.SetCells(len(tracts), vtk_ids)
	points = ns.numpy.vstack(tracts).astype(ns.get_vtk_to_numpy_typemap()[vtk.VTK_DOUBLE])
	points_array = ns.numpy_to_vtk(points, deep=True)

	poly_data = vtk.vtkPolyData()
	vtk_points = vtk.vtkPoints()
	vtk_points.SetData(points_array)
	poly_data.SetPoints(vtk_points)
	poly_data.SetLines(cell_array)
	poly_data.BuildCells()

	if filename.endswith('.xml') or filename.endswith('.vtp'):
		writer = vtk.vtkXMLPolyDataWriter()
		writer.SetDataModeToBinary()
	else:
		writer = vtk.vtkPolyDataWriter()
		writer.SetFileTypeToBinary()


	writer.SetFileName(filename)
	if hasattr(vtk, 'VTK_MAJOR_VERSION') and vtk.VTK_MAJOR_VERSION > 5:
		writer.SetInputData(poly_data)
	else:
		writer.SetInput(poly_data)
	writer.Write() 

def create_pointcloud_polydata(points, colors=None):
    """Creates a vtkPolyData object with the point cloud from numpy arrays
    
    points: numpy.ndarray
        pointcloud with shape (n,3)
    
    colors: numpy.ndarray
        uint8 array with colors for each point. shape is (n,3)

    Returns vtkPolyData object
    """
    import vtk
    vpoints = vtk.vtkPoints()
    vpoints.SetNumberOfPoints(points.shape[0])
    for i in range(points.shape[0]):
        vpoints.SetPoint(i, points[i])
    vpoly = vtk.vtkPolyData()
    vpoly.SetPoints(vpoints)
    
    if not colors is None:
        vcolors = vtk.vtkUnsignedCharArray()
        vcolors.SetNumberOfComponents(3)
        vcolors.SetName("Colors")
        vcolors.SetNumberOfTuples(points.shape[0])
        for i in range(points.shape[0]):
            vcolors.SetTuple3(i ,colors[i,0],colors[i,1], colors[i,2])
        vpoly.GetPointData().SetScalars(vcolors)

    vcells = vtk.vtkCellArray()
    
    for i in range(points.shape[0]):
        vcells.InsertNextCell(1)
        vcells.InsertCellPoint(i)
        
    vpoly.SetVerts(vcells)
    
    return vpoly 

def __init__(self, cells=None, n_cells=None, deep=False):
        """Initialize a vtkCellArray."""
        if cells is not None:
            self._set_cells(cells, n_cells, deep)