# Geometric Registration of High-genus Surfaces

This paper presents a method to obtain geometric registrations between high-genus (g≥ 1) surfaces. Surface registration between simple surfaces, such as simply-connected open surfaces, has been well studied. However, very few works have been carried out for the registration of high-genus surfaces. The high-genus topology of the surface poses great challenge for surface registration. A possible approach is to partition surfaces into simply-connected patches and registration is done patch by patch. Consistent cuts are required, which are usually difficult to obtain and prone to error. In this work, we propose an effective way to obtain geometric registration between high-genus surfaces without introducing consistent cuts. The key idea is to conformally parameterize the surface into its universal covering space, which is either the Euclidean plane or the hyperbolic disk embedded in R^2. Registration can then be done on the universal covering space by minimizing a shape mismatching energy measuring the geometric dissimilarity between the two surfaces. Our proposed algorithm effectively computes a smooth registration between high-genus surfaces that matches geometric information as much as possible. The algorithm can also be applied to find a smooth and bijective registration minimizing any general energy functionals. Numerical experiments on high-genus surface data show that our proposed method is effective for registering high-genus surfaces with geometric matching. We also applied the method to register anatomical structures for medical imaging, which demonstrates the usefulness of the proposed algorithm.

## Authors

• 4 publications
• 25 publications
• ### A Novel Stretch Energy Minimization Algorithm for Equiareal Parameterizations

Surface parameterizations have been widely applied to computer graphics ...
08/05/2017 ∙ by Mei-Heng Yueh, et al. ∙ 0

• ### Flipping Geometric Triangulations on Hyperbolic Surfaces

We consider geometric triangulations of surfaces, i.e., triangulations w...
12/10/2019 ∙ by Vincent Despré, et al. ∙ 0

• ### 3D Registration of Curves and Surfaces using Local Differential Information

This article presents for the first time a global method for registering...
04/02/2018 ∙ by Carolina Raposo, et al. ∙ 0

• ### Cortical surface registration using unsupervised learning

Non-rigid cortical registration is an important and challenging task due...
04/09/2020 ∙ by Jieyu Cheng, et al. ∙ 0

• ### Inconsistent Surface Registration via Optimization of Mapping Distortions

We address the problem of registering two surfaces, of which a natural b...
08/24/2019 ∙ by Di Qiu, et al. ∙ 0

• ### A Subpixel Registration Algorithm for Low PSNR Images

This paper presents a fast algorithm for obtaining high-accuracy subpixe...
03/31/2018 ∙ by Song Feng, et al. ∙ 0

• ### Surfaces Representation with Sharp Features Using Sqrt(3) and Loop Subdivision Schemes

This paper presents a hybrid algorithm that combines features form both ...
02/10/2014 ∙ by Yasser M. Abd El-Latif, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Registration refers to the process of finding an optimal one-to-one correspondence between images or surfaces. It has been extensively applied to different areas such as medical imaging, computer graphics and computer visions. For example, in medical imaging, registration is always needed for statistical shape analysis, morphometry and processing of signals on brain surfaces (e.g., denoising or filtering). While in computer graphics, surface registration is needed for texture mapping, which aligns each vertex to a position of the texture image, to improve the visualization of the surface mesh. Developing an effective algorithm for registration is therefore very important.

Surface registration between simple surfaces, such as simply-connected open surfaces or genus-0 closed surfaces, has been extensively studied. A lot of effective algorithms have been proposed. However, as far as we know, very few literatures have been reported on the registration between high-genus surfaces. The high-genus topology of the surfaces poses a great challenge to register the surfaces. A possible approach to cope with high-genus surface registration is by introducing cuts to partition the surface into several simply-connected patches. Registration can then be carried out in a patch-by-patch manner. As a result, consistent cuts are required, which is usually difficult to locate and prone to error. Motivated by this, we are interested in developing a geometric registration algorithm for high-genus surfaces, which does not involve the introduction of boundary cuts.

