Period Matrices of Polyhedral Surfaces

Finding a conformal parameterization for a surface and computing its period matrix is a classical problem which is useful in a lot of contexts, from statistical mechanics to computer graphics.

The 2D-Ising model [19, 8, 9] for example takes place on a cellular decomposition of a surface whose edges are decorated by interaction constants, understood as a discrete conformal structure. In certain configurations, called critical temperature, the model exhibits conformal invariance properties in the thermodynamical limit and certain statistical expectations become discrete holomorphic at the finite level. The computation of the period matrix of higher genus surfaces built from the rectangular and triangular lattices from discrete Riemann theory has been addressed in the cited papers by Costa-Santos and McCoy.

Global conformal parameterization of a surface is important in computer graphics [16, 12, 2, 26, 17, 27] in issues such as texture mapping of a flat picture onto a curved surface in $\mathbb{R}^{3}$. When the texture is recognized by the user as a natural texture known as featuring round grains, these features should be preserved when mapped on the surface, mainly because any shear of circles into ellipses is going to be wrongly interpreted as suggesting depth increase. Characterizing a surface by a few numbers is as well a desired feature in computer graphics, for problems like pattern recognition. Computing numerically the period matrix of a surface has been addressed in the cited papers by Gu and Yau.

This paper uses the general framework of discrete Riemann surfaces theory [14, 13, 19, 4] and the computation of period matrices within this framework (based on theorems and not only numerical analogies). We translate these theorems into algorithms. Some tests are performed to check the validity of the approach.

We start with some surfaces with known period matrices and compute numerically their discrete period matrices, at different level of refinement. In particular, some genus two surfaces made out of squares and the Wente torus are considered. We observe numerically good convergence properties. Moreover, we compute the yet unknown period matrix of the Lawson surface, recognized it numerically as one of the tested surfaces, which allowed us to conjecture their conformal equivalence, and finally prove it.

Consider a polyhedral surface in $\mathbb{R}^{3}$. It has a unique Delaunay tesselation, generically a triangulation [5]. That is to say each face is associated with a circumcircle drawn on the surface and this disk contains no other vertices than the ones on its boundary. Let’s call $\Gamma$ the graph of this cellular decomposition, $\Gamma_{0}$ its vertices, $\Gamma_{1}$ its edges and complete it into a cellular decomposition with $\Gamma_{2}$ the set of triangles. Each edge $(x,x^{\prime})=e\in\Gamma_{1}$ is adjacent to a pair of triangles, associated with two circumcenters $y,y^{\prime}$. The ratio of the (intrinsic) distances between the circumcenters and the length of the (orthogonal) edge $e$ is called $\rho(e)$.

Following [19], we call this data of a graph $\Gamma$, whose edges are equipped with a positive real number a discrete conformal structure. A discrete Riemann surface is a conformal equivalence class of surfaces with the same discrete conformal structure. It leads to a theory of discrete Riemann surfaces and discrete analytic functions, developed in [14, 13, 19, 20, 21, 4, 11], that shares a lot of features with the continuous theory and these features are recovered in a proper refinement limit. We are going to summarize these results.

In our examples, the extrinsic triangulations are Delaunay. That is to say the triangulations come from the embedding in $\mathbb{R}^{3}$ and the edges $(x,x^{\prime})\in\Gamma_{1}$ of the triangulation are the edges of the polyhedral surface in $\mathbb{R}^{3}$. On the contrary, the geodesic connecting the circumcenters $y$ and $y^{\prime}$ on the surface is not an interval of a straight line and its length is generically greater than the distance $||y-y^{\prime}||$ in $\mathbb{R}^{3}$. The later gives a more naive definition of $\rho$, easier to compute that we will call extrinsinc, because contrarily to the intrinsic $\rho$, they depend on the embedding of the surface in $\mathbb{R}^{3}$ and are not preserved by isometries. For surfaces which are refined and flat enough, the difference between extrinsic and intrinsic distances is not large. We compared numerically the two ways to compute $\rho$. The conclusion is that, in the examples we tested, the intrinsic distance is marginally better, see Sec. 4.2.

