Shape Matching
|
|||||||||||||||||
Shape Matching algorithm:This example shows the shape matching simulation method for deformable solids introduced by Müller et al. [MHTG05, BMM17].
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.2. Compute center of massIn 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 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_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 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. 5. Determine goal positionsThe goal positions are computed as: $$ \mathbf g_i = \mathbf R_i \mathbf q_i + \mathbf t_i. $$ 6. 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 user-defined stiffness parameter. A body is rigid if $\alpha=1$.References
|