def testVisually(self):
        '''blocks selected visually.'''
#        if self.comm.rank == 0:
        g2z,zvec = PETSc.Scatter().toZero(self.bnd.gindBlockWBand)
        g2z.scatter(self.bnd.gindBlockWBand,zvec, PETSc.InsertMode.INSERT)
        x = self.bnd.BlockSub2CenterCarWithoutBand(\
                                                   self.bnd.BlockInd2SubWithoutBand(zvec.getArray()) )
        lx = self.bnd.BlockSub2CenterCarWithoutBand(\
                                                    self.bnd.BlockInd2SubWithoutBand(self.bnd.gindBlockWBand.getArray()))
        try:
            try:
                from mayavi import mlab
            except ImportError:
                from enthought.mayavi import mlab

            if self.comm.rank == 0:
                mlab.figure()
                mlab.points3d(x[:,0],x[:,1],x[:,2])
            mlab.figure()
            mlab.points3d(lx[:,0],lx[:,1],lx[:,2])
            mlab.show()
            #fig.add(pts1)
            #fig.add(pts2)
        except ImportError:
            import pylab as pl
            from mpl_toolkits.mplot3d import Axes3D #@UnusedImport
            fig = pl.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter3D(x[:,0],x[:,1],x[:,2],c='blue',marker='o')
            ax.scatter3D(lx[:,0],lx[:,1],lx[:,2],c='red',marker='D')
            pl.savefig('testVis{0}.png'.format(self.comm.rank))
            pl.show()

    def show_volume(self, bbox):
        """
        show the volume with the given bounding box
        """
        # load the data
        with h5py.File(self.h5fn, 'r') as f:
            seg1 = f["segmentation"]["labels"][bbox[2]:bbox[5],bbox[1]:bbox[4],bbox[0]:bbox[3]]
            raw = f["raw"]["volume"][bbox[2]:bbox[5],bbox[1]:bbox[4],bbox[0]:bbox[3]]
        
        print "Drawing volume"
        t0 = time.time()
        
        #draw everything
        fig1 = mlab.figure(1, size=(500,450))
        mlab.clf(fig1)
        visCell.drawImagePlane(fig1, raw, 'gist_ncar')
        visCell.drawVolumeWithoutReferenceCell(fig1, seg1, np.array((-1,)), (0,0,1),0.5)
        with h5py.File(self.h5fn, 'r') as f:
            visCell.drawLabels(fig1, f, seg1, bbox)

        fig2 = mlab.figure(2, size=(500,450))
        mlab.clf(fig2)

        visCell.draw2DView(fig2, raw[20:-20,:,:], seg1[20:-20,:,:], -1)
        
        t = time.time() - t0
        print "Time for drawing:",t

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 zoncaview(m):
    """
    m is a healpix sky map, such as provided by WMAP or Planck.
    """

    nside = hp.npix2nside(len(m))
    vmin = -1e3; vmax = 1e3

    # Set up some grids:
    xsize = ysize = 1000
    theta = np.linspace(np.pi, 0, ysize)
    phi   = np.linspace(-np.pi, np.pi, xsize)
    longitude = np.radians(np.linspace(-180, 180, xsize))
    latitude = np.radians(np.linspace(-90, 90, ysize))

    # Project the map to a rectangular matrix xsize x ysize:
    PHI, THETA = np.meshgrid(phi, theta)
    grid_pix = hp.ang2pix(nside, THETA, PHI)
    grid_map = m[grid_pix]

    # Create a sphere:
    r = 0.3
    x = r*np.sin(THETA)*np.cos(PHI)
    y = r*np.sin(THETA)*np.sin(PHI)
    z = r*np.cos(THETA)

    # The figure:
    mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
    mlab.clf()

    mlab.mesh(x, y, z, scalars=grid_map, colormap="jet", vmin=vmin, vmax=vmax)

    mlab.draw()

    return

