def NumpyCombinePixelValues(self,inputVolumeNode1, inputVolumeNode2):
    """This method gets Numpy array information from the input volumes and combines 
    the pixel information to give an output numpy array. All input volumes must 
    be the same size"""

    # Print to Slicer CLI
    print('Combining input pixels...'),
    start_time = time.time()

    # Get Image Data for Input Volumes
    imdata1 = inputVolumeNode1.GetImageData()
    imdata2 = inputVolumeNode2.GetImageData()

    # Get Dimensions of First Input (all inputs should match)
    x,y,z   = imdata1.GetDimensions()

    # Get scalar data for all inputs
    scalars1 = imdata1.GetPointData().GetScalars()
    scalars2 = imdata2.GetPointData().GetScalars()

    ## Make Numpy Arrays from Scalar Data
    array1 = numpy_support.vtk_to_numpy(scalars1)
    array2 = numpy_support.vtk_to_numpy(scalars2)

    # Combine Arrays (must divide before summing or will add to more than 256 and wrap around values)
    #outputNumpyarray = array1/2+array2/2
    outputNumpyarray = np.around(27*np.true_divide(array1,array2)) # For Normalization

    # Print to Slicer CLI
    end_time = time.time()
    print('done (%0.2f s)') % float(end_time-start_time)

    return outputNumpyarray

def read_data_and_build_snapshot_matrix(x_min,x_max,y_min,y_max,z_min,z_max,snapshot_dir,file_list,\
						num_snapshots,num_points,num_components,var_name, A):
	reader = vtk.vtkUnstructuredGridReader()
	reader.ReadAllScalarsOn()
	reader.ReadAllVectorsOn()
	extract_region = "false"
	if ( (x_min<x_max) and (y_min<y_max) and (z_min<z_max) ):
		extract_region = "true"
		extract = vtk.vtkExtractUnstructuredGrid()
		extract.SetExtent(x_min, x_max, y_min, y_max, z_min, z_max)
		extract.MergingOn()
	u_temp_read = np.array(np.zeros((num_points,3), dtype=np.float64))
	print '\n   Reading data ...'
	for j in range(0,num_snapshots+1):
		print '      Reading file ', file_list[j].strip(), 'file number ', j, ' of ', num_snapshots
		reader.SetFileName(snapshot_dir + file_list[j].strip())
		reader.Update()
		if (extract_region=="true"):
			extract.SetInput(reader.GetOutput())
			extract.Update()
			u_temp_read = VN.vtk_to_numpy(extract.GetOutput().GetPointData().GetVectors(var_name))
		else:
			u_temp_read = VN.vtk_to_numpy(reader.GetOutput().GetPointData().GetVectors(var_name))
		for k in range(0,num_components):
			A[k*num_points:(k+1)*num_points,j] = u_temp_read[:,k]

def get_normals(polydata):
    #normals = polydata.GetPointData().GetArray("Normals")
    normals = polydata.GetPointData().GetArray("Normals")
    print numpy_support.vtk_to_numpy(normals)
    

    return numpy_support.vtk_to_numpy(normals)

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 _extractRectGridByBounds(vtrObj,boundObj):
    '''
    Function that extracts cell from a rectilinear grid (vtr) using bounds.

    Should be signifacantly faster the extractBounds method.

    '''
    import numpy as np, SimPEG as simpeg, vtk
    import vtk.util.numpy_support as npsup

    bO = boundObj.GetBounds()
    xC = npsup.vtk_to_numpy(vtrObj.GetXCoordinates())
    yC = npsup.vtk_to_numpy(vtrObj.GetYCoordinates())
    zC = npsup.vtk_to_numpy(vtrObj.GetZCoordinates())
    iLT = np.where(xC <= bO[0])[0]
    iL = 0 if iLT.shape[0] == 0 else iLT[-1]
    iUT = np.where(xC >= bO[1])[0]
    iU = len(xC) if iUT.shape[0] == 0 else iUT[0]
    jLT = np.where(yC <= bO[2])[0]
    jL = 0 if jLT.shape[0] == 0 else jLT[-1]
    jUT = np.where(yC >= bO[3])[0]
    jU = len(yC) if jUT.shape[0] == 0 else jUT[0]
    kLT = np.where(zC <= bO[4])[0]
    kL = 0 if kLT.shape[0] == 0 else kLT[-1]
    kUT = np.where(zC >= bO[5])[0]
    kU = len(zC) if kUT.shape[0] == 0 else kUT[0]
    extRect = vtk.vtkExtractRectilinearGrid()
    extRect.SetInputData(vtrObj)
    extRect.SetVOI((iL,iU,jL,jU,kL,kU))
    extRect.Update()
    return extRect

