def PlotHorizon3d(tss):
    """
    Plot a list of horizons.

    Parameters
    ----------

    tss : list of trappedsurface
        All the trapped surfaces to visualize.
    """
    from mayavi import mlab
    cmaps = ['bone', 'jet', 'hot', 'cool', 'spring', 'summer', 'winter']
    assert len(cmaps) > len(tss)
    extents = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    for ts, cm in zip(tss, cmaps):
        mlab.mesh(ts.X, ts.Y, ts.Z, colormap=cm, opacity=0.4)
        extents[0] = min(extents[0], np.min(ts.X))
        extents[1] = max(extents[1], np.max(ts.X))
        extents[2] = min(extents[2], np.min(ts.Y))
        extents[3] = max(extents[3], np.max(ts.Y))
        extents[4] = min(extents[4], np.min(ts.Z))
        extents[5] = max(extents[5], np.max(ts.Z))
    mlab.axes(extent=extents)
    mlab.outline(extent=extents)
    mlab.show()

def plot_predicted_labels(points, labels):
    print '[plot_points] Plotting points!'
    xs = np.array([int(point._x) for point in points])
    ys = np.array([int(point._y) for point in points])
    zs = np.array([int(point._z) for point in points])
    mlab.points3d(xs, ys, zs, labels, scale_factor = .4, mode='cube')
    mlab.show()

def draw_min_cost_surface(graph):

	graph = sim.graph

	X = []
	Y = []
	Z = []

	for day,sizes in graph.items():
		for size,days_const in sizes.items():
			days_const, node = sorted(days_const.items())[0]
			X.append(day)
			Y.append(size/10.)
			Z.append(node.min_cost)

	pts = mlab.points3d(X, Y, Z, Z)

	mesh = mlab.pipeline.delaunay2d(pts)

	# Remove the point representation from the plot
	pts.remove()

	# Draw a surface based on the triangulation
	surf = mlab.pipeline.surface(mesh)

	mlab.xlabel("time")
	mlab.ylabel("size")
	mlab.zlabel("cost")
	mlab.show()

def main(hdf5_path):

    with h5py.File(hdf5_path, 'r') as f:
        all_verts = np.array(f.get('all_verts'))
        faces = np.array(f.get('faces'))

    fig = mlab.figure(1, bgcolor=(1, 1, 1))

    @mlab.animate(delay=1000, ui=True)
    def animation():
        for i in count():
            frame = i % all_verts.shape[2]
            verts = all_verts[:, :, frame].T
            mlab.clf()
            mlab.triangular_mesh(
                verts[:, 0],
                verts[:, 1],
                verts[:, 2],
                faces,
                color=(.9, .7, .7))
            fig.scene.z_minus_view()
            mlab.view(azimuth=180)
            mlab.title('mesh %d' % i, size=0.5, height=0, color=(0, 0, 0))
            yield

    a = animation()
    mlab.show()

def vis():

    def vis_mesh(mesh, include_wireframe=False, **kwargs):
        from util3d.mayavi_vis import vis_mesh as vm
        v, f = (np.array(mesh[k]) for k in ('vertices', 'faces'))
        vm(v, f, include_wireframe=include_wireframe, **kwargs)

    example_ids = list(get_example_ids(cat_id, 'eval'))
    random.shuffle(example_ids)

    with all_ds:
        for example_id in example_ids:
            print(example_id)
            image, gt_mesh, predictions = all_ds[example_id]
            meshes = top_k_mesh_fn(
                *(np.array(predictions[k]) for k in ('probs', 'dp')))
            plt.imshow(image)
            mlab.figure()
            vis_mesh(gt_mesh, color=(0, 0, 1))
            for mesh in meshes:
                v, f, ov = (mesh[k] for k in
                            ('vertices', 'faces', 'original_vertices'))
                mlab.figure()
                vis_mesh({'vertices': v, 'faces': f}, color=(0, 1, 0))
                mlab.figure()
                vis_mesh({'vertices': ov, 'faces': f}, color=(1, 0, 0))

            plt.show(block=False)
            mlab.show()
            plt.close()

