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 Merge5 ():
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName('test10.vtk')

    reader2 = vtk.vtkPolyDataReader()
    reader2.SetFileName('test12.vtk')

    reader3 = vtk.vtkPolyDataReader()
    reader3.SetFileName('test13.vtk')

    app = vtk.vtkAppendPolyData()

    app.AddInputConnection(reader.GetOutputPort())

    tr = vtk.vtkTransform()
    tr.Translate(-50.916666, -1.083333, 0)
    tr.RotateZ(90)

    tp = vtk.vtkTransformPolyDataFilter()
    tp.SetTransform(tr)
    tp.SetInputConnection(reader2.GetOutputPort())

    app.AddInputConnection(tp.GetOutputPort())

    tr2 = vtk.vtkTransform()
    tr2.Translate(-50.916666, -1.083333, 0)

    tp2 = vtk.vtkTransformPolyDataFilter()
    tp2.SetTransform(tr2)
    tp2.SetInputConnection(reader3.GetOutputPort())

    app.AddInputConnection(tp2.GetOutputPort())

    return app

def test_airway_particles():
  try:
    # Get the path to the this test so that we can reference the test data
    this_dir = os.path.dirname(os.path.realpath(__file__))

    # Set up the inputs to AirwayParticles
    input_ct = this_dir + '/../../../Testing/Data/Input/airwaygauss.nrrd'
    input_mask = this_dir + '/../../../Testing/Data/Input/airwaygauss_mask.nrrd'
      
    #tmp_dir = this_dir + '/../../../Testing/tmp/'
    tmp_dir = tempfile.mkdtemp()
    output_particles = os.path.join(tmp_dir,'airway_particles.vtk')
    print tmp_dir

    max_scale = 6.0
    live_th = 40.0
    seed_th = 30.0
    scale_samples = 5
    down_sample_rate = 1.0
    min_intensity = -1100
    max_intensity = -400

    # Generate the airway particles
    ap = AirwayParticles(input_ct, output_particles, tmp_dir, input_mask,
                         max_scale, live_th, seed_th, scale_samples,
                         down_sample_rate, min_intensity, max_intensity)
    ap.execute()
    
    # Read in the reference data set for comparison
    ref_reader = vtk.vtkPolyDataReader()
    ref_reader.SetFileName(this_dir +
                    '/../../../Testing/Data/Input/airway_particles.vtk')
    ref_reader.Update()

    # Now read in the output data set
    test_reader = vtk.vtkPolyDataReader()
    test_reader.SetFileName(output_particles)
    test_reader.Update()
    
    pm = ParticleMetrics(ref_reader.GetOutput(),
                         test_reader.GetOutput(), 'airway')
      
    assert pm.get_particles_dice() > 0.97, \
        "Airway particle Dice score lower than expected"

  finally:
    #Clear particles cache
    ap._clean_tmp_dir=True
    ap.clean_tmp_dir()
    shutil.rmtree(tmp_dir)

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 __init__(self):
        self._extension_map = {}
        self._all_readers = []
        self._reader = vtk.vtkPolyDataReader()

        # register file types
        self.registerFileTypes()

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 get_num_fibers(filename):
    #reader = vtk.vtkDataSetReader()
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.Update()
    nlines = reader.GetOutput().GetNumberOfLines()
    return nlines

def load_streamlines_poyldata(file_name):
    # get file extension (type)
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(file_name)
    reader.Update()
    line_data = reader.GetOutput()
    return line_data

def vtk2vtp(inputfile, outputfile):
    """Read a vtk polydata and save as vtp."""
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(inputfile)
    reader.Update()
    surface = reader.GetOutput()
    writevtp(surface,outputfile)

    def CanReadFile(self, filename):
        '''
        Static methods for I/O
        @param filename: str 
        @rtype: unsigned int
        '''
        if self.IsMeshExtension(os.path.splitext(filename)[1]):
            return self.FILE_IS_VTK
        if self.IsOBJExtension(os.path.splitext(filename)[1]):
            reader = vtk.vtkOBJReader()
            reader.SetFileName(filename)
            try:
                reader.Update()
            except:
#                del reader
                return 0
#            del reader
            return self.FILE_IS_OBJ
        if not self.IsVtkExtension(os.path.splitext(filename)[1]):
            return 0
        try:
            reader = vtk.vtkPolyDataReader()
            reader.SetFileName(filename)
            if reader.IsFilePolyData():