In this paper, we propose an effective way to obtain registrations between high-genus surfaces without introducing any cuts, which matches the geometry as much as possible. The key idea is to conformally parameterize the surface into its universal covering, which is either the 2D Euclidean plane or the hyperbolic disk , using the discrete Ricci flow method. Registration can then be done on the universal covering space by minimizing a shape mismatching energy measuring the geometric dissimilarity between the surfaces. Our proposed algorithm effectively computes a smooth registration between high-genus surfaces that matches geometric information as much as possible. To test the performance of the proposed method, numerical experiments have been done on synthetic high-genus surface data. Results show that our proposed algorithm is effective in registering high-genus surfaces with complete geometric matching. The proposed method has also been applied to register anatomical structures for medical imaging, which demonstrates the usefulness of the proposed algorithm.

The rest of the paper is organized as follows. In section 2, we describe some previous works closely related to our paper. In section 3, we introduce some basic mathematical concepts. The proposed algorithm for high-genus surface registration is explained in detail in section 4. The detailed numerical implementation of the algorithm will be described in section 5. In section 6, we show the numerical experimental results. Conclusion and future works are described in section 7.

## 2 Previous works

In this section, we will describe some previous works closely related to our works.

Our proposed algorithm requires surface parameterization of the high-genus surface onto its universal covering space. Surface parameterization has been extensively studied, for which different representations of bijective surface maps have been proposed. Conformal registration, which minimizes angular distortion, has been widely used to obtain a smooth 1-1 correspondence between surfaces [4, 5, 7, 1, 13, 6]. For example, Hurdal et al. [13] proposed to compute the conformal parameterizations using circle packing and applied it to registration of human brains. Gu et al. [5, 7, 6] proposed to compute the conformal parameterizations of human brain surfaces for registration using harmonic energy minimization and holomorphic 1-forms. Conformal registration is advantageous for it preserves the local geometry well.

Surface registration, which aims to find an optimal one-to-one correspondence between surfaces, has also been extensive studied. Various algorithms have been proposed by different research groups. Landmark-free registration has been proposed to obtain 1-1 correspondences between shapes without feature landmarks. Different algorithms have been proposed to obtain registrations based on the shape information (such as curvatures) defined on the surfaces. Lyttelton et al. [2] computed surface parameterizations with surface curvature matching. Fischl et al. [4] proposed an algorithm for brain registration that better aligns cortical folding patterns, by minimizing the mean squared difference between the convexity of the surface and the average convexity across a set of subjects. Lord et al. [14] proposed to match surfaces by minimizing the deviation of the registration from isometry. Yeo et al. [15] proposed the spherical demons method, which adopted the diffeomorphic demons algorithm [17], to drive surfaces into correspondence based on the mean curvature and average convexity. Quasi-conformal mappings have been applied to obtain surface registration with bounded conformality distortion [24, 26, 25, 27]. For example, Lui et al. [26] proposed to compute quasi-conformal registration between hippocampal surfaces based on the holomorphic Beltrami flow method, which matches geometric quantities (such as curvatures) and minimizes the conformality distortion [24]. Most of the above registration algorithms cannot match feature landmarks, such as sulcal landmarks on the human brains, consistently. To alleviate this issue, landmark-matching registration algorithms are proposed by various research groups. Bookstein et al. [16] proposed to obtain a registration that matches landmarks as much as possible using a thin-plate spline regularization (or biharmonic regularization). Tosun et al. [30] proposed to combine iterative closest point registration, parametric relaxation and inverse stereographic projection to align cortical sulci across brain surfaces. These diffeomorphisms obtained can better match landmark features, although not perfectly. Wang et al. [18, 21]

proposed to compute the optimized conformal parameterizations of brain surfaces by minimizing a compounded energy. All of the above algorithms represent surface maps with their 3D coordinate functions. Special attention is required to ensure the bijectivity of the resulting registration. Besides, smooth vector field has also been proposed to represent surface maps. Lui et al.

