In this Python 2 notebook we use MOSEK Fusion to solve some geometric optimization problems on surfaces. The code demonstrates setting up a linear optimization problem, extending the Model class, reoptimization with different parameters and shows the practical difference between interior-point and basic solution.

`plyfile` is a package for reading the descriptions of surfaces in PLY format. It can be installed with pip install plyfile.

You will also need some data files with examples used in this presentation. The complete code and additional files are available from GitHub.

In [15]:

#Fetching examples predefined for this presentationdefgetExample(name):ifname=='mug.ply':fname='http://people.sc.fsu.edu/~jburkardt/data/ply/mug.ply'else:fname='https://raw.githubusercontent.com/MOSEK/Tutorials/master/Fusion/surfacecycles/'+namereturnplyfile.PlyData.read(urllib2.build_opener().open(urllib2.Request(fname)))defgetCycleExample(name):returnnp.load(np.DataSource().open('https://raw.githubusercontent.com/MOSEK/Tutorials/master/Fusion/surfacecycles/'+name))

One of the models used in this tutorial comes from John Burkardt's library of PLY files and the other ones were generated using CGAL. A comprehensive survey of topological optimization problems can be found in the survey paper by Jeff Erickson. The model used in this tutorial comes from the paper Optimal Homologous Cycles, Total Unimodularity, and Linear Programming by Dey, Hirani and Krishnamoorthy. The experts in the field are kindly asked to fill in some vague statements in the following text with the precise algebro-topological language on their own.

A triangulated surface is a mesh of 2D triangles in 3D that fit together to form a surface. The vertices and edges of a triangulation form a graph. The cycles in that graph can be used to detect interesting topological features. For an example of what that means, see the three loops on the surface below:

Now imagine that each loop is a rubber band that can be stretched in a continuous way and slide on the surface. Then the violet loop is not very interesting - it lies in a small patch of the surface and if we start shrinking it, it will collapse to a point (in other words, it is contractible). The yellow loop is different - it circles around a hole in the surface, and no matter how much we deform it, it will not disappear. We say this loop detects a feature - in this case the feature is one of the holes. The red loop has the same property. Moreover, the yellow and red loops are homotopic - each of them can be continuosly moulded into the second one - i.e. they detect the same feature (hole).

In some applications it is useful to have small representations of topological features, in this case the shortest possible loops in a given homotopy class. There are at least two reasonable ways of measuring the length of a path on a surface: with the Euclidean metric on the edges, or just counting the number of edges.

Suppose the surface has $v$ vertices, $e$ edges and $t$ triangles. It is convenient to represent it by boundary matrices, which are just variants of the incidence matrix of a graph. Specifically, let $V(i)$, $E(ij)$ and $T(ijk)$ denote the vertex, edge and triangle with the given vertex set. Then we have a linear map $D_1:\mathbb{R}^e\to\mathbb{R}^v$ given on the basis vectors by:

$$D_1(\mathbf{e}_{E(ij)})=\mathbf{e}_{V(j)}-\mathbf{e}_{V(i)}$$

and we identify $D_1\in \mathbb{R}^{v\times e}$ with the matrix of this map. Similarly, we can record the incidence between triangles and edges by $D_2:\mathbb{R}^t\to\mathbb{R}^e$:

and let $D_2\in \mathbb{R}^{e\times t}$ be its matrix. We refer to other sources for more details. Note that matrices $D_1,D_2$ are extremely sparse - they have just two, resp. three, nonzeros in each column.

A chain (more precisely, 1-chain) is just a vector $c\in \mathbb{R}^e$, that is an assignment of a real coefficient to every edge. If $c$ is an actual path or cycle in the graph, then these coefficients are $\pm1$ for edges in the cycle, and $0$ otherwise, but more general chains are allowed. We will mostly concentrate on cycles.

We can now formulate the following problem: given a cycle $c\in\mathbb{R}^e$, find the shortest cycle that detects the same topological features as $c$ (strictly speaking, the shortest cycle homologous to $c$). Suppose $w_1,\ldots,w_e$ are the weights (lengths) of the edges.

This optimization problem can easily be formulated as a linear program. Below is the implementation of this program in a class Surface, which extends a Mosek Fusion Model.

In [17]:

