How to replicate any face in FaceGen

To replicate any face (a mesh of vertices in three dimensions) in FaceGen, first you need to:

1. Align the mesh vertices with those of FaceGen, using key points (eyes, nose, lips...) and rigid registration algorithms.
2. Find the correspondence between the vertices of the mesh and those of FaceGen, using non-rigid registration algorithms.
3. Adjust the position, orientation, and scale of the mesh to make it easier to replicate with FaceGen.

Now what? Recall that FaceGen generates a face through a linear combination of modes (deviations):

\mathbf{v}_i = \mathbf{\bar v}_i + \sum_{j=1}^{N_s} s_j \mathbf{v}_i^{(j)} + \sum_{j=1}^{N_g} a_j \mathbf{u}_i^{(j)}

The approach is now reversed: how do we determine the symmetric ($s_j$) and asymmetric ($a_j$) coefficients given the vertices we want to obtain? First, we pack all the vertices into column matrices of size $3n \times 1$; that is:

$$ \begin{align} v &= \begin{bmatrix} v_{1x} & v_{1y} & v_{1z} & \cdots & v_{nx} & v_{ny} & v_{nz}\end{bmatrix}^\top \ \bar v &= \begin{bmatrix} \bar v_{1x} & \bar v_{1y} & \bar v_{1z} & \cdots & \bar v_{nx} & \bar v_{ny} & \bar v_{nz}\end{bmatrix}^\top \ D^{(j)} &= \begin{bmatrix} v_{1x}^{(j)} & v_{1y}^{(j)} & v_{1z}^{(j)} & \cdots & v_{nx}^{(j)} & v_{ny}^{(j)} & v_{nz}^{(j)}\end{bmatrix}^\top ~~\forall j \in {1,\ldots, N_s, \ldots, N_s + N_a} \end{align} $$

By stacking all the deviations $D^{(j)}$ into a single matrix, we obtain a matrix $D$ of size $3n \times m$, with $m = N_s$ if we only use symmetric modes or $m = N_s + N_a$ if we also use asymmetric modes. With all this, we can rewrite the first equation as a linear system:

$$ v - \bar{v} = D s $$

Where $s$ is the vector (column matrix) that contains all the symmetric (and asymmetric) coefficients. In general, this linear system has no exact solution. But we can find a quasi-solution, the one that minimizes the quadratic error given by:

$$ E = (v- \bar v-Ds)^\top (v -\bar v -Ds) $$

Whose optimum is:

$$ s_\min = (D^\intercal D)^{-1} D^\intercal (v - \bar v) $$

Let $D^\dagger = (D^\top D)^{-1} D^\top$ be the pseudoinverse of $D$. The minimum error is therefore:

$$ E = (v-\bar v)^\top (I - DD^\dagger)^\top(I - DD^\dagger) (v-\bar v) $$

Let $P = D D^\dagger$ be the (symmetric) projection matrix onto the subspace of faces generatable by FaceGen. Thus:

$$ E = (v - \bar v)^\intercal (I-P)^2 (v - \bar v) = (v - \bar v)^\intercal (I-P) (v - \bar v) = (v - \bar v)^\intercal N (v - \bar v) $$

Where $N = I-P$ is the orthogonal projection. Of course, if the target face was generated by the same FaceGen model error would be zero.

Setting the mesh for replication

Actually, that is not the best possible solution, but the best as the target mesh is currently aligned, and it was aligned to find the correspondence between vertices, not to minimize this error. For this reason, it is necessary to adjust position, orientation, and scale of the mesh to improve the solution and obtain the best possible replica. To do this, we must first decompose the matrix product into blocks of three elements:

$$ E = \sum_{ij} (v_i - \bar v_i)^\top N_{ij} (v_j - \bar v_j) $$

Here, $v_i$ are the (column) coordinate vectors of each vertex, and $N_{ij}$ is the $3\times 3$ block corresponding to vertices $i, j$, which satisfies $N_{ij}^\top = N_{ji}$ due to the symmetry of the projection. In short, we have unpacked the vertices to be able to apply the necessary transformation. This transformation must be rigid so as not to change the shape of the target face; that is:

$$ v _i^\prime = \lambda R v_i + t $$

We obtain the new mesh vertices $v_i^\prime$ from the original vertices by applying a scale change $\lambda$, a rotation $R$, and a translation $t$ (the translation, in 3 dimensions). The new quadratic error will then be:

$$ E = \sum_{ij} (\lambda R v_i + t - \bar v_i)^\intercal N_{ij} (\lambda R v_j + t - \bar v_j)\ $$