def m2screenshot(mayavi_fig=None, mpl_axes=None, autocrop=True):
    """ Capture a screeshot of the Mayavi figure and display it in the
        matplotlib axes.
    """
    import pylab as pl
    # Late import to avoid triggering wx imports before needed.
    try:
        from mayavi import mlab
    except ImportError:
        # Try out old install of Mayavi, with namespace packages
        from enthought.mayavi import mlab

    if mayavi_fig is None:
        mayavi_fig = mlab.gcf()
    else:
        mlab.figure(mayavi_fig)
    if mpl_axes is not None:
        pl.axes(mpl_axes)

    filename = tempfile.mktemp('.png')
    mlab.savefig(filename, figure=mayavi_fig)
    image3d = pl.imread(filename)
    if autocrop:
        bg_color = mayavi_fig.scene.background
        image3d = autocrop_img(image3d, bg_color)
    pl.imshow(image3d)
    pl.axis('off')
    os.unlink(filename)

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 plot3D(name, X, Y, Z, zlabel):
    """
    Plots a 3d surface plot of Z using the mayavi mlab.mesh function.

    Parameters
    ----------
    name: string
        The name of the figure.
    X: 2d ndarray
        The x-axis data.
    Y: 2d ndarray
        The y-axis data.
    Z: 2d nd array
        The z-axis data.
    zlabel: The title that appears on the z-axis.
    """
    mlab.figure(name)
    mlab.clf()
    plotData = mlab.mesh(X/(np.max(X) - np.min(X)),
                         Y/(np.max(Y) - np.min(Y)),
                         Z/(np.max(Z) - np.min(Z)))
    mlab.outline(plotData)
    mlab.axes(plotData, ranges=[np.min(X), np.max(X),
                                np.min(Y), np.max(Y),
                                np.min(Z), np.max(Z)])
    mlab.xlabel('Space ($x$)')
    mlab.ylabel('Time ($t$)')
    mlab.zlabel(zlabel)

def Main(outputFolder, isData, normalized = True):
    assert isinstance(outputFolder, str)
    min, max = 0., 1.
    if os.path.splitext(args.file)[1] == '.arff':
        datasets, targets, targetMap = loadFromArff(args.file)
    elif os.path.splitext(args.file)[1] == '.npy':
        datasets = numpy.load(args.file)
        min = -1.
    else:
        assert False

    datasets = (lambda x : histogramEqualization(x, min=min, max=max) if normalized else x)(datasets)
    if normalized : assert (datasets.min(), datasets.max()) == (min, max)

    if not os.path.isdir("%s/Pictures" % outputFolder):
        os.makedirs("%s/Pictures" % outputFolder)

    global listIndex
    if listIndex is None or (len(listIndex) >= len(datasets)):
        listIndex = xrange(len(datasets))

    for index in listIndex:
        assert 0 <= index < len(datasets)
        mlab.figure("Index : %d" % index, bgcolor=(1,1,1))
        showArray(datasets[index], isData)
        mlab.savefig("%s/Pictures/Index_%d.png" % (outputFolder, index))
        if isData:
            saveData('%s/Pictures/%s_Index_%d.txt' % (outputFolder, targetMap.reverse_mapping[targets[index]], index), datasets[index])
        else:
            saveData('%s/Pictures/Index_%d.txt' % (outputFolder, index), datasets[index])
        mlab.close()

def short_branches():
    """
    Visualization of short branches of the skeleton.
    
    """
    data1_sk = glob.glob('/backup/yuliya/vsi05/skeletons_largdom/*.h5')
    data1_sk.sort()

    for i,j, k in zip(d[1][37:47], data1_sk[46:56], ell[1][37:47]):
        g = nx.read_gpickle(i)
        dat = tb.openFile(j)
        skel = np.copy(dat.root.skel)
        bra = np.copy(dat.root.branches)
        mask = np.zeros_like(skel)    
        dat.close()
    
        length = nx.get_edge_attributes(g, 'length')
        number = nx.get_edge_attributes(g, 'number')
        num_dict = {}
        for m in number:
            for v in number[m]:
                num_dict.setdefault(v, []).append(m)
        find_br = ndimage.find_objects(bra)
        for l in list(length.keys()):
            if length[l]<0.5*k: #Criteria
                for b in number[l]:
                    mask[find_br[b-1]] = bra[find_br[b-1]]==b
        mlab.figure(bgcolor=(1,1,1), size=(1200,1200))
        mlab.contour3d(skel, colormap='hot')
        mlab.contour3d(mask)
        mlab.savefig('/backup/yuliya/vsi05/skeletons/short_bran/'+ i[42:-10] + '.png')
        mlab.close()