#A model of a surface (or a more general 2D complex in 3D)classSurface(mosek.fusion.Model):#Initialize from plydata using edge weights given by metricdef__init__(self,plydata,metric):super(Surface,self).__init__()### Initialize basic dataself.v=plydata.elements[0].count# number of verticesself.t=plydata.elements[1].count# number of trianglesself.coo=plydata['vertex']# coordinates of vertices in 3Dself.xc,self.yc,self.zc=zip(*self.coo)self.tri=[tuple(sorted(f[0]))forfinplydata['face']]# list of triangles, as sorted triples of verticesself.edg=list(set([tuple(sorted([f[i],f[j]]))# list of edges, as sorted pairs of verticesforfinself.trifori,jin[(0,1),(0,2),(1,2)]]))self.e=len(self.edg)# number of edgesself.wgh=[metric(self.coo[f[0]],self.coo[f[1]])# weights of edgesforfinself.edg]### Homological algebra# Construct the sparse matrix D2d1=[(self.edg[i][k],i,sgn)foriinrange(self.e)fork,sgnin[(1,+1),(0,-1)]]s1,s2,val=zip(*d1)self.D1=Matrix.sparse(self.v,self.e,list(s1),list(s2),list(val))# Construct the sparse matrix D3d2=[(self.edg.index((self.tri[i][k],self.tri[i][l])),i,sgn)foriinrange(self.t)fork,l,sgnin[(1,2,+1),(0,2,-1),(0,1,+1)]]s1,s2,val=zip(*d2)self.D2=Matrix.sparse(self.e,self.t,list(s1),list(s2),list(val))### Set up a MOSEK Fusion model for solving shortest homologous cyclesx=self.variable(self.e,Domain.unbounded())# the chain we are looking forxabs=self.variable(self.e,Domain.unbounded())# absolute value of xy=self.variable(self.t,Domain.unbounded())# a 2-chain whose boundary is x-CC=self.variable(self.e,Domain.unbounded())# the input chain# -xabs <= x <= xabsself.constraint(Expr.add(xabs,x),Domain.greaterThan(0))self.constraint(Expr.sub(xabs,x),Domain.greaterThan(0))# x - c = D2*y, so x-c is a boundaryself.constraint(Expr.sub(Expr.sub(x,C),Expr.mul(self.D2,y)),Domain.equalsTo(0))# min weighted 1-norm of xself.objective(ObjectiveSense.Minimize,Expr.dot(self.wgh,xabs))# Save the model for later cloningself.x,self.c,self.C=x,self.constraint(C,Domain.equalsTo(0)),0#Plot the surface and a number of chainsdefplot(self,chains=[],colors=[]):fig=plt.figure()ax=fig.gca(projection='3d')ax.plot_trisurf(self.xc,self.yc,self.zc,triangles=self.tri,linewidth=0.5,alpha=1.0,edgecolor='green')i=0forCinchains:forjinrange(self.e):ifabs(C[j])>0.0001:# Extract the edge and plot with intensity depending on its coefficient in Cxc,yc,zc=zip(*[self.coo[self.edg[j][0]],self.coo[self.edg[j][1]]])ax.plot(xc,yc,zc,color=colors[i],alpha=min(1,abs(C[j])))i+=1plt.axis('off')plt.show()#Find the shortest chain homologous to a given chain Cdefshort(self,C):# We add a constraint c=C to the model# This is a little hack, because we must remember the old Cself.c.add(self.C-C)self.C=Cself.solve()returnself.x.level()

The inititlization part processes the PLY description of the surface, initializes the sparse boundary matrices and sets up a linear program with the parameter $c$ as a variable. In subsequent optimizations (the method short) $c$ is set to the actual input cycle.

OK, but what if we don't have any predefined cycle, and we would simply like to visualize all topological features on a given input surface using shortest possible loops? In general, this is a hard problem, for which some algorithms are known on surfaces. However, we can use our model to try a more pedestrian approach, which is not perfect, but often good enough for visualization.

Without going into too much algebra, the space of features on the surface can be identified with the nullspace of the matrix $$A=\left[\begin{array}{l} D_1 \\ D_2^T\end{array}\right]$$ (for the experts, we are now computing first homology with real coefficients). This can be easily found using singular value decomposition:

In [19]:

#Find the basis for the nullspace of a matrixdefnullspace(A):u,s,vh=np.linalg.svd(A)nnz=(s>=1e-13).sum()ns=vh[nnz:].conj().Treturnns,ns.shape[1]defhomology(S):returnnullspace(np.reshape(np.concatenate((S.D1.getDataAsArray(),S.D2.transpose().getDataAsArray())),[S.v+S.t,S.e]))

For an orientable surface the dimension of the nullspace of $A$ equals $2g$, where $g$ is the genus, or the number of holes. For example, the torus has one hole and two essential features: the parallel and the meridian cycle.

Since the SVD decomposition is oblivious to the underlying geometry, the chains it returns are typically dense vectors with real coordinates, not very useful for visualization. However, we can replace them with their shortest equivalents and see if they are any better. It turns out in practice they often are. To make this easier to see we also normalize so that the maximal coefficient of each chain is 1.

In [20]:

#Normalize a chain so that l_\infty norm is 1defnormalize(C):returnC/max(abs(C))S=Surface(getExample('torus.ply'),euclideanM)generators,betti1=homology(S)G=[normalize(generators[:,i])foriinrange(betti1)]GS=[normalize(S.short(g))forginG]foriinrange(betti1):S.plot([G[i],GS[i]],['yellow','red'])