Shape Matching algorithm:

This example shows the shape matching simulation method for deformable solids introduced by Müller et al. [MHTG05, BMM17].

1. time integration to predict particle positions
2. compute center of mass of current particle configuration
3. compute moment matrix
4. extract rotation
5. determine goal positions
6. update positions and velocities

Note that shape matching is a meshless method which is geometrically motivated. So it can only generate visually plausible results. The general idea is to match the initial configuration of the particles to the current deformed configuration to get goal positions. Then the particles are pulled towards these goal positions to simulate elasticity.

1. Time integration

In the first step the particle positions are advected using a symplectic Euler method to obtain predicted positions $\mathbf x^*$: $$\begin{align*} \mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a_\text{ext}(t) \\ \mathbf x^* &= \mathbf x(t) + \Delta t \mathbf v^*, \end{align*}$$ where $\mathbf a_\text{ext}$ are accelerations due to external forces.

2. Compute center of mass

In order to determine the goal positions for all particles we search for a rotation matrix $\mathbf R$ and translation vectors $\mathbf t, \mathbf t_0$ which minimize the following equation: $$\begin{equation*} \sum_i m_i (\mathbf R(\mathbf X_i - \mathbf t_0) + \mathbf t - \mathbf x^*_i)^2, \end{equation*} $$ where $m$ is the mass, $\mathbf X$ defines the initial configuration and $\mathbf x^*$ the predicted deformed configuration.

The optimal translation vectors are the centers of mass of the initial configuration and the deformed configuration: $$ \begin{eqnarray*} \mathbf t_0 &=& \frac{\sum_i m_i \mathbf X_i}{\sum_i m_i} \\ \mathbf t &=& \frac{\sum_i m_i \mathbf x^*_i}{\sum_i m_i} \end{eqnarray*} $$

3. Compute moment matrix

To compute the optimal rotation we move the both configurations so that their centers of mass lie in the origin: $$\begin{align*} \mathbf q_i &= \mathbf X_i- \mathbf t^0 \\ \mathbf p_i &= \mathbf x^*_i- \mathbf t \end{align*} $$ Moreover, to simplify the problem we search for the optimal linear transformation $\mathbf A$ instead of the optimal rotation $\mathbf R$. This means we have to minimize: $$ \begin{equation*} \sum_i m_i (\mathbf A \mathbf q_i -\mathbf p_i)^2 \end{equation*} $$ Setting the derivatives with respect to the coefficients of A to zero gives us the optimal transformation: $$\begin{equation*} \mathbf A = \left (\sum_i m_i \mathbf p_i \mathbf q_i^T \right ) \left (\sum_i m_i \mathbf q_i \mathbf q_i^T \right )^{-1} = \mathbf A_{pq} \mathbf A_{qq}. \end{equation*} $$

4. Extract rotation

The term $\mathbf A_{qq}$ is symmetric and therefore contains no rotation. Thus, the optimal rotation is the rotational part of $\mathbf A_{pq}$ which can be determined by a polar decomposition: $$\mathbf A_{pq} = \mathbf R \mathbf S,$$ where $\mathbf R$ is an orthogonal matrix and $\mathbf S$ a symmetric matrix. If the configuration is not inverted, $\mathbf R$ is a rotation matrix. Otherwise $\mathbf R$ contains a reflection.

5. Determine goal positions

The goal positions are computed as: $$ \mathbf g_i = \mathbf R_i \mathbf q_i + \mathbf t_i. $$

6. Update positions and velocities

Finally, we update the positions and velocity by the following modified time integration scheme: $$ \begin{eqnarray*} \mathbf v_i(t+h) &=& \mathbf v_i(t) + \alpha \frac{\mathbf g_i(t) - \mathbf x_i(t)}{h} + \Delta t \frac{\mathbf f_{ext}(t)}{m} \\ \mathbf x_i(t+h) &=& \mathbf x_i(t) + \Delta t \mathbf v_i(t+h), \end{eqnarray*} $$ where $\alpha \in [0,1]$ is a user-defined stiffness parameter. A body is rigid if $\alpha=1$.

References

• [MHTG05] Matthias Müller, Bruno Heidelberger, Matthias Teschner, Markus Gross. Meshless Deformations Based on Shape Matching. ACM Transactions on Graphics 24, 3, 2005
• [BMM17] Jan Bender, Matthias Müller, Miles Macklin. A Survey on Position Based Dynamics, 2017. Eurographics Tutorials, 2017