[20] proposed the use of vector fields to represent surface maps and reconstruct them through integral flow equations. They obtained shape-based landmark matching harmonic maps by looking for the best vector fields minimizing a shape energy. The use of vector fields to represent surface maps makes optimization easier, but they cannot describe all surface maps. Time dependent vector fields can be used to represent the set of all surface maps. For example, Joshi et al. [8] proposed the generation of large deformation diffeomorphisms for landmark point matching, where the registrations are generated as solutions to the transport equation of time dependent vector fields. The time dependent vector fields facilitate the optimization procedure, although it may not be a good representation of surface maps since it requires more memory. Later, Lin et al. [34] propose a unified variational approach for registration of gene expression data to neuroanatomical mouse atlas in two dimensions that matches feature landmarks. Again, landmarks cannot be exactly matched. Note that inexact landmark-matching registrations are sometimes beneficial. In the case when landmark points/curves are not entirely accurate, this method is more tolerant of errors in labeling landmarks and gives better parameterization. Most of the above algorithms deal with the registration problem between simply-connected open or closed surfaces.

## 3 Mathematical background

In this section, we describe some basic mathematical concepts related to our algorithms. For details, we refer the readers to [3][28][29].

A surface with a Riemannian metric is called a Riemann surface. Given two Riemann surfaces and , a map is conformal if it preserves the surface metric up to a multiplicative factor called the conformal factor. An immediate consequence is that every conformal map preserves angles. With the angle-preserving property, a conformal map effectively preserves the local geometry of the surface structure.

According to the Riemann uniformization theorem, every Riemann surface admits a conformal Riemannian metric of constant Gaussian curvature. Such metric is called the uniformization metric. The uniformization metric for a genus surface induces 0 Gaussian curvature, whereas a genus surface induces Gaussian curvature, which is called the hyperbolic metric of the surface.

Given a high-genus surface (with genus ), is associated with a universal covering space . A universal covering space is a simply-connected space with a continuous surjective conformal map satisfying the following: for any , there exists an open neighborhood of such that is a disjoint union of open sets in . When , is equal to the whole plane . When , is the unit disk equipped with the hyperbolic metric, which is called the Poincarè disk . The Poincarè disk is a unit disk with metric defined as follows:

 ds2=4dzd¯z(1−z¯z)2 (1)

The distance between two points and on Poincare disk is given by:

 d(z,z0)=tanh−1|z−z01−z¯z0| (2)

All rigid motions on Poincarè disk are Mobiüs transformations:

 z→eiθz−z01−z¯z0,z0∈D,θ∈[0,2π] (3)

A generalization of conformal maps is the quasi-conformal maps, which are orientation preserving homeomorphisms between Riemann surfaces with bounded conformality distortion, in the sense that their first order approximations takes small circles to small ellipses of bounded eccentricity [3]. Thus, a conformal homeomorphism that maps a small circle to a small circle can also be regarded as quasi-conformal. Surface registrations and parameterizations can be considered as quasi-conformal maps. Mathematically, is quasi-conformal provided that it satisfies the Beltrami equation:

 ∂f∂¯¯¯z=μ(z)∂f∂z. (4)

for some complex valued Lebesgue measurable satisfying . is called the Beltrami coefficient, which is a measure of non-conformality. In particular, the map is conformal around a small neighborhood of when . From , we can determine the angles of the directions of maximal magnification and shrinking and the amount of them as well. Specifically, the angle of maximal magnification is with magnifying factor ; The angle of maximal shrinking is the orthogonal angle with shrinking factor . The distortion or dilation is given by:

 K=1+|μ(p)|1−|μ(p)|. (5)

Thus, the Beltrami coefficient gives us all the information about the properties of the map.

Given a Beltrami coefficient with . There is always a quasiconformal mapping from onto itself which satisfies the Beltrami equation in the distribution sense [3].

## 4 Algorithms

In this section, we explain our algorithm for registering high-genus surfaces in details. The basic idea is to embed the surfaces into their universal covering spaces and register them on the universal covering space. Our proposed algorithm can be divided into three main stages:

1. Embedding of the high-genus surface into the universal covering space: The high-genus surfaces are first conformally parameterized into its universal covering spaces in , which is the Euclidean plane for and hyperbolic disk for , using the discrete Ricci flow method.

2. Computing the initial registration Harmonic registration between the fundamental domains is computed as an initial registration.

3. Shape matching registration: A surface registration which matches the geometry is obtained by minimizing a shape mismatching energy on the universal covering space.

In the following, we will describe each stages in detail.

### 4.1 Embedding of the high-genus surface into the universal covering space

In this work, the surface registration is computed on the universal covering space in of the high-genus surfaces. This simplifies the calculation, since all computations can be done in the two dimensional space.