def k_COV_plots(k_arr, COV_arr):
    sig_max_arr = np.zeros((len(k_arr), len(COV_arr)))
    wmax_arr = np.zeros((len(k_arr), len(COV_arr)))
    mu_r, mu_tau = 0.01, 0.1
    for i, k in enumerate(k_arr):
        for j, cov in enumerate(COV_arr):
            Vf = Vf_k(k)
            loc_r = mu_r * (1 - cov * np.sqrt(3.0))
            scale_r = cov * 2 * np.sqrt(3.0) * mu_r
            loc_tau = mu_tau * (1 - cov * np.sqrt(3.0))
            scale_tau = cov * 2 * np.sqrt(3.0) * mu_tau
            reinf = ContinuousFibers(r=RV('uniform', loc=loc_r, scale=scale_r),
                                  tau=RV('uniform', loc=loc_tau, scale=scale_tau),
                                  xi=WeibullFibers(shape=7.0, sV0=0.003),
                                  V_f=Vf, E_f=200e3, n_int=100)
            ccb_view.model.reinforcement_lst = [reinf]
            sig_max, wmax = ccb_view.sigma_c_max
            sig_max_arr[i, j] = sig_max/Vf
            wmax_arr[i, j] = wmax
    ctrl_vars = orthogonalize([np.arange(len(k_arr)), np.arange(len(COV_arr))])
    print sig_max_arr
    print wmax_arr
    mlab.surf(ctrl_vars[0], ctrl_vars[1], sig_max_arr / np.max(sig_max_arr))
    mlab.surf(ctrl_vars[0], ctrl_vars[1], wmax_arr / np.max(wmax_arr))
    mlab.show()

def showDsweep():
    k_s_sweep = [10**(x-5) for x in range(10)]
    sl_div_diam_sweep = [5.0*(x+1)/1000 for x in range(20)]
    vol_ratio_sweep = [5.0*(x+1)/100 for x in range(10)]    
    
    Dmsd = numpy.load(data_dir + "/D_numpy.npy")
    kDa = 1.660538921e-30;
    mass = 40.0*kDa; 
    viscosity = 8.9e-4; 
    diameter = 5e-9; 
    T = 300.0; 
    Dbase = k_b*T/(3.0*numpy.pi*viscosity*diameter);
    Dmsd = Dmsd/Dbase
    
    mlab.figure(1, size=(800, 800), fgcolor=(1, 1, 1),
                                    bgcolor=(0.5, 0.5, 0.5))
    mlab.clf()
    contours = numpy.arange(0.01,2,0.2).tolist()
    obj = mlab.contour3d(Dmsd,contours=contours,transparent=True,vmin=contours[0],vmax=contours[-1])
    outline = mlab.outline(color=(.7, .7, .7),extent=(0,10,0,20,0,10))
    axes = mlab.axes(outline, color=(.7, .7, .7),
            nb_labels = 5,
            ranges=(k_s_sweep[0], k_s_sweep[-1], sl_div_diam_sweep[0], sl_div_diam_sweep[-1], vol_ratio_sweep[0], vol_ratio_sweep[-1]), 
            xlabel='spring stiffness', 
            ylabel='step length',
            zlabel='volume ratio')
    mlab.colorbar(obj,title='D',nb_labels=5)

    mlab.show()