These transformations depend on 7 parameters $$(x, y, z, \alpha, \beta, \gamma, \eta)$$:

• $t = x e_x + y e_y + z e_z$
• $R = R_z(\gamma) R_y(\beta) R_x(\alpha) = \exp(\gamma S_z)\exp(\beta S_y)\exp(\alpha S_x)$
• $\lambda = 1 + \eta$

These are the displacements $x, y, z$, the Euler rotation angles $\alpha, \beta, \gamma$, and the scale variation $\eta$. Additionally, $e_x, e_y, e_z$ are the unit vectors (column matrices) of each axis, while $S_x, S_y, S_z$ are the matrices that generate rotations about the respective axes:

\begin{bmatrix}
e_x & e_y & e_z
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}
\quad
S_x =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & -1 \\
0 & 1 & 0
\end{bmatrix}
\quad
S_y =
\begin{bmatrix}
0 & 0 & 1 \\
0 & 0 & 0 \\
-1 & 0 & 0
\end{bmatrix}
\quad
S_z =
\begin{bmatrix}
0 & -1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}

Since the vertices had been previously aligned, it is to be expected that the optimal transformations will be small, as will the value of these parameters, so a second-order approximation of the error will suffice; that is:

$$ \begin{align} E \simeq \sum_{ij} &[(v_i-\bar v_i)^\top + x e_x^\top + y e_y^\top + z e_z^\top + u_i^\top(\eta I + \alpha S_x^\top + \beta S_y^\top + \gamma S_z^\top) \ &+v_i^\top(\eta\alpha S_x^\top + \eta \beta S_y^\top +\eta\gamma S_z^\top +\gamma\beta S_y^\top S_z^\top + \gamma\alpha S_x^\top S_z^\top + \beta\alpha S_x^\top S_y^\top) \ &+\tfrac{1}{2}v_i^\top(\alpha^2 (S_x^2)^\top +\beta^2 (S_y^2)^\top +\gamma^2 (S_z^2)^\top)]\ &N_{ij}\ &[(v_j-\bar v_j) + x e_x + y e_y + z e_z +(\eta I + \alpha S_x + \beta S_y + \gamma S_z)v_j \ &+(\eta\alpha S_x + \eta \beta S_y +\eta\gamma S_z +\gamma\beta S_z S_y + \gamma\alpha S_z S_x + \beta\alpha S_y S_x)v_j \ &+\tfrac{1}{2}(\alpha^2 S_x^2 +\beta^2S_y^2 +\gamma^2 S_z^2)v_j]\ \end{align} $$

Let $\Delta_i = v_i - \bar v_i$, the difference between the vertices of the original face and the reference. Taking advantage of the symmetry of the projection and the fact that the indices $$i,j$$ are dummy, the Hessian matrix (second derivatives) of this quadratic function is:

