def main(argv):
    if ( len(argv) != 1 ):
        print ("Usage: python visualizeField.py <datafile.json>")
        return

    # Read input file
    try:
        infile = open( argv[0], 'r' )
    except:
        print ("Error when opening file %s"%(argv[0]))
        return
    data = json.load( infile )
    infile.close()

    y = np.array( data["points"]["y"] )
    z = np.array( data["points"]["z"] )
    x = np.zeros( len(y) )
    Ex = np.array( data["field"]["x"] )
    pts = mlab.points3d( y, z, x, Ex, scale_mode="scalar", scale_factor=0.0, mode="point")
    mesh = mlab.pipeline.delaunay2d( pts )
    surf = mlab.pipeline.surface( mesh )
    
    fname = "Figures/fieldyzPlane.png"
    mlab.view(0.0,0.0,1.0, (0.0,0.0,0.0))
    mlab.scalarbar()
    mlab.xlabel( "y" )
    mlab.ylabel( "z" )
    mlab.show()

 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 setup():
    sc = mlab.figure(0)
    mlab.view(90., -90., 3.)
    mlab.roll(90.)
    arm = Arm()
    neuron = Neuron(arm)
    return neuron, arm

	def _plotbutton1_fired(self):
		mlab.clf()
		self.loaddata()
		self.sregion[np.where(self.sregion<self.datamin)] = self.datamin
		self.sregion[np.where(self.sregion>self.datamax)] = self.datamax

		# The following codes from: http://docs.enthought.com/mayavi/mayavi/auto/example_atomic_orbital.html#example-atomic-orbital
		field = mlab.pipeline.scalar_field(self.sregion)     # Generate a scalar field
		colored = self.sregion
		vol=self.sregion.shape
		for v in range(0,vol[2]-1):
			colored[:,:,v] = self.extent[4] + v*(-1)*abs(self.hdr['cdelt3'])
		new = field.image_data.point_data.add_array(colored.T.ravel())
		field.image_data.point_data.get_array(new).name = 'color'
		field.image_data.point_data.update()

		field2 = mlab.pipeline.set_active_attribute(field, point_scalars='scalar')
		contour = mlab.pipeline.contour(field2)
		contour2 = mlab.pipeline.set_active_attribute(contour, point_scalars='color')

		mlab.pipeline.surface(contour2, colormap='jet', opacity=self.opacity)
		
		## Insert a continuum plot
		##im = pyfits.open('g28_SMA1.cont.image.fits')
		##dat = im[0].data
		##dat0 = dat[0]
		##channel = dat0[0]
		##region = np.swapaxes(channel[self.xstart:self.xend,self.ystart:self.yend]*1000.,0,1)
		##field = mlab.contour3d(region, colormap='gist_ncar')
		##field.contour.minimum_contour = 5

		self.field = field
		self.labels()
		mlab.view(azimuth=0, elevation=0, distance='auto')
		mlab.show()

    def plot_points(self, label, X, color=None, size=None, mode=None):

        mlab.figure(self.figure.name)

        if color == None:
            color = (1, 0, 0)

        if size == None and mode == None or size == 0:
            size = 1
            mode = "point"
        if size == None:
            size = 1
        if mode == None:
            mode = "sphere"

        if isinstance(X, list):
            X = scipy.array(X)

        if len(X.shape) == 1:
            X = scipy.array([X])

        mlab_obj = self.plots.get(label)
        if mlab_obj == None:
            if isinstance(color, tuple):
                self.plots[label] = mlab.points3d(X[:, 0], X[:, 1], X[:, 2], color=color, scale_factor=size, mode=mode)
            else:
                self.plots[label] = mlab.points3d(
                    X[:, 0], X[:, 1], X[:, 2], color, scale_factor=size, scale_mode="none", mode=mode
                )

        else:
            self.figure.scene.disable_render = True
            view = mlab.view()

            ### Commented out since VTK gives an error when using mlab_source.set
            # ~ if X.shape[0] == mlab_obj.mlab_source.x.shape[0]:
            # ~ if isinstance(color, tuple):
            # ~ mlab_obj.mlab_source.set(x=X[:,0], y=X[:,1], z=X[:,2])
            # ~ mlab_obj.actor.property.color = color
            # ~ else:
            # ~ mlab_obj.mlab_source.set(x=X[:,0], y=X[:,1], z=X[:,2], scalars=color)
            # ~
            # ~
            # ~ else:
            # ~ self.clear(label)
            # ~ if isinstance(color, tuple):
            # ~ self.plots[label] = mlab.points3d(X[:,0], X[:,1], X[:,2], color=color, scale_factor=size, mode=mode)
            # ~ else:
            # ~ self.plots[label] = mlab.points3d(X[:,0], X[:,1], X[:,2], color, scale_factor=size, scale_mode='none', mode=mode)

            self.clear(label)
            if isinstance(color, tuple):
                self.plots[label] = mlab.points3d(X[:, 0], X[:, 1], X[:, 2], color=color, scale_factor=size, mode=mode)
            else:
                self.plots[label] = mlab.points3d(
                    X[:, 0], X[:, 1], X[:, 2], color, scale_factor=size, scale_mode="none", mode=mode
                )

            mlab.view(*view)
            self.figure.scene.disable_render = False

