Computing Harmonic Maps and Conformal Maps on Point Clouds

Comparing to methods for triangle meshes, there are much fewer works on meshless methods of computing conformal and harmonic maps, especially for higher genus surfaces. Guo et al. [17] computed global conformal parameterizations of surfaces by computing holomorphic 1-forms on point clouds. Li et al. [18] computed harmonic volumetric maps by grids discretizations. Meng-Lui [19] developed the theory of computational quasiconformal geometry on point clouds. Using approximations of the differential operators on point clouds, Liang et al. [20] and Choi et al. [21] computed the spherical conformal parameterizations of genus-0 closed surfaces, and Meng et al. [22] computed quasiconformal maps on topological disks. Li-Shi-Sun [23] computed quasiconformal maps from surfaces to planar domain, using the so-called point integral method for discretizing integral equations for point clouds.

There is an extensive literature on computing conformal maps for triangle meshes. Gu-Yau [24][25] developed the method of computing conformal structures of surfaces by computing the discrete holomorphic one-forms. Pinkall-Polthier [13] proposed a method of conformal parameterization by computing a pair of conjugate harmonic functions. Lévy et al. [26] and Lipman [27] and Lui et al. [28] computed conformal or quasiconformal maps by minimizing or controlling the conformal distortion. There is also a big family of methods based on various notions of discrete conformality for triangle meshes, such as circle patterns [29][30][31][32], and inversive distances [33][34], and vertex scalings [35][36][32], and modified vertex scalings [37] [38] [39][40] allowing diagonal switches. Some related convergence results for discrete conformality can be found in [29][41][42][43][44], and other mathematical analysis can be found in [45][46][47][48][49]. Other works on computing conformal maps on triangle meshes include [50][51][52][53][26][54][55][56][57][58].

For computing harmonic maps on triangle meshes, Gaster-Loustau-Monsaingeon [59][60] give detailed discussions on the mathematical analysis and algorithms, and produced a computer software. Notions of discrete Dirichlet energy have been discussed and used to compute discrete harmonic maps, not only in Gaster-Loustau-Monsaingeon’s work [59][60], but also extensively in other works such as [61][13][62][63][64][65][66][67][68].

To the best of the authors’ knowledge, our method is the first meshless method in computing harmonic and conformal maps to surfaces of genus greater than 1.

Existing meshless methods always use approximating differential operators on point clouds, to compute harmonic or conformal maps. Our method provides a novel type of approach, by jumping out of the point clouds to cubic lattices. It also has the following nice properties.

• •

  The idea of lattice approximation is conceptually simple, and could be implemented for any topological types of surfaces.

• •

  One can simply tune the density of the lattices to get satisfactory accuracy, within the ability of the computing powers.

• •

  Numerical experiments indicate that sparse lattices already work well.

• •

  The lattice structure should be suitable for parallel computation, and multi-grid methods.

We proposed algorithms of computing harmonic maps and conformal maps for any closed surfaces and topological disks, embedded in $\mathbb{R}^{3}$, to constant curvature surfaces.

Given a point, or equivalently a vector $x$ in $\mathbb{R}^{3}$, denote

as the $l^{2}$-norm of $x$.

Given two points $x,y\in\mathbb{R}^{3}$, denote

as the distance between $x$ and $y$.

Given a subset $A\subset\mathbb{R}^{3}$ and a point $x\in\mathbb{R}^{3}$, denote

as the distance from $x$ to $A$.

Given a subset $A\subset\mathbb{R}^{3}$ and $r>0$, denote

as the $r$-neighborhood of the subset $A$.

Given a closed surface $M$ embedded in $\mathbb{R}^{3}$, we say a finite subset $P\subset\mathbb{R}^{3}$ is a $\delta$-point cloud of $M$ if

If $G=(V,E)$ is an undirected connected simple graph, and $w\in\mathbb{R}^{E}_{>0}$ is an edge weight, then $L=L_{G,w}$ denotes the discrete Laplacian such that for any $f:V\rightarrow\mathbb{R}^{k}$,

