RegionBased Shape Matching


RegionBased Shape Matching algorithm:This example shows the regionbased shape matching simulation method for deformable solids [BMM17].
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. However, the original shape matching algorithm allows only for small deviations from the initial shape. The regionbased method performs shape matching in overlapping regions and averages the resulting goal positions for each particle in order to enable the simulation of large deformations. 1. Time integrationIn 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. Moreover, we initialize the goal positions as $\mathbf g_i = \mathbf 0$.2. Loop over all regionsShape matching is performed for overlapping regions of particles. The resulting goal positions of a particle (one for each region that contains the particle) are averaged. In contrast to performing shape matching for all particles in one single large region, this approach enables the simulation of large deformations. Note that in this example we defined one region per particle. 3. Compute center of massIn order to determine the goal positions for all particles $i$ 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 for region $i$ are the centers of mass of the initial configuration and the deformed configuration: $$ \begin{eqnarray*} \mathbf t_i^0 &=& \frac{1}{\tilde M_i} \sum_{j \in \Re_i} \tilde m_j \mathbf X_j \\ \mathbf t_i &=& \frac{1}{\tilde M_i} \sum_{j \in \Re_i} \tilde m_j \mathbf x^*_j, \end{eqnarray*} $$ where $\tilde M_i = \sum_{j \in \Re_i} \tilde m_j$. Here we use modified particle masses $\tilde m_i = m_i /  \Re_i$ to consider the fact that a particle $i$ is part of $ \Re_i$ regions, where $\Re_i$ is the set of regions containing particle $i$. 4. Compute moment matrixTo 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_{j \in \Re_i} \tilde m_j \mathbf p_j \mathbf q_j^T \right ) \left (\sum_{j \in \Re_i} \tilde m_j \mathbf q_j \mathbf q_j^T \right )^{1} = \mathbf A_{pq} \mathbf A_{qq}. \end{equation*} $$ 5. Extract rotationThe 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. 6. Determine goal positionsIn this step we compute the goal position for each particle of the current region $r$. Since multiple regions can contain the same particle $i$, we sum up all contributions and divide by the number of regions which contain the particle: $$ \mathbf g_i := \mathbf g_i + \frac{1}{\Re_i} (\mathbf R \mathbf q_i + \mathbf t). $$ 7. Update positions and velocitiesFinally, 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 userdefined stiffness parameter. A body is rigid if $\alpha=1$.References
