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)