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