The remaining of the paper is organized as follows. We review the basic mathematical theory in Section 2, and then introduce the lattice approximation of a surface in Section 3. The algorithms, as well as numerical examples, are given in Section 4, 5, 6, 7, for the cases of spheres, rectangles, flat tori, and hyperbolic surfaces respectively.

This work is partially supported by Center of Mathematical Sciences and Applications at Harvard University, and Yau Mathematical Sciences Center at Tsinghua University, and NSF 1760471.

We give a formal definition of a conformal map as follows.

In this paper we always assume that a conformal map is a diffeomorphism between two surfaces. By the celebrated uniformization theorem, it is well known that any closed orientable surface is conformally equivalent to a surface of constant Gaussian curvature.

The Riemannian surface $(N,h)$ above is called a uniformization of $(M,g)$ and could be viewed as a canonical representative in the conformal equivalence class. A closed surface $(N,h)$ of constant curvature $-1$ or $0$ can be naturally represented as a planar domain, after a few cuts.

1. (a)

   If $(N,h)$ has constant curvature 0, then it is isometric to the flat complex plane $\mathbb{C}$ modulo a lattice $\mathbb{Z}\oplus\mathbb{Z}\tau$ where $\tau\in\mathbb{C}$ and $Im(\tau)>0$. If we properly cut two loops on $N$, then $(N,h)$ can be displayed isometrically as a planar domain in $\mathbb{C}$.

2. (b)

   if $(N,h)$ has constant curvature -1, then it is isometric to the hyperbolic plane $\mathbb{H}^{2}$ modulo a discrete subgroup $\Gamma$ of the orientation-preserving isometry group $Isom^{+}(\mathbb{H}^{2})$. In the Poincaré disk model, the hyperbolic plane $\mathbb{H}^{2}$ is identified as the unit disk $\{z\in\mathbb{C}:|z|<1\}$ with the complete Riemannian metric

   $$\frac{4|dz|^{2}}{(1-|z|^{2})^{2}}=\frac{4(dx^{2}+dy^{2})}{(1-x^{2}-y^{2})^{2}}.$$

   Similar to the flat case, if we properly cut a few loops, $(N,h)$ can be displayed, isometrically in the hyperbolic sense, as a domain in the unit disk $\{|z|<1\}$.

So for genus 0 surfaces, there is essentially a unique conformal equivalence class. For genus 1 surfaces, the conformal equivalence classes can be parameterized by a complex number $\tau$ with $Im(\tau)>0$. For a surface $M$ with genus $>1$, the conformal equivalence classes can be parameterized by the generators of the discrete subgroup $\Gamma$ of $Isom^{+}(\mathbb{H}^{2})$. As shown above, a surface of constant curvature $\pm 1$ or $0$ can always been identified as a planar domain, by a stereographic projection or cutting loops. So computing diffeomorphisms to such surfaces gives rise to global parameterizations and surface flattening, which have fundamental applications in computational graphics.

For nonclosed surfaces, the most important case is the topological disk. Assume $(M,g)$ is a Riemannian surface homeomorphic to a closed disk, then again by the uniformization theorem it is conformally equivalent to a closed unit disk. However, for our convenience in utilizing the method of harmonic maps, we consider conformal maps and harmonic maps to rectangles instead of disks. The following uniformization-type theorem for rectangles is well-known.

Let $(M,g)$ and $(N,h)$ be two smooth Riemannian surfaces. The Dirichlet energy, i.e., the stretching energy of a smooth map $f:M\rightarrow N$ is formally defined as

where $dv_{g}$ denotes the volume element on $(M,g)$, and $\|df\|$ is the norm of the differential of $f$, with respect to the induced metric on $T^{*}M\otimes f^{*}(TN)$.

Smooth harmonic maps have been extensively studied. See [69][70] for examples. Here we focus on harmonic diffeomorphisms to surfaces of constant curvature $\pm 1$ or $0$, and harmonic diffeomorphisms to rectangles with given boundary correspondences. Some relevant well-known facts are summarized as the following 2 theorems.