The circumcenters and their adjacencies define a 3-valent abstract (locally planar) graph, dual to the graph of the surface, that we call $\Gamma^{*}$, with vertices $\Gamma^{*}_{0}=\Gamma_{2}$, edges $\Gamma^{*}_{1}\simeq\Gamma_{1}$. We equip the edge $(y,y^{\prime})=e^{*}\in\Gamma^{*}_{1}$, dual to the primal edge $e\in\Gamma_{1}$, with the positive real constant $\rho(e^{*})=1/\rho(e)$. We define $\Lambda:=\Gamma\oplus\Gamma^{*}$ the double graph, with vertices $\Lambda_{0}=\Gamma_{0}\sqcup\Gamma_{0}^{*}$ and edges $\Lambda_{1}=\Gamma_{1}\sqcup\Gamma_{1}^{*}$. Each pair of dual edges $e,e^{*}\in\Lambda_{1}$, $e=(x,x^{\prime})\in\Gamma_{1}$, $e^{*}=(y,y^{\prime})\in\Gamma^{*}_{1}$, are seen as the diagonals of a quadrilateral $(x,y,x^{\prime},y^{\prime})$, composing a quad-graph $\lozenge$, with vertices $\lozenge_{0}=\Lambda_{0}$, edges $\lozenge_{1}$ composed of couples $(x,y)$ and faces $\lozenge_{2}$ composed of quadrilaterals $(x,y,x^{\prime},y^{\prime})$.

The Hodge star, which in the continuous theory is defined by $*(f\,dx+g\,dy)=-g\,dx+f\,dy$, is in the discrete case the duality transformation multiplied by the conformal structure:

A $1$-form $\alpha\in C^{1}(\Lambda)$ is of type $(1,0)$ if and only if, for each quadrilateral $(x,y,x^{\prime},y^{\prime})\in\lozenge_{2}$, $\int_{(y,y^{\prime})}{\alpha}=i\,\rho(x,x^{\prime})\int_{(x,x^{\prime})}{\alpha}$, that is to say if $*\alpha=-i\,\alpha$. Similarly forms of type $(0,1)$ are defined by $*\alpha=+i\,\alpha$. A form is holomorphic, resp. anti-holomorphic, if it is closed and of type $(1,0)$, resp. of type $(0,1)$. A function $f:\Lambda_{0}\to\mathbb{C}$ is holomorphic iff $d_{\Lambda}f$ is.

We define a wedge product for $1$-forms living whether on edges $\lozenge_{1}$ or on their diagonals $\Lambda_{1}$, as a $2$-form living on faces $\lozenge_{2}$. The formula for the latter is:

The exterior derivative $d$ is a derivation for the wedge product, for functions $f,g$ and a $1$-form $\alpha$:

Together with the Hodge star, they give rise, in the compact case, to the usual scalar product on $1$-forms:

The adjoint $d^{*}=-*\,d\,*$ of the coboundary $d$ allows to define the discrete Laplacian $\Delta=d^{*}\,d+d\,d^{*}$, whose kernel are the harmonic forms and functions. It reads, for a function at a vertex $x\in\Lambda_{0}$ with neighbours $x^{\prime}\sim x$:

Hodge theorem: The two $\pm i$-eigenspaces decompose the space of $1$-forms, especially the space of harmonic forms, into an orthogonal direct sum. Types are interchanged by conjugation: $\alpha\in C^{(1,0)}(\Lambda)\iff\overline{\alpha}\in C^{(0,1)}(\Lambda)$ therefore

where the projections on $(1,0)$ and $(0,1)$ spaces are

The harmonic forms of type $(1,0)$ are the holomorphic forms, the harmonic forms of type $(0,1)$ are the anti-holomorphic forms.

