def view_thresholdedT(design, contrast, threshold, inequality=np.greater):
    """
    A mayavi isosurface view of thresholded t-statistics

    Parameters
    ----------
    design : {'block', 'event'}
    contrast : str
    threshold : float
    inequality : {np.greater, np.less}, optional
    """
    maska = np.asarray(MASK)
    tmap = np.array(load_image_fiac('group', design, contrast, 't.nii'))
    test = inequality(tmap, threshold)
    tval = np.zeros(tmap.shape)
    tval[test] = tmap[test]

    # XXX make the array axes agree with mayavi2
    avganata = np.array(AVGANAT)
    avganat_iso = mlab.contour3d(avganata * maska, opacity=0.3, contours=[3600],
                               color=(0.8,0.8,0.8))

    avganat_iso.actor.property.backface_culling = True
    avganat_iso.actor.property.ambient = 0.3

    tval_iso = mlab.contour3d(tval * MASK, color=(0.8,0.3,0.3),
                            contours=[threshold])
    return avganat_iso, tval_iso

def test3():
    import numpy as np
    from mayavi import mlab
    sdf = np.load('/Users/yue/data/segment/sdf_diff_3d_1/imseg_iter_100.npz')['sdf']
    mlab.contour3d(sdf[0], contours=[0])
    # mlab.savefig('/Users/yue/data/segment/sdf_diff_3d_1/imseg_iter_100.png')
    mlab.show()

def visual_test1():
    
    poly = polygon_gaussian_density(3)
    #s = mlab.contour3d(poly)
    #mlab.show()
    
    grid_dimesions = poly.shape
    grid_center = np.array(grid_dimesions) / 2.0
    radii = np.arange(2.0, 10.0)
    L_max = 25
    
    shg = SphHarmGrid(grid_dimesions, grid_center, radii, L_max)
    
    coeff = shg(poly)
    reconstructed_poly = shg.expand_sph_harm(coeff)
    
    print 'magnitudes'
    print 'real:', np.sum(np.abs(reconstructed_poly.real)), 'imag:', np.sum(np.abs(reconstructed_poly.imag))
    
    poly  /= poly.sum()
    rpoly  = np.abs(reconstructed_poly)
    rpoly /= rpoly.sum()
    
    err = np.sum( np.abs(poly - rpoly) / float(np.product(grid_dimesions)) )
    print "Reconstruction error:", err
    
    print 'showing model'
    s1 = mlab.contour3d(rpoly)
    mlab.show()
    
    print 'showing error'
    s2 = mlab.contour3d(rpoly - poly)
    mlab.show()
    
    return

def vis_breakups():
    """
    This function allows to visualize break-ups as the highlighted branches in
    the original skeleton.
    """
    data_sk = glob.glob('/backup/yuliya/vsi05/skeletons_largdom/*.h5')
    data_sk.sort()
    
    data_br = glob.glob('/backup/yuliya/vsi05/breakups_con_correction/dict/*.npy')
    data_br.sort()
    
    for i,j in zip(data_sk[27:67][::-1], data_br[19:][::-1]):
        d = tb.openFile(i, mode='r')
        bran1 = np.copy(d.root.branches)
        skel1 = np.copy(d.root.skel)
        d.close() 

        br = np.load(j).item()
        
        mlab.figure(1, bgcolor=(1,1,1), size=(1200,1200))
        surf_skel = mlab.contour3d(skel1, colormap='hot')
        surf_skel.actor.property.representation = 'points'
        surf_skel.actor.property.line_width = 2.3
        mask = np.zeros_like(skel1)
        for k in br:
            sl = br[k]            
            mask[sl] = (bran1[sl] == k)
        surf_mask = mlab.contour3d(mask.astype('uint8'))
        surf_mask.actor.property.representation = 'wireframe'
        surf_mask.actor.property.line_width = 5
        
        mlab.savefig('/backup/yuliya/vsi05/breakups_con_correction/video_hot/' + i[39:-3] + '.png')
        mlab.close()

 def showContour(self):
     """
     Display real space model using contour
     """
     mlab.contour3d(self.array, contours=8, opacity=0.3, transparent=True)
     mlab.outline()
     mlab.show()

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_LAM_contours(self, readfile):
		if self.LAMgauss == None:
	 		print "Reading file..."
	 		self.LAMgauss = self.load_data(readfile)

	 	print "Plotting..."
		mlab.contour3d(self.LAMgauss, contours = 75, opacity=0.10)
		mlab.show()

