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

## 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);
``````