def vtk_ensure_trilist(polydata):
    try:
        import vtk
        from vtk.util.numpy_support import vtk_to_numpy

        trilist = vtk_to_numpy(polydata.GetPolys().GetData())

        # 5 is the triangle type - if we have another type we need to
        # use a vtkTriangleFilter
        c = vtk.vtkCellTypes()
        polydata.GetCellTypes(c)

        if c.GetNumberOfTypes() != 1 or polydata.GetCellType(0) != 5:
            warnings.warn('Non-triangular mesh connectivity was detected - '
                          'this is currently unsupported and thus the '
                          'connectivity is being coerced into a triangular '
                          'mesh. This may have unintended consequences.')
            t_filter = vtk.vtkTriangleFilter()
            t_filter.SetInputData(polydata)
            t_filter.Update()
            trilist = vtk_to_numpy(t_filter.GetOutput().GetPolys().GetData())

        return trilist.reshape([-1, 4])[:, 1:]
    except Exception as e:
        warnings.warn(str(e))
        return None

  def execute(self):
    vessel=self._vessel
    array_v =dict();
            
    for ff in ["ChestRegion"]:
      tmp=vessel.GetPointData().GetArray(ff)
      if isinstance(tmp,vtk.vtkDataArray) == False:
        tmp=vessel.GetFieldData().GetArray(ff)
      array_v[ff]=vtk_to_numpy(tmp)

    xyz_arr=vtk_to_numpy(vessel.GetPoints().GetData())
    
    fig=plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_aspect('equal')
    X=xyz_arr[:,0]
    Y=xyz_arr[:,1]
    Z=xyz_arr[:,2]
    ax.scatter(X,Y,Z,s=1,c=array_v['ChestRegion'],marker='.',cmap=plt.cm.jet,linewidth=0)
    ax.grid(True)
    plt.xlabel('x')
    plt.ylabel('y')
    # Create cubic bounding box to simulate equal aspect ratio
    max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max()
    Xb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][0].flatten() + 0.5*(X.max()+X.min())
    Yb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][1].flatten() + 0.5*(Y.max()+Y.min())
    Zb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][2].flatten() + 0.5*(Z.max()+Z.min())
    # Comment or uncomment following both lines to test the fake bounding box:
    for xb, yb, zb in zip(Xb, Yb, Zb):
      ax.plot([xb], [yb], [zb], 'w')
    ax.view_init(elev=20.,azim=80)
    fig.savefig(self.output_prefix+'_vasculatureQualityControl.png',dpi=180)

