Category: tutorial week

By SGI Fellows Xinwen Ding and Miles Silberling-Cook

The authors of this article appreciate Derek Liu and Eris Zhang for their awesome lectures on the third day of SGI. The rest of this summary is based on what we learned from their lectures.

Introduction

One of the motivations for shape deformation is the desire of reproducing the motion of a given initial shape as time evolves. It is clear that modeling every piece of information associated with the initial shape can be hard and unnecessary. So, a simple but nice-looking deformation strategy is needed. During the lecture, SGI fellows were introduced two deformation strategies listed below.

• Energy-based Deformation: This method turns the deformation problem into an optimization problem, where possible objective functions can be Dirichlet Energy (\(E = \int_{M} ||\nabla f||^2 dx\)), Laplacian Energy (\(E = \int_{M} ||\Delta f||^2 dx\)), or Conformal Energy. It is worth mentioning that the nonlinear optimization problems are, in general, difficult to solve.
• Direct Method: We mainly refer to a method called Linear Blend Skinning (LBS). In this method, we predefine a sequence of movable handles. Furthermore, we assume the movement of every handle will contribute linearly to the final coordinate of our input vertex. Based on this assumption, we write down the equation for calculating the target position: \(u_i = \displaystyle\sum_{j = 1}^{\text{# of handles}} w_{ij} T_j v_i\), where \(v_i, u_i\) are the \(i^{th}\) input and output vertices, respectively, \(w_{ij}\) is the weight, and \(T_j\) is the tranformation matrix of the \(j^{th}\) handle recording rotation and transformation in terms of homogeneous coordinates.

The remaining part of this summary will focus on the Linear Blend Skinning (LBS) method.

Workflow

The process of a 2D LBS deformation can be divided into three steps:

• Step 1: Define the handles.
• Step 2: Compute the weights \(w_{ij}\) and store the tranformation matrix \(T_j = \begin{bmatrix}\cos\theta & -\sin\theta & \Delta x \\ \sin\theta & \cos\theta & \Delta y\end{bmatrix}\), where the first two columns of \(T_j\) denote the rotation of the \(j^{th}\) handle, and the last column denotes the shift. Notice that in this case the \(i^{th}\) vertex should be \(v_i = \begin{bmatrix} v_x \\ v_y \\ 1\end{bmatrix}\).
• Step 3: Put all parts together and compute \(u_i = \displaystyle\sum_{j = 1}^{\text{# of handles}} w_{ij} T_j v_i\).

A concrete example will be given in the final section.

Weight Calculation

Requirements on weights:

To minimize the artifacts in our deformation, we need to be smart when we generate the weights. Briefly speaking, we want our weights to satisfy the following criteria:

• Smooth: the weight function should be \(C^{1}\) differentiable. This can be achieved by introducting the Laplacian operator.
• Local: the effective area in a 2D domain should be within some neighborhood of the handle.
• Shape aware: the effective area in a 2D domain should not cross the boundary of our domain of interest.
• Sum up to 1.

Let \(\Omega\) be the domain of the interest, and \(w\) be the weight. Based on the aforementioned requirements, we can use the following minimization problem to characterize the weights.

\begin{equation}
\begin{aligned}
&\min_{w} \int_{\Omega} || \Delta w||^2 dx \\
&\textrm{s.t.} w(j^{th} \text{handle}) = e_j \hspace{20pt} \text{($e_j$ is a standard basis vector)}
\end{aligned}
\end{equation}

Discretize the Laplacian Operator

One may notice that in the objective function of the above optimization problem, we introduced a Laplacian operator. Luckily, we can play some math tricks to discretize the Laplacian operator, namely \(\Delta = M^{-1}C\), where \(M\) is the diagonal mass matrix, and \(C\) is the cotangent formula for the Laplacian.

With this knowledge in mind and with some knowledge from Finite Element Method, we can rewrite the objective function into a quadratic form.

\[
\int_{\Omega} || \Delta w||^2 dx
= ||M^{-1}Cw||_{M}^2
= (M^{-1} C w)^{T} M (M^{-1} C w)
= w^{T} C^{T} M^{-1} C w
\]

Let \(B = C^{T} M^{-1} C\), then we have \(\int_{\Omega} || \Delta w||^2 dx = w^{T} C^{T} M^{-1} C w = w^{T} B w\). Moreoever, we name the matrix \(B\) as bilaplacian or square Laplacian. Hence, the optimization problem is rewritten as follows:

\begin{equation}
\begin{aligned}
&\min_{w} w^{T} B w \\
&\textrm{s.t.} w(j^{th} \text{handle}) = e_j \hspace{20pt} \text{($e_j$ is a standard basis vector)}
\end{aligned}
\end{equation}

Solution Derivation to the Optimization Problem

Recall that \(w_{ij}\) is the degree to which vertex \(i\) is effected the motion of handle \(j\). To form our constraint, we pick a vertex for each handle that will only be affected by that handle. This constraint is much easier to express when the fixed weightings are separated from the free weightings, so we re-organize the weightings as \(w = [x, y]\). In this notation, \(x\) contains the free variables for the non-handle vertices, and \(y\) contains the fixed weighting for each handle. The boundary can now be concisely stated as \(y = \mathbb{I}\).

We now solve for the \(x\) that minimizes \([x \mid y]^T B [x \mid y]\). Following Derek Liu and Eris Zhang’s lead, we use subscripts to express the indices of \(B\) pertaining to either \(x\) or \(y\).

\begin{align}
E
&= w^T B w \\
&= [x \mid y]^T B [x \mid y] \\
&= [x \mid 0]^T \begin{bmatrix}B_{x,x} & 0 \\ 0 & 0 \end{bmatrix} [x \mid 0]
+ [x \mid 0]^T \begin{bmatrix}0 & B_{x,y} \\ 0 & 0 \end{bmatrix} [0 \mid y]
+ [0 \mid y]^T \begin{bmatrix}0 & 0 \\ B_{y,x} & 0 \end{bmatrix} [x \mid 0]
+ [0 \mid y]^T \begin{bmatrix}0 & 0 \\ 0 & B_{y,x} \end{bmatrix} [0 \mid y] \\
&= x^T B_{x,x} x + x^T B_{x,y} y + y^T B_{y,x} x + y^T B_{y,y} y \\
&=
\underbrace{ x^T B_{x,x} x }_\text{Quadratic}
+ \underbrace{ x^T \left(B_{x,y} + B_{y,x}^T \right) y }_\text{Linear}
+ \underbrace{ y^T B_{y,y} y }_\text{Constant}
\end{align}

As expected, this yields another quadratic equation. Differentiating gives:

\begin{equation}
\frac{\partial E}{\partial x}
= 2 A_{x,x,} x + \left(B_{x,y} + B_{y,x}^T \right) y
\end{equation}

and the extremizer is therefore \(x = \left( -2A\right)^{-1}\left(B_{x,y} + B_{y,x}^T \right) y\).

Code Implementation and Acceleration

Implementation

To calculate the solution in Matlab, see the following.

function W = compute_skinning_weight(V,F,b)
% V: Vertices
% F: Triangles
% b: Handle vertex indices

% Compute the Laplacian using the cotangent-matrix.
L = -cotmatrix(V,F);

% Compute the mass matrix.
M = massmatrix(V,F);

% Compute the BiLaplacian
B = L*(M\L);

% Create a list of non-handle verticies
a = setdiff(1:size(V,1), b);

% Initialize the fixed handle constraint
y = eye(length(b));

% Compute the minimum
A = Q(a,a);
B = (Q(a,b) + Q(b,a).') * y;
x = (-2*A) \ B;

% Assemble the final weight matrix
W(a,:) = x;
W(b,:) = y;
end

Once the weightings are done, we can propagate the transformation of the handles to the mesh (using homogeneous coordinates) as follows:

function [U] = linear_blend_skinning(V, T, W)
% V: Vertices
% T: Handle transforms
% W: Weighting matrix

U = zeros(size(V));

for i=1:size(V)

for j=1:size(T,3)
U(i,:) = U(i,:) + ( W(i,j) * [V(i,:) 1] * T(:,:,j).' );
end

end

end

Acceleration

Although the solution above guarantees correctness, this is a fairly inefficient implementation. Some of this computation can be done offline, and the nested loops are not required. The blending weights and the vertex positions can be compiled into a single matrix ahead of time, like so:

function C = compile_blending_weights(V, N, W)
% V: Vertices
% N: Number of handles
% W: Weighting matrix

for i=1:length(V)
v = repmat( [V(i,:), 1], 1, N)';
wi = reshape(repmat(W(i, :), 3, 1), [], 1);
C(:,i) = v.*wi;
end
end

Then propagating the transform is as simple as composing it with the precomputed weightings.

function [U] = linear_blend_skinning(T, C)
% T: Handle transforms
% C: Compiled blending weights

B = reshape(T,[], size(T, 3)*size(T, 2));
U = (B*C)';
end

Shout-out to SGI Fellow Dimitry Kachkovski for providing this extra-efficient implementation!

Example: 2D Gingerbread Man Deformation

In this example, we try to implement the whole method over a gingerbread man domain. Following the define handle – computing weights – combine and transform workflow, we provide a series of video examples.

Step 1: Define handles

Step 2: Computing Weights

Step 3: Combine and Transform