By part (b) and (e) of Theorem 2.5, for closed surfaces it makes good sense to compute the harmonic map within a fixed homotopy class. Such a topological constrain also often naturally arises in applications. Part (b) of Theorem 2.5 and part (b) of Theorem 2.6 somewhat justify the physical intuition that harmonic maps minimize the stretching energy. Minimizing the Dirichlet energy would also be an important idea for computing harmonic maps. Part (c) of Theorem 2.5 and part (c) of Theorem 2.6 connect the notions of conformal maps and harmonic maps. So conformal maps are really special harmonic maps, and the converse is true if the target surface is a unit sphere. Part (a) of Theorem 2.5 and Part (a) of Theorem 2.6 point out that conformal maps minimize the normalized Dirichlet energy $\mathcal{E}(f)/Area(N)$ among the harmonic maps. This optimality provides a method to compute the uniformization using harmonic maps.

1. (a)

   If $M$ is a topological sphere, a harmonic diffeomorphism from $(M,g)$ to a unit sphere is already a conformal map.

2. (b)

   If $M$ is a genus 1 closed orientable surface, one can first compute the harmonic map to a flat torus $(N,h)=\mathbb{C}/\mathbb{Z}\oplus\mathbb{Z}\tau$ for some parameter $\tau$, and then minimize the normalized Dirichlet energy, by tuning the parameter $\tau$ in the upper half plane.

3. (c)

   If $M$ is a closed orientable surface with genus $>1$, one can first compute the harmonic map to a hyperbolic surface $(N,h)=\mathbb{H}^{2}/\Gamma$ for some discrete subgroup $\Gamma$ of $Isom^{+}(\mathbb{H}^{2})$, and then minimize the normalized Dirichlet energy, by continuously deforming the subgroup $\Gamma$.

4. (d)

   If $M$ is homeomorphic to a closed disk, then one can compute the harmonic map to a rectangle $[0,a^{-1}]\times[0,a]$, and then minimize the Dirichlet energy, by tuning the parameter $a\in\mathbb{R}_{>0}$.

Since we will be computing discrete harmonic maps on lattice graphs, here we briefly introduce the theory of discrete harmonic maps, which has been well-studied. See [68][71][61] for references of the related definitions and theorems. In this subsection, we always assume that $G=(V,E)$ is an undirected connected simple graph, and $w\in\mathbb{R}^{E}_{>0}$ denotes the edge weight, and $(N,h)$ is a flat rectangle or a closed orientable surface of constant curvature $\pm 1$ or 0.

A discrete map $f$ from $G$ to $(N,h)$ maps each vertex $i\in V$ to a point in $N$, and each edge $ij\in E$ to a geodesic segment in $N$ with endpoints $f(i),f(j)$. Given a discrete map $f$, we denote $l_{ij}(f)$ as the length of $f(ij)$. Assuming the stretching energy on a single edge $ij$ is $\frac{1}{2}w_{ij}l_{ij}^{2}$, we have the following definition of the discrete Dirichlet energy.

The notion of discrete harmonic maps can be defined by the balanced conditions on vertices.

This definition consists with the continuous one in the following sense.

Here the assumption on $N$ and $l_{ij}(f)$ is to guarantee that locally $f$ can be parameterized by $[f(i)]_{i\in V}\subset N^{V}$. For the cases of non-positive curvatures, we have the following two well-known theorems partially analogue to Theorem 2.5 and 2.6.

Assume that $M$ is a closed smooth surface embedded in $\mathbb{R}^{3}$, and $P$ is a $\delta$-point cloud of $M$. If $\delta$ is sufficiently small, comparing to some constant $\epsilon>0$, then the $\epsilon$-neighborhood $B(P,\epsilon)$ of $P$ would be a good approximation of the $\epsilon$-neighborhood $B(M,\epsilon)$ of $M$. If $n$ is sufficiently large, comparing to $1/\epsilon$, then