def GetRawDICOMData(filenames,fileID):
  print filenames,fileID
  vtkRealDcmReader = vtk.vtkDICOMImageReader()
  vtkRealDcmReader.SetFileName("%s/%s"%(rootdir,filenames[0]) )
  vtkRealDcmReader.Update()
  vtkRealData = vtk.vtkImageCast()
  vtkRealData.SetOutputScalarTypeToFloat()
  vtkRealData.SetInput( vtkRealDcmReader.GetOutput() )
  vtkRealData.Update( )
  real_image = vtkRealData.GetOutput().GetPointData() 
  real_array = vtkNumPy.vtk_to_numpy(real_image.GetArray(0)) 

  vtkImagDcmReader = vtk.vtkDICOMImageReader()
  vtkImagDcmReader.SetFileName("%s/%s"%(rootdir,filenames[1]) )
  vtkImagDcmReader.Update()
  vtkImagData = vtk.vtkImageCast()
  vtkImagData.SetOutputScalarTypeToFloat()
  vtkImagData.SetInput( vtkImagDcmReader.GetOutput() )
  vtkImagData.Update( )
  imag_image = vtkImagData.GetOutput().GetPointData() 
  imag_array = vtkNumPy.vtk_to_numpy(imag_image.GetArray(0)) 

  vtkAppend = vtk.vtkImageAppendComponents()
  vtkAppend.SetInput( 0,vtkRealDcmReader.GetOutput() )
  vtkAppend.SetInput( 1,vtkImagDcmReader.GetOutput() )
  vtkAppend.Update( )

  vtkDcmWriter = vtk.vtkDataSetWriter()
  vtkDcmWriter.SetFileName("rawdata.%04d.vtk" % fileID )
  vtkDcmWriter.SetInput(vtkAppend.GetOutput())
  vtkDcmWriter.Update()

  return (real_array,imag_array)

    def read_vtk(self, data_dir='./data', file_name='separatrices.vtk'):
        """
        Read the separatrices from a vtk file.

        call signature:

            read_vtk(data_dir='./data', file_name='separatrices.vtk')

        Arguments:

        *data_dir*:
            Origin data directory.

        *file_name*:
            Origin file name.
        """

        reader = vtk.vtkUnstructuredGridReader()
        reader.SetFileName(os.path.join(data_dir, file_name))
        reader.Update()
        output = reader.GetOutput()

        # Read the separatrices.
        points = output.GetPoints()
        cells = output.GetCells()
        self.separatrices = []
        self.connectivity = []
        for separatrix in range(points.GetNumberOfPoints()):
            self.separatrices.append(points.GetPoint(separatrix))
        self.separatrices = np.array(self.separatrices)
        self.connectivity = np.array([VN.vtk_to_numpy(cells.GetData())[1::3],
                                      VN.vtk_to_numpy(cells.GetData())[2::3]])
        self.connectivity = self.connectivity.swapaxes(0, 1)

    def PlotData(self):

        # remove any measurement line that might exist
        self.removeMeasurementLine(False)

        if self._usebar is False:
            # handle 'pixel' vs 'mm' options:
            self._transform.Identity()
            self._transform.Scale(self._unit_scalings[self._unit], self._unit_scalings[
                                  self._unit], self._unit_scalings[self._unit])
            self._plotData.Update()
            points = numpy_support.vtk_to_numpy(
                self._plotData.GetPoints().GetData())
            points = (points - points[0]).transpose()
            self.xdata = sum(points * points) ** 0.5
            self.ydata = numpy_support.vtk_to_numpy(
                self._plotData.GetPointData().GetScalars()).astype('float32')
        else:
            inc = abs(self.xdata[1] - self.xdata[0])
            self.lower_panel.m_spinCtrlLower.SetIncrement(inc)
            self.lower_panel.m_spinCtrlUpper.SetIncrement(inc)
            self.lower_panel.m_spinCtrlLower.SetMin(self.xdata[0])
            self.lower_panel.m_spinCtrlLower.SetMax(self.xdata[-1])
            self.lower_panel.m_spinCtrlUpper.SetMin(self.xdata[0])
            self.lower_panel.m_spinCtrlUpper.SetMax(self.xdata[-1])

        self.draw_plot()
        self.canvas.draw()

 def test_pointdata(self):
     self.nrrdArray = ns.vtk_to_numpy(
         self.rnrrd.GetOutput().GetPointData().GetScalars()
     )
     self.itkArray = ns.vtk_to_numpy(
         self.ritk.GetOutput().GetPointData().GetScalars()
     )
     self.assertTrue(numpy.allclose(self.nrrdArray, self.itkArray))

def rdmmag(X1, X2):
    """Rdm and Mag calculation."""
    X1 = vtk_to_numpy(X1)
    X2 = vtk_to_numpy(X2)
    rdm = la.norm(X1 / la.norm(X1) - X2 / la.norm(X2))
    rmag = abs(1. - la.norm(X2) / la.norm(X1))
    print("  RDM = ", rdm, "\t  rMAG = ", rmag)
    return rdm, rmag