def view():
    from mayavi import mlab
    from mayavi.modules.glyph import Glyph
    from simphony_mayavi.sources.api import CUDSFileSource

    mayavi.new_scene()

    # Mayavi Source
    src = CUDSFileSource()
    src.initialize('lattices.cuds')

    # choose a dataset for display
    src.dataset = 'orthorhombic'

    mayavi.add_source(src)

    # customise the visualisation
    g = Glyph()
    gs = g.glyph.glyph_source
    gs.glyph_source = gs.glyph_dict['sphere_source']
    g.glyph.glyph.scale_factor = 0.05
    g.glyph.scale_mode = 'data_scaling_off'
    mayavi.add_module(g)

    # add legend
    module_manager = src.children[0]
    module_manager.scalar_lut_manager.show_scalar_bar = True
    module_manager.scalar_lut_manager.show_legend = True

    # customise the camera
    mlab.view(63., 38., 3., [5., 4., 0.])

    def set_parallel_view(self, forward=None, up=None, scale=None):
        """Set view to parallel projection

        Parameters
        ----------
        forward : scalar
            Move the view forward (mm).
        up : scalar
            Move the view upward (mm).
        scale : scalar
            Mayavi parallel_scale parameter. Default is 95 for the inflated
            surface, 75 otherwise. Smaller numbers correspond to zooming in.
        """
        if scale is True:
            surf = self.geo['rh' if self._hemi == 'rh' else 'lh'].surf
            if surf == 'inflated':
                scale = 95
            else:
                scale = 75  # was 65 for WX backend

        from mayavi import mlab

        for figs in self._figures:
            for fig in figs:
                if forward is not None or up is not None:
                    mlab.view(focalpoint=(0, forward or 0, up or 0),
                              figure=fig)
                if scale is not None:
                    fig.scene.camera.parallel_scale = scale
                fig.scene.camera.parallel_projection = True
                fig.render()

 def anim():
   global bondcolors
   while True:
     print mlab.view(), mlab.roll()
     bondcolors[:] = bondcolors*.5
     bondsrc.update()
     yield

    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 draw_volume(scalars):
    # filts out low data
    nmax = np.amax(scalars)
    nmin = np.amin(scalars)
    #nptp = nmax - nmin
    #lowBound = nmin + nptp*0
    #for i in np.nditer(scalars, op_flags=['readwrite']):
    #    if i<lowBound: i[...] = 0
        
    # Get the shape of axes
    shape = scalars.shape
    
    # draw
    fig = mlab.figure(bgcolor=(0,0,0), size=(800,800) )
    src = mlab.pipeline.scalar_field(scalars)
    src.update_image_data = True
    #src.spacing = [1, 1, 1]
    
    vol = mlab.pipeline.volume(src)
    
    change_volume_property(vol)
    #change_test(vol)
    
    #vol = mlab.pipeline.image_plane_widget(src,plane_orientation='z_axes', slice_index=shape[2]/2)
    ax = mlab.axes(nb_labels=5, ranges=(0,shape[0],0,shape[1],0,shape[2]))
    
    mlab.savefig("result.jpg")
    print "nmax =", nmax
    print "nmin =", nmin
    
    viewkeeper = mlab.view()
    print "View", mlab.view()
    draw_map(shape[0],shape[1])
    mlab.view(viewkeeper[0],viewkeeper[1],viewkeeper[2],viewkeeper[3])

 def plotvtk3D(self):
     """
     3D plot using the vtk libraries
     """
     from tvtk.api import tvtk
     from mayavi import mlab
     # Plot the data on a 3D grid
     xy = np.column_stack((self.grd.xp,self.grd.yp))
     dvp = self.F(xy)
     vertexag = 50.0
     
     points = np.column_stack((self.grd.xp,self.grd.yp,dvp*vertexag))
     tri_type = tvtk.Triangle().cell_type
     #tet_type = tvtk.Tetra().cell_type
     ug = tvtk.UnstructuredGrid(points=points)
     ug.set_cells(tri_type, self.grd.cells)
     
     ug.cell_data.scalars = self.grd.dv
     ug.cell_data.scalars.name = 'depths'
     
     f=mlab.gcf()
     f.scene.background = (0.,0.,0.)
     d = mlab.pipeline.add_dataset(ug)
     h=mlab.pipeline.surface(d,colormap='gist_earth')
     mlab.colorbar(object=h,orientation='vertical')
     mlab.view(0,0)
     
     outfile = self.suntanspath+'/depths.png'
     f.scene.save(outfile)      
     print 'Figure saved to %s.'%outfile

 def plot_text(self, label, X, text, size=1):
     view = mlab.view()
     roll = mlab.roll()
     self.figure.scene.disable_render = True
     
     scale = (size, size, size)
     mlab_objs = self.plots.get(label)
     
     if mlab_objs != None:
         if len(mlab_objs) != len(text):
             for obj in mlab_objs:
                 obj.remove()
         self.plots.pop(label)
     
     mlab_objs = self.plots.get(label)
     if mlab_objs == None:
         text_objs = []
         for x, t in zip(X, text):
             text_objs.append(mlab.text3d(x[0], x[1], x[2], str(t), scale=scale))
         self.plots[label] = text_objs
     elif len(mlab_objs) == len(text):
         for i, obj in enumerate(mlab_objs):
             obj.position = X[i,:]
             obj.text = str(text[i])
             obj.scale = scale
     else:
         print "HELP, I shouldn\'t be here!!!!!"
     
     self.figure.scene.disable_render = False
     mlab.view(*view)
     mlab.roll(roll)

    def reset_view(self):
        center = self.shape * self.spacing / 2. + (self.shape + 1) % 2 * self.spacing / 2.
        width = (self.shape * self.spacing)[:2]
        width = np.min(width) * 0.5

        self.scene.scene.background = (0, 0, 0)
        mlab.view(*([(0, 0), (90, 0), (0, 0)][self.axis]),
                  focalpoint=center,
                  figure=self.scene.mayavi_scene)
        self.scene.scene.parallel_projection = True
        self.scene.scene.camera.parallel_scale = width * 1.2
        self.scene.scene.interactor.interactor_style = tvtk.InteractorStyleImage()

        try: #WX window
            self.scene.scene_editor.control.SetFocusFromKbd
            def focusfunc(vtkobj, i):
                self.scene.scene_editor.control.SetFocusFromKbd()
        except AttributeError: #QT window
            self.scene.scene_editor.control.setFocus
            def focusfunc(vtkobj, i):
                self.scene.scene_editor.control.setFocus()

        self.scene.interactor.add_observer("MouseMoveEvent", focusfunc)
        self.scene.interactor.add_observer("KeyReleaseEvent", self.handle_keys)
        self._outline_color_changed()

    def setUp(self):
        # set up source
        sgrid = datasets.generateStructuredGrid()
        source = VTKDataSource(data=sgrid)

        self.engine = mlab.get_engine()

        # set up scene, first scene is empty
        # second scene has the settings we want to restore
        for _ in range(2):
            fig = mlab.figure()
            fig.scene.off_screen_rendering = True

        # add source
        self.engine.add_source(source)

        # add more modules
        self.engine.add_module(IsoSurface())
        self.engine.add_module(Text3D())
        self.modules = source.children[0].children

        # set camera
        self.view = (25., 14., 20., [0., 0., 2.5])
        mlab.view(*self.view)

        # save the visualisation
        self.temp_dir = tempfile.mkdtemp()
        self.filename = os.path.join(self.temp_dir, "test_vis.mv2")
        self.engine.save_visualization(self.filename)

        # save the scene as an image for comparison later
        self.ref_saved_filename = os.path.join(self.temp_dir, "ref_saved.png")
        mlab.savefig(self.ref_saved_filename)

