def filter_step(self, unfiltered_file, include_file, exclude_file):
        import shutil

        fiber_basename = self.get_unique_name()
        filtered_temp_file = '%s_filtered.vtp' % fiber_basename
        shutil.copy2(unfiltered_file, filtered_temp_file)        

        if include_file: 
            cmd='slicerFilterFibers.sh --pass 1 %s %s %s' % (include_file, filtered_temp_file, filtered_temp_file)
            exec_cmd(cmd)

        if exclude_file:
            cmd='slicerFilterFibers.sh --nopass 1 %s %s %s' % (exclude_file, filtered_temp_file, filtered_temp_file)
            exec_cmd(cmd)

        vreader = vtkXMLPolyDataReader()
        vreader.SetFileName(filtered_temp_file)
        vreader.Update()
        polydata = vreader.GetOutput()

        vwriter = vtkPolyDataWriter()
        output_file = '%s_filtered.vtk' % fiber_basename
        vwriter.SetFileName(output_file)
        vwriter.SetInput(polydata)
        vwriter.Write()

        return output_file 

def read_polydata(filename):
    """Read whole-brain tractography as vtkPolyData format."""

    if VERBOSE:
        print "Reading in data from", filename, "..."

    basename, extension = os.path.splitext(filename)

    if   (extension == '.vtk'):
        reader = vtk.vtkPolyDataReader()
    elif (extension == '.vtp'):
        reader = vtk.vtkXMLPolyDataReader()
    else:
        print 'Cannot recognize model file format'
        return None

    reader.SetFileName(filename)
    reader.Update()
    outpd = reader.GetOutput()
    del reader
    if VERBOSE:
        print "Done reading in data from", filename
        print "Number of lines found:", outpd.GetNumberOfLines()

    return outpd

def _read_vtp(file_name):
    '''Internal function to read vtp mesh files
    '''

    reader = vtk.vtkXMLPolyDataReader()
    reader.SetFileName(file_name)
    reader.Update()

    mesh_data = PanelMesh(file_name)
    mesh_data.orig_type = 'VTK Polydata (.vtp)'
    readerOut = reader.GetOutput()
    mesh_data.num_faces =    int(readerOut.GetNumberOfCells())
    mesh_data.num_points =   int(readerOut.GetNumberOfPoints())

    for i in xrange(mesh_data.num_points):
        mesh_data.points.append(readerOut.GetPoint(i))
    mesh_data.points = np.array(mesh_data.points)

    for i in xrange(mesh_data.num_faces):
        c = readerOut.GetCell(i)
        numCellPoints = int(c.GetNumberOfPoints())
        idsTemp = []
        for i in xrange(numCellPoints):
            idsTemp.append(int(c.GetPointId(i)))
        mesh_data.faces.append(np.array(idsTemp))
    mesh_data.faces = np.array(mesh_data.faces)


    return mesh_data

