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

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

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

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

    mlab.title(subject)

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

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

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

        surf.module_manager.scalar_lut_manager.lut.table = lut

    mlab.draw()

    return f

    def draw_conns(self,new_edges=None):
        try:
            self.thres.set(lower_threshold=self.ds.thresval)
            lo=self.thres.lower_threshold; hi=self.thres.upper_threshold

            set_lut(self.vectors,self.ds.opts.activation_map)

            if new_edges is not None:
                new_starts=self.ds.lab_pos[new_edges[:,0]]
                new_vecs=self.ds.lab_pos[new_edges[:,1]] - new_starts

                self.vectors.mlab_source.reset(
                    x=new_starts[:,0],y=new_starts[:,1],z=new_starts[:,2],
                    u=new_vecs[:,0],v=new_vecs[:,1],w=new_vecs[:,2])

            if self.ds.curr_node is not None:
                self.vectors.actor.property.opacity=.75
                self.txt.set(text='  %s'%self.ds.labnam[self.ds.curr_node])
            else:
                self.vectors.actor.property.opacity=(
                    .5 if self.ds.opts.tube_conns else .3)
                self.txt.set(text='')

            mlab.draw()

        # In case the user changes the threshold while there are no connections
        # present and so the VTK objects have not been created yet
        except AttributeError:
            pass

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 draw(self,**kwargs):
    x = np.linspace(self._base_square_x[0] + 0.5*(self._base_square_x[1]-self._base_square_x[0])/self._Nl,
                    self._base_square_x[1] - 0.5*(self._base_square_x[1]-self._base_square_x[0])/self._Nl,
                    self._Nl)
    y = np.linspace(self._base_square_y[0] + 0.5*(self._base_square_y[1]-self._base_square_y[0])/self._Nw,
                    self._base_square_y[1] - 0.5*(self._base_square_y[1]-self._base_square_y[0])/self._Nw,
                    self._Nw)
    x,y = np.meshgrid(x,y)
    z = 0.0*x
    for trans in self._transforms:
      p = np.concatenate((x[:,:,None],
                          y[:,:,None],
                          z[:,:,None]),
                          axis=-1)
      p = trans(p)
      pflat = np.reshape(p,(-1,3))
      x,y,z = pflat[:,0],pflat[:,1],pflat[:,2]
      value = self._f(pflat,*self._f_args,**self._f_kwargs)
      u,v,w = value[:,0],value[:,1],value[:,2]
      m = mlab.quiver3d(x,y,z,u,v,w,mode='arrow',color=(1.0,1.0,1.0),scale_factor=self.scale_units,
                        resolution=20,**kwargs)
      mlab.draw()
 

      if self._plots is None:
        self._plots = (m,trans),

      else:
        self._plots += (m,trans),

    return [i[0] for i in self._plots]

def draw_scene():
    s = mlab.pipeline.triangular_mesh_source(x, y, z, triIndices)
    s.data.cell_data.scalars = np.cos(phaseAngle)
    surf = mlab.pipeline.surface(s)
    surf.contour.filled_contours = True
    surf.contour.minimum_contour = 0.0
    surf.contour.maximum_contour = 1.0
    surf.module_manager.scalar_lut_manager.data_range = (0,1)
    mlab.plot3d(xSun_plt, ySun_plt, zSun_plt, tube_radius=fStretch/1000, color=(1,1,0))
    mlab.plot3d(xSC_plt, ySC_plt, zSC_plt, tube_radius=fStretch/1000, color=(0,0,1))

    ball_x = []
    ball_y = []
    ball_z = []
    for i in range(nPixelsX):
        for j in range(nPixelsY):
            p = np.dot(R, pVectors[:,i,j])
            p_tan = np.dot(rCG, p) * p + rSC

            xVIR_plt, yVIR_plt, zVIR_plt = plt_coords(rSC,  1.1*fStretch*p)
            mlab.plot3d(xVIR_plt, yVIR_plt, zVIR_plt, tube_radius=fStretch/5000, color=(0,0,0))
            ball_x.append(p_tan[0])
            ball_y.append(p_tan[1])
            ball_z.append(p_tan[2])

    mlab.points3d(ball_x, ball_y, ball_z, np.ones(len(ball_x)), scale_factor=150,
                  color=(1,0.7,0.1))

    mlab.draw()

