History: Horta's Mestrado Thesis

In this work several proposals for modeling non-rigid objects were analyzed and the approach suggested by Terzopoulos et al. is implemented in C language with use of the numerical library Meschach. The relaxation method and, instead of the Choleski decomposition which requires the matrix being a positive definite one, the LU factorization are used for solving sparse linear systems. For faciliting the input of modeling parameter values, a simple description language was defined and its parser was implemented with help of lex and yacc. For producing nice pictures of the deforming surfaces, PoVRay and SIPP are used for rendering the generated results.

An Elastically Deformable Model

The deformation dynamics are ruled by the equation of motion in its Lagrangean formulation:


\begin{displaymath} \mu\frac{{\partial}^{2}\vec{r}}{{\partial}{t}^{2}}+ \gamma\f... ...elta\varepsilon(\vec{r})}{\delta\vec{r}}= \vec{f}(\vec{r},t)~. \end{displaymath} (1)

In equation (1), $\mu$ is the mass density and $\gamma$ is the dumping constant at a point $\vec{r}$. The vector $\vec{f}$ denotes the total contribution of external forces at $\vec{r}$ in an instant $t$. The term corresponding to the internal energies accumulated due to elastical deformation $\varepsilon(\vec{r})$ is estimated from the following empirical consideration:


\begin{displaymath} \epsilon({\bf\vec{r}})=\int_{\Omega}\Vert{\bf {G}}-{\bf {G}}... ...\Vert{\bf {B}}-{\bf {B}}^{0}\Vert^{2}_{\beta} da_{1} da_{2}, \end{displaymath} (2)

which takes, for the metric and curvature tensors, a weighted norm of the difference between the deformed state and rest state values. That measure can reasonably estimate the elastical energy of a surface, attaching the amount of energy to the variations in the surface's geometric shape. In other words, that norm is a measure of the energy needed to displace the surface's points, defined over a region $\Omega$, from their rest state.

By applying the weighted norms of equation (2), the following simplified deformation energy:


\begin{displaymath} \epsilon({\bf\vec{r}})=\int_{\Omega}\sum_{i,j}(\eta_{ij}(G_{... ...G^{0}_{ij})^{2}+\xi_{ij}(B_{ij}-B^{0}_{ij})^{2}) da_{1} da_{2} \end{displaymath} (3)

The weights $\eta_{ij}$ and $\xi_{ij}$ are denominated elasticity parameters.

From equation (3) an approximation for the internal force ${\delta \epsilon({\bf\vec{r}})}/{\delta {\bf\vec{r}}}$ is suggested:

\begin{displaymath} {\bf\vec{e}}({\bf\vec r})=\sum_{i,j}-\frac{\partial}{\partia... ...{2}{\bf\vec r}}{\partial {{a}_{i}}\partial {{a}_{j}}}\right), \end{displaymath} (4)

where the variables $\alpha_{ij}$ and $\beta_{ij}$ are defined as:
$\displaystyle \alpha_{ij}(\vec{a},\vec{r})~=~\eta_{ij}(\vec{a})(G_{ij}~-~G_{ij}^{0})$      
$\displaystyle \beta_{ij}(\vec{a},\vec{r})~=~\xi_{ij}(\vec{a})(B_{ij}~-~B_{ij}^{0})~.$     (5)

Since quantities $G_{ij}$ are related to surface stretching, while the values for $B_{ij}$ are related to curvature, the measures of deformation follow from these quantities and the surface's behavior of resistance to external forces will be as much effective as greater are the values assigned to the elasticity parameters.

Discretization

The discretization turns the partial differential equation of motion into a system of coupled ordinary differential equations.

The continuous space $\Omega$ is discretized to a MxN-node mesh, where each node $(m,n)$ represents a discrete point (or a nodal variable) ${\bf\vec r(m,n)}$ in 3D space. To the set of nodal variables ${\bf\vec r(m,n)}$ defined for MN nodes we call function mesh and we denote it ${\bf\vec r [m,n]}$.

Equations (4) and (5) are respectively discretized to

\begin{displaymath} e_{ij}[m,n]=\sum^{2}_{i,j}-D^{-}_{i}({\bf\vec p})[m,n]+D^{(-)}_{ij}({\bf\vec q})[m,n] \end{displaymath} (6)