def load_polydata(file_name):
    # get file extension (type)
    file_extension = file_name.split(".")[-1].lower()

    # todo better generic load
    if file_extension == "vtk":
        reader = vtk.vtkPolyDataReader()
    elif file_extension == "vtp":
        reader = vtk.vtkPolyDataReader()
    elif file_extension == "fib":
        reader = vtk.vtkPolyDataReader()
    elif file_extension == "ply":
        reader = vtk.vtkPLYReader()
    elif file_extension == "stl":
        reader = vtk.vtkSTLReader()
    elif file_extension == "xml":
        reader = vtk.vtkXMLPolyDataReader()
    elif file_extension == "obj":
        reader = vtk.vtkOBJReader()
        #try:  # try to read as a normal obj
        #    reader = vtk.vtkOBJReader()
        #except:  # than try load a MNI obj format
        #    reader = vtk.vtkMNIObjectReader()
    else:
        raise "polydata " + file_extension + " is not suported"

    reader.SetFileName(file_name)
    reader.Update()
    print file_name, " Mesh ", file_extension, "Loaded"
    return reader.GetOutput()

    def CreateSurfaceFromFile(self, filename):
        scalar = False
        if filename.lower().endswith('.stl'):
            reader = vtk.vtkSTLReader()
        elif filename.lower().endswith('.ply'):
            reader = vtk.vtkPLYReader()
        elif filename.lower().endswith('.obj'):
            reader = vtk.vtkOBJReader()
        elif filename.lower().endswith('.vtp'):
            reader = vtk.vtkXMLPolyDataReader()
            scalar = True
        else:
            wx.MessageBox(_("File format not reconized by InVesalius"), _("Import surface error"))
            return

        if _has_win32api:
            reader.SetFileName(win32api.GetShortPathName(filename).encode(const.FS_ENCODE))
        else:
            reader.SetFileName(filename.encode(const.FS_ENCODE))

        reader.Update()
        polydata = reader.GetOutput()

        if polydata.GetNumberOfPoints() == 0:
            wx.MessageBox(_("InVesalius was not able to import this surface"), _("Import surface error"))
        else:
            name = os.path.splitext(os.path.split(filename)[-1])[0]
            self.CreateSurfaceFromPolydata(polydata, name=name, scalar=scalar)

def pointToCellData():
    
    print "Using angle: " + str(angle)
    atpFiles = glob.glob(str(angle) + '/*.vtp')
    if not atpFiles:
        exit("No atp files found")
        
    if not os.path.exists(str(angle) + '_new'):
        os.makedirs(str(angle) + '_new')
    
    for inputFile in atpFiles: 

        print 'Reading', inputFile
        atpReader = vtk.vtkXMLPolyDataReader()
        atpReader.SetFileName(inputFile)
        atpReader.Update()
    
        atpDataset = atpReader.GetOutput()
        
        pointToCell = vtk.vtkPointDataToCellData()
        pointToCell.SetInputData(atpDataset)
        pointToCell.PassPointDataOn()
        pointToCell.Update()
        
        newArray = pointToCell.GetOutput()
        
        # Remove the point data arrays, they exist as cell data
        newArray.GetPointData().RemoveArray('ATP')
        newArray.GetPointData().RemoveArray('tau')
        newArray.GetCellData().RemoveArray('p')
    
        atpWriter = vtk.vtkXMLPolyDataWriter()
        atpWriter.SetInputData(newArray)
        atpWriter.SetFileName(str(angle) + '_new/' + inputFile[3:])
        atpWriter.Update()

def trim_vtp(fn = "c.vtp", fnout = "t.vtk", s = (1.0,1.0,1.0)):
    if not os.path.exists(fn):
        print "%s does not exists, skip" % (fn,)
        return
    
    # Read vtp
    rd = vtk.vtkXMLPolyDataReader()
    rd.SetFileName(fn)
    rd.Update()

    # the poly data
    pd = vtk.vtkPolyData()
    pd = rd.GetOutput()

    # times by 10.0 (CM -> MM)
    MM = [10.0 for each in range(0,len(s))]
    pd = scale_vtkPolyData(pd, MM)
    # then divides the spacing
    ss = [1/each for each in s]
    pd = scale_vtkPolyData(pd, ss)

    # Write vtk
    wt = vtk.vtkPolyDataWriter()
    wt.SetInputData(pd)
    wt.SetFileName(fnout)
    wt.Write()    