def surfcf(gridx, gridy, phase, modulus, colormap=None):
    r"""Plot the modulus of a complex valued function :math:`f:R^2 -> C`
    together with its phase in a color coded fashion.

    :param gridx: The grid nodes along the :math:`x` axis of the real domain :math:`R^2`
    :param gridy: The grid nodes along the :math:`y` axis of the real domain :math:`R^2`
    :param phase: The phase of the complex domain result f(grid)
    :param modulus: The modulus of the complex domain result f(grid)
    :param colormap: The colormap to use, if none is given, compute the 'default' QM colormap.
    """
    if colormap is None:
        colormap = compute_color_map()

    # The real(.) is necessary just to get an array with dtype real
    mesh = mlab.mesh(gridx, gridy, real(modulus), scalars=phase)

    # Set the custom color map
    mesh.module_manager.scalar_lut_manager.use_default_range = False
    mesh.module_manager.scalar_lut_manager.data_range = [-pi, pi]
    lut = mesh.module_manager.scalar_lut_manager.lut.table.to_array()
    lut[:, 0:3] = colormap.copy()
    mesh.module_manager.scalar_lut_manager.lut.table = lut

    # Update the figure
    mlab.draw()

    return mesh

def plot_isosurface(crystal):
    filename = './output/potentialfield.txt'
    data = np.genfromtxt(filename, delimiter='\t')
    size = np.round((len(data))**(1/3))
    X = np.reshape(data[:,0], (size,size,size))
    Y = np.reshape(data[:,1], (size,size,size))
    Z = np.reshape(data[:,2], (size,size,size))
    DeltaU = np.reshape(data[:,3], (size,size,size))
    average = np.average(crystal.coordinates[:,0])
    start = average - crystal.a
    end = average + crystal.a
    coords1 = np.array([[start, start, start]])
    coords2 = np.array([[end, end, end]])
    array1 = np.repeat(coords1,len(crystal.coordinates),axis=0)
    array2 = np.repeat(coords2,len(crystal.coordinates),axis=0)
    basefilter1 = np.greater(crystal.coordinates,array1)
    basefilter2 = np.less(crystal.coordinates,array2)
    basefilter = np.nonzero(np.all(basefilter1*basefilter2, axis=1))
    base = crystal.coordinates[basefilter]   

    mlab.figure(bgcolor=(1, 1, 1), fgcolor=(1, 1, 1), size=(2048,2048))
    dataset = mlab.contour3d(X, Y, Z, DeltaU, contours=[3.50],color=(1,0.25,0))
    scatter = mlab.points3d(base[:,0], 
                            base[:,1], 
                            base[:,2],
                            color=(0.255,0.647,0.88),
                            resolution=24, 
                            scale_factor=1.0, 
                            opacity=0.40)
    mlab.view(azimuth=17, elevation=90, distance=10, focalpoint=[average,average-0.2,average])
    mlab.draw()
    savename = './output/3Dpotential.png'
    mlab.savefig(savename, size=(2048,2048))
    mlab.show()