def getDataArray(vtkObj,name,arrType='Cell'):
    """Function that returns the cell/point data array. """
    return npsup.vtk_to_numpy(vtkObj.GetCellData().GetArray(name))
    if arrType == 'Cell':
        return npsup.vtk_to_numpy(vtkObj.GetCellData().GetArray(name))
    elif arrType == 'Point':
        return npsup.vtk_to_numpy(vtkObj.GetPointData().GetArray(name))
    else:
        raise Exception('Not a support arrType')

 def __init__(self, dataDir = 'data', streamFile = 'stream.vtk'):
     """
     Read the initial streamlines.
     
     call signature:
     
       readStream(dataDir = 'data', streamFile = 'stream.vtk')
       
     Keyword arguments:
      *dataDir*:
         Data directory.
         
      *streamFile*:
         Read the initial streamline from this file.
     """
 
     # load the data
     reader = vtk.vtkPolyDataReader()
     reader.SetFileName(dataDir + '/' + streamFile)
     reader.Update()
     output = reader.GetOutput()
     
     # get the fields
     field = output.GetFieldData()
     nArrays = field.GetNumberOfArrays()
     class params: pass
     p = params()
     for i in range(nArrays):            
         arrayName = field.GetArrayName(i)
         if any(arrayName == np.array(['l', 'sl'])):
             setattr(self, arrayName, VN.vtk_to_numpy(field.GetArray(arrayName)))
         elif any(arrayName == np.array(['hMin', 'hMax', 'lMax', 'tol', 'iterMax', 'nt'])):
             setattr(self, arrayName, VN.vtk_to_numpy(field.GetArray(arrayName))[0])
         else:
             # change this if parameters can have more than one entry
             setattr(p, arrayName, VN.vtk_to_numpy(field.GetArray(arrayName))[0])
     setattr(self, 'p', p)
     
     # get the points
     points = output.GetPoints()
     pointsData = points.GetData()
     data = VN.vtk_to_numpy(pointsData)
     #data = np.swapaxes(data, 0, 1)
     print self.nt
     print self.sl
     print data.shape
     tracers = np.zeros([self.nt, np.max(self.sl), 3], dtype = data.dtype)
     sl = 0
     for i in range(self.nt):
         #if (i > 0):
             #sl = self.sl[i-1]
         #else:
             #sl = 0
         print sl, self.sl[i]
         tracers[i,:self.sl[i],:] = data[sl:sl+self.sl[i],:]
         sl += self.sl[i]
     setattr(self, 'tracers', tracers)

def read_paths(fold):
    print fold
    fname = fold+'/pathsq.vtk'
    print fname
    npArray = fold+'/pathsq'    # fname = '/home/florian/MEGAsync/calcul/LCS_tractor/data/paths.vtk'
    if os.path.isfile(npArray + '.npy'):
        print 'loading ', npArray
        t0 = time.time()
        xx = np.load(npArray + '.npy')
        print 'already existing array loaded in', time.time() - t0, 's'
    else:
        reader = vtk.vtkPolyDataReader()
        reader.SetFileName(fname)
        reader.Update()
        data = reader.GetOutput()
        ll = data.GetLines()
        n_pts = ll.GetSize()  # nb of points
        n_lines = ll.GetNumberOfCells()  # nb of lines
        idList = vtk.vtkIdList()
        idList.SetNumberOfIds(n_lines)
        print idList
        idList.SetId(0, 0)
        abscissaArray = vtk.vtkFloatArray()
        traj = []  # list of pos
        vel = []  # list of vel
        vort = []  # list of vort
        k = 0
        for i in xrange(n_lines):
            cell = data.GetCell(i)
            abscissa = 0.0
            previousPoint = None
            if i % np.int(n_lines / 100) == 0:
                print k, '% read'
                k += 1


            llength = cell.GetNumberOfPoints()
            pos = []  # np.empty([llength, 3])  # pos, u, vort
            u = []  # np.empty([llength, 3])  # pos, u, vort
            vorti = []  # np.empty([llength])  # pos, u, vort
            for j in range(llength):
                pointId = cell.GetPointId(j)
                pos.append(data.GetPoint(pointId))
                u.append(vtk_to_numpy(data.GetPointData().GetArray("U"))[pointId])
                vorti.append(vtk_to_numpy(data.GetPointData().GetArray("Vorticity"))[pointId])

            traj.append(pos)
            vel.append(u)
            vort.append(vorti)

        # print traj
        # x = [traj, vel, vort]
        xx = np.array([traj,vel,vort])
        np.save(npArray, xx)
        print 'shit\'s read'
        print 'end of path lines reading'
    return xx