def view():
    from mayavi.modules.glyph import Glyph
    from simphony_mayavi.sources.api import EngineSource
    from mayavi import mlab

    # Define EngineSource, choose dataset
    src = EngineSource(engine=engine_wrapper,
                       dataset="particles")

    # choose the CUBA attribute for display
    src.point_scalars_name = "TEMPERATURE"

    mayavi.add_source(src)

    # customise the visualisation
    g = Glyph()
    gs = g.glyph.glyph_source
    gs.glyph_source = gs.glyph_dict['sphere_source']
    g.glyph.glyph.scale_factor = 0.2
    g.glyph.scale_mode = 'data_scaling_off'
    mayavi.add_module(g)

    # add legend
    module_manager = src.children[0]
    module_manager.scalar_lut_manager.show_scalar_bar = True
    module_manager.scalar_lut_manager.show_legend = True

    # set camera
    mlab.view(-65., 60., 14., [1.5, 2., 2.5])

def mark_gt_box3d( lidar_dir, gt_boxes3d_dir, mark_dir):

    os.makedirs(mark_dir, exist_ok=True)
    fig   = mlab.figure(figure=None, bgcolor=(0,0,0), fgcolor=None, engine=None, size=(500, 500))
    dummy = np.zeros((10,10,3),dtype=np.uint8)

    for file in sorted(glob.glob(lidar_dir + '/*.npy')):
        name = os.path.basename(file).replace('.npy','')

        lidar_file   = lidar_dir     +'/'+name+'.npy'
        boxes3d_file = gt_boxes3d_dir+'/'+name+'.npy'
        lidar   = np.load(lidar_file)
        boxes3d = np.load(boxes3d_file)

        mlab.clf(fig)
        draw_didi_lidar(fig, lidar, is_grid=1, is_axis=1)
        if len(boxes3d)!=0:
            draw_didi_boxes3d(fig, boxes3d)

        azimuth,elevation,distance,focalpoint = MM_PER_VIEW1
        mlab.view(azimuth,elevation,distance,focalpoint)
        mlab.show(1)
        imshow('dummy',dummy)
        cv2.waitKey(1)

        mlab.savefig(mark_dir+'/'+name+'.png',figure=fig)

    def __view(self, viewargs=None, roll=None):
        """Wrapper for mlab.view()

        Parameters
        ----------
        viewargs: dict
            mapping with keys corresponding to mlab.view args
        roll: num
            int or float to set camera roll

        Returns
        -------
        camera settings: tuple
            view settings, roll setting

        """
        try:
            from mayavi import mlab
        except ImportError:
            from enthought.mayavi import mlab

        if viewargs:
            viewargs['reset_roll'] = True
            mlab.view(**viewargs)
        if not roll is None:
            mlab.roll(roll)
        return mlab.view(), mlab.roll()

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 curvature_normalization(data_dir, subj, close=True):
    """Normalize the curvature map and plot contour over fsaverage."""
    surf_dir = op.join(data_dir, subj, "surf")
    snap_dir = op.join(data_dir, subj, "snapshots")
    for hemi in ["lh", "rh"]:

        cmd = ["mri_surf2surf",
               "--srcsubject", subj,
               "--trgsubject", "fsaverage",
               "--hemi", hemi,
               "--sval", op.join(surf_dir, "%s.curv" % hemi),
               "--tval", op.join(surf_dir, "%s.curv.fsaverage.mgz" % hemi)]

        sub.check_output(cmd)

        b = Brain("fsaverage", hemi, "inflated",
                  config_opts=dict(background="white",
                                   width=700, height=500))
        curv = nib.load(op.join(surf_dir, "%s.curv.fsaverage.mgz" % hemi))
        curv = (curv.get_data() > 0).squeeze()
        b.add_contour_overlay(curv, min=0, max=1.5, n_contours=2, line_width=4)
        b.contour["colorbar"].visible = False
        for view in ["lat", "med"]:
            b.show_view(view)
            mlab.view(distance=330)
            png = op.join(snap_dir, "%s.surf_warp_%s.png" % (hemi, view))
            b.save_image(png)

        if close:
            b.close()

def savefig(filename,mlabview=[],magnification = 3):
    """
    Save mayavi figure

    Parameters
    ----------

    name : str
        name of the figure
    mlabview : [] | (x,y,z, np.array([ xroll,yroll,zroll]))
        specifyy angle of camera view ( see mayavi.view )
    magnification : int
        resolution of the generated image ( see mayavi.savefig)

    """
    import os


    path = os.path.dirname(filename)
    if not mlabview == []:
        mlab.view(mlabview)
    if os.path.exists(path):
        mlab.savefig(filename+'.png',magnification=magnification )
    else:
        os.mkdir(path)
        mlab.savefig(filename+'.png',magnification=magnification )
    mlab.close()