def process_launch():
    '''Procédure reliant une fenetre graphique et le coeur du programme'''
    global nb_etapesIV
    nb_etapes=nb_etapesIV.get()#On récupère le nombre d'étapes
    fig=mlab.figure(1)
    mlab.clf()#La fenêtre de dessin est initialisée
    mlab.draw(terrain([(0,1,2),(2,3,4),(4,5,6)],[(Point(0,0,0),Point(1,0,0)),(Point(1,0,0),Point(1,1,0)),(Point(0,0,0),Point(1,1,0)),(Point(1,1,0),Point(0,1,0)),(Point(0,0,0),Point(0,1,0)),(Point(0,0,0),Point(-1,1,0)),(Point(-1,1,0),Point(0,1,0))],nb_etapes))#On affiche le dessin

 def __init__(self,start,end,maxd=5000.,n=100):
     '''
     Constructor
     '''
     self.start=start
     self.end=end
     self.lopath=numpy.linspace(start[0], end[0], n)
     self.lapath=numpy.linspace(start[1], end[1], n)
     self.set_proj()
     self.tile=TiffReader(lon=self.lopath[0],lat=self.lapath[0])
     self.tile.readit()
     for i in range(n):
             if not self.tile==TiffReader(lon=self.lopath[i],lat=self.lapath[i]):
                     self.tile=TiffReader(lon=self.lopath[i],lat=self.lapath[i])
                     self.tile.readit()
             if not hasattr(self,'mesh'):
                     lo,la,z=self.tile.subset(rect=None, around=(self.lopath[i],self.lapath[i],maxd))
                     x,y=self.proj(lo,la)
                     x=x-x.mean()
                     y=y-y.mean()
                     self.mesh=mlab.mesh(x,y,z,scalars=z,vmax=1500.,vmin=0.)
                     mlab.view(180.,45.,maxd,numpy.array([x.max(),0,z.max()]))
             else:
                     lo,la,self.mesh.mlab_source.z=self.tile.subset(rect=None, around=(self.lopath[i],self.lapath[i],maxd))
                     self.mesh.mlab_source.scalars=self.mesh.mlab_source.z
                     mlab.view(180.,45.,5*x.max(),numpy.array([x.max(),0,self.mesh.mlab_source.z.max()]))
             mlab.draw()

def mlab_imshowColor(im, alpha=255, **kwargs):
    """
    Plot a color image with mayavi.mlab.imshow.
    im is a ndarray with dim (n, m, 3) and scale (0->255]
    alpha is a single number or a ndarray with dim (n*m) and scale (0->255]
    **kwargs is passed onto mayavi.mlab.imshow(..., **kwargs)
    """
    try:
        alpha[0]
    except:
        alpha = pl.ones(im.shape[0] * im.shape[1]) * alpha
    if len(alpha.shape) != 1:
        alpha = alpha.flatten()

    # The lut is a Nx4 array, with the columns representing RGBA
    # (red, green, blue, alpha) coded with integers going from 0 to 255,
    # we create it by stacking all the pixles (r,g,b,alpha) as rows.
    myLut = pl.c_[im.reshape(-1, 3), alpha]
    myLutLookupArray = pl.arange(im.shape[0] * im.shape[1]).reshape(im.shape[0], im.shape[1])

    #We can display an color image by using mlab.imshow, a lut color list and a lut lookup table.
    theImshow = mlab.imshow(myLutLookupArray, colormap='binary', **kwargs) #temporary colormap
    theImshow.module_manager.scalar_lut_manager.lut.table = myLut
    mlab.draw()

    return theImshow

    def generate_plots_3d(self):
        self.ax = mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(800, 600))
        self.clf = mlab.clf()

        minS, maxS = maxint, 0
        contour_plots = []
        for cond in self.conductors.itervalues():

            minS, maxS, face_data = self.generate_plot_data_for_faces_3d(cond, minS, maxS)
            for (x, y, z, s) in face_data:
                if isinstance(cond, conductor_type_3d['Unstructured']):
                    pts = mlab.points3d(x, y, z, s, scale_mode='none', scale_factor=0.002)
                    mesh = mlab.pipeline.delaunay3d(pts)
                    contour_plots.append(mlab.pipeline.surface(mesh, colormap='viridis'))
                else:
                    if np.min(s) < 0.0:
                        contour_plots.append(mlab.mesh(x, y, z, color=(0, 0, 0), colormap='viridis'))
                    else:
                        contour_plots.append(mlab.mesh(x, y, z, scalars=s, colormap='viridis'))

        for cp in contour_plots:
            cp.module_manager.scalar_lut_manager.trait_set(default_data_range=[minS * 0.95, maxS * 1.05])

        mlab.draw()
        mlab.colorbar(object=contour_plots[0], orientation='vertical')
        mlab.show()

  def draw(self):
    x = np.linspace(self._base_square_x[0],
                    self._base_square_x[1],self._Nl)
    y = np.linspace(self._base_square_y[0],
                    self._base_square_y[1],self._Nw)
    x,y = np.meshgrid(x,y)
    z = 0.0*x
    for trans in self._transforms:
      p = np.concatenate((x[:,:,None],
                          y[:,:,None],
                          z[:,:,None]),
                          axis=-1)
      p = trans(p)
      pflat = np.reshape(p,(-1,3))
      x,y,z = pflat[:,0],pflat[:,1],pflat[:,2]
      value = self._f(pflat,*self._f_args,**self._f_kwargs)
      u,v,w = value[:,0],value[:,1],value[:,2]
      m = mlab.quiver3d(x,y,z,u,v,w,mode='arrow',color=(0.4,0.4,0.4))
      mlab.draw()
 

      if self._plots is None:
        self._plots = (m,trans),

      else:
        self._plots += (m,trans),