def plot_matrix(connectmat_file, centers_file, threshold_pct=5, weight_edges=False,
                node_scale_factor=2, edge_radius=.5, resolution=8, name_scale_factor=1,
                names_file=None, node_indiv_colors=[], highlight_nodes=[], fliplr=False):
    """
    Given a connectivity matrix and a (x,y,z) centers file for each region, plot the 3D network
    """
    matrix = core.file_reader(connectmat_file)
    nodes = core.file_reader(centers_file)
    if names_file:
        names = core.file_reader(names_file,1)
    num_nodes = len(nodes)
    edge_thresh_pct = threshold_pct / 100.0
    matrix_flat = np.array(matrix).flatten()
    edge_thresh = np.sort(matrix_flat)[len(matrix_flat)-int(len(matrix_flat)*edge_thresh_pct)]
    
    matrix = core.file_reader(connectmat_file)
    ma = np.array(matrix)

    thresh = scipy.stats.scoreatpercentile(ma.ravel(),100-threshold_pct)
    ma_thresh = ma*(ma > thresh)
    
    if highlight_nodes:
        nr = ma.shape[0]
        subset_mat = np.zeros((nr, nr))
        for i in highlight_nodes:
            subset_mat[i,:] = 1
            subset_mat[:,i] = 1
        ma_thresh = ma_thresh * subset_mat
        
    if fliplr:
        new_nodes = []
        for node in nodes:
            new_nodes.append([45-node[0],node[1],node[2]]) # HACK
        nodes = new_nodes
    
    mlab.figure(bgcolor=(1, 1, 1), size=(400, 400))
    for count,(x,y,z) in enumerate(nodes):
        if node_indiv_colors:
            mlab.points3d(x,y,z, color=colors[node_indiv_colors[count]], scale_factor=node_scale_factor, resolution=resolution)
        else:
            mlab.points3d(x,y,z, color=(0,1,0), scale_factor=node_scale_factor, resolution=resolution)
        if names_file:
            width = .025*name_scale_factor*len(names[count])
            print width
            print names[count]
            mlab.text(x, y,names[count], z=z,width=.025*len(names[count]),color=(0,0,0))
    for i in range(num_nodes-1):    
        x0,y0,z0 = nodes[i]
        for j in range(i+1, num_nodes):
            #if matrix[i][j] > edge_thresh:
            if ma_thresh[i][j] > edge_thresh:
                x1,y1,z1 = nodes[j]
                if weight_edges:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=matrix[i][j]/matrix_flat.max(),
                            color=(1,1,1))
                else:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=edge_radius,
                            color=(1,1,1))

def test_normals_new5 ():

#pts1 = clouds.downsample(pts1, 0.02).astype('float64')
    
#     pts1 = np.array([[0.,0.,0.], [0.,0.5,0.], [0.,1.,0.], [1.,0.,0.0], [1.,0.5,0.0], [1.,1.,0.25]])
#     pts2 = np.array([[0.,0.,0.], [0.,0.5,0.], [0.,1.,0.], [1.,0.,1.], [1.,0.5,1.], [1.,1.,1.25]])
#     e1 = np.array([[1.,0.,0.], [1.,0.,0.], [1.,0.,0.], [-1.,0.,0.], [-1.,0.,0.], [-1.,0.,0.]])
#     e2 = np.array([[1.,0.,0.], [1.,0.,0.], [1.,0.,0.], [-1.,0.,0.], [-1.,0.,0.], [-1.,0.,0.]])
    
    pts1, pts2, e1, e2 = create_flap_points_normals(3.0,1,dim=3)
    f1 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, use_cvx=True)
    #f2 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, use_cvx=True)
    f2 = te.tps_eval(pts1, pts2, e1, e2, bend_coef=0.01, rot_coef=1e-5, wt_n=None, nwsize=0.15, delta=0.0001)
    #f2 = te.tps_fit_normals_cvx(pts1, pts2, e1, e2, bend_coef=0.1, rot_coef=1e-5, normal_coef=0.1, wt_n=None, nwsize=0.15, delta=0.0001)
    #f2 = te.tps_fit_normals_exact_cvx(pts1, pts2, e1, e2, bend_coef=0.1, rot_coef=1e-5, normal_coef = 0.1, wt_n=None, nwsize=0.15, delta=0.002)
    