def show_scatterer(scatterer, spacing=None):
    mlab = import_mayavi()
    if spacing == None:
        # if no spacing given, represent the object with 100 voxels
        # along each dimension
        spacing = [(b[1]-b[0])/100 for b in scatterer.bounds]
    vol = scatterer.voxelate(spacing, 0)
    mlab.contour3d(vol)

 def plotParcels(self, brain, label = None, contourVals = [], opacity = 0.5, cmap='autumn'):
     ''' plot an isosurface using Mayavi, almost the same as skull plotting '''
     
     if not(label):
         label = self.getAutoLabel()        
     
     if contourVals == []:            
         self.isosurfacePlots[label] = mlab.contour3d(brain.parcels, opacity = opacity, colormap=cmap)
     else:
         self.isosurfacePlots[label] = mlab.contour3d(brain.parcels, opacity = opacity, contours = contourVals, colormap=cmap)

def plot(ply):
    '''
    Plot vertices and triangles from a PlyData instance. Assumptions:
        `ply' has a 'vertex' element with 'x', 'y', and 'z'
            properties;

        `ply' has a 'face' element with an integral list property
            'vertex_indices', all of whose elements have length 3.

    '''
    vertex = ply['vertex']

    (x, y, z) = (vertex[t] for t in ('x', 'y', 'z'))

    # mlab.points3d(x, y, z, color=(1, 1, 1), mode='point')

    if 'face' in ply:
        tri_idx = ply['face']['vertex_indices']
        print type(tri_idx)
        idx_dtype = tri_idx[0].dtype
        print "idx_dtype" , idx_dtype
        triangles = numpy.fromiter(tri_idx, [('data', idx_dtype, (3,))],
                                   count=len(tri_idx))['data']
        k = triangles.byteswap()
        # p.from_array(numpy.array([[1,2,3],[3,4,5]], dtype=numpy.float32))
        p.from_array((triangles.byteswap().newbyteorder().astype('<f4')))
        # p1.add_points_from_input_cloud()
        p1.set_input_cloud(p)
        fil = p.make_statistical_outlier_filter()
        # cloud_filtered = fil.filter()
        # fil = p1.make_statistical_outlier_filter()
        fil.set_mean_k(20)
        fil.set_std_dev_mul_thresh(0.1)
        pointx = triangles[:,0]

       

        pc_1 = pcl.PointCloud()
        pc_1 = p
        pc_2 = pcl.PointCloud()
        pc_2 = p
        kd = pcl.KdTreeFLANN(pc_1)
        print('pc_1:')
        print type(triangles)
        # file1 = open("values.txt",'w')
        # file1.write(str(triangles[:]))

        # mlab.mesh(x,y,z)
        print numpy.shape(p)
        print numpy.shape(z)
        print numpy.shape(y)
        print numpy.shape(x)
        mlab.contour3d(x, y, z, p)

 def plotIsosurface(self, brain, label = None, contourVals = [], opacity = 0.1, cmap='autumn'):
     ''' plot an isosurface using Mayavi, almost the same as skull plotting '''
     
     if not(label):
         label = self.getAutoLabel()        
     
     if contourVals == []:            
         s = mlab.contour3d(brain.iso, opacity = opacity, colormap=cmap)
     else:
         s = mlab.contour3d(brain.iso, opacity = opacity, contours = contourVals, colormap=cmap)
         
     # get the object for editing
     self.isosurfacePlots[label] = s

	def show_fitting(self, depth_image, offsets, model_transform):
		v, u = np.nonzero(depth_image != 0)
		depth = depth_image[(v,u)]/1000.0
		
		points = np.vstack([depth*(u+offsets[0]), depth*(v+offsets[1]), depth])
		
		model_points = np.dot(d.K_default_inv, points)
		model_points = np.vstack([model_points, np.ones_like(model_points[0])])
		transformed_points = np.dot(inv(model_transform), model_points)[0:3,:]
		
		mlab.contour3d(self.x/self.scale[0], self.y/self.scale[1], self.z/self.scale[2], self.sdf, contours=[0.05], opacity=0.5)
		mlab.points3d(transformed_points[0], transformed_points[1], transformed_points[2], color=(1,0,0), mode='point')
		mlab.show()

 def plotBackground(self, brain, label = None, contourVals = [3000,9000], opacity = 0.1, cmap='Spectral'):
     ''' plot the skull using Mayavi '''
     
     if not(label):
         label = self.getAutoLabel()        
     
     if contourVals == []:            
         s = mlab.contour3d(brain.background, opacity = opacity, colormap=cmap)
     else:
         s = mlab.contour3d(brain.background, opacity = opacity, contours = contourVals, colormap=cmap)
         
     # get the object for editing
     self.skullPlots[label] = s

    def plotSkull(self, brain, label = None, contourVals = [], opacity = 0.1, cmap='Spectral'):
        ''' plot the skull using Mayavi '''
      
        if not(label):
            label = self.getAutoLabel()        

        # remove old version        
        if not(label in self.skullPlots):
            self.skullPlots[label].remove()
        
        if contourVals == []:            
            self.skullPlots[label] = mlab.contour3d(brain.background, opacity = opacity, colormap=cmap)
        else:
            self.skullPlots[label] = mlab.contour3d(brain.background, opacity = opacity, contours = contourVals, colormap=cmap)

