def plot_element_values(n, nodes, values, resample_n=None,
        node_tuples=None, show_nodes=False):
    dims = len(nodes)

    orig_nodes = nodes
    orig_values = values

    if resample_n is not None:
        import modepy as mp
        basis = mp.simplex_onb(dims, n)
        fine_nodes = mp.equidistant_nodes(dims, resample_n)

        values = np.dot(mp.resampling_matrix(basis, fine_nodes, nodes), values)
        nodes = fine_nodes
        n = resample_n

    from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \
            as gnitstam

    if dims == 1:
        import matplotlib.pyplot as pt
        pt.plot(nodes[0], values)
        if show_nodes:
            pt.plot(orig_nodes[0], orig_values, "x")
        pt.show()
    elif dims == 2:
        import mayavi.mlab as mlab
        mlab.triangular_mesh(
                nodes[0], nodes[1], values, submesh(list(gnitstam(n, 2))))
        if show_nodes:
            mlab.points3d(orig_nodes[0], orig_nodes[1], orig_values,
                    scale_factor=0.05)
        mlab.show()
    else:
        raise RuntimeError("unsupported dimensionality %d" % dims)

 def plotMesh(
     cls, coordinates, triangles=None, values=None, axes=True, displacement=False, newfigure=True, scale=1, **kwargs
 ):
     # http://docs.enthought.com/mayavi/mayavi/auto/mlab_helper_functions.html#mayavi.mlab.triangular_mesh
     if newfigure:
         cls.figure()
     if isinstance(coordinates, Function) and triangles is None:
         coordinates, triangles, values = cls._function_data(coordinates)
     else:
         assert triangles is not None
     x = coordinates[:, 0]
     y = coordinates[:, 1]
     representation = "surface"
     scalars = None
     if displacement:
         representation = "wireframe"
         assert values.shape[1] == 2
         x = coordinates[:, 0] + scale * values[:, 0]
         y = coordinates[:, 1] + scale * values[:, 1]
         scalars = sqrt(sum(values * values, axis=1))
     else:
         assert values is None or len(values.shape) == 1
     if values is None or displacement:
         values = 0 * x
     mlab.triangular_mesh(x, y, values, triangles, representation=representation, scalars=scalars, **kwargs)
     if axes:
         cls.axes()

def plot3dformesh(x,cv,f):
    return
    cv = band.toZeroStatic(cv)
    if MPI.COMM_WORLD.rank == 0:
        v2 = cv.getArray()
        pl.figure(bgcolor=(1,1,1),fgcolor=(0.5,0.5,0.5))
        pl.triangular_mesh(x[:,0],x[:,1],x[:,2],f,scalars=v2)

def show_contrasts(subject, contrasts, side, threshold):

    x, y, z, triangles = get_geometry(subject, side, "inflated")   ## inflated or white
    curv = get_curvature_sign(subject, side)

    f = mlab.figure()
    mlab.clf()

    # anatomical mesh
    mlab.triangular_mesh(x, y, z, triangles, transparent=False,
                         opacity=1., name=subject,
        scalars=curv, colormap="bone", vmin=-1, vmax=2)

    mlab.title(subject)

    cmaps = [colormaps[c.split("-")[0]]['colormap'] for c in contrasts]

    for contrast, colormap in zip(contrasts, cmaps):
        # functional mesh
        data = get_contrast(subject, contrast, side)
        func_mesh = mlab.pipeline.triangular_mesh_source(x, y, z, triangles,
                                                     scalars=data)
            # threshold
        thresh = mlab.pipeline.threshold(func_mesh, low=threshold)

        surf = mlab.pipeline.surface(thresh, colormap='hot', transparent=True,
                          opacity=.8) # diminuer pour avoir plus de transparence
        lut = (np.array([colormap(v) for v in np.linspace(.25, 1., 256)]) * 255
                       ).astype(int)

        surf.module_manager.scalar_lut_manager.lut.table = lut

    mlab.draw()

    return f