#    import IPython
#    IPython.embed()    
    mlab.figure(1, bgcolor=(0,0,0))
    mayavi_utils.plot_warping(f1, pts1, pts2, fine=False, draw_plinks=False)
    _,f1e2 = te.transformed_normal_direction(pts1, e1, f1, delta=0.0001)#np.asarray([tu.tps_jacobian(f2, pt, 2).dot(nm) for pt,nm in zip(pts1,e1)])
    test_normals_pts(f1.transform_points(pts1), f1e2, wsize=0.15,delta=0.15)
    test_normals_pts(pts2, e2, wsize=0.15,delta=0.15)
    #mlab.show()
    mlab.figure(2,bgcolor=(0,0,0))
    #mlab.clf()
    mayavi_utils.plot_warping(f2, pts1, pts2, fine=False, draw_plinks=False)
    _,f2e2 = te.transformed_normal_direction(pts1, e1, f2, delta=0.0001)#np.asarray([tu.tps_jacobian(f2, pt, 2).dot(nm) for pt,nm in zip(pts1,e1)])
    test_normals_pts(f2.transform_points(pts1), f2e2, wsize=0.15,delta=0.15)
    test_normals_pts(pts2, e2, wsize=0.15,delta=0.15)
    mlab.show()

def test_normals_new4 (n=2,l=0.5,dim=2):

    pts1, pts2, e1, e2 = create_flap_points_normals(n,l,dim)

    delta = 1e-2
    f1 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, use_cvx=True)
    #f1 = te.tps_fit_normals_cvx(pts1, pts2, e1, e2, bend_coef=0.1, rot_coef=1e-5, normal_coef=0.1, wt_n=None, nwsize=0.15, delta=0.0001)
    #f2 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, use_cvx=True)
    #f2 = te.tps_eval(pts1, pts2, e1, e2, bend_coef=0.0, rot_coef=1e-5, wt_n=None, nwsize=0.15, delta=1e-8)
    #f2 = te.tps_fit_normals_cvx(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, normal_coef=10, wt_n=None, nwsize=0.15, delta=1e-6)
    f2 = te.tps_fit_normals_cvx(pts1, pts2, e1, e2, bend_coef=0.0, rot_coef=1e-5, normal_coef=1, wt_n=None, nwsize=0.15, delta=delta)

    mlab.figure(1, bgcolor=(0,0,0))
    mayavi_utils.plot_warping(f1, pts1, pts2, fine=False, draw_plinks=False)
    _,f1e2 = te.transformed_normal_direction(pts1, e1, f1, delta=delta)#np.asarray([tu.tps_jacobian(f2, pt, 2).dot(nm) for pt,nm in zip(pts1,e1)])
    test_normals_pts(np.c_[f1.transform_points(pts1),np.zeros((pts2.shape[0],1))], np.c_[f1e2,np.zeros((f1e2.shape[0],1))], wsize=0.15,delta=0.15)
    test_normals_pts(np.c_[pts2,np.zeros((pts2.shape[0],1))], np.c_[e2,np.zeros((e2.shape[0],1))], wsize=0.15,delta=0.15)
    #mlab.show()
    mlab.figure(2,bgcolor=(0,0,0))
    #mlab.clf()
    mayavi_utils.plot_warping(f2, pts1, pts2, fine=False, draw_plinks=False)
    _,f2e2 = te.transformed_normal_direction(pts1, e1, f2, delta=delta)
    test_normals_pts(np.c_[f2.transform_points(pts1),np.zeros((pts2.shape[0],1))], np.c_[f2e2,np.zeros((f2e2.shape[0],1))], wsize=0.15,delta=0.15)
    test_normals_pts(np.c_[pts2,np.zeros((pts2.shape[0],1))], np.c_[e2,np.zeros((e2.shape[0],1))],  wsize=0.15,delta=0.15)
    mlab.show()