def plot_volume(nifti_file,v_color=(.98,.63,.48),v_opacity=.1, fliplr=False, newfig=False):
    """
    Render a volume from a .nii file
    Use fliplr option if scan has radiological orientation
    """
    import nibabel as nib
    if newfig:
        mlab.figure(bgcolor=(1, 1, 1), size=(400, 400))
    input = nib.load(nifti_file)
    input_d = input.get_data()
    d_shape = input_d.shape
    if fliplr:
        input_d = input_d[range(d_shape[0]-1,-1,-1), :, :]
    mlab.contour3d(input_d, color=v_color, opacity=v_opacity) # good natural color is (.98,.63,.48)

def calc():
    """
    The function that implements break-ups detection algorithm and
    write the data in a file. Input skeletons and graphs must correspond each
    other.
    """
    
    #input graphs
    data_gr = glob.glob('/path/graphs_largedom/*_1.gpickle')
    data_gr.sort()

    #input skeletons
    data_sk = glob.glob('/path/skeletons_largedom/*.h5')
    data_sk.sort()

    for i, j, w in zip(data_gr[:-1][::-1], data_sk[:-1][::-1], data_sk[1:][::-1]):
        g = nx.read_gpickle(i)

        #We simplify the graph
        g1 = remove_terminal_br(g)
        g1 = remove_nodes(g1)
        
        #We use labeled branches arrays, skeletons, distans transfrom. All are of the same shape.
        d = tb.openFile(j, mode='r')
        bran1 = np.copy(d.root.branches1)
        skel1 = np.copy(d.root.skel) #this 'current' skeletons are used only for visualization.
        thick = np.copy(d.root.dist)
        d.close()
        
        #Skeletons shifted in time by one step.
        d = tb.openFile(w, mode='r')
        skel2 = np.copy(d.root.skel)
        d.close()
        
        b, pos = branches_map(g1, skel2, bran1, thick)
        
        #Visualization of break-ups if necessary.
        mlab.figure(1, bgcolor=(1,1,1), size=(1500,1500))
        mlab.contour3d(skel1, colormap='hot')
        mask = np.zeros_like(skel1)
        for k in b:
            sl = b[k]            
            mask[sl] = (bran1[sl] == k)
        mlab.contour3d(mask.astype('uint8'))
        
        np.save('/path/breakups/' + j[41:-3] + '_1.npy', b)  #j[41:-3] must be replaced according to the path
        np.save('/path/breakups/' + j[41:-3] + '_pos.npy', pos)
        mlab.savefig('/path/breakups/' + j[41:-3] + '_1.png')
        mlab.close()

def make3dHist(data):
    # n, m, l = data.shape
    # results = np.zeros((n*l,3))
    # for i in range(n):
    #     for j in range(m):
    #         results[i*l:(i+1)*l,j] = data[i,j,:]

    # H, edges = np.histogramdd(results, bins = (100,100,100))
    # print H.shape
    # data = earthSpherify(data)
    mlab.contour3d(data)
    data = data / np.max(data) * 3.74
    print "Still fine"
    mlab.pipeline.volume(mlab.pipeline.scalar_field(data))
    mlab.show()

    def plotIsosurface(self, brain, label=None, contourVals=[], opacity=0.1, cmap="autumn"):
        """ plot an isosurface using Mayavi, almost the same as skull plotting """

        if not (label):
            label = self.getAutoLabel()

        if label in self.isosurfacePlots:
            self.isosurfacePlots[label].remove()

        if contourVals == []:
            self.isosurfacePlots[label] = mlab.contour3d(brain.iso, opacity=opacity, colormap=cmap)
        else:
            self.isosurfacePlots[label] = mlab.contour3d(
                brain.iso, opacity=opacity, contours=contourVals, colormap=cmap
            )

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 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()