This is a cute little example we didn't have time to put into our upcoming SIGGRAPH paper Linear Subspace Design for Real-Time Shape Deformation.
In this example I first decimate the high resolution (blue) octopus to create the low-res orange octopus. I'm careful to use a decimate that keeps vertices at (or at least near) their original positions. Then I compute "biharmonic coordinates" for vertices of the blue mesh with respect to the vertices of the coarse orange mesh. Notice that the orange mesh does not need to be an enclosing cage. Rather its points just need to lie in or on the blue octopus (there're even ways to deal with points outside the blue octopus, but right now I'm skipping that). These weights are computed in gptoolbox with:
b = knnsearch(orange_V,blue_tet_V);
W = biharmonic_coordinates(blue_tet_V,blue_tet_T,b);
Then I run a sort of co-rotational linear-elasticity simulation on the orange mesh with collisions activated for the floor. The blue mesh tracks the orange mesh by simply interpolating the orange vertex positions via the coordinates:
blue_V = W * orange_V;
Simple as that. Smooth physics "skinning" with a matrix-matrix multiplication. This is due to the affine precision of our coordinates and their squared Laplacian energy minimization. The "standard" approach to physics "skinning" would be to create the orange tet-mesh so that it envelops the blue mesh and then use piecewise-linear (i.e. non-smooth) interpolation to move the blue mesh vertices.
Update: I should make it super clear that this is not a subspace (aka model reduction) simulation. The simulation of the orange mesh is not aware of the blue mesh at all. Sometimes subspace physics would be more appropriate, though other times pushing the subspace through all parts of the simulation pipeline is impractical/impossible. Hence this prototypical result.
For posterity here's the entire example code:
% load hi-res mesh
[uV,uF] = load_mesh('~/Dropbox/models/octopus-med.ply');
% load low-res mesh
[~,CV,CF] = qslim(uV,uF,900,'QSlimFlags','-O 0');
[CV,CF] = meshfix(CV,CF);
% Find closest points
b = knnsearch(uV,CV);
CV = uV(b,:);
% Tet-mesh fine resolution mesh and compute weights
[dV,dT,dF] = tetgen(uV,uF,'Flags','-q2');
W = biharmonic_coordinates(dV,dT,b);
W = W(1:size(uV,1),:);
[VV,VT,VF] = tetgen(CV,CF,'Flags','-q100 -Y');
off = mean(VV);
VV=bsxfun(@minus,VV,off);
UU = bsxfun(@plus,VV,[0 0 1]);
vel = zeros(size(UU));
data = [];
M = massmatrix(VV,VT);
% This normalization is very important
M = M/max(M(:));
dt = 1e-2;
t = tsurf(VF,UU, ...
'FaceColor',[0.8 0.5 0.2], ...
'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
%t.Visible = 'off';
l = light('Style','local','Position',[5 0 20]);
l2 = light('Style','infinite','Position',[-1 -1 2]);
hold on;
up = @(UU) bsxfun(@plus,[1 0 0],W*UU(1:numel(b),:));
u = tsurf(uF,up(UU),'EdgeColor','none', ...
'FaceVertexCData',repmat([0.3 0.4 1.0],size(uV,1),1), ...
fphong,'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
shadows = false;
if shadows
e = -1e-4;
QV = [-0.8,-0.5 e; -0.8 0.5 e; 1.5 0.5 e;1.5,-0.5 e];
QQ = [1 2 3 4];
trisurf(QQ,QV(:,1),QV(:,2),QV(:,3),'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);
st = add_shadow(t,l,'Ground',[0 0 -1 0]);
su = add_shadow(u,l,'Ground',[0 0 -1 0]);
st.FaceColor = [0.5 0.5 0.5];
su.FaceColor = [0.5 0.5 0.5];
set(gca,'Visible','off');
set(gcf,'Color','w');
end
hold off;
axis equal;
axis([-0.8 1.8 -0.5 0.5 0 1.5]);
axis manual;
camproj('persp');
T = get(gca,'tightinset');
set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);
while true
dur = 0;
UU0 = UU;
[UU,data] = arap(VV,VT,[],[], ...
'Energy','elements', ...
'V0',UU0, ...
'Data',data, ...
'Dynamic',M*repmat([0 0 -9.8],size(VV,1),1), ...
'MaxIter',100, ...
'TimeStep',dt, ...
'Vm1',UU-vel);
vel = UU-UU0;
C = UU(:,3)<0;
UU(C,3) = -UU(C,3);
vel(C,3) = -vel(C,3);
dur = dur+dt;
t.Vertices = UU;
u.Vertices = up(UU);
if shadows
light_pos = [l.Position strcmp(l.Style,'local')];
ground = [0 0 -1 0];
d = ground * light_pos';
shadow_mat = d*eye(4) - light_pos'*ground;
H = [t.Vertices ones(size(t.Vertices,1),1)]*shadow_mat';
st.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
H = [u.Vertices ones(size(u.Vertices,1),1)]*shadow_mat';
su.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
end
drawnow;
end