def test_normals_new3 ():
    #pts1 = clouds.downsample(pts1, 0.02).astype('float64')
    
    pts1 = gen_circle_points(0.5, 30)
    pts2 = gen_circle_points_pulled_in(0.5,30,6,0.4)#gen_circle_points(0.5, 30) + np.array([0.1,0.1])
    wt_n = None#np.linalg.norm(pts1-pts2,axis=1)*2+1
    

    f1 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=wt_n, use_cvx=True)
    #f2 = fit_ThinPlateSpline(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, use_cvx=True)
    #f2 = te.tps_eval(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, wt_n=None, nwsize=0.15, delta=0.0001)
    #f2 = te.tps_fit_normals_cvx(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, normal_coef=10, wt_n=wt_n, nwsize=0.15, delta=0.0001)
    #f2 = te.tps_fit_normals_cvx(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, normal_coef=0.1, wt_n=None, nwsize=0.15, delta=0.0001)
    f2 = te.tps_fit_normals_exact_cvx(pts1, pts2, bend_coef=0.1, rot_coef=1e-5, normal_coef = 1, wt_n=None, nwsize=0.15, delta=0.0001)    
    mlab.figure(1, bgcolor=(0,0,0))
    mayavi_utils.plot_warping(f1, pts1, pts2, fine=False, draw_plinks=True)
    test_normals_pts(np.c_[f1.transform_points(pts1),np.zeros((pts2.shape[0],1))], wsize=0.15,delta=0.15)
    test_normals_pts(np.c_[pts2,np.zeros((pts2.shape[0],1))], wsize=0.15,delta=0.15)
    #mlab.show()
    mlab.figure(2,bgcolor=(0,0,0))
    #mlab.clf()
    mayavi_utils.plot_warping(f2, pts1, pts2, fine=False, draw_plinks=True)
    test_normals_pts(np.c_[f2.transform_points(pts1),np.zeros((pts2.shape[0],1))], wsize=0.15,delta=0.15)
    test_normals_pts(np.c_[pts2,np.zeros((pts2.shape[0],1))], wsize=0.15,delta=0.15)
    mlab.show()

def test_base_line2 (pts1, pts2):
    pts1 = clouds.downsample(pts1, 0.02)
    pts2 = clouds.downsample(pts2, 0.02)
    print pts1.shape
    print pts2.shape

    #plotter = PlotterInit()

    def plot_cb(src, targ, xtarg_nd, corr, wt_n, f):
        plot_requests = plot_warping(f.transform_points, src, targ, fine=False)
        for req in plot_requests:
            plotter.request(req)

    f1,_ = tps_rpm_bij(pts1, pts2, reg_init=10, reg_final=1, rot_reg=np.r_[1e-3,1e-3,1e-1], n_iter=50, plot_cb=plot_cb, plotting=0)
    #raw_input("Done with tps_rpm_bij")
    #plotter.request(gen_mlab_request(mlab.clf))
    f2,_ = tps_rpm_bij_normals(pts1, pts2, reg_init=10, reg_final=01, n_iter=50, rot_reg=np.r_[1e-3,1e-3,1e-1], normal_coeff = 0.01,  
                                    nwsize = 0.07, plot_cb=plot_cb, plotting =0)
    #raw_input('abcd')

    from tn_rapprentice import tps
    #print tps.tps_cost(f1.lin_ag, f1.trans_g, f1.w_ng, pts1, pts2, 1)
    #print tps.tps_cost(f2.lin_ag, f2.trans_g, f2.w_ng, pts1, pts2, 1)
    #plotter.request(gen_mlab_request(mlab.clf))
    mlab.figure(1)
    mayavi_utils.plot_warping(f1, pts1, pts2, fine=False, draw_plinks=True)
    #mlab.show()
    mlab.figure(2)
    #mlab.clf()
    mayavi_utils.plot_warping(f2, pts1, pts2, fine=False, draw_plinks=True)
    mlab.show()

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 vis_mesh(model_id, edge_length_threshold, shuffle=False, wireframe=False):
    import numpy as np
    from mayavi import mlab
    from shapenet.core.meshes import get_mesh_dataset
    from shapenet.core import cat_desc_to_id
    from util3d.mayavi_vis import vis_mesh
    from template_ffd.inference.meshes import get_inferred_mesh_dataset
    from template_ffd.model import load_params
    import random

    def vis(mesh, **kwargs):
        v, f = (np.array(mesh[k]) for k in ('vertices', 'faces'))
        vis_mesh(v, f, include_wireframe=wireframe, **kwargs)

    cat_id = cat_desc_to_id(load_params(model_id)['cat_desc'])
    inf_mesh_dataset = get_inferred_mesh_dataset(
        model_id, edge_length_threshold)
    with inf_mesh_dataset:
        with get_mesh_dataset(cat_id) as gt_mesh_dataset:
            example_ids = list(inf_mesh_dataset.keys())
            if shuffle:
                random.shuffle(example_ids)

            for example_id in example_ids:
                inf = inf_mesh_dataset[example_id]
                gt = gt_mesh_dataset[example_id]
                mlab.figure()
                vis(inf, color=(0, 1, 0), opacity=0.2)
                mlab.figure()
                vis(gt, opacity=0.2)
                mlab.show()