def ply_importer(filepath, asset=None, texture_resolver=None, **kwargs):
    """Allows importing Wavefront (OBJ) files.

    Uses VTK.

    Parameters
    ----------
    asset : `object`, optional
        An optional asset that may help with loading. This is unused for this
        implementation.
    texture_resolver : `callable`, optional
        A callable that recieves the mesh filepath and returns a single
        path to the texture to load.
    \**kwargs : `dict`, optional
        Any other keyword arguments.

    Returns
    -------
    shape : :map:`PointCloud` or subclass
        The correct shape for the given inputs.
    """
    import vtk
    from vtk.util.numpy_support import vtk_to_numpy

    ply_importer = vtk.vtkPLYReader()
    ply_importer.SetFileName(str(filepath))

    ply_importer.Update()

    # Get the output
    polydata = ply_importer.GetOutput()

    # We must have point data!
    points = vtk_to_numpy(polydata.GetPoints().GetData()).astype(np.float)

    trilist = np.require(vtk_ensure_trilist(polydata), requirements=['C'])

    texture = None
    if texture_resolver is not None:
        texture_path = texture_resolver(filepath)
        if texture_path is not None and texture_path.exists():
            texture = mio.import_image(texture_path)

    tcoords = None
    if texture is not None:
        try:
            tcoords = vtk_to_numpy(polydata.GetPointData().GetTCoords())
        except Exception:
            pass

        if isinstance(tcoords, np.ndarray) and tcoords.size == 0:
            tcoords = None

    colour_per_vertex = None
    return _construct_shape_type(points, trilist, tcoords, texture,
                                 colour_per_vertex)

 def feature_extractor(self,in_vtk):
   points = vtk_to_numpy(in_vtk.GetPoints().GetData())
   vec= vtk_to_numpy(in_vtk.GetPointData().GetArray("hevec0"))
   #features =np.concatenate((points, vec[:,0:1]),axis=1)
   features = points
   #features=StandardScaler().fit_transform(features)
   pca = PCA(n_components=4)
   pca.fit(features)
   features_t=pca.transform(features)
   return features_t[:,[0,1]]

def get_points_normals_from(polydata):

    nodes_vtk_array = polydata.GetPoints().GetData()
    array = vtk_to_numpy(nodes_vtk_array)

    numberArrays = polydata.GetPointData().GetNumberOfArrays()
    for i in range(numberArrays):
        if polydata.GetPointData().GetArrayName(i) == "Normals":
            array_normals = vtk_to_numpy(polydata.GetPointData().GetScalars("Normals"))
            return array, array_normals
    else:
        raise Exception("No available Normals")

  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 PrepareMesh(self, nclus, subratio, verbose):
     """ Initialize mesh parameters, to include subdividing mesh if necessary"""
     
     # Check if the mesh needs to be subdivided
     if self.mode == 'face':
         n = self.origmesh.GetPolys().GetNumberOfCells()
     else:
         n = self.origmesh.GetPoints().GetNumberOfPoints()
     
     # Required number of points/faces
     nreq = nclus*subratio
     
     if nreq > n:
         
         # Determine the number of subdivisions given that each subdivision generates 4 triangles
         # per single triangle input or two points per triangle
         if self.mode == 'face':
             nsub = int(ceil(log(float(nreq)/n, 4))) # solved given nreq = n*4**nsub
         else:
             nsub = int(ceil(log(float(nreq)/n, 3))) # solved given nreq = n*2**nsub
         
         if verbose:
             print('Subdividing mesh with {:d} subdivision(s)'.format(nsub)); sys.stdout.flush()
         
         # Perform subdivision
         self.mesh = SubdivideMesh(self.origmesh, nsub)
     else:
         self.mesh = self.origmesh
     
     # Extract vertices and faces of mesh
     f = VN.vtk_to_numpy(self.mesh.GetPolys().GetData()).reshape((-1, 4)).astype(np.int32)
     v = VN.vtk_to_numpy(self.mesh.GetPoints().GetData()).astype(np.float)     
     
     # Store faces without padding
     self.fc = f[:, 1:]
     
     
     if self.mode == 'face':
         # Get face areas and centroids
         self.area = GetTriangleArea(v, f)
         self.cent = v.take(self.fc, axis=0).mean(1)
         
     else:
         # Compute the mean of area of the triangles adjcent to a point
         self.area = GetPointWeight(v, f)
         
         # Centers are simply the surface points
         self.cent = v
         
     # Develop face/edge associations or get the points at either edge of unique edges
     self.edgeID = GetEdgeID(self.fc, self.mode)
     self.nclus = nclus