def surfacePlot(
    Hmap,
    nrows,
    ncols,
    xyspacing,
    zscale,
    name,
    hRange,
    file_path,
    lutfromfile,
    lut,
    lut_file_path,
    colorbar_on,
    save,
    show,
):

    # Create a grid of the x and y coordinates corresponding to each pixel in the height matrix
    x, y = np.mgrid[0 : ncols * xyspacing : xyspacing, 0 : nrows * xyspacing : xyspacing]

    # Create a new figure
    mlab.figure(size=(1000, 1000))

    # Set the background color if desired
    # bgcolor=(0.16, 0.28, 0.46)

    # Create the surface plot of the reconstructed data
    plot = mlab.surf(x, y, Hmap, warp_scale=zscale, vmin=hRange[0], vmax=hRange[1], colormap=lut)

    # Import the LUT from a file if necessary
    if lutfromfile:
        plot.module_manager.scalar_lut_manager.load_lut_from_file(lut_file_path)

    # Draw the figure with the new LUT if necessary
    mlab.draw()

    # Zoom in to fill the entire window
    f = mlab.gcf()
    f.scene.camera.zoom(1.05)

    # Change the view to a top-down perspective
    mlab.view(270, 0)

    # Add a colorbar if indicated by colorbar_on (=True)
    if colorbar_on:
        # mlab.colorbar(title='Height (nm)', orientation='vertical')
        mlab.colorbar(orientation="vertical")

    # Save the figure if indicated by save (=True)
    if save:
        mlab.savefig(file_path, size=(1000, 1000))
        if show == False:
            mlab.close()

    # Keep the figure open if indicated by show (=True)
    if show:
        mlab.show()

def force_render( figure=None ):
    from mayavi import mlab
    figure.scene.render()
    mlab.draw(figure=figure)
    from pyface.api import GUI
    _gui = GUI()
    orig_val = _gui.busy
    _gui.set_busy(busy=True)
    _gui.process_events()
    _gui.set_busy(busy=orig_val)
    _gui.process_events()

def draw_nucleus(v):
    mlab.clf()
    s = mlab.pipeline.triangular_mesh_source(x, y, z, triIndices)
    s.data.cell_data.scalars = np.cos(phaseAngle)
    surf = mlab.pipeline.surface(s)
    surf.contour.filled_contours = True
    surf.contour.minimum_contour = 0.0
    surf.contour.maximum_contour = 1.0
    surf.module_manager.scalar_lut_manager.data_range = (0,1)
    mlab.view(v)
    mlab.draw()