In this work, the embedding of into its universal covering space is computed using the Ricci flow method introduced by Gu et al. [32][33]. Ricci flow is the process to conformally deform the surface metric according to its induced Gaussian curvature . The process is similar to heat flow on manifolds:

 dgij(t)dt=−2(K(t)−¯K)gij(t) (6)

where () or () is target curvature. Convergence of this process is guaranteed by Hamilton’s theorem. is the desired uniformization metric.

To obtain the embedding, the surface is firstly sliced along the cut graph to get the fundamental domain of , denoted as . Let be a base point on the surface . Then, there are many closed loops based at . Two loops and are said to be equivalent if one can be deformed into the other without breaking. Mathematically, there exists a homotopy such that and . All equivalent closed loops equivalent to each other form an equivalence class. The set of all equivalence classes form a group, which is called fundamental group, , of . Suppose is a basis of . Slicing along the basis, the high-genus surface will become a simply connected open surface, which is called fundamental domain.

The fundamental group basis is called canonical if any two loops intersect only at the base point . From algebraic topology, the boundary of the fundamental domain with respect to the canonical loops or cuts is given by

 a1b1a−11b−11a2b2a−12b−12⋯agbga−1gb−1g (7)

In this paper, we apply the greedy approach proposed in [31] to compute the homotopic basis. Each canonical cut is the chosen to be the shortest path in its equivalent class.

With the uniformization metric, the fundamental domain can be embedded onto a region in its universal covering space. For genus one closed surface, the universal covering is the Euclidean plane (See Figure 1(A)). For genus greater than one, the universal covering is the Poincarè disk (See Figure 1(B)). Let , where is the index set. belongs to one of the pieces, . Also, and intersect at the boundaries only if . By glueing all ’s together, the universal covering can be obtained.

Note that the canonical cuts are introduced to obtain the universal covering only. During the registration process, the canonical cuts on the source surface are allowed to move freely on the target surface, since the whole process will be done on the universal coverings. In other words, the correspondences between canonical cuts are not required. It avoids the issue of finding a consistent cuts to obtain the accurate registration.

### 4.2 Initial registration between fundamental domains

We first compute an initial surface registration between two high-genus surfaces and of genus . An initial map can be chosen as the harmonic map by fixing the correspondence of the boundary cuts.

Let and be the canonical fundamental domains of and respectively, computed in the first stage as described in the previous subsection. In this subsection, the metric used is always chosen to be the Euclidean metric if and the hyperbolic metric if .

The canonical polygon has vertices and hence edges. Vertices of the canonical polygon corresponding to the single base point on the surface. All edges are geodesics.

Let and be the base point of and respectively, corresponds to . With these base points, and can be conformally mapped to their fundamental domains and in their universal covering spaces (See Figure 1). We denote the conformal parameteriztions by and . To obtain an initial registration between and , we compute a mapping between and . gives us an initial mapping between and .

Here, we assume corresponding edges between two canonical polygons can be matched. In other words, we assume the boundary condition is given, through the arc-length parameterization. Note that the boundary cuts on each surfaces might not exactly correspond to each others. However, since canonical cuts are chosen, the edges corresponds to the shortest loops on the surfaces. As a result, the initial boundary correspondence is a reasonable guess for the initial registration. With the boundary correspondence , a unique harmonic map between the two canonical polygons can be computed.

The harmonic map can be computed by minimizing harmonic energy

 E(g)=∫Ω1|∇g|2Ω2,  given   g|∂Ω1=h. (8)

Minimizing the above energy functional is equivalent to solve following PDE

 Δg=0  subject to   g|∂Ω1=h. (9)

where is Laplace-Beltrami operator under the uniformization metric.

This initial registration provide us with a smooth mapping between and . Note that there are also other choices of initial maps, such as patch by patch registration or landmark matching registration.

Now, to obtain a geometric matching registration, we propose to refine the registration from the initial registration to match a geometric mismatching energy. The detailed numerical implementation of the initial registration will be described in section 4.3.

### 4.3 Shape matching registration