$$ H_E = \sum_{ij}\begin{bmatrix} e_x^\top N_{ij} e_x & e_x^\top N_{ij} e_y & e_x^\top N_{ij} e_z & e_x^\top N_{ij} S_x v_j & e_x^\top N_{ij} S_y v_j & e_x^\top N_{ij} S_z v_j & e_x^\top N_{ij} v_j \ e_y^\top N_{ij} e_x & e_y^\top N_{ij} e_y & e_y^\top N_{ij} e_z & e_y^\top N_{ij} S_x v_j & e_y^\top N_{ij} S_y v_j & e_y^\top N_{ij} S_z v_j & e_y^\top N_{ij} v_j \ e_z^\top N_{ij} e_x & e_z^\top N_{ij} e_y & e_z^\top N_{ij} e_z & e_z^\top N_{ij} S_x v_j & e_z^\top N_{ij} S_y v_j & e_z^\top N_{ij} S_z v_j & e_z^\top N_{ij} v_j \ v_i^\top S_x^\top N_{ij} e_x & v_i^\top S_x^\top N_{ij} e_y & v_i^\top S_x^\top N_{ij} e_z & v_i^\top S_x^\top N_{ij} S_x v_j + v_i^\top (S_x^2)^\top N_{ij} \Delta_j & v_i^\top S_x^\top N_{ij} S_y v_j + v_i^\top S_x^\top S_y^\top N_{ij} \Delta_j & v_i^\top S_x^\top N_{ij} S_z v_j + v_i^\top S_x^\top S_z^\top N_{ij} \Delta_j & v_i^\top S_x^\top N_{ij} v_j + v_i^\top S_x^\top N_{ij} \Delta_j \ v_i^\top S_y^\top N_{ij} e_x & v_i^\top S_y^\top N_{ij} e_y & v_i^\top S_y^\top N_{ij} e_z & v_i^\top S_y^\top N_{ij} S_x v_j + v_i^\top S_x^\top S_y^\top N_{ij} \Delta_j & v_i^\top S_y^\top N_{ij} S_y v_j + v_i^\top (S_y^2)^\top N_{ij} \Delta_j & v_i^\top S_y^\top N_{ij} S_z v_j + v_i^\top S_y^\top S_z^\top N_{ij} \Delta_j & v_i^\top S_y^\top N_{ij} v_j + v_i^\top S_y^\top N_{ij} \Delta_j \ v_i^\top S_z^\top N_{ij} e_x & v_i^\top S_z^\top N_{ij} e_y & v_i^\top S_z^\top N_{ij} e_z & v_i^\top S_z^\top N_{ij} S_x v_j + v_i^\top S_x^\top S_z^\top N_{ij} \Delta_j & v_i^\top S_z^\top N_{ij} S_y v_j + v_i^\top S_y^\top S_z^\top N_{ij} \Delta_j & v_i^\top S_z^\top N_{ij} S_z v_j + v_i^\top (S_z^2)^\top N_{ij} \Delta_j & v_i^\top S_z^\top N_{ij} v_j + v_i^\top S_z^\top N_{ij} \Delta_j \ v_i^\top N_{ij} e_x & v_i^\top N_{ij} e_y & v_i^\top N_{ij} e_z & v_i^\top N_{ij} S_x v_j + v_i^\top S_x^\top N_{ij} \Delta_j & v_i^\top N_{ij} S_y v_j + v_i^\top S_y^\top N_{ij} \Delta_j & v_i^\top N_{ij} S_z v_j + v_i^\top S_z^\top N_{ij} \Delta_j & v_i^\top N_{ij} v_j \end{bmatrix} $$

Assuming that $\Delta_i$ are small, that is, once aligned there is little distance between the vertices of both faces, we can reduce the Hessian matrix to:

$$ H_E \simeq \sum_{ij}\begin{bmatrix} e_x^\top N_{ij} e_x & e_x^\top N_{ij} e_y & e_x^\top N_{ij} e_z & e_x^\top N_{ij} S_x v_j & e_x^\top N_{ij} S_y v_j & e_x^\top N_{ij} S_z v_j & e_x^\top N_{ij} v_j \ e_y^\top N_{ij} e_x & e_y^\top N_{ij} e_y & e_y^\top N_{ij} e_z & e_y^\top N_{ij} S_x v_j & e_y^\top N_{ij} S_y v_j & e_y^\top N_{ij} S_z v_j & e_y^\top N_{ij} v_j \ e_z^\top N_{ij} e_x & e_z^\top N_{ij} e_y & e_z^\top N_{ij} e_z & e_z^\top N_{ij} S_x v_j & e_z^\top N_{ij} S_y v_j & e_z^\top N_{ij} S_z v_j & e_z^\top N_{ij} v_j \ v_i^\top S_x^\top N_{ij} e_x & v_i^\top S_x^\top N_{ij} e_y & v_i^\top S_x^\top N_{ij} e_z & v_i^\top S_x^\top N_{ij} S_x v_j & v_i^\top S_x^\top N_{ij} S_y v_j & v_i^\top S_x^\top N_{ij} S_z v_j & v_i^\top S_x^\top N_{ij} v_j \ v_i^\top S_y^\top N_{ij} e_x & v_i^\top S_y^\top N_{ij} e_y & v_i^\top S_y^\top N_{ij} e_z & v_i^\top S_y^\top N_{ij} S_x v_j & v_i^\top S_y^\top N_{ij} S_y v_j & v_i^\top S_y^\top N_{ij} S_z v_j & v_i^\top S_y^\top N_{ij} v_j\ v_i^\top S_z^\top N_{ij} e_x & v_i^\top S_z^\top N_{ij} e_y & v_i^\top S_z^\top N_{ij} e_z & v_i^\top S_z^\top N_{ij} S_x v_j & v_i^\top S_z^\top N_{ij} S_y v_j & v_i^\top S_z^\top N_{ij} S_z v_j & v_i^\top S_z^\top N_{ij} v_j \ v_i^\top N_{ij} e_x & v_i^\top N_{ij} e_y & v_i^\top N_{ij} e_z & v_i^\top N_{ij} S_x v_j & v_i^\top N_{ij} S_y v_j & v_i^\top N_{ij} S_z v_j & v_i^\top N_{ij} v_j \end{bmatrix} $$