def display(fun_dir, fs_dir, contrast):
    mlab.figure(bgcolor=(1, 1, 1))
    left_mesh = freesurfer.read_geometry(os.path.join(fs_dir, 'lh.inflated'))
    right_mesh = freesurfer.read_geometry(os.path.join(fs_dir, 'rh.inflated'))
    left_curv = os.path.join(fs_dir, 'lh.curv')
    right_curv = os.path.join(fs_dir, 'rh.curv')
    meshes = [left_mesh, right_mesh]
    curves = [left_curv, right_curv]
    for hemisphere, mesh_file, curv_file in zip(['lh', 'rh'], meshes, curves):
        fun_file = os.path.join(fun_dir, '%s_z_map_%s.gii' % (
                contrast, hemisphere))
        coords, triangles = mesh_file
        x, y, z = coords.T

        if hemisphere == 'lh':
            x -= 50
        else:
            x += 50

        curv = freesurfer.read_morph_data(curv_file).astype(np.float)
        tex = np.array([darrays.data for darrays in 
                        read(fun_file).darrays]).ravel()
        print fun_file, tex.min(), tex.max()
        name = ''
        cmin = -1
        cmax = 1
        mlab.triangular_mesh(x, y, z, triangles, transparent=True, opacity=1.,
                             name=name, scalars=curv, colormap="bone",
                             vmin=cmin, vmax=cmax)
        func_mesh = mlab.pipeline.triangular_mesh_source(
            x, y, z, triangles, scalars=tex)
        thresh = mlab.pipeline.threshold(func_mesh, low=THRESHOLD)
        mlab.pipeline.surface(thresh, colormap="hot", vmin=THRESHOLD, vmax=7)

def display_collada(dae_file):
    """
    Display the DAE mesh. Requires pycollada.
    """
    print("Displaying %s" % dae_file)

    from collada import Collada, DaeUnsupportedError, DaeBrokenRefError
    mesh = Collada(dae_file, ignore=[DaeUnsupportedError, DaeBrokenRefError])

    for geometry in mesh.scene.objects('geometry'):
        for prim in geometry.primitives():
            # use primitive-specific ways to get triangles
            prim_type = type(prim).__name__
            if prim_type == 'BoundTriangleSet':
                triangles = prim
            elif prim_type == 'BoundPolylist':
                triangles = prim.triangleset()
            else:
                # Unsupported mesh type
                triangles = None

            if triangles is not None:
                x = triangles.vertex[:,0]
                y = triangles.vertex[:,1]
                z = triangles.vertex[:,2]

                mlab.triangular_mesh(x, y, z, triangles.vertex_index,
                                     color=(0, 0, 1))

def cluster_show(_3D_chan1, _3D_chan2, _3D_labels, n_clusters_, means, std, hulls):
    """

    :param _3D_chan1:
    :param _3D_chan2:
    :param _3D_labels:
    """
    red = (1.0, 0.0, 0.0)
    green = (0.0, 1.0, 0.1)

    mlab.pipeline.volume(mlab.pipeline.scalar_field(_3D_chan1), color=green)
    mlab.pipeline.volume(mlab.pipeline.scalar_field(_3D_chan2), color=red)
    mlab.show()

    cm = get_cmap('gist_rainbow')

    for i in range(0, n_clusters_):
        re_img = np.zeros(_3D_chan2.shape)
        re_img[_3D_labels == i] = _3D_chan2[_3D_labels == i]
        mlab.pipeline.volume(mlab.pipeline.scalar_field(re_img), color=tuple(cm(1.*i/n_clusters_)[:-1]))

    # mean_arr = np.zeros((n_clusters_, 3))
    # std_arr = np.zeros((n_clusters_))
    # for i in range(0, n_clusters_):
    #     mean_arr[i, :] = means[i]
    #     std_arr[i] = std[i]
    # x,y,z = mean_arr.T
    # mlab.points3d(x,y,z, std_arr)

    for hull in hulls:
        x,y,z = hull.points.T
        triangles = hull.simplices
        mlab.triangular_mesh(x, y, z, triangles, representation='wireframe', color=(0, 0, 0))

    mlab.show()

    def draw(self, test=None, u=None):
        """Draw all tetrahedra."""
        if test is not None:
            xs = 1./4.*(self.p[0, self.t[0, :]] +
                        self.p[0, self.t[1, :]] +
                        self.p[0, self.t[2, :]] +
                        self.p[0, self.t[3, :]])
            ys = 1./4.*(self.p[1, self.t[0, :]] +
                        self.p[1, self.t[1, :]] +
                        self.p[1, self.t[2, :]] +
                        self.p[1, self.t[3, :]])
            zs = 1./4.*(self.p[2, self.t[0, :]] +
                        self.p[2, self.t[1, :]] +
                        self.p[2, self.t[2, :]] +
                        self.p[2, self.t[3, :]])
            tset = np.nonzero(test(xs, ys, zs))[0]
        else:
            tset = range(self.t.shape[1])

        fset = np.unique(self.t2f[:, tset].flatten())

        if u is None:
            u = self.p[2, :]

        if OPT_MAYAVI:
            mlab.triangular_mesh(self.p[0, :], self.p[1, :], self.p[2, :],
                                 self.facets[:, fset].T, scalars=u)
            mlab.triangular_mesh(self.p[0, :], self.p[1, :], self.p[2, :],
                                 self.facets[:, fset].T,
                                 representation='wireframe', color=(0, 0, 0))
        else:
            raise ImportError("Mayavi not supported "
                              "by the host system!")