def loadVtk(filename):
    if 'vtp' in filename:
        vreader = vtk.vtkXMLPolyDataReader()
    else:
        vreader = vtk.vtkPolyDataReader()

    vreader.SetFileName(filename)
    vreader.Update()
    polydata = vreader.GetOutput()
    polydata.ReleaseDataFlagOn()

    streamlines = []
    verts = vtk_to_numpy(polydata.GetPoints().GetData())
    scalars = {}

    pointdata = polydata.GetPointData()
    for si in range(pointdata.GetNumberOfArrays()):
        sname =  pointdata.GetArrayName(si)
        scalars[sname] = vtk_to_numpy(pointdata.GetArray(si))

    for i in range(polydata.GetNumberOfCells()):
        pids =  polydata.GetCell(i).GetPointIds()
        ids = [ pids.GetId(p) for p in range(pids.GetNumberOfIds())]
        streamlines.append(ids)

    res = {'points':verts, 'values':scalars, 'streamlines':streamlines}
    return res

    def Read(self, m_forceRead=False):
        """
        Initial step, should be performed before everything else starts, considering to add
        this into the class constructor so that it is read automatically.
        :return:
        """

        # Skip redundant read if it is already done
        if self._IS_READ_FLAG and not m_forceRead:
            return

        if self.filename.split(".")[-1] == "vtp" or self.filename.split(".")[-1] == "vtk":
            m_reader = vtk.vtkXMLPolyDataReader()
        elif self.filename.split(".")[-1] == "stl":
            m_reader = vtk.vtkSTLReader()
        else:
            raise IOError("Input file for arm surface is of incorrect format")
        m_reader.SetFileName(self.filename)
        m_reader.Update()

        m_mapper = vtk.vtkPolyDataMapper()
        m_mapper.SetInputConnection(m_reader.GetOutputPort())

        m_actor = vtk.vtkActor()
        m_actor.SetMapper(m_mapper)

        self._reader = m_reader
        self._renderer.AddActor(m_actor)
        self._data = m_reader.GetOutput()
        self._IS_READ_FLAG = True
        pass

def read_vtk_file(filename):
    rd = vtk.vtkXMLPolyDataReader()
    rd.SetFileName(filename)
    rd.Update()
    data = rd.GetOutput()

    N = data.GetNumberOfPoints()
    x = [];
    y = [];
    z = [];
    itype = [];
    ilagr = [];
    exist = [];
    T = [];

    for i in range(N):
        pt = data.GetPoint(i)
        x.append(pt[0])
        y.append(pt[1])
        z.append(pt[2])
        itype.append( data.GetPointData().GetArray('ITYPE').GetTuple1(i) )
        ilagr.append( data.GetPointData().GetArray('ILAGR').GetTuple1(i) )
        exist.append( data.GetPointData().GetArray('EXIST').GetTuple1(i) )
        T.append( data.GetPointData().GetArray('T').GetTuple1(i) )

    df = pd.DataFrame({'COORX':x, 'COORY':y, 'COORZ':z, 'ITYPE':itype, 'ILAGR':ilagr, 'EXIST':exist, 'T':T })
    return df

def read_vtkPolyData(filename):
    r'''
    Reads a VTKPolyData file and outputs a tracts/tracts_data pair

    Parameters
    ----------
    filename : str
        VTKPolyData filename

    Returns
    -------
    tracts : list of float array N_ix3
        Each element of the list is a tract represented as point array,
        the length of the i-th tract is N_i
    tract_data : dict of <data name>= list of float array of N_ixM
        Each element in the list corresponds to a tract,
        N_i is the length of the i-th tract and M is the
        number of components of that data type.
    '''

    if filename.endswith('xml') or filename.endswith('vtp'):
        polydata_reader = vtk.vtkXMLPolyDataReader()
    else:
        polydata_reader = vtk.vtkPolyDataReader()

    polydata_reader.SetFileName(filename)
    polydata_reader.Update()

    polydata = polydata_reader.GetOutput()

    return vtkPolyData_to_tracts(polydata)

