Python mlab.figure函数代码示例

    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)