would be a good cubic lattice approximation of $B(P,\epsilon)\approx B(M,\epsilon)$ (see Figure 1). Connect two vertices $x,y$ in $V$ by an edge if $d(x,y)=1/n$, and then we obtain a lattice graph $G=(V,E)$. Since we are using standard cubic lattice, the weight function is set to be constant $w_{ij}=1$. Such weighted graph $(G,w)$ would be our lattice approximation of $M$.

In practice, if the pointcloud data $P$ is given, we need to choose a proper $\epsilon>0$ such that

1. (a)

   $\epsilon$ is sufficiently large such that $B(P,\epsilon)$ forms a connected domain in $\mathbb{R}^{3}$ and has uniform thickness, and

2. (b)

   $\epsilon$ is sufficiently small such that $B(P,\epsilon)\approx B(M,\epsilon)$ is a good approximation of $M$.

For our second parameter $n$, if we choose a larger integer $n$, the lattice graph $G$ should approximates the domain $B(P,\epsilon)\approx B(M,\epsilon)$ better. However, in practice we find that sparse lattices already work well.

Suppose $G=(V,E)$ is a lattice approximation of $P$ constructed as before, then given a function $f:V\rightarrow\mathbb{R}^{k}$, we can use the standard trilinear interpolation to extend $f$ to the point cloud $P$. Assume $p=(p_{1},p_{2},p_{3})\in P$ lies inside a cubic cell $[x_{0},x_{1}]\times[y_{0},y_{1}]\times[z_{0},z_{1}]$ of the lattice $G$. Denote $a_{ijk}=f(x_{i},y_{j},z_{k})$ where $i,j,k=\in\{0,1\}$, and

then

is the trilinear interpolation of $f$ on $P$.

Our method is smoothly adapted from Gu-Yau’s method [24] for computing harmonic maps to spheres. First we pick a lattice approximation $G=(V,E)$ as in Section 3, and then

1. (1)

   compute the (approximated) discrete harmonic map $f:G\rightarrow\mathbb{S}^{2}$, and

2. (2)

   extend $f$ to $P$ by trilinear interpolation, and then do a normalization $x\mapsto x/|x|_{2}$ and get an approximation of a harmonic map from $M$ to the unit sphere.

Step (2) is pretty clear and simple, so let us discuss the details of step (1). Recall that a discrete map $f$ from $G$ to $\mathbb{S}^{2}\subset\mathbb{R}^{3}$ is harmonic if it is a critical point of the energy

If the edge length $l_{ij}(f)$ is small, $l_{ij}(f)\approx|f(j)-f(i)|_{2}$, and the discrete energy is approximated by

So for simplicity we will minimize $\mathcal{E}_{0}(f)$ instead of $\mathcal{E}(f)$ to compute the approximation of the discrete harmonic map. Since

a critical point $f$ of $\mathcal{E}_{0}$ satisfies that for any $i\in V$

i.e.,

where $L^{\|}f(i)$ denotes the tangential component of $Lf(i)$, i.e., the orthogonal projection of $Lf(i)$ to the plane $f(i)^{\perp}$. We view a map $f$ satisfying equation (2) as an approximation of a discrete harmonic map, and wish to compute it by the following discrete heat flow on $\mathbb{S}^{2}$

We can simply solve the above ODE by the explicit Euler’s method, with a normalization after each iteration to keep that all the points are on the unit sphere. More specifically, if the step size is $\delta t$ we have that

where $\pi(x)=x/|x|_{2}$.

Two harmonic maps to a sphere can differ by a Möbius transformation. So a smooth harmonic map to a sphere is not unique, and would have large distortion if most of the mass is concentrated near a point of the sphere. To avoid such big distortion, we add a correction mapping $m$ in each iteration. Here the correction mapping $m$ first translate all the points $f(i)$’s in $\mathbb{R}^{3}$ to make the center of mass be at the origin, and then do a normalization $x\mapsto x/|x|_{2}$. Now the modified iteration is

and we have the following algorithm.

See the two numerical examples in Figure 2 and Figure 3.

First we construct a lattice approximation $G=(V,E)$ of $P$ as in Section 3, and let $V_{i}=\cup_{x\in P_{i}}V_{x}$ be the lattice approximation of $\gamma_{i}$, where $V_{x}$ consists of the 8 vertices of the cubic cell of $G$ containing $x\in P$. Then we would like to