def parse_box(box_data,indx):
  """
  from box_data produce a 3D image of surfaces and display it using an
  animation of a slice moving back and forth through the box
  NOT FUNCTIONAL AT PRESENT - NEED TO WORK OUT THE MAYAVI SYNTAX
  -can't get different slices to display in a nice fashion
  """

  s=box_data
  sp=1-box_data

  mlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),
                            plane_orientation='z_axes',
                            slice_index=0,colormap='black-white'
                        )
  mlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),
                            plane_orientation='x_axes',
                            slice_index=0,colormap='black-white'
                        )
  l=mlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),
                            plane_orientation='y_axes',
                            slice_index=0,colormap='black-white'
                        )


  #s=sp

  fig=mlab.gcf()
  ms=l.mlab_source
  print 'indx=',dir(ms)
  for i in range(90):
    #l=mlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),
    #                        plane_orientation='y_axes',
    #                        slice_index=i,colormap='black-white'
    #                    )
    if(i%2):
      ms.reset(scalars=sp)
    else:
      ms.reset(scalars=s)
      
    #print dir(ms)
    #print ms
    #camera=fig.scene.camera
    #camera.yaw(i)
    fig.scene.reset_zoom()
    mlab.draw()
    filename='volume_slice_frame'+str(i)+'.png'
    print filename
    #mlab.savefig(filename)
    
  mlab.show()
    
  return

	def draw_nodes(self): 
		nc=self.ds.scalar_display_settings.node_color
		ns=self.ds.scalar_display_settings.node_size

		lhn=self.nodes_lh; lhn_ix=self.ds.lhnodes
		rhn=self.nodes_rh; rhn_ix=self.ds.rhnodes

		for nodes,ixes in ((lhn,lhn_ix),(rhn,rhn_ix)):
			nr=len(ixes)
			#set node size
			if self.ds.display_mode=='scalar' and ns:
				nodes.glyph.scale_mode='scale_by_vector'
				nodes.glyph.glyph.scale_factor=8
				nodes.mlab_source.dataset.point_data.vectors=(
					np.tile(self.ds.node_scalars[ns][ixes],(3,1)).T)
			else:
				nodes.glyph.scale_mode='data_scaling_off'
				nodes.glyph.glyph.scale_factor=3

		#set node color -- we dont care about ds.node_colors for mayavi
			if (self.ds.display_mode=='normal' or
					(self.ds.display_mode=='scalar' and not nc)):
				set_lut(nodes,self.ds.opts.default_map)
				nodes.mlab_source.dataset.point_data.scalars=np.tile(.3,nr)

			elif self.ds.display_mode=='scalar': #and nc must be true
				set_lut(nodes,self.ds.opts.scalar_map)
				nodes.mlab_source.dataset.point_data.scalars=(
					self.ds.node_scalars[nc][ixes])

			elif self.ds.display_mode=='module_single':
				set_lut(nodes,self.ds.opts.default_map)
				new_colors=np.tile(.3,self.ds.nr_labels)
				new_colors[self.ds.get_module()]=.8
				nodes.mlab_source.dataset.point_data.scalars=new_colors[ixes]

			elif self.ds.display_mode=='module_multi':
				new_colors=np.array(self.ds.module_colors[:self.ds.nr_modules])
				manager=nodes.module_manager.scalar_lut_manager
				#set the mayavi object to the dummy cmap that we hide from user
				#so that when changed notifications will work correctly
				manager.lut_mode='black-white'
				#now adjust its LUT manually
				manager.number_of_colors=self.ds.nr_modules
				manager.lut.table=new_colors	

				#set the mayavi scalars to be fractions between 0 and 1
				import bct
				nodes.mlab_source.dataset.point_data.scalars=(bct.ls2ci(
					self.ds.modules,zeroindexed=True)/self.ds.nr_modules)[ixes]

		mlab.draw()

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 plot_shape(self, kind, spatialSize=None, renderSize=None,
                   features=None):
        """ outline=True  - box edges
            title=True - print kind of picture
            **kwargs # one of: vertA, faceA, features, x, y, z, spatialSize
        """
        f = mlab.figure(bgcolor=(1, 1, 1))
        if kind == 'solid':
            vertA, faceA = self.vertices, self.faces
            mlab.triangular_mesh(vertA[:, 0], vertA[:, 1], vertA[:, 2], faceA)
        elif kind == 'transparent':
            vertA, faceA = self.vertices, self.faces
            mlab.triangular_mesh(vertA[:,0], vertA[:, 1], vertA[:, 2], faceA,
                                 opacity=0.1)
        elif kind == 'wireframe':
            vertA, faceA = self.vertices, self.faces
            mlab.triangular_mesh(vertA[:, 0], vertA[:, 1], vertA[:, 2], faceA,
                                 representation='wireframe')
        elif kind == "Bool":
            x, y, z = map(np.array, self.voxelize(spatialSize, renderSize))
            assert len(x) == len(y) == len(z)
            N = len(x)
            scalars = np.arange(N) # Key point: set an integer for each point
            colors = np.zeros((N, 4), dtype=np.uint8)
            colors[:, -1] = 255  # No transparency
            if features is not None:
                features = features.ravel()
                colors[:, 0] = 255
                colors[:, 1] = (
                255 * (1 - features / np.max(features))).astype(np.uint8)
                colors[:, 2] = (
                255 * (1 - features / np.max(features))).astype(np.uint8)
            else:
                colors[:, 0] = 0
                colors[:, 1] = 255
                colors[:, 2] = 0

            pts = mlab.quiver3d(
                x-spatialSize/2,
                y-spatialSize/2+0.5,
                z-spatialSize/2+0.5,
                np.ones(N), np.zeros(N), np.zeros(N),
                scalars=scalars, mode='cube', scale_factor=0.7, line_width=10)
            pts.glyph.color_mode = 'color_by_scalar'
            try:
                pts.module_manager.scalar_lut_manager.lut.table = colors
            except:
                pass
            mlab.draw()

        return f