def show_sensors(pth, red=None, green=None):
    """ 
    show sensors stored in a h5. 
    :param red: a list of labels to be shown in red    
    """
    if red is None:
        red = []
    if green is None:
        green = []

    mlab.figure()
    sensors = h5py.File(pth)
    d = numpy.array(sensors['locations'])    
    labels = list(sensors['labels'])    
    highlight = numpy.ones(d.shape[0])

    for i, l in enumerate(labels):
        if l in red:
            highlight[i] = 5.0
        elif l in green:
            highlight[i] = 7.0
        else:
            highlight[i] = 2.0    

    mlab.points3d(d[:,0], d[:,1], d[:,2], highlight, scale_mode='none')
    mlab.axes()

def main():
    fileDir='./ascfiles2'
    mapFile = './maps/worldmap.png'
    colorFile = './colormaps/color1.txt'

    # color map used for test
    color1 = [[0,0,0,0,0],
        [1,255,255,0,10],
        [10,255,255,0,50],
        [100,255,100,0,150],
        [1000,255,0,0,255]]
    
    color2 = [[0,0,0,0,0],
        [100,0,255,255,10],
        [1000,0,255,255,50],
        [10000,0,100,255,150],
        [100000,0,0,255,255]]

    s = loadData(fileDir)
    mlab.figure(bgcolor=(0,0,0), size=(800,800) )
    vol = drawVolume(mlab,s)
    colormap = readColorFile(colorFile)
    changeVolumeColormap(vol,colormap)
    drawMap(mlab, s, mapFile)
    changeToBestView(mlab)

def draw_coordinate_system_axes(fig, coordinate_system, offset=0.0, scale=1.0, draw_labels=True):
    points, lengths = coordinate_system_arrows(coordinate_system, offset=offset, scale=scale)
    mlab.figure(fig, bgcolor=fig.scene.background)
    arrows = mlab.quiver3d(
        points[:, 0],
        points[:, 1],
        points[:, 2],
        lengths[0, :],
        lengths[1, :],
        lengths[2, :],
        scalars=np.array([3, 2, 1]),
        mode="arrow",
    )
    arrows.glyph.color_mode = "color_by_scalar"
    arrows.glyph.glyph.scale_factor = scale
    data = arrows.parent.parent
    data.name = coordinate_system.name
    glyph_scale = arrows.glyph.glyph.scale_factor * 1.1
    # label_col = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
    labels = []
    if draw_labels:
        for i in range(3):
            labels.append(
                mlab.text(
                    points[i, 0] + glyph_scale * coordinate_system.basis[i, 0],
                    points[i, 1] + glyph_scale * coordinate_system.basis[i, 1],
                    coordinate_system.labels[i],
                    z=points[i, 2] + glyph_scale * coordinate_system.basis[i, 2],
                    # color=label_col[i],
                    width=0.1 * scale,
                )
            )
    return arrows, labels

def plot_sphere_func(f, grid='Clenshaw-Curtis', theta=None, phi=None, colormap='jet', fignum=0):

    # Note: all grids except Clenshaw-Curtis have holes at the poles

    import matplotlib
    matplotlib.use('WxAgg')
    matplotlib.interactive(True)
    from mayavi import mlab

    if grid == 'Driscoll-Healy':
        b = f.shape[0] / 2
    elif grid == 'Clenshaw-Curtis':
        b = (f.shape[0] - 2) / 2
    elif grid == 'SOFT':
        b = f.shape[0] / 2
    elif grid == 'Gauss-Legendre':
        b = (f.shape[0] - 2) / 2

    if theta is None or phi is None:
        theta, phi = meshgrid(b=b, convention=grid)

    phi = np.r_[phi, phi[0, :][None, :]]
    theta = np.r_[theta, theta[0, :][None, :]]
    f = np.r_[f, f[0, :][None, :]]

    x = np.sin(theta) * np.cos(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(theta)

    mlab.figure(fignum, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 400))
    mlab.clf()
    mlab.mesh(x, y, z, scalars=f, colormap=colormap)

    # mlab.view(90, 70, 6.2, (-1.3, -2.9, 0.25))
    mlab.show()