1. (1)

   compute the (constrained) discrete harmonic map $f:G\rightarrow\mathbb{[}0,a^{-1}]\times[0,a]$, such that

   	$\displaystyle f(V_{1})\subset$	$\displaystyle[0,a^{-1}]\times\{0\},$
   	$\displaystyle f(V_{2})\subset$	$\displaystyle\{a^{-1}\}\times[0,a],$
   	$\displaystyle f(V_{3})\subset$	$\displaystyle[0,a^{-1}]\times\{a\},$
   	$\displaystyle f(V_{4})\subset$	$\displaystyle\{0\}\times[0,a],$

   and

2. (2)

   extend $f$ to $P$ by trilinear interpolation, and then $f|_{P}$ is our approximation of the desired harmonic map.

Step (2) is clear and simple. Computing $f=(f^{1},f^{2})$ in step (1) amounts to solve the following 2 systems of linear equations

These two equations can be solved efficiently by the preconditioned conjugate gradient method.

In a summary we have the following Algorithm 2.

Assume $f_{a}$ is the discrete harmonic map from the lattice approximation $G=(V,E)$ to the rectangle $[0,a^{-1}]\times[0,a]$. To compute conformal maps, we only need to find a proper edge length $a$ such that the computed discrete harmonic map $f_{a}$ approximates a conformal mapping. Inspired by part (a) of Theorem 2.6, we compute the conformal map by minimizing $\mathcal{E}(f_{a})$ over $a\in\mathbb{R}_{>0}$. If $\bar{a}$ is the minimizer, then $f_{\bar{a}}:V\rightarrow[0,\bar{a}^{-1}]\times[0,\bar{a}]$ would approximate a conformal map.

Our method is summarized as Algorithm 3.

Here are two numerical examples for harmonic maps and conformal maps respectively.

First we construct a lattice approximation $G=(V,E)$ as in Section 3, and then need to

(1) compute the discrete harmonic map $f$ from $G$ to $\mathbb{C}/\mathbb{Z}\oplus\mathbb{Z}\tau$, or equivalently, to a fundamental domain in $\mathbb{C}$, and

(2) extend $f$ to $P$ by trilinear interpolation, and then $f|_{P}$ represents an approximation of the desired harmonic map.

Step (2) is kind of clear and simple, so let us discuss the details of step (1). First we ”cut” the lattice $G$ along two ”loops” $\gamma_{1},\gamma_{2}$, so that we can represent a discrete map from $G$ to a torus as a map from $V$ to a planar domain in $\mathbb{C}$. Secondly we fix the homotopy class of discrete maps $f:G\rightarrow\mathbb{C}/\mathbb{Z}\oplus\mathbb{Z}\tau$, by requiring that $\gamma_{1}$ corresponds to the translation $z\mapsto z+1$, and $\gamma_{2}$ corresponds to the translation $z\mapsto z+\tau$. Now computing the discrete harmonic map in the fixed homotopy class amounts to solve a system of complex linear equations. For a vertex $i$ that is away from the loops $\gamma_{1}$ and $\gamma_{2}$, the balanced condition (1) gives us that

Things become a bit subtle when the vertex $i$ is near a loop. Assume the edge $ij$ pass through the loop $\gamma_{1}$, and $l_{ij}(f)$ is small, and all the other edges adjacent to vertex $i$ do not pass through $\gamma_{1}$ or $\gamma_{2}$. Then

where the sign depends on the orientation of $\gamma_{1}$ and the relative position between $\gamma_{1}$ and $i$. For such a vertex $i$, the corrected balanced condition should be

Similar corrections should be made for all the edges passing through $\gamma_{1}$ or $\gamma_{2}$. After making all the corrections, we arrived at a system of complex linear equations of the form

where $b(i)$ is computed by adding up all the correction terms. The equation (9) can be solved efficiently by the preconditioned conjugate gradient method.

In a summary we have the following Algorithm 4.

