Harmoic Map Tutorial with C++ Source
Smooth Case
Suppose $S$ is a surface embedded in the Euclidean space $R^3$. Furthermore, here we only consider the simplest topology: the surface $S$ is a topological disk.
A topological disk is a surface homeomorphic to the planar unit disk, namely, an orientable, genus zero surface with a single boundary.
The surface has an induced Euclidean metric $\mathbf{g}$.
Suppose $f$ is a function from the surface to $R$. The harmonic energy of the function is defined as
$$ E(f) = \int_S |\nabla_{\mathbf{g}} f |^2 dA $$. The minimizer of the harmonic energy is called a harmonic function. A harmonic function
satisfies the Laplace equation
\[
\Delta_{\mathbf{g}} f = 0.
\]
with Dirichlet boundary condition
\[
f|_{\partial S} = h,
\]
where $\partial S$ represents the boundary of the surface, $h$ is a given function.
Suppose we choose an isothermal coordiante system $(u,v)$, then the Riemannian metric is represented as
\[
\mathbf{g} = e^{2\lambda(u,v)} (du^2+dv^2),
\]
then the gradient operator is
\[
\nabla_{\mathbf{g}} = e^{-\lambda(u,v)} (\frac{\partial}{\partial u},\frac{\partial}{\partial v})^T,
\]
the area element is
\[
dA = e^{2\lambda(u,v)} du \wedge dv,
\]
and the Laplace-Beltrami operator is
\[
\Delta_{\mathbf{g}} = e^{-2\lambda(u,v)} (\frac{\partial^2}{\partial u^2} + \frac{\partial^2}{\partial v^2}).
\]
Suppose $\phi:S \to D$ is a mapping, such that $\phi(u,v)=(x(u,v),y(u,v))$, both $x(u,v)$ and $y(u,v)$ are harmonic functions, then $\phi$
is called a harmonic map. Intuitively, we can imaginge the surface $S$ is made of rubber sheet, then we flattned the surface onto the planar
unit disk, fixing the boundary. The mapping will naturally minimizes the membrane energy (stretching energy), and minizer is the harmonic map.
Harmonic map is very useful for engineering applications, one fundamental reason is because harmonic map gives diffeomrophism under some appropriate conditions.
Namely the mapping is bijective and smooth. This is formulated as Rado's theorem,
Rado's Theorem Suppose $\phi: S\to D$ is a harmonic map, $D$ is a planar convex domain. If the restriction of the map on the boundary
is a homeomorphism, then the interior mapping is a diffeomorphism.
Discrete Case
The surface is approximated by simplicial complex, namely triangle mesh $M$. A vertex is denoted as $v_i$; an oriented edge (halfedge) $[v_i,v_j]$,
where $v_i$ is the source, $v_j$ is the target; an oriented face is $[v_i,v_j,v_k]$, where $v_i$, $v_j$ and $v_k$ are sorted counter-clock-wisely.
A function $f: M\to R$ is defined on vertices. Suppose $p$ is a point inside a triangle $[v_i,v_j,v_k]$, then $p$ can be represented as the lineary combination of
three vertices,
\[
p = \lambda_i v_i + \lambda_j v_j + \lambda_k v_k,
\]
where $(\lambda_i,\lambda_j,\lambda_k)$ are called the barycentric coordinates of $p$, $\lambda_i + \lambda_j + \lambda_k = 1$. Barycentric coordinates can be computed
as follows
\[
\lambda_i = \frac{|[p,v_j,v_k]|}{|[v_i,v_j,v_k]|}
\]
where $|\cdot |$ represent the signed area of an oriented triangle. Then the function $f$ becomes a piece-wise linear function
\[
f(p) = \lambda_i f(v_i) + \lambda_j f(v_j) + \lambda_k f(v_k).
\]
By direct computation, it is easy to show that the harmonic energy of the function $f$ can be written as
\[
E(f) = \sum_{[v_i,v_j]\in M} w_{ij} (f(v_i) - f(v_j))^2.
\]
where $w_{ij}$ is called the cotangent weight on edge $[v_i,v_j]$. Suppose $[v_i,v_j]$ is adjacent to two faces $[v_i,v_j,v_k]$ and $[v_j,v_i,v_l]$, $\theta_k$ is the corner
angle in $[v_i,v_j,v_k]$ with $v_k$ as the apex, similalary $\theta_l$ is the corner angle in $[v_j,v_i,v_l]$ at vertex $v_l$, then the edge weight
\[
w_{ij} = \frac{1}{2}(\cot\theta_k + \cot \theta_l),
\]
if $[v_i,v_j]$ is an boundary edge, and only adjacent to one triangle $[v_i,v_j,v_k]$, then
\[
w_{ij} = \frac{1}{2} \cot \theta_k.
\]
The discrete Laplace equation becomes, for all interior vertex $v_i \not\in \partial M$,
\[
\sum_{[v_i,v_j]\in M} w_{ij} ( f(v_i) - f(v_j) ) = 0, \forall v_i\not\in \partial M,
\]
This is linar equation system.
Algorithm
The computational algorithm is as follows:
- Trace the mesh boundary, and get a sequence boundary vertices $\{v_0,v_1,\cdots, v_{n-1}\}$
- Set the image of boundary vertex $v_i$ as
\[
f(v_i) = (\cos\theta_i, \sin \theta_i), \theta_i = 2\pi \frac{\sum_{l=1}^i |v_{l}-v_{l-1}|}{\sum_{k=1}^n |v_{k}-v_{k-1}|},
\]
- For all edges $[v_i,v_j]$ compute the cotangent edge weight $w_{ij}$ using the formula above.
- For all invertior vertex $v_k$, go through all the surrounding vertices $v_j$, construct the discrete Laplace equation
\[
\sum_{[v_k,v_j]\in M} w_{kj}(f(v_k)-f(v_j)) = 0,
\]
- Solve the linear system.
The figure on top shows a discrete harmonic map result. Harmonic map is close to be conformal, so it preserves local shapes.
Source Code
The code is written in C++ using VisualStudio on Windows platform.
- Download the HarmonicMapper.zip file, unzip it.
- Enter Harmonicmapper/HarmonicMapper directory, double click on RiemannMapper.sln file. Compile the project using Visual Studio.
- Copy the output binary to HarmonicMapper/bin/HarmonicMapper.exe
- Go to HamronicMapper/demo/Alex directory, double click on HarmonicMap_test.bat
- The code will start running, and the results will be displayed simpleviewer.exe. The output file is alex.uv.m.
The triangle mesh viewer can be found here.
Contact
If you encounter any problem during the compilation and running the code, please contact the author David Gu at $gu@cs.stonybrook.edu$.
