# Chapter 3: Shape deformation¶

```
import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact
import os
root_folder = os.getcwd()
```

Modern mesh-based shape deformation methods satisfy user deformation
constraints at handles (selected vertices or regions on the mesh) and propagate
these handle deformations to the rest of the shape *smoothly* and *without removing
or distorting details*. Libigl provides implementations of a variety of
state-of-the-art deformation techniques, ranging from quadratic mesh-based
energy minimizers, to skinning methods, to non-linear elasticity-inspired
techniques.

## Biharmonic deformation¶

The period of research between 2000 and 2010 produced a collection of techniques that cast the problem of handle-based shape deformation as a quadratic energy minimization problem or equivalently the solution to a linear partial differential equation.

There are many flavors of these techniques, but a prototypical subset are those that consider solutions to the bi-Laplace equation, that is a biharmonic function (Botsch, 2004). This fourth-order PDE provides sufficient flexibility in boundary conditions to ensure C^1 continuity at handle constraints in the limit under refinement (Jacobson, 2010).

### Biharmonic surfaces¶

Let us first begin our discussion of biharmonic *deformation*, by considering
biharmonic *surfaces*. We will casually define biharmonic surfaces as surface
whose *position functions* are biharmonic with respect to some initial
parameterization:

\Delta^2 \mathbf{x}' = 0

and subject to some handle constraints, conceptualized as “boundary conditions”:

\mathbf{x}'_{b} = \mathbf{x}_{bc}.

where \mathbf{x}' is the unknown 3D position of a point on the surface. So we are asking that the bi-Laplacian of each of spatial coordinate function to be zero.

In libigl, one can solve a biharmonic problem with `harmonic`

and setting k=2 (*bi*-harmonic).

This produces a smooth surface that interpolates the handle constraints, but all
original details on the surface will be *smoothed away*. Most obviously, if the
original surface is not already biharmonic, then giving all handles the
identity deformation (keeping them at their rest positions) will **not**
reproduce the original surface. Rather, the result will be the biharmonic
surface that does interpolate those handle positions.

Thus, we may conclude that this is not an intuitive technique for shape deformation.

### Biharmonic deformation fields¶

Now we know that one useful property for a deformation technique is “rest pose reproduction”: applying no deformation to the handles should apply no deformation to the shape.

To guarantee this by construction we can work with *deformation fields* (ie.
displacements)
\mathbf{d} rather
than directly with positions \mathbf{x}. Then the deformed positions can be
recovered as

\mathbf{x}' = \mathbf{x}+\mathbf{d}.

A smooth deformation field \mathbf{d} which interpolates the deformation
fields of the handle constraints will impose a smooth deformed shape
\mathbf{x}'. Naturally, we consider *biharmonic deformation fields*:

\Delta^2 \mathbf{d} = 0

subject to the same handle constraints, but rewritten in terms of their implied deformation field at the boundary (handles).

\mathbf{d}_b = \mathbf{x}_{bc} - \mathbf{x}_b.

Again we can use `harmonic`

with k=2, but this time solve for the
deformation field and then recover the deformed positions:

```
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "decimated-max.obj"))
v[:,[0, 2]] = v[:,[2, 0]] # Swap X and Z axes
u = v.copy()
s = igl.read_dmat(os.path.join(root_folder, "data", "decimated-max-selection.dmat"))
b = np.array([[t[0] for t in [(i, s[i]) for i in range(0, v.shape[0])] if t[1] >= 0]]).T
## Boundary conditions directly on deformed positions
u_bc = np.zeros((b.shape[0], v.shape[1]))
v_bc = np.zeros((b.shape[0], v.shape[1]))
for bi in range(b.shape[0]):
v_bc[bi] = v[b[bi]]
if s[b[bi]] == 0: # Don't move handle 0
u_bc[bi] = v[b[bi]]
elif s[b[bi]] == 1: # Move handle 1 down
u_bc[bi] = v[b[bi]] + np.array([[0, -50, 0]])
else: # Move other handles forward
u_bc[bi] = v[b[bi]] + np.array([[-25, 0, 0]])
p = subplot(v, f, s, shading={"wireframe": False, "colormap": "tab10"}, s=[1, 4, 0])
for i in range(3):
u_bc_anim = v_bc + i*0.6 * (u_bc - v_bc)
d_bc = u_bc_anim - v_bc
d = igl.harmonic(v, f, b, d_bc, 2)
u = v + d
subplot(u, f, s, shading={"wireframe": False, "colormap": "tab10"}, s=[1, 4, i+1], data=p)
p
# @interact(deformation_field=True, step=(0.0, 2.0))
# def update(deformation_field, step=0.0):
# # Determine boundary conditions
# u_bc_anim = v_bc + step * (u_bc - v_bc)
# if deformation_field:
# d_bc = u_bc_anim - v_bc
# d = igl.harmonic(v, f, b, d_bc, 2)
# u = v + d
# else:
# u = igl.harmonic(v, f, b, u_bc_anim, 2)
# p.update_object(vertices=u)
```