Examples using pymstri

This module is probably offers the most flexibility for those hoping to use this software for data analysis and visualization. Support for python scripting and plugins are very common in large Visualization software packages like VTK , paraview, etc. , or even domain specific tools such as PyMol or Chimera. Thus, adding Morse-Smale complex based analysis capablity to these tools becomes very easy with the python modules.

The following are some basic python scripts to compute the Morse-Smale complex to do some analysis. The code is commented, so it should be quite easy to follow. The scripts and the data are available here. Unpack them in the install dir and execute the script.

  • Computing Homology

    #import necessary modules
    import pymstri
    
    #create the mscomplex object. Call help(msc) for info. 
    msc = pymstri.mscomplex()
    
    # compute the Morse-Smale complex.
    # Scalar function is given in the 4th field of the off file.
    msc.compute_off("3wde.off")
    
    # simplify maximally using persistence
    msc.simplify_pers(1.0,True,0,0)
    
    # get the infinite persistence cps
    minima  = [cp for cp in msc.cps() if msc.index(cp) == 0]
    saddles = [cp for cp in msc.cps() if msc.index(cp) == 1]
    maxima  = [cp for cp in msc.cps() if msc.index(cp) == 2]
    
    # print the homology
    print "# components  : ", len(minima)
    print "# tunnels     : ", len(saddles)
    print "# voids       : ", len(maxima)
    
  • Compute Extremum Graph

    # In this script, we will extract the an extremum graph from the 
    # combinatorial Morse-Smale complex.
    #
    # The extremum graph will consist of only maxima and 1-saddles. 
    #
    # The extrema that are retained are those that have persistence above 5% 
    #
    # We will also generate a persistence hierarchy to figure out 
    # persistence pairs. 
    #
    # We will retain only those 1-saddles a persistence pair with a maxima
    # and infinite persistent 1-saddles. 
    
    #import necessary modules
    import pymstri
    
    #create the mscomplex object. Call help(msc) for info. 
    msc = pymstri.mscomplex()
    
    # compute the Morse-Smale complex. 
    # Scalar function is given in the 4th field of the off file.
    msc.compute_off("3wde.off")
    
    # simplify using persistence upto 5% of the function range
    msc.simplify_pers(0.05,True,0,0)
      
    # Generate the persistence hierarchy
    msc.gen_pers_hierarchy()
    
    # get the set of persistence pairs and form a dictionary for lookups
    ppairs     = [msc.canc(i) for i in range(msc.num_canc())]
    ppairsDict = dict(ppairs + [(b,a) for a,b in ppairs])
    
    # get the trianulation object for vertex information
    tcc        = msc.get_tri_cc()
    
    # compute a list of cps that fulfill the above criteria
    cpset = [ cp for cp in msc.cps() 
    		if msc.index(cp) == 2 or 
    		(msc.index(cp) == 1 and 
    			(cp not in ppairsDict or 
    			msc.index(ppairsDict[cp]) == 2) )]
    
    # Write the extremum graph data to a file  
    fo = open("3wde_extremumGraph.txt","w")
    
    # First, write the number of critical points
    fo.write("#numcps\n")
    fo.write(str(len(cpset))+"\n")
    
    # Write the per cp information. Note: Each cp is identified by a unique id.
    fo.write("#    id       x         y         z         fn        index\n")
    for cp in msc.cps():
    	x,y,z = tcc.point(msc.vertid(cp))
    	fn,index = msc.fn(cp),msc.index(cp)
    	fo.write("%5d   %8.4f  %8.4f  %8.4f  %8.6f  %5d\n"%(cp,x,y,z,fn,index))
    
    # Write the combinatorial adjaceny list of each cp
    fo.write("#  id  NumAdj adjCp1 adjCp2 ..\n")
    for cp in msc.cps():
    	# filter connections to retain only relevent connections
    	conn = msc.asc(cp) if msc.index(cp) == 1 else msc.des(cp)    
    	conn = [i for i in conn if i in cpset]
    	
    	ln  = ""
    	ln += str(cp).ljust(10)
    	ln += (str(len(conn))).ljust(10)
    	ln += "".join(str(i).ljust(10) for i in conn)
    	ln += "\n"
    
    	fo.write(ln)
    
    #finished
    fo.write("# The end")
    
  • Segment Protrusions on Mol surface and Save to Vtk File

    # In this script, we will compute a segmentation of a molecular surface.
    #
    # Significant protrusions are identified and segmented. 
    #
    # The segments are then color tagged in a surface file that is 
    # visualizable in vtk/paraview. 
    
    #import necessary modules
    import pymstri,vtk,random
    
    #create the mscomplex object. 
    msc = pymstri.mscomplex()
    
    # compute the Morse-Smale complex. 
    # Scalar function is given in the 4th field of the off file.
    msc.compute_off("3wde.off")
    
    # simplify using persistence upto 5% of the function range
    msc.simplify_pers(0.05,True,0,0)
    
    # collect the msc geometry after simplification
    msc.collect_geom()
    
    # get the object representing the triangulation. 
    tcc = msc.get_tri_cc()
    
    #get number of cells of each type
    nv,ne,nt = tcc.num_dcells(0),tcc.num_dcells(1),tcc.num_dcells(2)
    
    #create and insert the points
    pts = vtk.vtkPoints()
    for cellid in range(0,nv):
    	pts.InsertNextPoint(tcc.point(cellid))
    	
    #create a vtk_poly_data
    pd = vtk.vtkPolyData()
    pd.Allocate()
    
    #set its pts
    pd.SetPoints(pts)
    
    #iterate over triangle cells  and insert the polys
    tcell = vtk.vtkIdList()
    tcell.SetNumberOfIds(3)
    
    for cellid in range(nv+ne,nv+ne+nt):
    	i,j,k = tcc.cell_vert(cellid)
    	tcell.SetId(0,i)
    	tcell.SetId(1,j)
    	tcell.SetId(2,k)	
    	pd.InsertNextCell(vtk.VTK_TRIANGLE,tcell)
    	
    
    #create a color array
    colors = vtk.vtkDoubleArray()
    colors.SetNumberOfComponents(3)
    colors.SetNumberOfTuples(nt)
    colors.SetName("Color")
    
    # get a list of surviving maxima
    survMax = [cp for cp in msc.cps() if msc.index(cp) == 2]
    
    # iterate over the des manifold of each survigin max. 
    for cp in survMax:
    	#pick a random color
    	color = random.random(),random.random(),random.random()
    
    	# The des mfold comprises of a set of triangles. 	
    	for cellid in msc.des_geom(cp):
    		#triangle cellids are offset by nv+ne in tcc
    		tid = cellid - nv - ne 
    		colors.SetTuple(tid,color)
    
    #Add the data as cell data
    pd.GetCellData().AddArray(colors)
    
    
    #write the data to a vtk file
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetFileName( "3wde_maxSeg.vtp" )
    writer.SetInput( pd)
    writer.Write()