def plot_potential(grid, potential, sparsify=1, along_axes=False, view=None, interactive=False, path='.'):
    """Plot the potential

    :param iom: An :py:class:`IOManager` instance providing the simulation data.
    """
    # The Grid
    u, v = grid.get_nodes(split=True, flat=False)
    u = real(u[::sparsify, ::sparsify])
    v = real(v[::sparsify, ::sparsify])

    # Create potential and evaluate eigenvalues
    potew = potential.evaluate_eigenvalues_at(grid)
    potew = [real(level).reshape(grid.get_number_nodes(overall=False))[::sparsify, ::sparsify] for level in potew]

    # Plot
    if not interactive:
        mlab.options.offscreen = True

    fig = mlab.figure(size=(800, 700))

    for level in potew:
        # The energy surfaces of the potential
        src = mlab.pipeline.grid_source(u, v, level)

        # Clip to given view
        if view is not None:
            geometry_filter = mlab.pipeline.user_defined(src, filter='GeometryFilter')
            geometry_filter.filter.extent_clipping = True
            geometry_filter.filter.extent = view
            src = mlab.pipeline.user_defined(geometry_filter, filter='CleanPolyData')

        # Plot the surface
        normals = mlab.pipeline.poly_data_normals(src)
        mlab.pipeline.surface(normals)

    mlab.axes()

    fig.scene.parallel_projection = True
    fig.scene.isometric_view()
    # fig.scene.show_axes = True

    mlab.draw()
    if interactive:
        mlab.show()
    else:
        mlab.savefig(os.path.join(path, "potential_3D_view.png"))
        mlab.close(fig)