# Surface Ricci Flow Tutorial with C++ Source

### Smooth Case

Suppose $S$ is a surface embedded in the Euclidean space $R^3$, it has a natural induced Euclidean metric $\mathbf{g}$. Suppose $K$ is the Gaussian curvature induced by the metric $\mathbf{g}$. Choose a local coordinates, the metric tensor is $\mathbf{g}=(g_{ij})$. Then Hamilton's surface Ricci flow is defined as $\frac{ dg_{ij}(t) }{dt} = -2K(t)g_{ij}(t).$ Ricci flow deforms the Riemannian metric in a conformal way, $\mathbf{g}(t) = e^{2\lambda(t)}\mathbf{g}(0),$ the flow is given by $\frac{d\lambda(t)}{dt} = -K(t).$ The normalized Ricci flow is given by $\frac{ dg_{ij}(t) }{dt} = 2(\bar{K} -K(t) ) g_{ij}(t),$ where $\bar{K} = \frac{2\pi \chi(S)}{A(0)},$ $\chi(S)$ is the Euler characteristic of the surface, $A(0)$ is the total area at time $0$. Normalized surface Ricci flow converges to the metric which induces constant Gaussian curvature $\bar{K}$.

### 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. The discrete Riemannian metric is the edge length. Let $l_{ij}$ be the edge length for edge $[v_i,v_j]$, which satisfies the triangle inequality $l_{ij} + l_{jk} > l_{ki}.$ The discrete Gaussian curvature is the angle deficit. Suppose $v_i$ is an interior vertex of the mesh, $K_i = 2\pi - \sum_{jk}\theta_{i}^{jk},$ where $\theta_i^{jk}$ is the corner angle at vertex $v_i$ in triangle $[v_i,v_j,v_k]$. If $v_i$ is on the boundary, then $K_i = \pi - \sum_{jk}\theta_{i}^{jk}.$ Surface Ricci flow conformally deforms the Riemannian metric, which transforms infinitesimal circles to infinitesimal circles. In discrete case, we replace infinitesimal circles by circles with finite radii. Each vertex $v_i$ is associated with a circle, $c(v_i,r_i)$. The edge length $l_{ij}^2 = \gamma_i^2 + \gamma_j^2 + 2I_{ij} \gamma_i\gamma_j,$ The discrete conformal metric deformation varies the radii $\gamma_i$'s, but preserves all the weights $I_{ij}$'s. Let $u_i = \log \gamma_i.$ Let the target curvature at vertex $v_i$ be $\bar{K}_i$. Denote $\mathbf{u}=(u_1,u_2,\cdots,u_n)$. Define an Energy $E(\mathbf{u}) = \int^\mathbf{u} \sum_{i=1}^n (\bar{K}_i - K_i) du_i.$ Then $\nabla E(\mathbf{u}) = (\bar{K}_1 - K_1, \bar{K}_2 - K_2, \cdots, \bar{K}_n - K_n )^T.$ The Hessian matrix of $E$ is positive definite on the sublinear space $\sum_{i=1}^n u_i = 0.$ Each triangle $[v_i,v_j,v_k]$ has a circle, which is orthogonal to three vertex circles $c(v_i,\gamma_i)$, $c(v_j,\gamma_j)$ and $c(v_k,\gamma_k)$, which is called the power circle. Suppose an edge $[v_i,v_j]$ attaches to two triangles $[v_i,v_j,v_k]$ and $[v_j,v_i,v_k]$, two power circle centers are $o_k$ and $o_l$. Define the edge weight $w_{ij}$ as $w_{ij} = \frac{|o_l - o_k| }{|v_i-v_j|}.$ Then we obtain $d\mathbf{K} = \Delta d\mathbf{u},$ where $\Delta=(\Delta_{ij})$, $\Delta_{ij} = \left\{ \begin{array}{cl} -w_{ij} & [v_i,v_j] \in M\\ \sum_k w_{ik} & i= j \\ 0 & otherwise \end{array} \right.$ Namely, for each vertex $v_i$, $dK_i = \sum_{[v_i,v_j]\in M} w_{ij} ( du_i - du_j).$

### Algorithm

The computational algorithm is to use Newton's method to optimize the convex energy. The figure on top shows a conformal mapping from the human face to a rectangle. On the boundary, four corner vertices are selected, denoted as $\{v_0,v_1,v_2,v_3\}$. The target curvature on these 4 vertices equal to $\frac{\pi}{2}$, and zeros for any other vertex.
1. Set the target curvature $\bar{K}_i$ for each vertex $v_i$, which satisfies the Gauss-Bonnet theorem $\sum_{v_i \in M} K_i = 2\pi \chi(M).$
2. For each vertex $v_i$, estimate the circle radius $\gamma_i$; for each edge, estimate the weight $I_{ij}$.
3. Compute the edge length $l_{ij}$ from $\gamma_i$, $\gamma_j$ and $I_{ij}$, compute corner angles using cosine law.
4. Compute current Gaussian curvature on each vertex $K_i$. If $\max_i |\bar{K}_i - K_i | < \delta$, return.
5. For each face, compute the power circle; for each edge, compute the edge weight $w_{ij}$. Construct the Hessain matrix $\Delta$.
6. Solve linear system, $\mathbf{\bar{K}} - \mathbf{K} = \Delta \delta\mathbf{u},$
7. Update the radius vector $\mathbf{u} = \mathbf{u}+\delta \mathbf{u}, \gamma_i = e^{u_i}.$
8. Repeat from step 3 through step 7.

### Source Code

The code is written in C++ using VisualStudio on Windows platform.
If you encounter any problem during the compilation and running the code, please contact the author David Gu at $gu@cs.stonybrook.edu$.