#                del reader
                return self.FILE_IS_VTK
#            del reader
        except:
            pass
        return 0

def readMeshFile(filename, verbose=False):
    """Read mesh file.
    The input format is determined by file name extension. Degenerate data gets
    removed and polygons get split into triangles to support varios restrictive
    output formats."""

    informat = path.splitext(options.infilename)[1].strip('.')
    # set reader based on filename extension
    if informat=='stl':
        reader = vtk.vtkSTLReader()
    elif informat=='vtk':
        reader = vtk.vtkPolyDataReader()
    elif informat=='obj':
        reader = vtk.vtkMNIObjectReader()
    #elif informat=='tag':
    #    reader = vtk.vtkMNITagPointReader()
    else:
        raise ValueError('cannot read input format' + informat)
    reader.SetFileName(filename)

    # merge duplicate points, and/or remove unused points and/or remove degenerate cells
    clean = vtk.vtkCleanPolyData()
    clean.SetInputConnection(reader.GetOutputPort())

    # convert input polygons and strips to triangles
    triangles = vtk.vtkTriangleFilter()
    triangles.SetInputConnection(clean.GetOutputPort())

    #triangles = reader.GetOutputPort()  # skipping above 'cleaning' doesn't work
    if verbose:
        print "read", filename

    return triangles

    def __init__(self, module_manager):
        """Constructor (initialiser) for the PD reader.

        This is almost standard code for most of the modules making use of
        the FilenameViewModuleMixin mixin.
        """

        # call the constructor in the "base"
        ModuleBase.__init__(self, module_manager)

        # setup necessary VTK objects
	self._reader = vtk.vtkPolyDataReader()

        # ctor for this specific mixin
        FilenameViewModuleMixin.__init__(
            self,
            'Select a filename',
            'VTK data (*.vtk)|*.vtk|All files (*)|*',
            {'vtkPolyDataReader': self._reader})
        
        module_utils.setup_vtk_object_progress(
            self, self._reader,
            'Reading vtk polydata')

        # set up some defaults
        self._config.filename = ''
        self.sync_module_logic_with_config()

def mris2vtk(input_name, output_name):
  """
  Converts a FreeSurfer surface to VTK file format
  """
  # convert surface to VTK format
  check_call(['mris_convert', input_name, output_name])
  # get surface RAS translation
  out = check_output(["mris_info", input_name], stderr=STDOUT)
  m = re.search("c_\(ras\)\s:\s\((-?\d+\.\d+),\s(-?\d+\.\d+),\s(-?\d+\.\d+)\)", out)
  if m is None: raise RuntimeError('Could not find c_(ras) coordinates in mris_info output!')
  tx = float(m.group(1))
  ty = float(m.group(2))
  tz = float(m.group(3))
  # transform vertex positions to scanner RAS of orig.mgz
  reader = vtkPolyDataReader()
  reader.SetFileName(output_name)
  reader.Update()
  surface = reader.GetOutput()
  points = surface.GetPoints()
  for i in range(points.GetNumberOfPoints()):
    x, y, z = points.GetPoint(i)
    points.SetPoint(i, x + tx, y + ty, z + tz)
  surface.SetPoints(points)
  writer = vtkPolyDataWriter()
  writer.SetFileName(output_name)
  writer.SetInput(surface)
  writer.Write()

def vtk2objfcn(vtkfile,objfile,decimFactor):
    import vtk
    import random
    from math import *
    import numpy
    obj = open(objfile,'w')
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(vtkfile)
    reader.Update()
    out = reader.GetOutput()
    verttotal = 0
    #print out.GetNumberOfLines()
    for i in xrange(0,out.GetNumberOfLines()-1,decimFactor):
        bounds = out.GetCell(i).GetBounds()
        for j in range(out.GetCell(i).GetNumberOfPoints()):
            voutstr = "v %7.3f %7.3f %7.3f\n" % (out.GetPoint(out.GetCell(i).GetPointId(j))[0]-90.,out.GetPoint(out.GetCell(i).GetPointId(j))[1]-126.,out.GetPoint(out.GetCell(i).GetPointId(j))[2]-72.)
            obj.write(voutstr)
        loutstr = 'l'
        for k in range(out.GetCell(i).GetNumberOfPoints()-1):
            loutstr += ' %7d' % (k+1+verttotal)
        loutstr += '\n'
        obj.write(loutstr)
        verttotal += out.GetCell(i).GetNumberOfPoints()
    obj.close()
    print('Done!')

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