In the previous subsection, we use hyperbolic harmonic map between two canonical polygons as initial registration. We assume that the boundary cuts are properly matched. Edges of canonical polygon are (hyperbolic) shortest loops on surface which pass the base point. However, the shortest loops depend on the uniformaization metric, which do not directly take the geometric information of the surfaces into consideration. Constraining the boundary cuts to be exactly matched often induces error in the final registration. To obtain a better geometric matching registration, we propose a variational approach, which minimizes a geometric mismatching energy, without fixing the correspondences of the boundary cuts.

Surface curvatures are important quantities to describe the surface geometry. We therefore consider an energy functional which measures the curvature mismatching under a registration . More specifically, we propose to find an optimal diffeomorphism which minimizes the following energy functional :

 E(f)=12∫S1|∇g1f|2S2+α22∫S1(H1−H2∘f)2+β22∫S1(K1−K2∘f)2 (10)

where are mean curvatures on and respectively, and are the Gauss curvatures on and respectively.

The first term, which is the harmonic energy, controls the smoothmess of the registration. The last two terms, which measure the mistmatching of surface curvatures, are used to match the surface geometry.

Solving the above variational problem (10) directly on the surfaces is challenging. To simplify the optimization process, we propose to solve the problem on the universal covering spaces of and .

#### 4.3.1 Optimization on universal covering spaces

Let and be the covering maps of and respectively. Suppose , where is the index set and intersects with at their boundaries if . Similarly, we let , where intersects with at their boundaries if . We then proceed to look for a diffeomorphism , which is the lifting of the optimal registration . In other words, we require that

 π−11|˜D1α∘g∗∘π2=f∗, for any α∈U (11)

Equation (11) ensures that satisfies the periodic condition on the covering spaces. In practice, suppose the canonical cuts on is given by , we require that

 φi(g∗(ai))=g∗(a−1i) and ϕi(g∗(bi))=g∗(b−1i) (12)

where and are the deck transformations.

Since is the lifting of , it minimizes the following energy functional:

 EH(g)=12∫˜S1|∇g|2+α22∫˜S1(~H1−~H2∘g)2+β22∫˜S1(~K1−~K2∘g)2 (13)

subject to the constraint that for all .

are called the Fuchsian group generators, which are the generators of the Deck transformation group of . When , and are just translations in . When , and are Mobiüs transformations of the unit disk, which can be computed explicitly. We will describe the computation of . The other Fuchsian group generators can be obtained in the same way. Suppose the starting point and ending points of are and , and the starting point and ending points of are and . We need to look for a Mobiüs transformation such that and . We first compute a Mobiüs transformation to map to the origin, which is given by: . Then, maps to a radial Euclidean line. Let the angle between and the real axis be , and let . Then, maps to the origin and to the real axis. Similarly, we can find Mobiüs transformation and such that maps to the origin and to the real axis. The deck transformation is then given by: .

To solve the optimization problem (13), we use a splitting method to minimize:

 EH(g,h)=12∫˜S1|∇g|2+μ22∫˜S1|g−h|2+α22∫˜S1(~H1−~H2∘h)2+β22∫˜S1(~K1−~K2∘h)2 (14)

Fixing , we first minimize :

 E1(h)=μ22∫˜S1|g−h|2+α22∫˜S1(~H1−~H2∘h)2+β22∫˜S1(~K1−~K2∘h)2 (15)

At each point , we consider the Taylor’s expansion of and about ,

 H2(h)(p)≈H2(g)(p)+∇H2(g)(p)⋅(h−g)(p)K2(h)(p)≈K2(g)(p)+∇K2(g)(p)⋅(h−g)(p) (16)

Plugging equations (16) into equation (15), we look for a small perturbation from to such that is minimized. It can be done by solving the following PDE:

 μ2(g−h)+α2(H1−H2(h))∇H2(h)+β2(K1−K2(h))∇K2(h)=0 (17)

In the discrete case, the above problem can be solved by the Guass-Newton method, which will be described in the next section.

Next, fixing , we minimize

 E2(g)=12∫˜S1|∇g|2+μ22∫˜S1|g−h|2 (18)

can be minimized by solving the elliptic PDE:

 Δg−μ2(g−h)=0 (19)

Recall that the registration computed should satisfy the constraint (11). Hence, we enforce this constraint when solving equation (19). In the discrete case, the above problem becomes a nonlinear system, which can be solved effectively using Newton’s method.