where
$\displaystyle {\bf {\vec p}}$ $\textstyle =$ $\displaystyle \alpha_{ij}[m,n]D^{+}_{j}({\bf {\vec r}})[m,n]$  
$\displaystyle {\bf {\vec q}}$ $\textstyle =$ $\displaystyle \beta_{ij}[m,n]D^{(+)}_{ij}({\bf {\vec r}})[m,n]$ (7)


\begin{displaymath} \alpha_{ij}[m,n]=\eta_{ij}[m,n](D^{+}_{i}({\bf\vec r})[m,n].D^{+}_{j}({\bf\vec r})[m,n]-G^{0}_{ij}) \end{displaymath} (8)

and
\begin{displaymath} \beta_{ij}[m,n]=\xi_{ij}[m,n]({\bf\vec n}[m,n].D^{(+)}_{ij}({\bf\vec r})[m,n]-B^{0}_{ij}). \end{displaymath} (9)

One can observe that the values for the difference operators are not determined for points laying at the boundaries of domain $\Omega$. Nevertheless, a natural condition of boundary can be simulated by assigning a $zero$ value to any difference operator of equation (7) that refers to points ${\bf\vec r(m,n)}$ not belonging to the set of points MN of the mesh.

If the nodal variables in function meshes ${\bf r[m,n]}$ and ${\bf e[m,n]}$ are grouped, respectively, into column matrices ${\cal R}$ e ${\cal E}$ of dimension MN, then equation (6) can be written in matrix form

\begin{displaymath} {{\cal E}}={\bf {K}}({\bf {\vec r}}){{\cal R}} \end{displaymath} (10)

where ${\bf {K}}$ is known as strength matrix.

The discrete form of the equation of motion can then be expressed by the following coupled system of differential equations:

\begin{displaymath} {\bf {M}}\frac {\partial^{2}{\vec{\cal R}}}{\partial {\bf t}... ...rtial {\bf t}}+ {\bf {K}}({\bf {\vec r}}){{\cal R}}={{\cal F}} \end{displaymath} (11)

where
  • ${\bf {M}}$ is the diagonal matrix formed by the mass density of each element,
  • ${\bf {C}}$, the diagonal matrix formed by the dumping density of each element and
  • ${{\cal F}}$, the column matrix containing the external force applied to each element, calculated from ${\bf\vec f}({\bf\vec r},t)$.

To simulate the dynamics of a non-rigid object, the system of differential equations (11) must be integrated through time. Those equations will be integrated using a step-by-step process, which now converts a system of non-linear differential equations into a sequence of linear systems.

The time interval from t=0 to t=T is subdivided into smaller time intervals of same duration $\Delta$ and the integration process carries out the calculations for the sequence of approximated solutions for instants t, t+$\Delta t$, t+2$\Delta t$, ... , T. Computing ${{\cal E}}$ in $t+\Delta t$ and ${{\cal F}}$ in t, and substituting the discrete-time approximations

\begin{displaymath} \frac {\partial ^{2} {{\cal R}}}{\partial t} = ({{\cal R}}_{t+\Delta t}-2{{\cal R}}+{{\cal R}}_{t-\Delta t})/\Delta t^{2} \end{displaymath} (12)


\begin{displaymath} \frac {\partial {{\cal R}}}{\partial t} = ({{\cal R}}_{t+\Delta t}-{{\cal R}}_{t-\Delta t})/2\Delta t \end{displaymath} (13)

in equation (11) it is obtained
\begin{displaymath} {\bf A}_{t}{{\cal R}}_{t+\Delta t}={{\cal G}}, \end{displaymath} (14)

where
\begin{displaymath} {\bf A}_{t}={\bf K}({\bf\vec r}_{t})+\left(\frac{1}{\Delta t ^{2}}{\bf M}+\frac{1}{2\Delta t}{\bf C}\right) \end{displaymath} (15)

and
\begin{displaymath} {{\cal G}}_{t} = {{\cal F}}_{t}+ \left(\frac{1}{\Delta t ^{2... ...{1}{\Delta t}{\bf M}-\frac{1}{2}{\bf C}\right) {{\cal V}}_{t}. \end{displaymath} (16)

The column matrix of speed ${{\cal V}}_t$ is given by:

\begin{displaymath} {{\cal V}}_{t}=({{\cal R}}_{t}-{{\cal R}}_{t-\Delta t})/\Delta t. \end{displaymath} (17)

An Animation

An animation of a cloth was generated for the Sibgrapi 94 Video Festival.