def plot(ply):
    '''
    Plot vertices and triangles from a PlyData instance. Assumptions:
        `ply' has a 'vertex' element with 'x', 'y', and 'z'
            properties;

        `ply' has a 'face' element with an integral list property
            'vertex_indices', all of whose elements have length 3.

    '''
    vertex = ply['vertex']

    (x, y, z) = (vertex[t] for t in ('x', 'y', 'z'))

    mlab.points3d(x, y, z, color=(1, 1, 1), mode='point')

    if 'face' in ply:
        tri_idx = ply['face']['vertex_indices']
        idx_dtype = tri_idx[0].dtype

        triangles = numpy.fromiter(tri_idx, [('data', idx_dtype, (3,))],
                                   count=len(tri_idx))['data']

        mlab.triangular_mesh(x, y, z, triangles,
                             color=(1, 0, 0.4), opacity=0.5)

def plot_inverse_operator(subject, fnout_img=None, verbose=None):
    """"Reads in and plots the inverse solution."""

    # get name of inverse solution-file
    subjects_dir = os.environ.get('SUBJECTS_DIR')
    fname_dir_inv = os.path.join(subjects_dir,
                                 subject)

    for files in os.listdir(fname_dir_inv):
        if files.endswith(",cleaned_epochs_avg-7-src-meg-inv.fif"):
            fname_inv = os.path.join(fname_dir_inv,
                                     files)
            break

    try:
        fname_inv
    except NameError:
        print "ERROR: No forward solution found!"
        sys.exit()


    # read inverse solution
    inv = min_norm.read_inverse_operator(fname_inv)

    # print some information if desired
    if verbose is not None:
        print "Method: %s" % inv['methods']
        print "fMRI prior: %s" % inv['fmri_prior']
        print "Number of sources: %s" % inv['nsource']
        print "Number of channels: %s" % inv['nchan']


    # show 3D source space
    lh_points = inv['src'][0]['rr']
    lh_faces = inv['src'][0]['tris']
    rh_points = inv['src'][1]['rr']
    rh_faces = inv['src'][1]['tris']

    # create figure and plot results
    fig_inverse = mlab.figure(size=(1200, 1200),
                              bgcolor=(0, 0, 0))

    mlab.triangular_mesh(lh_points[:, 0],
                         lh_points[:, 1],
                         lh_points[:, 2],
                         lh_faces)

    mlab.triangular_mesh(rh_points[:, 0],
                         rh_points[:, 1],
                         rh_points[:, 2],
                         rh_faces)

    # save result
    if fnout_img is not None:
        mlab.savefig(fnout_img,
                     figure=fig_inverse,
                     size=(1200, 1200))

    mlab.close(all=True)