The $L^{2}$ norm of the $1$-form $df$, called the Dirichlet energy of the function $f$, is the average of the usual Dirichlet energies on each independent graph

The conformal energy of a map measures its conformality defect, relating these two harmonic functions. A conformal map fulfills the Cauchy-Riemann equation

Therefore a quadratic energy whose null functions are the holomorphic ones is

It is related to the Dirichlet energy through the same formula as in the continuous case:

where the area of the image of the application $f$ in the complex plane has the same formulae (the second one meaningful on a simply connected domain)

as in the continuous case. For a face $(x,y,x^{\prime},y^{\prime})\in\lozenge_{2}$, the algebraic area of the oriented quadrilateral $\Bigl{(}f(x),f(x^{\prime}),f(y),f(y^{\prime})\Bigr{)}$ is given by

When a holomorphic reference map $z:\Lambda_{0}\to\mathbb{C}$ is chosen, a holomorphic (resp. anti-holomorphic) $1$-form $df$ is, locally on each pair of dual diagonals, proportional to $dz$, resp. $d\bar{z}$, so that the decomposition of the exterior derivative into holomorphic and anti-holomorphic parts yields $df\wedge\overline{df}=\left(|\partial f|^{2}+|\bar{\partial}f|^{2}\right)dz\wedge d\bar{z}$ where the derivatives naturally live on faces and are not be confused with the boundary operator $\partial$.

The theory described above is straightforward to implement. The most sensitive part is based on a minimizer procedure which finds the minimum of the Dirichlet energy for a discrete Riemann surface, given some boundary conditions. Here is the crude algorithm that we are going to detail. A normalized homotopy basis $\aleph$ of a quad-graph $\lozenge(S)$ is a set of loops $\aleph_{k}\in\ker\partial$

Finding a normalized homotopy basis of a connected cellular decomposition is performed by several well known algorithms. The way we did it is to select a root vertex and grow from there a spanning tree, by computing the vertices at combinatorial distance $d$ from the root and linking each one of them to a unique vertex at distance $d-1$, already in the tree. Repeat until no vertices are left.

Then we inflate this tree into a polygonal fundamental domain by adding faces one by one to the domain, keeping it simply connected: We recursively add all the faces which have only one edge not in the domain. We stop when all the remaining faces have at least two edges not in the domain.

Then we pick one edge (one of the closest to the root) as defining the first element of our homotopy basis: adding this edge to the fundamental domain yields a non simply connected cellular decomposition and the spanning tree gives us a rooted cycle of this homotopy type going down the tree to the root. It is (one of) the combinatorially shortest in its (rooted) homology class. We add faces recursively in a similar way until we can no go further, we then choose a new homotopy basis element, and so on until every face is closed. At the end we have a homotopy basis. We compute later on the intersection numbers in order to normalize it.

\begin{picture}(3229.0,3376.0)(247.0,-2615.0)\end{picture}\begin{picture}(3229.0,3376.0)(247.0,-2615.0)\end{picture}\begin{picture}(3229.0,3376.0)(247.0,-2615.0)\end{picture}

We compute the unique real harmonic form $\eta$ associated with each cycle $\aleph$ such that $\oint_{\gamma}\eta=\gamma\circ\aleph$. This is done by a minimizing procedure which finds the unique harmonic function $f$ on the graph $\Gamma$, split along $\aleph$, whose vertices are duplicated, which is zero at the root and increases by one when going across $\aleph$. This is done by linking the values at the duplicated vertices, in effect yielding a harmonic function on the universal cover of $\Gamma$. The harmonic $1$-form $d\,f$ doesn’t depend on the chosen root nor on the representative $\aleph$ in its homology class.

We begun with testing discrete surfaces of known moduli in order to investigate the quality of the numerics and the robustness of the method. We purposely chose to stick with raw double $15$-digits numbers and a linear algebra library which is fast but not particularly accurate. In order to be able to compare period matrices, we used a Siegel reduction algorithm [10] to map them by a modular transformation to the same canonical form.