def make_figures(coor, fun, interp_results, error):
    # Figure of harmoinc function on sphere in fine cordinates
    # Points3d showing interpolation training points coloured to their value
    mlab.figure()
    vmax, vmin = np.max(fun.fine), np.min(fun.fine)
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=fun.fine, vmax=vmax, vmin=vmin)
    mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, fun.coarse,
                  scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)
    mlab.colorbar(title='Spherical Harmonic', orientation='vertical')
    # mlab.savefig('interppointssphere.png')

    # Figure showing results of rbf interpolation
    mlab.figure()
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=interp_results, vmax=vmax, vmin=vmin)
    mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, fun.coarse,
                  scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)
    mlab.colorbar(title='Interpolation', orientation='vertical')
    # mlab.savefig('interpsphere.png')

    mlab.figure()
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=error.errors, vmax=error.max, vmin=-error.max)
    mlab.colorbar(title='Error', orientation='vertical')
    # mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, scalars, scale_factor=0.1, scale_mode='none',vmax=vmax, vmin=vmin)
    # mlab.savefig('interpsphere.png')

    mlab.show()

    def plot_mcontour(self, ndim0, ndim1, z, show_mode):
        "use mayavi.mlab to plot contour."
        if not mayavi_installed:
            self.__logger.info("Mayavi is not installed on your device.")
            return
        #do 2d interpolation
        #get slice object
        s = np.s_[0:ndim0:1, 0:ndim1:1]
        x, y = np.ogrid[s]
        mx, my = np.mgrid[s]
        #use cubic 2d interpolation
        interpfunc = interp2d(x, y, z, kind='cubic')
        newx = np.linspace(0, ndim0, 600)
        newy = np.linspace(0, ndim1, 600)
        newz = interpfunc(newx, newy)
        #mlab
        face = mlab.surf(newx, newy, newz, warp_scale=2)
        mlab.axes(xlabel='x', ylabel='y', zlabel='z')
        mlab.outline(face)
        #save or show
        if show_mode == 'show':
            mlab.show()
        elif show_mode == 'save':
            mlab.savefig('mlab_contour3d.png')
        else:
            raise ValueError('Unrecognized show mode parameter : ' +
                             show_mode)

        return

def test_surface_normals(plot=False, skip_asserts=False,
                         write_reference=False):
    "Test the surface normals of a horseshoe mesh"
    sim = openmodes.Simulation()
    mesh = sim.load_mesh(osp.join(mesh_dir, 'horseshoe_rect.msh'))
    part = sim.place_part(mesh)
    basis = sim.basis_container[part]

    r, rho = basis.integration_points(mesh.nodes, triangle_centres)
    normals = mesh.surface_normals
    r = r.reshape((-1, 3))

    if write_reference:
        write_2d_real(osp.join(reference_dir, 'surface_r.txt'), r)
        write_2d_real(osp.join(reference_dir, 'surface_normals.txt'), normals)

    r_ref = read_2d_real(osp.join(reference_dir, 'surface_r.txt'))
    normals_ref = read_2d_real(osp.join(reference_dir, 'surface_normals.txt'))

    if not skip_asserts:
        assert_allclose(r, r_ref)
        assert_allclose(normals, normals_ref)

    if plot:
        from mayavi import mlab
        mlab.figure()
        mlab.quiver3d(r[:, 0], r[:, 1], r[:, 2],
                      normals[:, 0], normals[:, 1], normals[:, 2],
                      mode='cone')
        mlab.view(distance='auto')
        mlab.show()

def test3():
    import numpy as np
    from mayavi import mlab
    sdf = np.load('/Users/yue/data/segment/sdf_diff_3d_1/imseg_iter_100.npz')['sdf']
    mlab.contour3d(sdf[0], contours=[0])
    # mlab.savefig('/Users/yue/data/segment/sdf_diff_3d_1/imseg_iter_100.png')
    mlab.show()

