def action(u, x, xv, y, yv, t, n):
        #print 'action, t=',t,'\nu=',u, '\Nx=',x, '\Ny=', y
        if plot == 1:
            mesh(xv, yv, u, title='t=%g' %t[n])
            time.sleep(0.2) # pause between frames

        elif plot == 2:
            # mayavi plotting
            mlab.clf()
            extent1 = (0, 20, 0, 20,-2, 2)
            s = mlab.surf(x , y, u, colormap='Blues', warp_scale=5,extent=extent1)
            mlab.axes(s, color=(.7, .7, .7), extent=extent1,
            ranges=(0, 10, 0, 10, -1, 1), xlabel='', ylabel='',
            zlabel='',
            x_axis_visibility=False, z_axis_visibility=False)
            mlab.outline(s, color=(0.7, .7, .7), extent=extent1)
            mlab.text(2, -2.5, '', z=-4, width=0.1)
            mlab.colorbar(object=None, title=None, orientation='horizontal', nb_labels=None, nb_colors=None, label_fmt=None)
            mlab.title('test 1D t=%g' % t[n])

            mlab.view(142, -72, 50)
            f = mlab.gcf()
            camera = f.scene.camera
            camera.yaw(0)
        
        
        if plot > 0:
            path = 'Figures_wave2D'
            time.sleep(0) # pause between frames
            if save_plot and plot != 2:
                filename = '%s/%08d.png' % (path, n)
                savefig(filename)  # time consuming!
            elif save_plot and plot == 2:
                filename = '%s/%08d.png' % (path,n)
                mlab.savefig(filename)  # time consuming!

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 plot_matrix_metric(connectmat_file,centers_file,threshold_pct,grp_metrics=None,node_metric='bc',
                       weight_edges=0,node_scale_factor=2,edge_radius=.5,resolution=8,name_scale_factor=1,names_file=0,
                       red_nodes=None):
    """
    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)]
    
    if grp_metrics: # regional metrics caclulated elsewhere, loaded in
        node_colors = {}     # define colors for each metric
        node_colors['s'] = (1,  0.733,  0)
        node_colors['cc'] = (0.53, 0.81, 0.98)
        node_colors['bc'] = (0.5,  1,  0)
        node_colors['eloc'] = (1,  0 ,  1)
        node_colors['ereg'] = (1,  1,  0)
        
        node_metrics={}
        metrics = np.array(core.file_reader(grp_metrics))
        cols = np.shape(metrics)[1]
        for i in range(cols):
            colmean = np.mean(metrics[:,i])
            colscale = 2 / colmean
            metrics[:,i] = metrics[:,i] * colscale # make node mean size 2
        node_metrics['s'] = metrics[:,0] # strength
        node_metrics['cc'] = metrics[:,1] # clustering coefficient
        node_metrics['bc'] = metrics[:,2] # betweenness centrality
        node_metrics['eloc'] = metrics[:,3] # local efficiency
        node_metrics['ereg'] = metrics[:,4] # regional efficiency
    
    mlab.figure(bgcolor=(1, 1, 1), size=(800, 800))
    for count,(x,y,z) in enumerate(nodes):
        if grp_metrics:
            mlab.points3d(x,y,z, color=node_colors[node_metric],
                          scale_factor=node_metrics[node_metric][count], resolution=resolution)
        else:
            mlab.points3d(x,y,z, color=(0,1,0), scale_factor=node_scale_factor, resolution=resolution)
        if names_file:
            mlab.text(x, y,names[count], z=z,width=.02*name_scale_factor*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:
                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 disp_pts(pts1, pts2, color1=(1,0,0), color2=(0,1,0)):
    figure = mlab.gcf()
    mlab.clf()
    figure.scene.disable_render = True

    pts1_glyphs   = mlab.points3d(pts1[:,0], pts1[:,1], pts1[:,2], color=color1, resolution=20, scale_factor=0.001)
    pts2_glyphs   = mlab.points3d(pts2[:,0], pts2[:,1], pts2[:,2], color=color2, resolution=20, scale_factor=0.001)
    glyph_points1 = pts1_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()
    glyph_points2 = pts2_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()

    dd = 0.001

    outline1 = mlab.outline(pts1_glyphs, line_width=3)
    outline1.outline_mode = 'full'
    p1x, p1y, p1z = pts1[0,:]
    outline1.bounds = (p1x-dd, p1x+dd,
                       p1y-dd, p1y+dd,
                       p1z-dd, p1z+dd)

    pt_id1 = mlab.text(0.8, 0.2, '0 .', width=0.1, color=color1)

    outline2 = mlab.outline(pts2_glyphs, line_width=3)
    outline2.outline_mode = 'full'
    p2x, p2y, p2z = pts2[0,:]
    outline2.bounds = (p2x-dd, p2x+dd,
                       p2y-dd, p2y+dd,
                       p2z-dd, p2z+dd)  
    pt_id2 = mlab.text(0.8, 0.01, '0 .', width=0.1, color=color2)
    
    figure.scene.disable_render = False


    def picker_callback(picker):
        """ Picker callback: this gets called during pick events.
        """
        if picker.actor in pts1_glyphs.actor.actors:
            point_id = picker.point_id/glyph_points1.shape[0]
            if point_id != -1:
                ### show the point id
                pt_id1.text = '%d .'%point_id
                #mlab.title('%d'%point_id)
                x, y, z = pts1[point_id,:]
                outline1.bounds = (x-dd, x+dd,
                                   y-dd, y+dd,
                                   z-dd, z+dd)
        elif picker.actor in pts2_glyphs.actor.actors:
            point_id = picker.point_id/glyph_points2.shape[0]
            if point_id != -1:
                ### show the point id
                pt_id2.text = '%d .'%point_id
                x, y, z = pts2[point_id,:]
                outline2.bounds = (x-dd, x+dd,
                                   y-dd, y+dd,
                                   z-dd, z+dd)


    picker = figure.on_mouse_pick(picker_callback)
    picker.tolerance = dd/2.
    mlab.show()

    def plot_u(u, x, xv, y, yv, t, n):
        """User action function for plotting."""
        if t[n] == 0:
            time.sleep(2)
        if plot_method == 1:
            # Works well with Gnuplot backend, not with Matplotlib
            st.mesh(x, y, u, title='t=%g' % t[n], zlim=[-1,1],
                    caxis=[-1,1])
        elif plot_method == 2:
            # Works well with Gnuplot backend, not with Matplotlib
            st.surfc(xv, yv, u, title='t=%g' % t[n], zlim=[-1, 1],
                  colorbar=True, colormap=st.hot(), caxis=[-1,1],
                  shading='flat')
        elif plot_method == 3:
            print 'Experimental 3D matplotlib...under development...'
            # Probably too slow
            #plt.clf()
            ax = fig.add_subplot(111, projection='3d')
            u_surf = ax.plot_surface(xv, yv, u, alpha=0.3)
            #ax.contourf(xv, yv, u, zdir='z', offset=-100, cmap=cm.coolwarm)
            #ax.set_zlim(-1, 1)
            # Remove old surface before drawing
            if u_surf is not None:
                ax.collections.remove(u_surf)
            plt.draw()
            time.sleep(1)
        elif plot_method == 4:
	    # Mayavi visualization
            mlab.clf()
            extent1 = (0, 20, 0, 20,-2, 2)
            s = mlab.surf(x , y, u,
                          colormap='Blues',
                          warp_scale=5,extent=extent1)
            mlab.axes(s, color=(.7, .7, .7), extent=extent1,
                      ranges=(0, 10, 0, 10, -1, 1),
                      xlabel='', ylabel='', zlabel='',
                      x_axis_visibility=False,
                      z_axis_visibility=False)
            mlab.outline(s, color=(0.7, .7, .7), extent=extent1)
            mlab.text(6, -2.5, '', z=-4, width=0.14)
            mlab.colorbar(object=None, title=None,
                          orientation='horizontal',
                          nb_labels=None, nb_colors=None,
                          label_fmt=None)
            mlab.title('Gaussian t=%g' % t[n])
            mlab.view(142, -72, 50)
            f = mlab.gcf()
            camera = f.scene.camera
            camera.yaw(0)

        if plot_method > 0:
            time.sleep(0) # pause between frames
            if save_plot:
                filename = 'tmp_%04d.png' % n
		if plot_method == 4:
                    mlab.savefig(filename)  # time consuming!
		elif plot_method in (1,2):
                    st.savefig(filename)  # time consuming!

 def test_text(self):
     """ Test the text module.
     """
     data = np.random.random((3, 3, 3))
     src = mlab.pipeline.scalar_field(data)
     # Some smoke testing
     mlab.text(0.1, 0.9, 'foo')
     mlab.text(3, 3, 'foo', z=3)
     mlab.title('foo')
     # Check that specifying 2D positions larger than 1 raises an
     # error
     self.assertRaises(ValueError, mlab.text, 0, 1.5, 'test')

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 show_variable_timecourse(self, var, time_point, start_value, end_value):
        """Show an animation of all the section that have 
        the recorded variable among time"""

        # Getting the new scalar
        new_scalar = self.get_var_data(var, time_point)

        d = self.dataset.point_data.get_array("diameter")
        if len(d) != len(new_scalar):
            message = (
                "ERROR! MISMATCH on the Vector Length. \
            If you assign the new vectors it will not work \
            Diameter length: %s New Scalar length: %s var: %s"
                % (len(d), len(new_scalar), var)
            )
            logger.error(message)
        # ReEnable the rendering
        self.mayavi.visualization.scene.disable_render = True

        self.redraw_color(new_scalar, var)

        if not self.colorbar:
            self.colorbar = mlab.colorbar(orientation="vertical")
            self.timelabel = mlab.text(0.05, 0.05, str(time_point), width=0.05)

        self.colorbar.data_range = [start_value, end_value]
        time = self.manager.groups["t"][time_point]
        self.timelabel.text = str(round(time, 3))

        self.mayavi.visualization.scene.disable_render = False

    def surface_timeseries(surface, data, step=1):
        """
        
        """
        fig = mlab.figure(figure="surface_timeseries", fgcolor=(0.5, 0.5, 0.5))
        #Plot an initial surface and colourbar #TODO: Change to use plot_surface function, see below.
        surf_mesh = mlab.triangular_mesh(surface.vertices[:, 0],
                                         surface.vertices[:, 1],
                                         surface.vertices[:, 2],
                                         surface.triangles,
                                         scalars=data[0, :],
                                         vmin=data.min(), vmax=data.max(),
                                         figure=fig)
        mlab.colorbar(object=surf_mesh, orientation="vertical")

        #Handle for the surface object and figure
        surf = surf_mesh.mlab_source

        #Time #TODO: Make actual time rather than points, where/if possible.
        tpts = data.shape[0]
        time_step = mlab.text(0.85, 0.125, ("0 of %s" % str(tpts)),
                              width=0.0625, color=(1, 1, 1), figure=fig,
                              name="counter")

        #Movie
        k = 0
        while 1:
            if abs(k) >= tpts:
                k = 0
            surf.set(scalars=data[k, :])
            time_step.set(text=("%s of %s" % (str(k), str(tpts))))
            k += step
            yield
        mlab.show(stop=True)

def label_init():
  from mayavi.mlab import points3d
  obj = points3d(0., 0., 0., scale_factor=0.01, color=(0., 0., 0.))
  global actlbl
  actlbl = text(.4, .1, 'odorname', width=.15)
  actlbl.text = ''
  obj.remove()

 def draw_axis_name(self, axis_name):
     global_axis_name = self.controller.get_global_axis_name(axis_name)
     if global_axis_name == 'w':
         width = 0.06
     else:
         width = 0.04
     t = mlab.text(0.01, 0.9, global_axis_name, width=width, color=(1, 0, 0))
     return t

def plot_matrix_path(connectmat_file,centers_file,paths_file,path_num=0,threshold_pct=5,weight_edges=False,
                     node_scale_factor=2,edge_radius=.5,resolution=8,name_scale_factor=1,names_file=0):
    """
    Given a connectivity matrix and a (x,y,z) centers file for each region, plot the 3D network;
    paths file contains columns listing nodes contained in path
    """
    matrix=core.file_reader(connectmat_file)
    nodes=core.file_reader(centers_file)
    paths=zip(*core.file_reader(paths_file))
    path=paths[path_num]
    path_pairs=zip(path[0:len(path)-1],path[1:])
    print path_pairs
    if names_file:
        names=core.file_reader(names_file)
    num_nodes=len(nodes)
    edge_thresh_pct=(float(threshold_pct)/100)
    matrix_flat=np.array(matrix).flatten()
    edge_thresh=np.sort(matrix_flat)[len(matrix_flat)-int(len(matrix_flat)*edge_thresh_pct)]
    mlab.figure(bgcolor=(1, 1, 1), size=(400, 400))
    for count,(x,y,z) in enumerate(nodes):
        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])
            mlab.text(x, y,names[count], z=z,width=.025*name_scale_factor,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:
                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))
    for n1,n2 in path_pairs:
        n1=int(n1)
        n2=int(n2)
        x0,y0,z0=nodes[n1]
        x1,y1,z1=nodes[n2]
        mlab.plot3d( [x0,x1], [y0,y1], [z0,z1],
                    tube_radius=1,
                    color=(0,0,1))

 def plot_u(u, x, xv, y, yv, t, n):
     print " ploting"
     if t[n] == 0:
         time.sleep(2)
     if plot_method == 1:
         mesh(x, y, u, title='t=%g' % t[n], zlim=[-1,1],
              caxis=[-1,1])
     elif plot_method == 2:
         surfc(xv, yv, u, title='t=%g' % t[n], zlim=[-1, 1],
               colorbar=True, colormap=hot(), caxis=[-1,1],
               shading='flat')
     elif plot_method == 3:
         # mayavi plotting
         mlab.clf()
         extent1 = (0, 20, 0, 20,-2, 2)
         s = mlab.surf(x , y, u, colormap='Blues', warp_scale=5,
                      extent=extent1)
         mlab.axes(s, color=(.7, .7, .7), extent=extent1,
         ranges=(0, 10, 0, 10, -1, 1), xlabel='', ylabel='',
         zlabel='',
         x_axis_visibility=False, z_axis_visibility=False)
         mlab.outline(s, color=(0.7, .7, .7), extent=extent1)
         mlab.text(2, -2.5, '', z=-4, width=0.1)
         mlab.colorbar(object=None, title=None, orientation='horizontal',
                       nb_labels=None, nb_colors=None, label_fmt=None)
         mlab.title('t=%g' % t[n])
         mlab.view(142, -72, 50)
         f = mlab.gcf()
         camera = f.scene.camera
         camera.yaw(0)
         #g = mlab.figure()
         #g.scene.movie_maker.record = True
         
     
     if plot_method > 0:
         if not os.path.exists(path):
             os.makedirs(path)
         time.sleep(0) # pause between frames
         if save_plot and plot_method == 3:
             filename = '%s/%08d.png'%(path,n)
             mlab.savefig(filename)  # time consuming!
         elif save_plot and plot_method != 3:
             filename = '%s/%08d.png'%(path,n)
             savefig(filename)  # time consuming!

def do_test(lines, scalars, show=False, txt=""):
    viscid.logger.info('--> ' + txt)
    title = txt + '\n' + "\n".join(textwrap.wrap("scalars = {0}".format(scalars),
                                                 width=50))

    try:
        from viscid.plot import vpyplot as vlt
        from matplotlib import pyplot as plt

        vlt.clf()
        vlt.plot_lines(lines, scalars=scalars)
        plt.title(title)
        vlt.savefig(next_plot_fname(__file__, series='q2'))
        if show:
            vlt.show()
    except ImportError:
        pass

    try:
        from mayavi import mlab
        from viscid.plot import vlab

        try:
            fig = _global_ns['figure']
            vlab.clf()
        except KeyError:
            fig = vlab.figure(size=[1200, 800], offscreen=not show,
                              bgcolor=(1, 1, 1), fgcolor=(0, 0, 0))
            _global_ns['figure'] = fig

        vlab.clf()
        vlab.plot_lines3d(lines, scalars=scalars)
        vlab.fancy_axes()
        mlab.text(0.05, 0.05, title)
        vlab.savefig(next_plot_fname(__file__, series='q3'))
        if show:
            vlab.show(stop=True)
    except ImportError:
        pass

def anim():
    global vessels, colorGradient, timeText, showingNode
    timeText = None
    while True:
        if timeText != None:
            timeText.remove()            
        timeText = mlab.text(0.01, 0.01, 'Time: ' + str(g.globals.time), width=0.3)
        for (top, bottom, host) in vessels:
            lutT = top.module_manager.scalar_lut_manager.lut.table.to_array()
            lutB = bottom.module_manager.scalar_lut_manager.lut.table.to_array()
            count = host.getBacteriaCount()
            colorNdx = int(float(count) / parameters.cell_count_color_mapping * (len(colorGradient) - 1))
            if colorNdx >= len(colorGradient):
                colorNdx = len(colorGradient) - 1
            color = colorGradient[colorNdx]
            assert len(color) == 6
            ones = np.ones(np.shape(lutT[:, 0]))
            R = int(color[0:2], 16)
            G = int(color[2:4], 16)
            B = int(color[4:6], 16)
            #R
            lutT[:,0] = ones * R
            lutB[:,0] = ones * R
            #G
            lutT[:,1] = ones * G
            lutB[:,1] = ones * G
            #B
            lutT[:,2] = ones * B
            lutB[:,2] = ones * B            
            top.module_manager.scalar_lut_manager.lut.table = lutT
            bottom.module_manager.scalar_lut_manager.lut.table = lutB
        
        if showingNode is not None:
            (ax1, ax2) = plots
            ax1.cla()
            history = copy.deepcopy(showingNode.getCellCountHistory())
            x = np.linspace(0, len(history) * parameters.cell_count_history_interval, len(history))
            y = history
            ax1.plot(x, y, 'r')
            ax1.set_title('Bacteria Count history')            

            ax2.cla()
            history = copy.deepcopy(showingNode.getFlowHistory())
            x = np.linspace(0, len(history) * parameters.cell_count_history_interval, len(history))
            y = history
            ax2.plot(x, y, 'b')
            ax2.set_title('Blood flow history')  
        mlab.draw()
        if parameters.verbose:      
            print("updating graph")
        yield

def draw_graph3d(
    graph,
    graph_colormap="winter",
    bgcolor=(1, 1, 1),
    node_size=0.03,
    edge_color=(0.8, 0.8, 0.8),
    edge_size=0.002,
    text_size=0.008,
    text_color=(0, 0, 0),
):

    H = nx.Graph()

    # add edges
    for node, edges in graph.items():
        for edge, val in edges.items():
            if val == 1:
                H.add_edge(node, edge)

    G = nx.convert_node_labels_to_integers(H)

    graph_pos = nx.spring_layout(G, dim=3)

    # numpy array of x,y,z positions in sorted node order
    xyz = np.array([graph_pos[v] for v in sorted(G)])

    # scalar colors
    scalars = np.array(G.nodes()) + 5
    mlab.figure(1, bgcolor=bgcolor)
    mlab.clf()

    pts = mlab.points3d(
        xyz[:, 0],
        xyz[:, 1],
        xyz[:, 2],
        scalars,
        scale_factor=node_size,
        scale_mode="none",
        colormap=graph_colormap,
        resolution=20,
    )

    for i, (x, y, z) in enumerate(xyz):
        label = mlab.text(x, y, str(i), z=z, width=text_size, name=str(i), color=text_color)
        label.property.shadow = True

    pts.mlab_source.dataset.lines = np.array(G.edges())
    tube = mlab.pipeline.tube(pts, tube_radius=edge_size)
    mlab.pipeline.surface(tube, color=edge_color)

    mlab.show()  # interactive window

def add_cbar_label(cbar,title):
    position = cbar.scalar_bar_representation.position
    position2 = cbar.scalar_bar_representation.position2
    
    x = position[0] + position2[0] #+ 0.002
    y = position[1]
    
    text = mlab.text(x,y,title)
    text.property.font_size = 24
    set_text(text.property)
    
    text.property.orientation = 90.
    text.width = 0.205
    
    return text

def draw_graph3d(graph, graph_colormap='summer', bgcolor = (0, 0, 0),
                 node_size=0.05,
                 edge_color=(0.5, 0.5, 0.5), edge_size=0.0005,
                 text_size=0.004, text_color=(1, 1, 1)):

    H=nx.Graph()

    # add edges
    for node, edges in graph.items():
        for edge, val in edges.items():
            if val == 1:
                H.add_edge(node, edge)

    G=nx.convert_node_labels_to_integers(H)
    #G=H

    graph_pos=nx.random_layout(G, dim=3)

    # numpy array of x,y,z positions in sorted node order
    xyz=np.array([graph_pos[v] for v in sorted(G)])

    # scalar colors
    scalars=np.array(G.nodes())+5
    mlab.figure(1, bgcolor=bgcolor)
    mlab.clf()

    pts = mlab.points3d(xyz[:,0], xyz[:,1], xyz[:,2],
                        scalars,
                        scale_factor=node_size,
                        scale_mode='none',
                        colormap=graph_colormap,
                        resolution=20)

    mlab.figure(1).scene.disable_render = True
    for i, (x, y, z) in enumerate(xyz):
        label = mlab.text(x, y, graph.keys()[i], z=z,
                          width=len(graph.keys()[i])*text_size, name=str(i), color=text_color)
        label.property.shadow = False
    

    pts.mlab_source.dataset.lines = np.array(G.edges())
    tube = mlab.pipeline.tube(pts, tube_radius=edge_size)
    mlab.pipeline.surface(tube, color=edge_color)
    
    mlab.figure(1).scene.disable_render = False

    mlab.show() # interactive window

def draw_graph3d(graph_colormap='winter', bgcolor = (1, 1, 1),
                 node_size=0.03,
                 edge_color=(0.8, 0.8, 0.8), edge_size=0.002,
                 text_size=0.008, text_color=(0, 0, 0)): 
    
    H = nx.Graph()
    
    db = db_utils.getDBInstance()
    nodes = db.webpages.find()
    edges = []
    for node in nodes:
        iurls = node["incoming_links"]
        for iurl in iurls:
            H.add_edge(iurl, node["url"])

    G= nx.convert_node_labels_to_integers(H)

    graph_pos=nx.spring_layout(G, dim=3)

    # numpy array of x,y,z positions in sorted node order
    xyz=np.array([graph_pos[v] for v in sorted(G)])

    # scalar colors
    scalars=np.array(G.nodes())+5
    mlab.figure(1, bgcolor=bgcolor)
    mlab.clf()

    pts = mlab.points3d(xyz[:,0], xyz[:,1], xyz[:,2],
                        scalars,
                        scale_factor=node_size,
                        scale_mode='none',
                        colormap=graph_colormap,
                        resolution=20)

    for i, (x, y, z) in enumerate(xyz):
        label = mlab.text(x, y, str(i), z=z,
                          width=text_size, name=str(i), color=text_color)
        label.property.shadow = True

    pts.mlab_source.dataset.lines = np.array(G.edges())
    tube = mlab.pipeline.tube(pts, tube_radius=edge_size)
    mlab.pipeline.surface(tube, color=edge_color)

    mlab.show() # interactive window

def picker_callback(picker):
    global vessels, lastSelected, bound, showingNode, plots
    if 'lastSelected' not in globals():
        lastSelected = None
    if 'bound' not in globals():
        bound = None
    picked = picker.actors
    for (top, bottom, node) in vessels:
        if top.actor.actor._vtk_obj in [o._vtk_obj for o in picked] or bottom.actor.actor._vtk_obj in [o._vtk_obj for o in picked]:
            if lastSelected != None:
                lastSelected.remove()
            lastSelected = mlab.text(node.start.x, node.start.y, node.name, width=0.2, z=node.start.z)
            if bound != None:
                (tBound, bBound) = bound
                tBound.remove()
                bBound.remove()
            tBound = mlab.outline(top)
            bBound = mlab.outline(bottom)
            bound = (tBound, bBound)
            if plots is None:
                f, plots = plt.subplots(2, 1)
            (ax1, ax2) = plots
            ax1.cla()
            history = copy.deepcopy(node.getCellCountHistory())
            x = np.linspace(0, len(history) * parameters.cell_count_history_interval, len(history))
            y = history
            ax1.plot(x, y, 'r')
            ax1.set_title('Bacteria Count history')            

            ax2.cla()
            history = copy.deepcopy(node.getFlowHistory())
            x = np.linspace(0, len(history) * parameters.cell_count_history_interval, len(history))
            y = history
            ax2.plot(x, y, 'b')
            ax2.set_title('Blood flow history')            

            if showingNode is None:
                plt.show()

            showingNode = node

            break