As-rigid-as-possible (ARAP) skinning

Now that we've talked about Laplacian surface editing, it's a good time to introduce another technique for surface deformation. We're going to talk about As-Rigid-As-Possible Surface Modeling (Sorkine et al., 2007), then we'll look at how to make use of the technique alongside traditional skeletal skinning.

Let's start with the definition of a cell. A cell $C_i$ contains the vertex $i$ as well as $i$'s one-ring neighborhood $j \in N(i)$. Rigidly transforming $C_i$ to produce $C_i'$ means there exists a rotation matrix $R_i$ such that:

$$
\begin{equation}
p_i' - p_j'= R_i (p_i - p_j), \forall j \in N(i).
\end{equation}
\label{eq:rigidrotation}
$$

But that would mean every neighbor of $i$ gets the exact same transformation, so the deformation will probably not look smooth. Instead, what we can do is look to minimize the error functional:

$$
\begin{equation}
E(C_i, C_i') = \sum_{j \in N(i)} w_{ij} ||(p_i' - p_j') - R_i (p_i - p_j) ||^2
\end{equation}
\label{eq:araperror}
$$

We can figure out a suitable weighting method for $w_{ij}$, but we're still dealing with two unknowns: the updated vertex positions $p_i'$ and the rotation matrices $R_i$. We can reform Equation$~\ref{eq:araperror}$ to determine the optimal $R_i$ for a given set of $p_i',$ which involves finding the singular value decomposition of a covariance matrix built out of $i$'s edges (see the paper for full details on this because it's really neat), but I can also just take a page out of the elastic implicit skinning book and cheat by computing $R_i$ for every $i$ by using skinning weights and skeleton transforms. I'm more concerned with finding the optimal $p_i'$ given $R_i$, so let's look at that!

We can compute the overall energy of a deformed mesh $M'$ by iterating over Equation$~\ref{eq:araperror}$ for all $i$:

$$
\begin{equation}
E(M') = \sum_{i=1}^{n} w_i E(C_i, C_i') = \
\sum_{i=1}^{n} w_i \sum_{j \in N(i)} w_{ij} ||(p_i' - p_j') - R_i (p_i - p_j) ||^2
\end{equation}
\label{eq:mesherror}
$$

The value of Equation$~\ref{eq:araperror}$ is an integrated quantity, which means the cell $C_i$'s energy is proportional to its area, and we can get away with settings $w_i=1$. Alternately, we could set $w_i=A_i$, the Voronoi area of $C_i$, and $w_{ij}'=\frac{1}{A_i}w_{ij}$. But then the cell area would cancel out, so why bother? Choice of $w_{ij}$ is somewhat flexible, but we want the per-edge weights to be useful for any kind of mesh we encounter, so the standard here is to use cotangent weights:

$$
\begin{equation}
w_{ij} = \frac{1}{2}(\cot \alpha_{ij} + \cot \beta_{ij}),
\end{equation}
\label{eq:cotangentweights}
$$

where $\alpha_{ij}$ and $\beta_{ij}$ are the angles opposite of the edge $(i, j) \in M$. So $\alpha_{ij}$ comes from one triangle and $\beta_{ij}$ comes from another.

Now that we have weights defined, we can determine how to compute new vertex positions $p_i'$ using by taking the partial derivative of $E(M')$ with respect to $p_i'$:

$$
\begin{equation}
\frac{\partial E(M')}{\partial p_i'} = \frac{\partial}{\partial p_i'} \Big( \sum_{j \in N(i)} w_{ij} || (p_i' - p_j') - R_i (p_i - p_j) ||^2 + \
\sum_{j \in N(i)} w_{ji} || (p_j' - p_i') - R_j (p_j - p_i) ||^2 \Big) = \
\sum_{j \in N(i)} 2w_{ij} ((p_i' - p_j') - R_i (p_i - p_j)) + \
\sum_{j \in N(i)} -2w_{ji} ((p_j' - p_i') - R_j (p_j - p_i)) = \
\sum_{j \in N(i)} 4w_{ij} \big( (p_i' - p_j') - \frac{1}{2}(R_i + R_j) (p_i - p_j) \big)
\end{equation}
\label{eq:energygradient}
$$

The last line of that equation is possible because $w_{ij} = w_{ji}$. If we set $\frac{\partial E(M')}{\partial p_i'} = 0,$ we can separate terms and get:

$$
\begin{equation}
\sum_{j \in N(i)} w_{ij} (p_i' - p_j') = \sum_{j \in N(i)} \frac{w_{ij}}{2} (R_i + R_j)(p_i - p_j)
\end{equation}
\label{eq:energygradient2}
$$

The left-hand side of Equation$~\ref{eq:energygradient2}$ is simply $Lp'$, where $L$ is the Laplace-Beltrami operator of $M$ and $p'$ are the new vertex positions. The right-hand side can be computed for each vertex $i$ given the rotation matrices for each vertex (which we can obtain from linear blend skinning) and stored in a vector $b$, giving us a very familiar-looking sparse linear system:

$$
\begin{equation}
Lp' = b
\end{equation}
\label{eq:lbp}
$$

To enforce constraints (static and handle anchors as discussed in the LSE post), we first specify the indices of constrained vertices as $k \in F$. Then we erase each $k$ row and column from $L$, and set the $k$-th entries of $b$ to the desired position. $L$ is symmetric positive definite, and the original authors use a Sparse Cholesky factorization with fill-reducing reordering to perform the solve. In my own implementation, I used Eigen's SPQR solver, which is just QR factorization accelerated with the SuiteSparse library.

ARAP and skinning

So how does ARAP fit into skinning? As you may recall, the implicit skinning folks use it to formalize their tangential relaxation energy, which is what helps give the mesh such a smooth global deformation - when a character twists its upper body, the skin on the lower body gently deforms to accommodate. Now that's snazzy!

We can apply ARAP to a mesh after any skinning method as a post-processing step. In my code, the skinning vertex shaders run once for the mesh vertices (using glDrawElements(GL_POINTS, numPoints, ...)) and the output (position, rotation, normal) gets saved to a transform feedback buffer. Having precomputed $L$ and primed a solver, I just have to update the $b$ vector and solve to get new vertex positions.

But that's not all! We could formulate a different energy functional for ARAP and minimize to a different goal. Incremental rotation skinning operates by 1) figuring out how surface vertices attach to an underling skeleton, and 2) rotating vertices about their attachment point on the skeleton in order to smoothly deform the surface. If a vertex is right in the middle of a bone (not close to the joints on either side), it will probably just receive the same rotation as its bone. If it's within a joint's radius, however, then it will probably receive additional rotation to make the deformation smooth.

This approach to skinning leaves us with a scale direction $s_i$ for every vertex. We might then rewrite Equation$~\ref{eq:energygradient2}$ to be:

$$
\begin{equation}
\sum_{j \in N(i)} w_{ij} (s_i' - s_j') = \sum_{j \in N(i)} \frac{w_{ij}}{2} (R_i + R_j)(s_i - s_j),
\end{equation}
\label{eq:sdenergygradient}
$$

leaving us with a new sparse system to solve:

$$
\begin{equation}
Ls' = b
\end{equation}
\label{eq:sd_lbp}
$$