def main_plot3():
    __import__("matplotlib").rcParams.update({'axes.labelsize': 20,
                                              'axes.titlesize': 20})
    from pylab import plot, show, imshow, figure, colorbar, xlabel, ylabel
    from pylab import legend, title, savefig, close, grid
    from mpl_toolkits.mplot3d import axes3d
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from mayavi import mlab

    names = ['res5.json', 'res24.json', 'res26.json', 'res28.json']
    ps = _load_jsons(names, 'ps')
    ns = _load_jsons(names, 'ns')
    print(len(ps[0]))

    T_fine = r_[0:len(ps) - 1]
    N_fine = r_[:101]
    T_fine, N_fine = np.meshgrid(T_fine, N_fine)

    P_all = array([(array(p[:len(p) // 2]) + p[len(p) // 2:]) for p in ps])
    P_fine = P_all[T_fine, N_fine]
    surf = mlab.mesh(N_fine / (len(N_fine) - 1), P_fine, T_fine / (len(ps) - 1),
                     colormap='blue-red')
    surf.module_manager.scalar_lut_manager.reverse_lut = True

    ax = mlab.axes(xlabel='State (n)', ylabel="Population", zlabel="Time",
                   nb_labels=6, extent=[0, 1, 0, 1, 0, 1],
                   ranges=[0, len(N_fine) - 1, 0, 1, 0, 1])
    ax.label_text_property.font_size = 5

    mlab.outline(surf, color=(.7, .7, .7),
                 extent=[0, 1, 0, 1, 0, 1])
    mlab.show()

 def plot_3D( self ):
     x = self.compute3D_plot[0]
     y = self.compute3D_plot[1]
     z = self.compute3D_plot[2]
     # print x_axis, y_axis, z_axis
     if self.autowarp_bool:
         x = x / x[-1]
         y = y / y[-1]
         z = z / z[-1] * self.z_scale
     mlab.surf( x, y , z, representation = 'wireframe' )
     engine = Engine()
     engine.start()
     if len( engine.scenes ) == 75.5:
         engine.new_scene()
         surface = engine.scenes[0].children[0].children[0].children[0].children[0].children[0]
         surface.actor.mapper.scalar_range = np.array( [ 6.97602671, 8.8533387 ] )
         surface.actor.mapper.scalar_visibility = False
         scene = engine.scenes[0]
         scene.scene.background = ( 1.0, 1.0, 1.0 )
         surface.actor.property.specular_color = ( 0.0, 0.0, 0.0 )
         surface.actor.property.diffuse_color = ( 0.0, 0.0, 0.0 )
         surface.actor.property.ambient_color = ( 0.0, 0.0, 0.0 )
         surface.actor.property.color = ( 0.0, 0.0, 0.0 )
         surface.actor.property.line_width = 1.
         scene.scene.isometric_view()
         mlab.xlabel( self.x_name3D )
         mlab.ylabel( self.y_name3D )
     mlab.outline()
     mlab.show()

 def plot(self):
   """Plots the geometry using the package Mayavi."""
   print('\nPlot the three-dimensional geometry ...')
   from mayavi import mlab
   x_init = self.gather_coordinate('x', position='initial')
   y_init = self.gather_coordinate('y', position='initial')
   z_init = self.gather_coordinate('z', position='initial')
   x = self.gather_coordinate('x')
   y = self.gather_coordinate('y')
   z = self.gather_coordinate('z')
   figure = mlab.figure('body', size=(600, 600))
   figure.scene.disable_render = False
   same = (numpy.allclose(x, x_init, rtol=1.0E-06) and
           numpy.allclose(y, y_init, rtol=1.0E-06) and
           numpy.allclose(z, z_init, rtol=1.0E-06))
   if not same:
     mlab.points3d(x_init, y_init, z_init, name='initial',
                   scale_factor=0.01, color=(0, 0, 1))
   mlab.points3d(x, y, z, name='current',
                 scale_factor=0.01, color=(1, 0, 0))
   mlab.axes()
   mlab.orientation_axes()
   mlab.outline()
   figure.scene.disable_render = True
   mlab.show()

    def surface_parcellation(cortex_boundaries, colouring, mapping_colours, colour_rgb, interaction=False):
        """
        """
        number_of_regions = len(cortex_boundaries.region_neighbours)
        alpha = 255
        lut = numpy.zeros((number_of_regions, 4), dtype=numpy.uint8)
        for k in range(number_of_regions):
            lut[k] = numpy.hstack((colour_rgb[mapping_colours[colouring[k]]], alpha))

        fig = mlab.figure(figure="surface parcellation", bgcolor=(0.0, 0.0, 0.0), fgcolor=(0.5, 0.5, 0.5))
        surf_mesh = mlab.triangular_mesh(cortex_boundaries.cortex.vertices[:, 0],
                                         cortex_boundaries.cortex.vertices[:, 1],
                                         cortex_boundaries.cortex.vertices[:, 2],
                                         cortex_boundaries.cortex.triangles,
                                         scalars=cortex_boundaries.cortex.region_mapping,
                                         figure=fig)
        surf_mesh.module_manager.scalar_lut_manager.lut.number_of_colors = number_of_regions
        surf_mesh.module_manager.scalar_lut_manager.lut.table = lut

        #TODO: can't get region labels to associate with colorbar...
        #mlab.colorbar(object=surf_mesh, orientation="vertical")

        x = cortex_boundaries.boundary[:, 0]
        y = cortex_boundaries.boundary[:, 1]
        z = cortex_boundaries.boundary[:, 2]
        bpts = mlab.points3d(x, y, z, color=(0.25, 0.05, 0.05), scale_factor=1)
        mlab.show(stop=interaction)
        return surf_mesh, bpts

    def plot_field(self, **kwargs):
        "plot scalar field for elf data"
        if not mayavi_installed:
            self.__logger.warning("Mayavi is not installed on your device.")
            return
        # set parameters
        vmin = kwargs['vmin'] if 'vmin' in kwargs else 0.0
        vmax = kwargs['vmax'] if 'vmax' in kwargs else 1.0
        axis_cut = kwargs['axis_cut'] if 'axis_cut' in kwargs else 'z'
        nct = kwargs['nct'] if 'nct' in kwargs else 5
        widths = kwargs['widths'] if 'widths' in kwargs else (1, 1, 1)
        elf_data, grid = self.expand_data(self.elf_data, self.grid, widths)
        #create pipeline
        field = mlab.pipeline.scalar_field(elf_data)  # data source
        mlab.pipeline.volume(field, vmin=vmin, vmax=vmax)  # put data into volumn to visualize
        #cut plane
        if axis_cut in ['Z', 'z']:
            plane_orientation = 'z_axes'
        elif axis_cut in ['Y', 'y']:
            plane_orientation = 'y_axes'
        elif axis_cut in ['X', 'x']:
            plane_orientation = 'x_axes'
        cut = mlab.pipeline.scalar_cut_plane(
            field.children[0], plane_orientation=plane_orientation)
        cut.enable_contours = True  # 开启等值线显示
        cut.contour.number_of_contours = nct
        mlab.show()
        #mlab.savefig('field.png', size=(2000, 2000))

        return

    def surface_pattern(surface, vertex_colours):
        """
        Plot a surface and colour it based on a vector of length number of 
        vertices (vertex_colours).

        * How to obtain a pretty picture (from Mayavi's gui): 
        
          - set surf_mesh color to rgb(237, 217, 221)
          - add a surface module derived from surf_mesh; set 'Actor' 
            representation to wireframe; colour 'gray'.
          - enable contours of scalar_surf  
        """
        #surf_mesh = plot_surface(surface, name="surface pattern")
        fig = mlab.figure(figure="surface pattern", fgcolor=(0.5, 0.5, 0.5))
        surf_mesh = mlab.triangular_mesh(surface.vertices[:, 0],
                                         surface.vertices[:, 1],
                                         surface.vertices[:, 2],
                                         surface.triangles,
                                         figure=fig)
        sm_obj = surf_mesh.mlab_source
        scalar_data = surf_mesh.mlab_source.dataset.point_data
        scalar_data.scalars = vertex_colours
        scalar_data.scalars.name = 'Scalar data'
        scalar_data.update()
        scalar_mesh = mlab.pipeline.set_active_attribute(surf_mesh, point_scalars='Scalar data')
        scalar_surf = mlab.pipeline.surface(scalar_mesh)
        mlab.show(stop=True)
        return sm_obj

def plotSpherical(dataset, vertices, triangles, ptitle="", tsize=0.4, theight=0.95):
    """Plot the spherical data given a data set, triangle set and vertex set.

    The vertex set defines the direction cosines of the individual samples.
    The triangle set defines how the surfrace must be structured between the samples.
    The data set defines, for each direction cosine, the length of the vector.

    Args:
        | dataset(numpy.array(double)): array of data set values
        | vertices(numpy.array([])): array of direction cosine vertices as [x y z]
        | triangles(numpy.array([])): array of triangles as []
        | ptitle(string): title or header for this display
        | tsize(double): title size (units not quite clear)
        | theight(double): title height (y value) (units not quite clear)

    Returns:
        | provides and mlab figure.

    Raises:
        | No exception is raised.
"""

    # calculate a (x,y,z) data set from the direction vectors
    x = dataset * vertices[:, 0]
    y = dataset * vertices[:, 1]
    z = dataset * vertices[:, 2]

    mlab.figure(1, fgcolor=(0, 0, 0), bgcolor=(1, 1, 1))

    # Visualize the points
    pts = mlab.triangular_mesh(x, y, z, triangles)  # z, scale_mode='none', scale_factor=0.2)
    mlab.title(ptitle, size=tsize, height=theight)
    mlab.show()

def vis_voxels(model_id, edge_length_threshold, filled, shuffle=False):
    from mayavi import mlab
    from util3d.mayavi_vis import vis_voxels
    from shapenet.core import cat_desc_to_id
    from template_ffd.inference.voxels import get_voxel_dataset
    from template_ffd.data.voxels import get_gt_voxel_dataset
    from template_ffd.model import load_params
    from template_ffd.data.ids import get_example_ids
    cat_id = cat_desc_to_id(load_params(model_id)['cat_desc'])
    gt_ds = get_gt_voxel_dataset(cat_id, filled)
    inf_ds = get_voxel_dataset(model_id, edge_length_threshold)
    example_ids = get_example_ids(cat_id, 'eval')
    if shuffle:
        example_ids = list(example_ids)
        example_ids.shuffle

    with gt_ds:
        with inf_ds:
            for example_id in example_ids:
                gt = gt_ds[example_id].data
                inf = inf_ds[example_id].data
                vis_voxels(gt, color=(0, 0, 1))
                mlab.figure()
                vis_voxels(inf, color=(0, 1, 0))
                mlab.show()

def make_figures(coor, fun, interp_results, error):
    '''Produce MayaVi figures for interpolation results'''
    mlab.figure()
    vmax, vmin = np.max(fun.fine), np.min(fun.fine)
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=fun.fine, vmax=vmax, vmin=vmin)
    mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, fun.coarse,
                  scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)
    mlab.colorbar(title='Spherical Harmonic', orientation='vertical')
    mlab.savefig('Figures/poorfunctionsphere.png')

    # Figure showing results of rbf interpolation
    mlab.figure()
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=interp_results, vmax=vmax, vmin=vmin)
    mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, fun.coarse,
                  scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)
    mlab.colorbar(title='Interpolation', orientation='vertical')
    mlab.savefig('Figures/poorinterpsphere.png')

    mlab.figure()
    mlab.mesh(coor.fine.x, coor.fine.y, coor.fine.z,
              scalars=error.errors, vmax=error.max, vmin=-error.max)
    mlab.colorbar(title='Error', orientation='vertical')
    # mlab.points3d(coor.coarse.x, coor.coarse.y, coor.coarse.z, scalars, scale_factor=0.1, scale_mode='none',vmax=vmax, vmin=vmin)
    mlab.savefig('Figures/poorerrorsphere.png')

    mlab.show()