Thus, it can be factorized as:

$$ H_E \simeq \sum_{ij}\begin{bmatrix} e_x & e_y & e_z & S_x v_i & S_y v_i & S_z v_i & v_i \end{bmatrix}^\top N_{ij} \begin{bmatrix} e_x & e_y & e_z & S_x v_j & S_y v_j & S_z v_j & v_j \end{bmatrix} $$

And in this way, we can repack the vertices:

$$ H_E \simeq \begin{bmatrix} e_x & e_y & e_z & S_x v_1 & S_y v_1 & S_z v_1 & v_1 \ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \ e_x & e_y & e_z & S_x v_n & S_y v_n & S_z v_n & v_n \ \end{bmatrix}^\top N \begin{bmatrix} e_x & e_y & e_z & S_x v_1 & S_y v_1 & S_z v_1 & v_1 \ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \ e_x & e_y & e_z & S_x v_n & S_y v_n & S_z v_n & v_n \ \end{bmatrix} $$

As for the gradient at point 0 (i.e., without transformation), it can also be factorized:

$$ \nabla_E = \sum_{ij} \begin{bmatrix} e_x^\top N_{ij}\Delta_j \ e_y^\top N_{ij}\Delta_j \ e_z^\top N_{ij}\Delta_j \ v_i^\top L_x^\top N_{ij}\Delta_j \ v_i^\top L_y^\top N_{ij}\Delta_j \ v_i^\top L_z^\top N_{ij}\Delta_j \ v_i^\top N_{ij}\Delta_j\ \end{bmatrix} = \sum_{ij}\begin{bmatrix} e_x & e_y & e_z & S_x v_i & S_y v_i & S_z v_i & v_i \end{bmatrix}^\top N_{ij} \Delta_j $$

As before, we pack to obtain:

$$ \nabla_E = \begin{bmatrix} e_x & e_y & e_z & S_x v_1 & S_y v_1 & S_z v_1 & v_1 \ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \ e_x & e_y & e_z & S_x v_n & S_y v_n & S_z v_n & v_n \ \end{bmatrix}^\top N \begin{bmatrix} \Delta_1\ \vdots \ \Delta_n \end{bmatrix} $$

Finally, finding the combination $(x, y, z, \alpha, \beta, \gamma, \eta)$ that minimizes the error (approximately) consists of solving the linear system:

$$ H_E \begin{bmatrix} x & y & z & \alpha & \beta & \gamma & \eta \ \end{bmatrix}^\top = -\nabla_E $$

If the Hessian matrix is invertible, the only possible solution is:

$$ \begin{bmatrix} x & y & z & \alpha & \beta & \gamma & \eta \end{bmatrix}^\top_\min = -H_E^{-1}\nabla_E $$

Once the solution is obtained, the associated transformations are applied and, if desired, since the solution is only approximate, it is possible to repeat the procedure with the new vertices until a satisfactory solution is obtained. In practice, one or two iterations are sufficient, unless the target face is an aberration.

Variants

In the definition of quadratic error, it is implicitly assumed that all points are equally important. This does not necessarily have to be the case; we can make some more important than others, such as the key points mentioned earlier. To achieve this, we express the quadratic error per vertex and assign distinct weights to each one:

$$ E = \sum_i w_i^2 |v_i - \bar v_i - \sum_j D_{ij} s_j|^2 = \sum_i |w_iv_i - w_i\bar v_i - \sum_j w_i D_{ij} s_j|^2 $$

We see that error is now equivalent to the original one after multiplying all vertices and deviations by these weights. Furthermore, when introducing rigid transformations:

$$ E = \sum_{ij} (\lambda R w_iv_i + w_it - w_i\bar v_i)^\intercal N_{ij} (\lambda R w_jv_j + w_jt - w_j\bar v_j) $$

We observe that the weights also affect the translation. However:

$$w_i t = x(w_i e_x) + y(w_i e_y) + z(w_i e_z)$$

reveals that it is sufficient to apply these weights to the unit vectors of the axes as well, to make the problem equivalent to the original one. In short, after making these changes, the problem is solved exactly the same way. Another possibility is that we seek not the most exact solution but the one that deviates least from the mean. That is, minimize:

$$E^\prime = s^\top s = (v-\bar v)^\top D^{\dagger\top}D^{\dagger}(v-\bar v)$$

By defining $N^\prime = D^{\dagger\top} D^\dagger$, we recover the same problem, and once again, the way to solve it remains the same.