def plot(ply):
    '''
    Plot vertices and triangles from a PlyData instance. Assumptions:
        `ply' has a 'vertex' element with 'x', 'y', and 'z'
            properties;

        `ply' has a 'face' element with an integral list property
            'vertex_indices', all of whose elements have length 3.

    '''
    vertex = ply['vertex']

    (x, y, z) = (vertex[t] for t in ('x', 'y', 'z'))

    # mlab.points3d(x, y, z, color=(1, 1, 1), mode='point')

    if 'face' in ply:
        tri_idx = ply['face']['vertex_indices']
        print type(tri_idx)
        idx_dtype = tri_idx[0].dtype
        print "idx_dtype" , idx_dtype
        triangles = numpy.fromiter(tri_idx, [('data', idx_dtype, (3,))],
                                   count=len(tri_idx))['data']
        k = triangles.byteswap()
        # p.from_array(numpy.array([[1,2,3],[3,4,5]], dtype=numpy.float32))
        p.from_array((triangles.byteswap().newbyteorder().astype('<f4')))
        # p1.add_points_from_input_cloud()
        p1.set_input_cloud(p)
        fil = p.make_moving_least_squares()

        # cloud_filtered = fil.filter()
        # fil = p1.make_statistical_outlier_filter()
        fil.set_mean_k(20)
        fil.set_std_dev_mul_thresh(0.1)
        # pointx = triangles[:,0]

        # k = pcl.PassThroughFilter(p)
        
        # print k.set_filter_field_name ("x");
        # print k.set_filter_limits(3.0,10.0)
        # print k.filter()

        pc_1 = pcl.PointCloud()
        pc_1 = p
        pc_2 = pcl.PointCloud()
        pc_2 = p
        kd = pcl.KdTreeFLANN(pc_1)
        print('pc_1:')
        print type(triangles)
        # file1 = open("values.txt",'w')
        # file1.write(str(triangles[:]))

        # mlab.mesh(x,y,z)
        print numpy.shape(p)
        print numpy.shape(z)
        print numpy.shape(y)
        print numpy.shape(x)
        # mlab.contour3d(x, y, z, p)
        mlab.triangular_mesh(x, y, z, fil,color=(1, 0.3, 0.9), opacity=0.5,colormap='cool')

 def _render_mesh(self):
     import mayavi.mlab as mlab
     mlab.triangular_mesh(self.points[:, 0],
                          self.points[:, 1],
                          self.points[:, 2],
                          self.trilist,
                          color=(0.5, 0.5, 0.5),
                          figure=self.figure)

 def draw_mesh(self, mesh, color=WHITE, opacity=0.7):
     x, y, z, triangles = self._get_triangular_mesh(mesh)
     mlab.triangular_mesh(x, y, z, triangles, color=color,
                          opacity=opacity, representation='surface')
     #mlab.triangular_mesh(x, y, z, triangles, colormap='bone',
     #                     opacity=0.8, representation='surface')
     mlab.triangular_mesh(x, y, z, triangles, color=(0, 0, 0),
                          line_width=1.0, representation='wireframe')

def plot_big(MH, MHexam, points_proj, points_projal, skel_projal, ViewP, proj_type, plot_type="normal", fnum=1):
  
    # Plotting
    print "Plot figure z-axis reversed (because z=0 is top of scan)"
    mlab.figure(fnum, bgcolor=(1,1,1), size=(800,800))
    mlab.clf()
    #airway plot    
    mlab.triangular_mesh(MH.vertices[:,0], MH.vertices[:,1], MH.vertices[:,2]*(-1), 
                         MH.faces, colormap='Blues', representation='wireframe', 
                         line_width=0.5)
    try:
        mlab.triangular_mesh(MHexam.vertices[:,0], MHexam.vertices[:,1], 
                             MHexam.vertices[:,2]*(-1), MHexam.faces,
                             colormap='Greens', representation='surface',
                             opacity=0.2)
    except:
        print "Example mesh not plotted"
    #airway axes
    #mlab.axes('off') #extent=([-200, 200, -200,200, 200,-200])) 
    #points on the airway
    mlab.points3d(MH.silvert1[:,0], MH.silvert1[:,1], MH.silvert1[:,2]*(-1), 
                  scale_factor=1, color=(0,0,0))
