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