def convert_mris(input_name, output_name):
  """
  Converts a FreeSurfer surface to VTK file format
  """
  # convert surface to VTK format
  if output_name.endswith('.vtp'): temp_name = output_name[0:-1] + 'k'
  else:                            temp_name = output_name
  if not temp_name.endswith('.vtk'):
    raise RuntimeError('Output file name extension must be either .vtk or .vtp')
  check_call(['mris_convert', input_name, temp_name])
  # get surface RAS translation
  out = check_output(["mris_info", input_name], stderr=STDOUT)
  m = re.search("c_\(ras\)\s:\s\((-?\d+\.\d+),\s(-?\d+\.\d+),\s(-?\d+\.\d+)\)", out)
  if m is None: raise RuntimeError('Could not find c_(ras) coordinates in mris_info output!')
  tx = float(m.group(1))
  ty = float(m.group(2))
  tz = float(m.group(3))
  # transform vertex positions to scanner RAS of orig.mgz
  reader = vtkPolyDataReader()
  reader.SetFileName(temp_name)
  reader.Update()
  surface = reader.GetOutput()
  points = surface.GetPoints()
  for i in range(points.GetNumberOfPoints()):
    x, y, z = points.GetPoint(i)
    points.SetPoint(i, x + tx, y + ty, z + tz)
  surface.SetPoints(points)
  if output_name.endswith('.vtp'): writer = vtkXMLPolyDataWriter()
  else:                            writer = vtkPolyDataWriter()
  writer.SetFileName(output_name)
  writer.SetInput(surface)
  writer.Write()
  if temp_name != output_name:
    remove(temp_name)

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_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 PolyDataReader(self, currentElement):
     reader = vtk.vtkPolyDataReader()
     try:
         reader.SetFileName(os.path.join(self.basedir, currentElement.get('SetFileName')))
     except:
         self.logger.error('  .. <PolyDataReader> failed to SetFileName')
     return reader

def MeshToVolume(Filename):
	reader = vtk.vtkPolyDataReader()
	pol2stenc = vtk.vtkPolyDataToImageStencil()
	imgstenc = vtk.vtkImageStencil()

	reader.SetFileName(os.path.join(subjects_dir,subject_dir,Filename))
	reader.Update()

	ref_mesh = reader.GetOutput()
	ref_volume = vtk.vtkImageData()

	# define output volume dimension
	spacing = (0.5,0.5,0.5)

	ref_volume.SetSpacing(spacing)

	bounds = ref_mesh.GetBounds()

	dim = [math.ceil(bounds[ii*2+1] - bounds[ii*2] / spacing[ii]) for ii in range(0,3)]
	origin = [bounds[ii*2] + spacing[ii] / 2 for ii in range(0,3)]
	extent = (0,dim[0] - 1,0,dim[1] -1 ,0,dim[2]-1)

	ref_volume.SetOrigin(origin)
	ref_volume.SetDimensions(dim)
	ref_volume.SetExtent(extent)

	ref_volume.SetScalarTypeToUnsignedChar()
	ref_volume.AllocateScalars()

	# Fill the image with white voxels
	for i in range(0,ref_volume.GetNumberOfPoints()):
		ref_volume.GetPointData().GetScalars().SetTuple1(i,255)

	print ref_volume.GetNumberOfPoints()

	pol2stenc.SetInput(ref_mesh)

	pol2stenc.SetOutputOrigin(origin)
	pol2stenc.SetOutputSpacing(spacing)
	pol2stenc.SetOutputWholeExtent(ref_volume.GetExtent())

	pol2stenc.Update()

	imgstenc.SetInput(ref_volume)
	imgstenc.SetStencil(pol2stenc.GetOutput())

	imgstenc.ReverseStencilOff()
	imgstenc.SetBackgroundValue(0)
	imgstenc.Update()
	tmp = imgstenc.GetOutput()

	writer = vtk.vtkImageWriter()
	writer.SetFileName('prova.nii.gz')
	writer.SetInput(ref_volume)
	writer.Update()

	out = v2n(tmp.GetPointData().GetScalars())

	return np.reshape(out, (dim[0],dim[1],dim[2]))