#    #projected points
    mlab.points3d(points_proj[:,0], points_proj[:,1], points_proj[:,2]*(-1), 
                  scale_factor=1, color=(0,0,0))  
    
    if plot_type is "normal":        
#    #alignment points
        mlab.points3d(MH.points[:,0], MH.points[:,1], MH.points[:,2]*(-1), 
                      scale_factor=2, color=(0,1,0))
        
        #skeleton
        mlab.points3d(MH.skel[:,0], MH.skel[:,1], MH.skel[:,2]*(-1), 
                      scale_factor=0.5, color=(0.5,0.5,1))
    #    #alignment points
        mlab.points3d(points_projal[:,0], points_projal[:,1], points_projal[:,2]*(-1), 
                      scale_factor=2, color=(0,1,0))
        
    try:
        mlab.points3d(skel_projal[:,0], skel_projal[:,1], skel_projal[:,2]*(-1), 
                      scale_factor=0.5, color=(0.5,0.5,1))
    except:
        print "skeleton not plotted"
    #view direction
    #if proj_type=="Lodox":
    #    ViewL=np.array([ViewP, ViewP])
    #    ViewL[0,2]=30; ViewL[1,2]=-30        
    #    
    #    mlab.plot3d(ViewL[:,0]/10, ViewL[:,1], ViewL[:,2]*(-1), 
    #                color=(0,0,0), 
    #                tube_radius=1, opacity=1)
    
    #elif proj_type=="normal":
    #    mlab.points3d(ViewP[0]/10, ViewP[1], ViewP[2]*(-1), 
    #                  scale_factor=4, color=(1,0,0))
    
    mlab.view(elevation=80,azimuth=20)

    def plot_beam(self,Z=1,xrot=None,yrot=None,zrot=None): #modified from
        n = BEAMN
        t = np.linspace(-np.pi, np.pi, n)
        z = np.exp(1j * t)
        x = z.real.copy()*BEAMR
        y = z.imag.copy()*BEAMR
        z = np.zeros_like(x)

        triangles = [(0, i, i + 1) for i in range(1, n)]
        x = np.r_[0, x]
        y = np.r_[0, y]
        z = np.r_[Z, z] - Z
        xcap = np.copy(x)
        ycap = np.copy(y)
        zcap = np.zeros_like(x)-Z
        t = np.r_[0, t]

        if xrot != 0.0 or yrot != 0.0 or zrot != 0.0:
            pts = np.matrix(np.vstack((x,y,z)))
            R = np.identity(3)
            if xrot != 0.0:
                c = cos(xrot)
                s = sin(xrot)
                R = np.matrix([[1,0,0],[0,c,-s],[0,s,c]])*R
            if yrot != 0.0:
                c = cos(yrot)
                s = sin(yrot)
                R = np.matrix([[c,0,s],[0,1,0],[-s,0,c]])*R
            if zrot != 0.0:
                c = cos(zrot)
                s = sin(zrot)
                R = np.matrix([[c,-s,0],[s,c,0],[0,0,1]])*R
            



            for i in range(len(x)):
                x[i],y[i],z[i] = R*pts[:,i]
            #x,y,z = R*pts
            #x = np.array(x).flatten()
            #y = np.array(y).flatten()
            #z = np.array(z).flatten()
            #print R*pts
            #print np.shape(np.array(x).flatten()),y,z
            #raise SystemExit
            ptscap = np.matrix(np.vstack((xcap,ycap,zcap)))
            for i in range(len(x)):
                xcap[i],ycap[i],zcap[i] = R*ptscap[:,i]
            #xcap,ycap,zcap = R*ptscap
            #xcap = np.array(xcap).flatten()
            #ycap = np.array(ycap).flatten()
            #zcap = np.array(zcap).flatten()
        

        beam = mlab.triangular_mesh(x+self.xp, y+self.yp, z+self.zp, triangles,color=COLOR_BEAM)
        cap = mlab.triangular_mesh(xcap+self.xp, ycap+self.yp, zcap+self.zp, triangles,color=COLOR_BEAM)
        return beam,cap

    def plot_mesh(self,):
        '''
        Plot the mesh that has been created

        SRH: 12 July2013
        '''
        from mayavi import mlab
        mlab.triangular_mesh(self.vertices[:,0], self.vertices[:,1], self.vertices[:,2],self.faces)
        mlab.triangular_mesh(self.vertices[:,0], self.vertices[:,1], self.vertices[:,2],self.faces,representation='wireframe',color=(0,0,0))