def readPData(
        filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readPData: "+filename+" ***")

    assert (os.path.isfile(filename)), "Wrong filename (\""+filename+"\"). Aborting."

    if ('vtk' in filename):
        pdata_reader = vtk.vtkPolyDataReader()
    elif ('vtp' in filename):
        pdata_reader = vtk.vtkXMLPolyDataReader()
    else:
        assert 0, "File must be .vtk or .vtp. Aborting."

    pdata_reader.SetFileName(filename)
    pdata_reader.Update()
    pdata = pdata_reader.GetOutput()

    myVTK.myPrint(verbose-1, "n_points = "+str(pdata.GetNumberOfPoints()))
    myVTK.myPrint(verbose-1, "n_verts = "+str(pdata.GetNumberOfVerts()))
    myVTK.myPrint(verbose-1, "n_lines = "+str(pdata.GetNumberOfLines()))
    myVTK.myPrint(verbose-1, "n_polys = "+str(pdata.GetNumberOfPolys()))
    myVTK.myPrint(verbose-1, "n_strips = "+str(pdata.GetNumberOfStrips()))

    return pdata

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkXMLPolyDataReader(), 'Reading vtkXMLPolyData.',
         (), ('vtkXMLPolyData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

def Import(filename):
    reader = vtk.vtkXMLPolyDataReader()
    try:
        reader.SetFileName(filename.encode(wx.GetDefaultPyEncoding()))
    except AttributeError:
        reader.SetFileName(filename)
    reader.Update()
    return reader.GetOutput()

def rescaleATPToRange():
    # This is where the data is for testing purposes.
    print "Current working directory:", os.getcwd()

    print 'Reading', inputFile
    atpReader = vtk.vtkXMLPolyDataReader()
    atpReader.SetFileName(inputFile)
    atpReader.Update()

    atpDataset = atpReader.GetOutput()
    atp = atpDataset.GetPointData().GetArray('ATP')

    inMin, inMax = atp.GetRange()

    atpR = vtk.vtkFloatArray()
    atpR.SetName('initialATP')    
    
    for i in range(atp.GetNumberOfTuples()):
        val = atp.GetTuple(i)[0]
        atpR.InsertNextValue(rescale(val, inMin, inMax))

    atpDataset.GetCellData().AddArray(atpR)
    
    atpDataset.GetPointData().RemoveArray('ATP')
    atpDataset.GetPointData().RemoveArray('tau_w')
    atpDataset.GetCellData().RemoveArray('initialJPLC')

    print 'Writing', outputFile
    atpWriter = vtk.vtkXMLPolyDataWriter()
    atpWriter.SetInput(atpDataset)
    atpWriter.SetFileName(outputFile)
    atpWriter.Update()

    print 'Reading', inputSurfaceFile
    surfaceReader = vtk.vtkXMLPolyDataReader()
    surfaceReader.SetFileName(inputSurfaceFile)
    surfaceReader.Update()

    surface = surfaceReader.GetOutput()
    surface.GetCellData().AddArray(atpR)

    print 'Writing', outputSurfaceFile
    surfaceWriter = vtk.vtkXMLPolyDataWriter()
    surfaceWriter.SetInput(surface)
    surfaceWriter.SetFileName(outputSurfaceFile)
    surfaceWriter.Update()

 def ReadVTKXMLSurfaceFile(self):
     if (self.InputFileName == ''):
         self.PrintError('Error: no InputFileName.')
     self.PrintLog('Reading VTK XML surface file.')
     reader = vtk.vtkXMLPolyDataReader()
     reader.SetFileName(self.InputFileName)
     reader.Update()
     self.Surface = reader.GetOutput()

    def CentreLinePolyData(self):
        """Compute centrelines based on the profile, reusing our
        memoed copy or reading from the cache file if possible.
        """
        if (os.path.exists(self.CentreLineFile ) and 
            os.path.getmtime(self.CentreLineFile ) > os.path.getmtime(self.StlFile) and
            os.path.getmtime(self.CentreLineFile ) > os.path.getmtime(self.FileName)):
            # Cached!
            reader = vtk.vtkXMLPolyDataReader()
            reader.SetFileName(self.CentreLineFile)
            reader.Update()
            return reader.GetOutput()

        # Have to compute it
        
        # Read the STL file
        reader = vtk.vtkSTLReader()
        reader.SetFileName(profile.StlFile)
        
        # Find the seed points for centreline calculation
        # Use points one iolet radius back along the normal.
        
        outletPts = []
        def scale(iolet):
            pt = (iolet.Centre - iolet.Radius * iolet.Normal)
            pt = pt / self.LengthUnit
            return pt.magnitude
        
        for iolet in self.Iolets:
            if isinstance(iolet._iolet, Inlet):
                inletPt = scale(iolet)
            else:
                outletPts.append(scale(iolet))
                pass
            continue
        
        srcPts, tgtPts = FindSeeds(reader, inletPt, outletPts)
        
        # Lazy import since it's so slow!
        from vmtk import vtkvmtk
        centreliner = vtkvmtk.vtkvmtkPolyDataCenterlines()
        centreliner.SetInputConnection(reader.GetOutputPort())
    
        centreliner.SetSourceSeedIds(srcPts)
        centreliner.SetTargetSeedIds(tgtPts)
        centreliner.SetRadiusArrayName("radius")
        centreliner.SetCostFunction("1/R")

        cleaner = vtk.vtkCleanPolyData()
        cleaner.SetInputConnection(centreliner.GetOutputPort())
        
        writer = vtk.vtkXMLPolyDataWriter()
        writer.SetInputConnection(cleaner.GetOutputPort())
        writer.SetFileName(self.CentreLineFile)
        writer.Write()
        return cleaner.GetOutput()

 def __read_vtk_mesh(self):
     """
     Import vtk mesh
     """
     reader = vtk.vtkXMLPolyDataReader()
     reader.SetFileName(self.__vtk_mesh)
     reader.Update()
     self.__vtk_model = reader.GetOutput()
     self.__nb_vtk_cells = self.__vtk_model.GetNumberOfCells()
     self.__nb_vtk_points = self.__vtk_model.GetNumberOfPoints()

    def  vtp_reader(filename):
        reader = vtk.vtkXMLPolyDataReader()
        reader.SetFileName(filename)
        reader.Update()

        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInput(reader.GetOutput())

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

def HDF5toVTKCells():

    input_meshes = []

    # Read input meshes.
    for in_file in input_mesh_files[output]:
        reader = vtk.vtkXMLPolyDataReader()
        reader.SetFileName(in_file)
        reader.Update()

        input_meshes += [reader.GetOutput()]
        
        # Only add parent mesh for a tube (non-bifurcation)
        if branches == 1:
            break   

    for time_step in range(args.start, args.end + 1):
        print("Time:", time_step)
        
        append_filter = vtk.vtkAppendFilter()
        
        for branch in range(branches):
            mesh = vtk.vtkPolyData()
            mesh.DeepCopy(input_meshes[branch])

            # The base input h5 filename given the branch and from which writer it came on said branch.
            h5_file_base = base_names[output] + str(time_step) + '_b_' + str(branch + 1) + '_' + 'x' + '.h5'
            print("Processing file", h5_file_base)
            
            # Group all datasets of a branch at a specific time point given
            # the number of writers the data was split into.
            species_array = append_datasets(writers, h5_file_base, "data")

            # Loop through all attirbutes and append them to a new array in the 
            # correct order given the quad to task ratio.            
            for attribute in attributes[output]:
                reordered_array = vtk.vtkDoubleArray()
                reordered_array.SetName(attribute)
                reordered_array.SetNumberOfValues(numCells[output][0] * numCells[output][1] * circQuads * axialQuads)
                
                reorder_species(species_array[attribute], reordered_array, output)
                mesh.GetCellData().AddArray(reordered_array)
               
            append_filter.AddInputData(mesh)
            
        append_filter.Update()

        # Write the result.
        vtu_file = base_names[output] + str(time_step) + '.vtu'
        print("Writing file", os.path.abspath(vtu_file))

        writer = vtk.vtkXMLUnstructuredGridWriter()
        writer.SetFileName(vtu_file)
        writer.SetInputData(append_filter.GetOutput())
        writer.Update()