Assume that $f_{\tau}$ is the discrete harmonic map from the lattice approximation $G=(V,E)$ to the flat torus $\mathbb{C}/\mathbb{Z}\oplus\mathbb{Z}\tau$. Inspired by part (a) of Theorem 2.5, we compute the conformal map by minimizing

over the parameter $\tau$. If $\bar{\tau}$ is the minimizer, then $f_{\bar{\tau}}:V\rightarrow\mathbb{C}/\mathbb{Z}\oplus\mathbb{Z}{\bar{\tau}}$ would approximate a conformal map.

Our method is summarized as Algorithm 5.

Here are two numerical examples for harmonic maps and conformal maps respectively.

First we construct a lattice approximation $G=(V,E)$ as in Section 3, and then need to

(1) compute the discrete harmonic map $f$ from $G$ to $\mathbb{H}^{2}/\Gamma$, or equivalently, to a fundamental domain in the unit disk representation $\mathbb{H}^{2}=\{|z|<1\}$, and

(2) extend $f$ to $P$ by trilinear interpolation, and then $f|_{P}$ represents an approximation of the desired harmonic map.

Step (2) is kind of clear and simple, so let us discuss the details of step (1). First we need to ”cut” the lattice $G$ along $2\cdot genus(M)$ ”loops” $\gamma_{i}$, so that we can represent a discrete map from $G$ to a hyperbolic surface as a map from $V$ to a planar domain in the unit disk. Secondly we fix the homotopy class of discrete maps $f:G\rightarrow\mathbb{H}^{2}/\Gamma$, by requiring that each loop $\gamma_{i}$ corresponds to a certain deck transformation $\alpha_{i}$ in $\Gamma$.

Now computing the discrete harmonic map in the fixed homotopy class amounts to solve a system of nonlinear equations. For a vertex $i$ that is away from any loop $\gamma_{i}$, the balanced condition (1) gives us that

where $\vec{e}(x,y)$ denotes the unit tangent vector in $T_{x}\mathbb{H}^{2}$ indicating the direction to point $y$. Things become subtle when the vertex $i$ is near some loop $\gamma_{i}$. Assume the edge $ij$ pass through the loop $\gamma_{k}$, and $l_{ij}(f)$ is small. Then

where $\alpha_{ij}$ is a deck transformation in $\Gamma$, which is determined by

1. (1)

   the choice of the loops $\gamma_{i}$, and

2. (2)

   the choice of $\alpha_{i}$’s, and

3. (3)

   the relative position between $i$ and $\gamma_{k}$.

Properly choosing $\alpha_{i}$’s and computing $\alpha_{ij}$’s are involved, particularly for higher genus surfaces. One need to use polygonal representations for higher genus surfaces, and parameterizations of the deformation spaces of $\Gamma$. See [72][73] for references. Once we computed all the $\alpha_{ij}$’s, it remains to solve the following system of nonlinear equations.

The existence and uniqueness of the solution is guaranteed by Theorem 2.10. Notice that equation (10) indicates that $f(i)$ is the hyperbolic center of mass of $\alpha_{ij}(f(j))$’s where $j$ goes over all the neighbors of $i$. So one can use the so-called center of mass method to solve equation (10). This is an iterating method where in each iteration $f(i)$ is replaced by the mass-center of its neighbors $\alpha_{ij}(f(j))$, i.e., the $(n+1)$-th function $f_{n+1}(i)$ is determined by

It has been proved by Gaster-Loustau-Monsaingeon [59] that the center of mass method will converge to the unique discrete harmonic map. But the disadvantage is that there is no direct way to compute a hyperbolic mass-center and this iteration is slow. In our numerical experiments, we use the so-called cosh-center of mass method introduced also in Gaster-Loustau-Monsaingeon’s paper [59]. Instead of minimizing the exact discrete Dirichlet energy, we minimize

which is a good approximation of the discrete Dirichlet energy if all the edges lengths $l_{ij}$ are small. The balanced condition for critical points of the new energy ${\mathcal{E}}_{0}$ is