Implementing "Linear Subspace Design for Real-Time Shape Deformation"

Alec Jacobson

May 12, 2015

weblog/

Here's a way to implement the linearly precise smoothness energy using gptoolbox in manner consistent with the paper.

First construct the usual cotangent laplacian L:

C = cotangent(V,F);
L = sparse( ...
  F(:,[2 3 1 3 1 2 2 3 1 3 1 2]), ...
  F(:,[3 1 2 2 3 1 2 3 1 3 1 2]), ...
  [C C -C -C], ...
  size(V,1),size(V,1));

This throws the cotangent of the angle across from each half-edge (i,j) to the off-diagonal entries L(i,j) and L(j,i) and minus that value to the diagonal entries L(i,i) and L(j,j).

You could also just use:

L = cotmatrix(V,F);

To compute the normal derivative matrix N, first find all of the half-edges on the boundary and the subscript indices into F so that if on_b_f(1) = f and on_b_c(1) = c then we know the half-edge across from vertex F(f,c) is on the boundary.

[~,on_b] = on_boundary(F);
[on_b_f,on_b_c] = find(on_b);

Now build lists of triplets (p,i,q) so that the half-edge i-->p is on the boundary and q is opposite it:

P = F(sub2ind(size(F),on_b_f,mod(on_b_c+1,3)+1));
I = F(sub2ind(size(F),on_b_f,mod(on_b_c,3)+1));
Q = F(sub2ind(size(F),on_b_f,on_b_c));

Collect the cotangents across from vertices i and j and throw accordingly at entries N(i,j), N(i,k), N(i,i), etc.:

CP = C(sub2ind(size(F),on_b_f,mod(on_b_c+1,3)+1));
CI = C(sub2ind(size(F),on_b_f,mod(on_b_c,3)+1));
N = sparse( ...
  [P P P P I I I I], ...
  [Q P Q I Q I Q P], ...
  [-CI CI -CP CP -CP CP -CI CI], ...
  size(V,1),size(V,1));

Now we can compute the "linearly precise Laplacian" K:

K = L+N;

Build the per-vertex mass matrix, and square our Laplacian to result in the bilaplacian smoothness energies quadratic form A:

M = massmatrix(V,F);
A = K' * (M\ K);