In this way, we can minimize alternatively over and . More specifically, suppose is obtained at the n-th iteration, we fix to obtain by solving equation (17). We then fix to obtain by solving equation(19).

#### 4.3.2 Preservation of bijectivity

One crucial issue in computing the surface registration is to preserve its bijectivity. In this work, we propose to enforce the bijectivity using the Beltrami coefficient of the surface map.

Let be the mapping between the universal coverings of and . We need to ensure that is bijective. Every mapping is associated with a Beltrami coefficient, , which is a complex-valued function defined on . is bijective if and only if its Jacobian everywhere. Simple checking gives

 Jg=|∂g∂z|2(1−|μ(g)|2) (20)

Hence, is bijective if and only if everywhere.

Motivated by the above observation, we propose to enforce in each iterations during the optimization process described in the last subsection. This can be done as follows. Suppose is obtained at the n-th iteration. Let be a small parameter. We first compute:

 νn={max{|μn|,1−ϵ}μn|μn|,if |νn|≠00,if |νn|=0 (21)

We then smooth by minimizing the following energy functional:

 ∫H2|∇ν|2+λ2∫H2|ν−νn|2 (22)

The above minimization problem is equivalent to solving the following PDEs:

 Δν+λ(ν−νn)2=0 (23)

subject to the constraint that for every for all and for all .

Once a smooth Beltrami coefficient is obtained, we need to find a quasi-conformal map whose Beltrami coefficient closely resemble to . Suppose with Beltrami coefficient . We can write and as linear combinations of and ,

 −vy=α1ux+α2uy;vx=α2ux+α3uy. (24)

where ; ; .

Similarly,

 −uy=α1vx+α2vy;ux=α2vx+α3vy. (25)

Since , we obtain

 ∇⋅(D(uxuy))=0  and  ∇⋅(D(vxvy))=0 (26)

where .

Therefore, to construct , we let and solve equation (26) subject to the constraint that for all . The details of the numerical implementation will be explained in the next section.

We summarize our proposed high-genus surface registration algorithm as follows.

(High-genus surface registration)
High-genus surface and
Geometric matching surface registration

1. Compute the conformal parameterizations and of and respectively;

2. Compute the initial mapping ; Let ;

3. Given at n-th iteration, obtain by fixing and solving equation 17; Fixing , obtain by solving equation 19;

4. Compute the Beltrami coefficient of ; obtain a smooth Beltrami coefficient by solving equations 21 and 23;

5. Obtain a quasi-conformal map from by solving equation 26;

6. If , continue. Otherwise, stop the iteration.

## 5 Numerical implementation

In this section, we describe in details the numerical implementation of our proposed algorithm. All our computations are done on the universal covering space, which is for genus one surfaces and for high genus surfaces. The universal covering space consists of infinite copies of fundamental domains, which are unique up to deck transformations. Many important operators are identical on each fundamental domains. For example, the Laplace-Beltrami operator, which is crucial in our model, are identical on each fundamental domain, since it is invariant to rigid motions. Based on this observation, the numerical implementation can be done on one piece of the fundamental domain, while allowing its boundary to be mapped freely onto the universal covering space of the target surface. In other words, boundary correspondences between the canonical cuts of the two surfaces are not enforced.

### 5.1 Poisson’s equation on the universal covering

Laplace-Beltrami operator plays a crucial role in our proposed algorithm. Most key steps involve solving a Poisson equation on the universal covering space. In this subsection, we describe how to discretize Laplace-Beltrami operator on the fundamental domain, which can be lifted to the universal covering space. Poisson’s equation can then be solved on the universal covering space.

On triangle mesh, Laplace-Beltrami operator can be discretized by the cotangent formula:

 ΔMf(zi)=∑j∈N1(i)wij(f(zj)−f(zi))

where is the set of vertex indices of one-ring neighbors of vertex ; where and are the two angles facing the edge . We use denote both the vertex and it’s complex coordinate. Then we obtain the Poisson equation in matrix form:

 Af=b;

where is a square matrix, .

The computation on each vertex only uses its one-ring neighbor vertices. So we can discretize the Laplace-Beltrami operator on each vertex of the fundamental domain. But on boundary of fundamental domain, the discretization will use vertices outside the fundamental domain (see Figure 2).

Every vertex outside the fundamental domain has a unique copy inside the fundamental domain up to a rigid motion. Whenever computation involves vertices outside the fundamental domain, we will refer to its inside copy. Note that this is also valid for base points, so we have a valid discretization on base points. We always fix base points by letting if , for a base point . For a vertex on boundary except base points, the discretization will be

 ΔMf(zi)=∑j∈N1(i)wij(f(zj)−f(zi))+∑j∈˜N1(i)wij(f(zj)−f(zi)) (27)

where is the set of vertex indices of one-ring neighbors of vertex on the fundamental domain, is the set of vertex indices of one-ring neighbors of vertex outside the fundamental domain, while denotes the set of vertex indices of one-ring neighbors of vertex inside the fundamental domain. For example, in figure LABEL:fig:Poincare_boundary_illustration, , , .

For simplicity, we let . The Laplace-Beltrami operator becomes

 ΔM~zi=∑j∈N1(i)wij(~zj−~zi)+∑j∈˜N1(i)wij(~zj−~zi) (28)

Suppose is outside the fundamental domain. We denote the the inside copy of vertex by . Let be the deck transformation that moves to , that is, , we have

 ΔM~zi =∑j∈N1(i)wij(~zj−~zi)+∑j∈˜N1(i)wij(~zj−~zi) =∑j∈N1(i)wij(~zj−~zi)+∑j′∈∘N1(i′)wi′j′(~zj′−~zi′) =∑j∈N1(i)wij(~zj−~zi)+∑j∈˜N1(i)wij(φ(~zj)−φ(~zi))

The second equality uses the fact that Laplace-Beltrami operator is rigid-motion invariant.

With this discretization, the Poison’s equation can be rewriten in a matrix form:

 A~z+Q(~z)=b (29)

where is the matrix representation of the Laplace-Beltrami operator, and is a deck transformation that transforms outside neighbor of vertex to its inside copy and is zero elsewhere.

For genus one surfaces, deck transformations are linear translation and so is a linear operator. Combining into , equation 29 becomes a linear system and can be solved efficiently.

However, for higher genus surfaces, deck transformations are Mobiüs transformation, which is nonlinear. Equation 29 becomes a nonlinear system. It can be solved by Newton’s method efficiently. Let . is the gradient of , where is computed element-wisely. The problem can then be solved using standard Newton’s method:

1. initialize by , which is current position;

2. compute , if , stop the process;

3. compute , solve from equation ; if , stop the process; Otherwise, let and go to step 2.

The linear equation in step 3 can be solved by LU factorization, which turns out to be quite efficient. In our numerical computation, we observe that the Newton’s method converges very quickly: usually two or three iterations will achieve accuracy.

### 5.2 Solving energy minimizing problem

We use an alternating approach to minimize the proposed energy function. In each iteration, we first minimize to get , then minimize to get .

We first discuss the minimization of . With the linear approximation (16), we have

 E1(h) =μ22∫˜S1|g−h|2 +α22∫˜S1(~H1−~H2(g)−∇~H2(g)⋅(h−g))2 +β22∫˜S1(~K1−~K2(g)−∇~K2(g)⋅(h−g))2 =12∫˜S1∥∥ ∥ ∥∥⎛⎜ ⎜⎝cα(~H1−~H2(g))β(~K1−~K2(g))0⎞⎟ ⎟⎠−⎛⎜ ⎜⎝α∇~H2(g)β∇~K2(g)μI2×2⎞⎟ ⎟⎠(h−g)∥∥ ∥ ∥∥2

Then the minimization problem can be solved individually for each vertex in least square sense:

 ⎛⎜ ⎜⎝α∇~H2(g)β∇~K2(g)μI2×2⎞⎟ ⎟⎠(h−g)(p)=⎛⎜ ⎜⎝α(~H1−~H2(g))(p)β(~K1−~K2(g))(p)0⎞⎟ ⎟⎠

Let

 S=⎛⎜ ⎜⎝α∇~H2(g)β∇~K2(g)μI2×2⎞⎟ ⎟⎠,d=⎛⎜ ⎜⎝α(~H1−~H2(g))β(~K1−~K2(g))0⎞⎟ ⎟⎠

we have

 h(p)=g(p)+(STS)−1⋅(STd)

In computation, the inversion can be obtained by Sherman-Morrison formula. Let , , then , . Apply Sherman-Morrison formula twice, we have

 (STS)−1=1μ2(I−uuT+vvT+(uT⋅v⊥)2I1+uTu+vTv+(uT⋅v⊥)2)

Hence we have a simple solution for . If we consider either mean curvature or Gaussian curvature , i.e., or , the expression of can be further simplified. For example, if , we have

 (STS)−1=1μ2(I−uuT1+uTu)

hence,

 h(p)=g(p)+(~H1−~H2(g))∇~H2(g)μ2α2+∇~H2(g)T∇~H2(g)(p)

The minimization of is obtained by solving equation (19):

 Δg−μ2(g−h)=0 (30)

where Laplace operator is discretized by cotangent formula.

Since the operation discussed in section 5.1 will not affect identity matrix, it can be applied to this equation. So we have

 (A−μ2I)g+Q(g)=−μ2h (31)

This nonlinear equation is then solved by Newton’s method. For genus one surfaces, it is still linear, we can solve it directly.

### 5.3 Solving Beltrami equation

To ensure bijectivity, a smoothing operation on Beltrami coefficient is applied. Then we reconstruct the mapping from smoothed Beltrami coefficient by solving the Beltrami equation.
The Beltrami equation is in fact a Poisson equation with a generalized Laplace-Beltrami operator (see equation 26). We can solve the equation as described in section 5.1. Since we have a generalized Laplace-Beltrami operator, cotangent formula can’t be used. We use discretization scheme proposed in [27], which also uses one-ring neighborhood to discretize the generalized Laplace-Beltrami operator. Hence, the method described in section 5.1 can still be applied.

More specifically, the gradient operator can be discretized by linear approximation. For a triangle , the coordinates of three vertices, let , we have

 ∇Tfi=14aT(fiti+fjtj+fktk)

where is the area of the triangle, .
Then we obtain the discrete gradient operator at vertex :

 ∇fi=∑T∈Ni14aT(fiti+fjtj+fktk)

where be the collection of neighborhood faces attached to vertex . Note that in the summation we omit the superscripts on and to avoid confusion.
Similarly, the discretization of divergence operator for a vector :

 ∇⋅Fi=∑T∈Ni14aT(Fi⋅ti+Fj⋅tj+Fk⋅tk)

The discretization of equation (26) can be obtained by applying above two formulas.
Following the discussion in section 5.1, the Beltrami equation can be formulated as

 A~z+Q(~z)=b (32)

The above equation is linear in the case of genus one surfaces and is nonlinear in the case of higher genus surfaces. By solving the equation, we will get reconstructed quasi-conformal map associated to the smoothed Beltrami coefficient.

## 6 Experimental results

To test the efficacy of the proposed algorithm, experiments have been carried out on synthetic high-genus surface data together with real medical data (vertebrae bone and vestibular system).

### 6.1 Synthetic surface data

We first test our algorithm on synthetic surface data.

##### Example 1

In our first examples, we test the proposed method on a standard torus of genus one. Figure 3(A) and (B) show two genus-1 torus, denoted by and respectively, with different intensity functions defined on each of them. The two surfaces are parameterized onto their universal covering spaces, and registration between the two surfaces is computed on the 2D parameter domains. The intensity functions on each surfaces are plotted on their universal covering spaces, which are shown in (C) and (D). Figure 4(A) shows the registration result that matches the intensity functions. The intensity function defined on is mapped to using the obtained registration. (B) shows the registration result on the universal covering spaces. The intensity functions are perfectly matched under the obtained registration (compared with Figure 3(D)). Note that the boundary cuts are not fixed. They move freely on the universal covering space, which satisfy the periodic conditions. Figure 5 shows the curvature mismatching energy, harmonic energy and total energy versus iterations. All of them decrease monotonically as iteration increases. It demonstrates that our algorithm computes the optimized harmonic map between the genus-1 surfaces that matches the intensity functions as much as possible.