def plotc(cs,before=None,representation='wireframe',color=(0.5,1,1)):
    mlab.clf()
    if None is not before:
        before()
    mlab.triangular_mesh(*cs.get_triangles(),color=color,representation=representation)
    plot_ale(cs)
    if cs.ale:
        plot_edge(cs.ale[-1])
        plot_edge(cs.ale[-1].succ,color=(1,0,1))

 def plotSurfaces(self, X, T, scalars=None, color=None, rep='surface'):
     
     if self.autoRemove: self.removeSurfaces()
     
     if color==None: color=(1,0,0)
     if scalars==None:
         self.surfaces.append(mlab.triangular_mesh(X[:,0], X[:,1], X[:,2], T, color=color, representation=rep))
     else:
         self.surfaces.append(mlab.triangular_mesh(X[:,0], X[:,1], X[:,2], T, scalars=scalars))

    def show(self, using='maya', fit_ellipse=False, show_hull=False):
        """
        A simple scatter-plot to represent the aggregate - either using mpl
        or mayavi
        """

        if fit_ellipse:

            (center, radii, rotation) = self.fit_ellipse()

            u = np.linspace(0.0, 2.0 * np.pi, 100)
            v = np.linspace(0.0, np.pi, 100)

            # cartesian coordinates that correspond to the spherical angles:
            x = radii[0] * np.outer(np.cos(u), np.sin(v))
            y = radii[1] * np.outer(np.sin(u), np.sin(v))
            z = radii[2] * np.outer(np.ones_like(u), np.cos(v))
            # rotate accordingly
            for i in range(len(x)):
                for j in range(len(x)):
                    [x[i,j],y[i,j],z[i,j]] = np.dot([x[i,j],y[i,j],z[i,j]], rotation) + center

        if show_hull:
            hull = self.chull()
            hull_x = hull.points[:,0]
            hull_y = hull.points[:,1]
            hull_z = hull.points[:,2]

        if using=='mpl':

            from mpl_toolkits.mplot3d import Axes3D
            fig = plt.figure()
            # ax = fig.add_subplot(111, projection='3d'
            ax = Axes3D(fig) 
            # ax.set_aspect('equal')
            h = ax.scatter(self.pos[:,0], self.pos[:,1], self.pos[:,2], s=100.)

            if fit_ellipse:
                ax.plot_wireframe(x, y, z,  rstride=4, cstride=4, color='k', alpha=0.2)
            plt.show()

        elif using=='maya':
            import mayavi.mlab as mlab
            fig = mlab.figure(bgcolor=(1, 1, 1), fgcolor=(0, 0, 0))
            h = mlab.points3d(self.pos[:,0], self.pos[:,1], self.pos[:,2], self.radius, scale_factor=2, resolution=16)

            if fit_ellipse:
                mlab.mesh(x,y,z, opacity=0.25, color=(1,1,1))

            if show_hull:
                mlab.triangular_mesh(hull_x, hull_y, hull_z, hull.simplices, representation='wireframe', color=(1,1,1))

        else:
            print('ERROR: using= must ne either mpl or maya')

        return h

 def draw_edges(self):
     """Draw all edges in a wireframe representation."""
     # use mayavi
     if OPT_MAYAVI:
         mlab.triangular_mesh(self.p[0, :], self.p[1, :], self.p[2, :],
                              self.facets.T, representation='wireframe',
                              color=(0, 0, 0))
     else:
         raise ImportError("Mayavi not supported "
                           "by the host system!")