Implementing Baraff & Witkin's Cloth Simulation

David Pritchard

Date: May 2003, with minor updates April 2012

1 Introduction

In April 2002, I implemented a cloth simulator based on Baraff & Witkin's paper [BW98] for a graduate animation course taught by Michiel van de Panne at the University of British Columbia. I later released that simulator as an open-source project called Freecloth. Over the following year, I improved the performance of the simulator and fixed some major bugs.

The initial version of Freecloth had some major bugs, and did not implement all parts of the original paper. Newer versions have eliminated the bugs in the simulator, but do not yet include all parts of the original paper.

I will assume that you are familiar with Baraff & Witkin's paper, and I will not try to re-explain concepts which are clearly presented there. For my own work, I relied heavily upon a paper by Dean Macri [Mac00], who attempted to interpret and implement Baraff & Witkin's paper. I will make reference to Macri's notation and some of his formulas, but his paper is not necessary reading. I rederived all of the formulas found in Macri's paper, and verified my derivation against Macri's. Initially, I found a few minor typographical errors in his paper, but my formulas were generally the same.

As my implementation progressed, I found more major errors in Macri's paper. The largest error occurred when he calculated the derivatives of the bend condition. I give a corrected version of this equation in (54) and (57).

In this paper, I will attempt to explain obstacles that I encountered while implementing Baraff & Witkin's paper.

2 Notation

In this paper, I will use the symbol ${\bf I}$ to denote an identity matrix. The symbol ${\bf I}_s$ refers to a vector containing the $s$ th column of the identity matrix. To denote the skew-symmetric matrix form of a vector cross product, I will use the notation

$$S({\bf v}) = \begin{pmatrix} 0 & -{\bf v}_z & {\bf v}_y \\ {\bf v}_z & 0 & -{\bf v}_x \\ -{\bf v}_y & {\bf v}_x & 0 \end{pmatrix}$$

I will also define the vector $S_s({\bf v})$ as the transpose of the $s$ th row of $S({\bf v})$ . $S_s$ is a column vector, but represents a row in the original $S$ matrix. I will use a hat (e.g. ${\bf {\hat w}}_u$ ) to refer to a normalised vector of unit length, and omit the hat (e.g. ${\bf w}_u$ ) for an unnormalised vector. This is the opposite convention to that used in my original project report.

I will retain Baraff & Witkin's notation and use ${\bf x}$ to refer to the positions of points, not ${\bf p}$ as in Macri's paper. Furthermore, I will refer to the points of a triangle as ${\bf x}_0$ , ${\bf x}_1$ and ${\bf x}_2$ , and not use the $i$ , $j$ and $k$ subscripts used by the other two papers. For the points of an edge, I will likewise use subscripts of 0, 1, 2, and 3 rather than the $h$ , $i$ , $j$ and $k$ subscripts used by Macri.

When referring to points, I will use ${\bf x}_m$ and ${\bf x}_n$ to refer to two possibly distinct points. I will use the subscript $m_s$ to refer to the to the $s$ th component of ${\bf x}_m$ , and likewise the subscript $n_t$ for the $t$ th component of ${\bf x}_n$ . I will also sometimes use subsubscripts of $x$ , $y$ and $z$ to refer to the components.

I like to think of the derivatives in this paper in terms of their components. Both Baraff & Witkin and Macri differentiate with respect to a vector in their papers, giving quantities like

$$\frac{\partial C}{\partial {\bf x}_m} = \begin{pmatrix}\frac{\partial C}{\partial {\bf x}_{m_x}} \\ \fr... ...rtial {\bf x}_{m_y}} \\ \frac{\partial C}{\partial {\bf x}_{m_z}} \end{pmatrix}$$

or

$$\frac{\partial {\bf v}}{\partial {\bf x}_m} = \begin{pmatrix}\frac{\partial {\bf v}_x}{\partial {\bf x}_{m_x}... ...\bf x}_{m_z}} & \frac{\partial {\bf v}_z}{\partial {\bf x}_{m_z}} \end{pmatrix}$$

This type of formula is common in some areas of graphics; for example, the Jacobian matrix has this form. However, what happens when we take the second derivative of ${\bf v}$ ? It will yield a third-order tensor, a $3\times3\times3$ linear algebraic quantity. However, this is unnecessarily complicated, especially for people (like me) who never covered tensors in their undergraduate education. In this paper, I will try to avoid differentiation with respect to a vector. Instead, I will do most of the work on a component-by-component basis, using only differentiation with respect to a scalar. Using this approach, I will only ever need to use vectors and scalars, not matrices or higher-order tensors. This makes operations such as cross-products and dot-products straightforward.

Some of Baraff & Witkin's formulas involve taking a second derivative with respect to two vectors. For these, just remember the layout of the resulting matrix. For example,

$$\frac{\partial^2 C}{\partial {\bf x}_m \partial {\bf x}_n} = \begin{pmatrix}\frac{\partial C}{\partial {\bf x}_{m_x}} {{\bf ... ...n_z}} & \frac{\partial C}{\partial {\bf x}_{m_z}} {{\bf x}_{n_z}} \end{pmatrix}$$

3 Condition Functions

Baraff & Witkin used condition functions as an indirect way of defining energy and forces. Their condition functions are almost the same as energy formulas, since $E=\frac{k}{2}{\bf C}^T{\bf C}$ , where ${\bf C}$ is the condition function. For the scalar condition functions (shear and bend), $E=\frac{k}{2}C^2$ , and the condition function is hence basically just the squart-root of the energy, except that it can be positive or negative.

The condition functions were defined explicitly in the paper. Energy and forces can also be easily calculated, if the derivatives of the condition functions are known. However, the paper didn't provide explicit equations for the derivatives or second derivatives.

Macri made a good attempt at deriving the equations for the derivatives and second derivatives. I found two major shortcomings in his approach, however. First, he made reference to tensors and other mathematical constructs that aren't really necessary for the derivation. Second, I found an error in his derivation for the bend condition.

3.1 Per-Triangle Quantities

To begin, I will define a few quantities that will help with explanation of the stretch and shear conditions. These quantities are defined for a triangle using points ${\bf x}_0$ , ${\bf x}_1$ , and ${\bf x}_2$ . The triangle's area in $u$ /$v$ space is given by

$$a = \frac{1}{2} \begin{Vmatrix}\begin{pmatrix}\Delta u_1 \\ \Delta ... ...\times \begin{pmatrix}\Delta u_2 \\ \Delta v_2 \\ 0 \end{pmatrix} \end{Vmatrix}$$ (1)

The ${\bf {\hat w}}_u$ and ${\bf {\hat w}}_v$ vectors represent the directions of the $u$ and $v$ axes in world space, as seen within a single triangle. They are defined by

$${\bf {\hat w}}_u = \frac{ {\bf w}_u}{\begin{Vmatrix}{\bf w}_u\end{Vmatrix}} {\bf {\hat w}}_v = \frac{ {\bf w}_v}{\begin{Vmatrix}{\bf w}_v\end{Vmatrix}}$$ (2)

$${\bf w}_u = \frac{ ({\bf x}_1 - {\bf x}_0) \Delta v_2 - ({\bf x}_2 - {\bf x}_0) \Delta v_1} {\Delta u_1 \Delta v_2 - \Delta u_2 \Delta v_1}$$ (3)

$${\bf w}_v = \frac{ -({\bf x}_1 - {\bf x}_0) \Delta u_2 + ({\bf x}_2 - {\bf x}_0) \Delta u_1} {\Delta u_1 \Delta v_2 - \Delta u_2 \Delta v_1}$$ (4)

The first derivatives of the unnormalised vectors ${\bf w}_u$ and ${\bf w}_v$ are given by

$$\frac{\partial {\bf w}_u}{\partial {\bf x}_m} = \frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{m_x}} {\bf I}$$ (5)

$$\frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{0_x}} = \frac{\Delta v_1 - \Delta v_2}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (6)

$$\frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{1_x}} = \frac{\Delta v_2}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (7)

$$\frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{2_x}} = \frac{-\Delta v_1}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (8)

$$\frac{\partial {\bf w}_v}{\partial {\bf x}_m} = \frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{m_x}} {\bf I}$$ (9)

$$\frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{0_x}} = \frac{\Delta u_2 - \Delta u_1}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (10)

$$\frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{1_x}} = \frac{-\Delta u_2}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (11)

$$\frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{2_x}} = \frac{\Delta u_1}{\Delta u_1\Delta v_2 - \Delta u_2\Delta v_1}$$ (12)

The second derivatives of ${\bf w}_u$ and ${\bf w}_v$ are all zero.

3.2 Stretch Condition

Baraff & Witkin define the stretch condition function as a vector ${\bf C}$ . Two special parameters are defined for stretch, $b_u$ and $b_v$ . These represent the desired stretchiness of the cloth in the $u$ and $v$ directions. A $b_u$ value below one will result in the cloth trying to compress in the $u$ direction, while a $b_u$ value above one will result in the cloth trying to stretch in the $u$ direction.

$${\bf C}= \begin{pmatrix}C_u\\ C_v\end{pmatrix} = \alpha \begin{pmatrix}\begin{Vmatrix}{\bf w}_u\end{Vmatrix} - b_u \\ \begin{Vmatrix}{\bf w}_v\end{Vmatrix} - b_v \end{pmatrix}$$ (13)

Baraff & Witkin defined $\alpha=a$ , the area of the triangle. We define $\alpha=a^\frac{3}{4}$ , for reasons described later.

The first and second derivatives of the stretch condition function are

$$\frac{\partial C_u}{\partial {\bf x}_m} = \alpha \frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{m_x}} {\bf {\hat w}}_u$$ (14)

$$\frac{\partial^2 C_u}{\partial {\bf x}_m \partial {\bf x}_n} = \frac{\alpha}{\begin{Vmatrix}{\bf w}_u\end{Vmatrix}} \frac{\par... ...rtial {\bf x}_{n_x}} \left( {\bf I}- {\bf {\hat w}}_u{\bf {\hat w}}_u^T \right)$$ (15)

$$\frac{\partial C_v}{\partial {\bf x}_m} = \alpha \frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{m_x}} {\bf {\hat w}}_v$$ (16)

$$\frac{\partial^2 C_v}{\partial {\bf x}_m \partial {\bf x}_n} = \frac{\alpha}{\begin{Vmatrix}{\bf w}_v\end{Vmatrix}} \frac{\par... ...rtial {\bf x}_{n_x}} \left( {\bf I}- {\bf {\hat w}}_v{\bf {\hat w}}_v^T \right)$$ (17)

These second derivatives exhibit some symmetry, since $\frac{\partial^2 C_u}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \frac{\partial^2 C_u}{\partial {\bf x}_{n_t} \partial {\bf x}_{m_s}}$ . Consequently, the $3\times3$ matrix $\frac{\partial^2 C_u}{\partial {\bf x}_m \partial {\bf x}_n} = \frac{\partial^2 C_u}{\partial {\bf x}_n \partial {\bf x}_m} ^T$ . To avoid any tensors when applying Baraff & Witkin's equation (8) to the shear equations, remember that

$$\frac{\partial^2 {\bf C}({\bf x})}{\partial {\bf x}_m \partial {\bf x}_n} {\bf C}({\bf x}) = \frac{\partial^2 C_u}{\partial {\bf x}_m \partial {\bf x}_n} C_u+ \frac{\partial^2 C_v}{\partial {\bf x}_m \partial {\bf x}_n} C_v$$ (18)

3.3 Shear Condition

Baraff & Witkin define the shear condition function as a scalar. This scalar is essentially the dot product of the $u$ axis with the $v$ axis in world space. If no shear is occurring, the axes are perpendicular and the condition function is zero. If shear is occurring, the condition function is equivalent to the cosine of the angle between them, weighted by the triangle's area in $u$ / $v$ space.

$$C = \alpha {\bf w}_u\cdot {\bf w}_v$$ (19)

Baraff & Witkin's condition assumes that ${\bf w}_u$ and ${\bf w}_v$ are approximately unit vectors. If $b_u$ and $b_v$ are both equal to one, then this approximation is a fair one. If, however, $b_u$ or $b_v$ are not equal to one, the approximation may give rise to some strange behaviour. I have not had a chance to investigate the effects of varying $b_u$ or $b_v$ .

Since we know the derivatives of ${\bf w}_u$ and ${\bf w}_v$ , taking the derivative of this condition function is trivial.

$$\frac{\partial C}{\partial {\bf x}_{m_s}} = \alpha \left( \frac{\partial {\bf w}_u}{\partial {\bf x}_{m_s}}... ... w}_v+ {\bf w}_u\cdot \frac{\partial {\bf w}_v}{\partial {\bf x}_{m_s}} \right)$$ (20)

$$= \alpha \left( \frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{m... ...s}+ {\bf w}_{u_s} \frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{m_x}} \right)$$ (21)

$$\frac{\partial^2 C}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$= \,\alpha \biggl( \frac{\partial^2 {\bf w}_u}{\partial {\bf x}_{m_... ...artial {\bf x}_{m_s}} \cdot \frac{\partial {\bf w}_v}{\partial {\bf x}_{n_t}} +$$

$$\frac{\partial {\bf w}_u}{\partial {\bf x}_{n_t}} \cdot \frac{\pa... ...ac{\partial^2 {\bf w}_u}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} \biggr)$$ (22)

$$=$$ (23)

These second derivatives exhibit even more symmetry than the stretch condition. Like the stretch condition, $\frac{\partial^2 C}{\partial {\bf x}_m \partial {\bf x}_n} = \frac{\partial^2 C}{\partial {\bf x}_n \partial {\bf x}_m} ^T$ but beyond that, the $\frac{\partial^2 C}{\partial {\bf x}_m \partial {\bf x}_n}$ is actually just an identity matrix multiplied by a scalar.

3.4 Bend Condition

Baraff & Witkin define the bend condition in terms of two adjoining triangles. I have labelled the unit normals of these triangles as ${\bf {\hat n}}^A$ and ${\bf {\hat n}}^B$ , rather than Baraff & Witkin's ${\bf n}_1$ or ${\bf n}_2$ . The unit vector parallel to their common edge is labelled ${\bf {\hat e}}$ . Quantities associated with the normals or the common edge receive a superscript $A$ , $B$ or $e$ .

$${\bf {\hat n}}^A = \frac{{\bf n}^A}{\begin{Vmatrix}{\bf n}^A\end{Vmatrix}} {\bf {\hat n}}^B = \frac{{\bf n}^B}{\begin{Vmatrix}{\bf n}^B\end{Vmatrix}} {\bf {\hat e}} = \frac{{\bf e}}{\begin{Vmatrix}{\bf e}\end{Vmatrix}}$$ (24)

$${\bf n}^A = ({\bf x}_2 - {\bf x}_0) \times ({\bf x}_1 - {\bf x}_0)$$ (25)

$${\bf n}^B = ({\bf x}_1 - {\bf x}_3) \times ({\bf x}_2 - {\bf x}_3)$$ (26)

$${\bf e} = {\bf x}_1 - {\bf x}_2$$ (27)

The scalar bend condition $C$ is defined quite simply as the angle between the two edges, $\theta$ .

$$\cos \theta = {\bf {\hat n}}^A \cdot {\bf {\hat n}}^B$$ (28)

$$\sin \theta = ({\bf {\hat n}}^A \times {\bf {\hat n}}^B) \cdot {\bf {\hat e}}$$ (29)

$$C = \theta = \arctan \frac{\sin \theta}{\cos \theta}$$ (30)

Macri defined $C$ using the $\arccos$ function, which is incorrect and will only yield positive values for $\theta$ . My definition of $C$ will yield both positive and negative values if implemented using the atan2 function.

In order to find the derivatives of $C$ , we need the derivatives of the quantities seen so far. The first and second derivatives of ${\bf n}^A$ , ${\bf n}^B$ and ${\bf e}$ are shown below, using auxiliary variables ${\bf q}_m$ . The derivatives of ${\bf q}_m$ are straightforward to derive.

$$\frac{\partial {\bf n}^A}{\partial {\bf x}_{m_s}} = S_s({\bf q}^A_m)$$ (31)

$$\frac{\partial {\bf n}^B}{\partial {\bf x}_{m_s}} = S_s({\bf q}^B_m)$$ (32)

$$\frac{\partial {\bf e}}{\partial {\bf x}_{m_s}} = q^e_m {\bf I}_s$$ (33)

$${\bf q}^A = \{ {\bf x}_2 - {\bf x}_1, {\bf x}_0 - {\bf x}_2, {\bf x}_1 - {\bf x}_0, {\bf0} \}$$ (34)

$${\bf q}^B = \{ {\bf0}, {\bf x}_2 - {\bf x}_3, {\bf x}_3 - {\bf x}_1, {\bf x}_1 - {\bf x}_2 \}$$ (35)

$$q^e = \{ 0, 1, -1, 0 \}$$ (36)

$$\frac{\partial^2 {\bf n}^A}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = S_s\left( \frac{\partial {\bf q}^A_m}{\partial {\bf x}_{n_t}} \right)$$ (37)

$$\frac{\partial^2 {\bf n}^B}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = S_s\left( \frac{\partial {\bf q}^B_m}{\partial {\bf x}_{n_t}} \right)$$ (38)

$$\frac{\partial^2 {\bf e}}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = {\bf0}$$ (39)

Using Baraff & Witkin's assumption that the normal and edge vectors have constant magnitude,

$$\frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}} = \frac{1}{\begin{Vmatrix}{\bf n}^A\end{Vmatrix}} \frac{\partial {\bf n}^A}{\partial {\bf x}_{m_s}}$$ (40)

$$\frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{m_s}} = \frac{1}{\begin{Vmatrix}{\bf n}^B\end{Vmatrix}} \frac{\partial {\bf n}^B}{\partial {\bf x}_{m_s}}$$ (41)

$$\frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{m_s}} = \frac{1}{\begin{Vmatrix}{\bf e}\end{Vmatrix}} \frac{\partial {\bf e}}{\partial {\bf x}_{m_s}}$$ (42)

$$\frac{\partial^2 {\bf {\hat n}}^A}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \frac{1}{\begin{Vmatrix}{\bf n}^A\end{Vmatrix}} \frac{\partial^2 {\bf n}^A}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$ (43)

$$\frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{m_s}} {{\bf x}_{n_t}} = \frac{1}{\begin{Vmatrix}{\bf n}^B\end{Vmatrix}} \frac{\partial^2 {\bf n}^B}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$ (44)

$$\frac{\partial^2 {\bf {\hat e}}}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = {\bf0}$$ (45)

Given the derivatives of the normal and edge vectors, the first and second derivatives of $\cos \theta$ and $\sin \theta$ are trivial to calculate.

$$\frac{\partial \cos \theta }{\partial {\bf x}_{m_s}} = \frac{\partial {\bf {\hat n}}^A \cdot {\bf {\hat n}}^B }{\partial {\bf x}_{m_s}}$$ (46)

$$= \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}} \cdot {\... ...{\bf {\hat n}}^A \cdot \frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{m_s}}$$ (47)

$$\frac{\partial \sin \theta }{\partial {\bf x}_{m_s}} = \frac{\partial ({\bf {\hat n}}^A \times {\bf {\hat n}}^B) \cdot {\bf {\hat e}}}{\partial {\bf x}_{m_s}}$$ (48)

$$= \left( \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}} \... ...\partial {\bf {\hat n}}^B}{\partial {\bf x}_{m_s}} \right) \cdot {\bf {\hat e}}$$

$$+ ({\bf {\hat n}}^A \times {\bf {\hat n}}^B) \cdot \frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{m_s}}$$ (49)

$$\frac{\partial^2 \cos \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \frac{\partial^2 {\bf {\hat n}}^A}{\partial {\bf x}_{m_s} \partia... ...l {\bf x}_{n_t}} \cdot \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}}$$

$$+ \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{n_t}} \cdot ... ...rac{\partial^2 {\bf {\hat n}}^B}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$ (50)

$$\frac{\partial^2 \sin \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$= \bigg( \frac{\partial^2 {\bf {\hat n}}^A}{\partial {\bf x}_{m_s} ... ... {\bf x}_{m_s}} \times \frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{n_t}}$$

$$+ \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{n_t}} \times... ...}^B}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} \bigg) \cdot {\bf {\hat e}}$$

$$+ \left( \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}}... ... x}_{m_s}} \right) \cdot \frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{n_t}}$$

$$+ \left( \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{n_t}}... ... x}_{n_t}} \right) \cdot \frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{m_s}}$$

$$+ ({\bf {\hat n}}^A \times {\bf {\hat n}}^B) \frac{\partial^2 {\bf {\hat e}}}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$ (51)

Finally, we can differentiate equation (30) to obtain the derivatives of $C$ .

$$\frac{\partial C}{\partial {\bf x}_{m_s}} = \frac{\partial \theta}{\partial {\bf x}_{m_s}}$$ (52)

$$= \frac{1}{\sin^2 \theta + \cos^2 \theta}$$

$$\times \left(\cos \theta \frac{\partial \sin \theta}{\partial {\b... ...m_s}} - \sin \theta \frac{\partial \cos \theta}{\partial {\bf x}_{m_s}} \right)$$ (53)

$$= \cos \theta \frac{\partial \sin \theta}{\partial {\bf x}_{m_s}} - \sin \theta \frac{\partial \cos \theta}{\partial {\bf x}_{m_s}}$$ (54)

Note that the second derivative formula below starts from (53), including the denominator term.

$$\frac{\partial^2 C}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \frac{\partial^2 \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$ (55)

$$= \frac{1}{(\sin^2 \theta + \cos^2 \theta)^2}$$

$$\times \Bigg[ (\sin^2 \theta + \cos^2 \theta) \bigg(\cos \theta \frac{\partial^2 \sin \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$- \sin \theta \frac{\partial^2 \cos \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} \bigg)$$

$$+ (\sin^2 \theta - \cos^2 \theta) \left( \frac{\partial \sin \the... ...tial {\bf x}_{m_s}} \frac{\partial \sin \theta}{\partial {\bf x}_{n_t}} \right)$$

$$+ 2 \sin\theta \cos\theta \left( \frac{\partial \cos \theta}{\par... ...bf x}_{m_s}} \frac{\partial \sin \theta}{\partial {\bf x}_{n_t}} \right) \Bigg]$$ (56)

$$= \cos \theta \frac{\partial^2 \sin \theta}{\partial {\bf x}_{m_s... ...ta \frac{\partial^2 \cos \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$+ (\sin^2 \theta - \cos^2 \theta) \left( \frac{\partial \sin \the... ...tial {\bf x}_{m_s}} \frac{\partial \sin \theta}{\partial {\bf x}_{n_t}} \right)$$

$$+ 2 \sin\theta \cos\theta \left( \frac{\partial \cos \theta}{\par... ...tial {\bf x}_{m_s}} \frac{\partial \sin \theta}{\partial {\bf x}_{n_t}} \right)$$ (57)

In a manner similar to the stretch condition, $\frac{\partial^2 C}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \frac{\partial^2 C}{\partial {\bf x}_{n_t} \partial {\bf x}_{m_s}}$ and the $3\times3$ matrix is symmetric. The second derivatives of $\cos \theta$ and $\sin \theta$ are likewise symmetric.

4 Verification

Macri didn't describe any way of verifying the correctness of his derivations or implementation. Once I had initially implemented the simulator, I found bugs but had no way of determining where an error was occurring. To address this, I used numeric derivatives.

For example, suppose that we need to verify the calculation of $\frac{\partial f(x)}{\partial x}$ . Our program calculates the derivative using an analytic formula, which we'll call $f'_a(x)$ . The derivative can also be calculated numerically:

$$f'_n(x) = \frac{f(x+h)-f(x)}{h}$$

If $f'_a(x) = f'_n(x)$ , then the implementation and analytic formula give the correct results.

Verification is a little bit trickier for the bend condition. Here, we can't just use a normal derivative, since we have to differentiate under the assumption that the vectors have constant magnitude. Suppose that we want to verify the derivative $\frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}}$ . Let ${\bf n}^A({\bf x}_{m_s})$ represent the normal with point $m$ at position ${\bf x}_{m_s}$ , and ${\bf n}^A({\bf x}_{m_s}+ h{\bf I}_s)$ represent the normal once point $m$ 's $s$ th component is shifted by $h$ . The unit normals are then given by

$${\bf {\hat n}}^A({\bf x}_{m_s}) = \frac{{\bf n}^A({\bf x}_{m_s})}{\begin{Vmatrix}{\bf n}^A({\bf x}_{m_s})\end{Vmatrix}}$$

$${\bf {\hat n}}^A({\bf x}_{m_s}+ h{\bf I}_s) = \frac{{\bf n}^A({\bf x}_{m_s}+ h{\bf I}_s)}{\begin{Vmatrix}{\bf n}^A({\bf x}_{m_s})\end{Vmatrix}}$$

And the numeric derivative is given by

$$\left({{\bf {\hat n}}^{A}}_n\right)' = \frac{{\bf {\hat n}}^A({\bf x}_{m_s}+ h{\bf I}_s) - {\bf {\hat n}}^A({\bf x}_{m_s})}{h}$$

The main trick is to normalise the shifted edge and normal vectors using the magnitude of the unshifted vector.

To be honest, I don't know why numeric derivatives can't be used for everything. I suspect that they may be less robust than analytic derivatives, and there may be a performance tradeoff involved, but I haven't confirmed these suspicions.

5 Discussion

5.1 Macri's Error

Equations (54) and (57) are very different from the equations obtained by Macri. Macri differentiated equations (28) and (29) and then rearranged to obtain two different formulas for $\frac{\partial C}{\partial {\bf x}_m}$ :

$$\frac{\partial C}{\partial {\bf x}_m} = -\frac{1}{\sin \theta} \frac{\partial \cos \theta}{\partial {\bf x}_m}$$
and

$$\frac{\partial C}{\partial {\bf x}_m} = \frac{1}{\cos \theta} \frac{\partial \sin \theta}{\partial {\bf x}_m}$$

Likewise, Macri obtained two different formulas for $\frac{\partial^2 C}{\partial {\bf x}_m \partial {\bf x}_n}$ . Macri suggested choosing the first set of equations if $\lvert\sin \theta\rvert > \lvert\cos \theta\rvert$ . The rationale for this choice is presumably to avoid numerical problems and division by zero when $\sin \theta=0$ . In order for the transition from one regime to the other to be smooth, they should be equal at some changeover point. However, it is not evident that there exists any changeover point where the two formulas are equal. Additionally, neither of these equations matches the numeric derivative of $C$ . My equations, however, give a result identical to the numeric derivative.

Macri's approach fails because he doesn't consider the transition between his two bend formulas. My approach succeeds in these cases because it relies on a single formula for all bend angles, instead of switching between formulas at critical bend points.

I derived equation (54) by the following logic. Define $y=\frac{\sin\theta}{\cos \theta}$ , and $\theta = \arctan y$ . Then, $\frac{\partial \theta}{\partial {\bf x}_{m_s}} =\cos^2\theta \frac{\partial y}{\partial {\bf x}_{m_s}}$ . Differentiating $y$ , we obtain

$$\frac{\partial y}{\partial {\bf x}_{m_s}} =\frac{\cos\theta \fra... ...\sin\theta \frac{\partial \cos\theta}{\partial {\bf x}_{m_s}} }{\cos^2\theta}$$

and thus,

$$\frac{\partial \theta}{\partial {\bf x}_{m_s}} = \cos\theta \fr... ...\bf x}_{m_s}} -\sin\theta \frac{\partial \cos\theta}{\partial {\bf x}_{m_s}}$$

The derivation of (57) follows directly from differentation of this equation. In some ways, this definition feels a little circular due to the introduction of the artifical $\arctan$ term. I believe that the derivation is nevertheless correct. One area of serious concern, however, is the fact that the second derivative does not exhibit typical symmetry: $\frac{\partial^2 \theta}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} \neq \frac{\partial^2 \theta}{\partial {\bf x}_{n_t} \partial {\bf x}_{m_s}}$ in general.

5.2 Baraff & Witkin's Errors

I believe that there is a very small error in equation (14) of Baraff & Witkin's paper: the ${\bf z}$ term should be scaled by $h$ .

There is an error in their definition of $\frac{\partial C}{\partial t}$ , just above equation (11). They define it as

$$\frac{\partial \bf C}{\partial t} = \left( \frac{\partial \bf C}{\partial {\bf x}} \right)^T \frac{\partial {\bf x}}{\partial t}$$

However, the condition function ${\bf C}$ is a function of the position of several particles ${\bf x}$ , and the paper does not specify which particle to choose. For my early implementations, I defined a different $\frac{\partial \bf C}{\partial t}$ for each particle, using

$$\frac{\partial {\bf C}_m}{\partial t} = \left( \frac{\partial \bf C}{\partial {\bf x}_m} \right)^T \frac{\partial {\bf x}_m}{\partial t} .$$

However, this yielded incorrect results, especially when calculating the damping forces for the bend condition, resulting in frequent instability. Later, I realised that correct application of the chain rule should yield a single $\frac{\partial \bf C}{\partial t}$ for all particles, of the form

$$\frac{\partial \bf C}{\partial t} = \sum_m\left( \frac{\partial \bf C}{\partial {\bf x}_m} \right)^T \frac{\partial {\bf x}_m}{\partial t} .$$ (58)

The sum should be over the particles that influence ${\bf C}$ . Visually, this yields vastly better results.

The final error in their paper is related to scale-invariance. Ideally, the parameters used for the simulation should be independent of the tesselation of the surface. For example, if you run the simulation with a $10\times10$ mesh and a $20\times20$ mesh and the same parameters, it should look basically the same. However, using Baraff & Witkin's choice of $\alpha=a$ , this does not occur. I looked at the math for the specific case where the number of particles is doubled in $u$ and $v$ (e.g., from $10\times10$ to $20\times20$ ). In this case, the gravity force on each original particle decreases by a factor of four (since the mass of the surrounding triangles decreases by four). However, the (strong) stretch and shear forces decrease by a factor of eight. These forces are responsible for most of the resistance against gravity, and consequently the cloth sags substantially as the particle density is increased.

However, by using $\alpha=a^\frac{3}{4}$ , the stretch and shear forces also decrease by a factor of four, as expected. The formula is not as intuitively obvious as $\alpha=a$ , and I have not done tests with irregularly triangulated meshes, but it does work much better for regularly triangulated meshes.

5.3 Unit Magnitude Approximation

During my original development, I chased down one wrong alley while looking for the error in my (and Macri's) derivations. I was highly suspicious of one of Baraff & Witkin's approximations, where they assume that vectors' magnitudes are constant. This occurs in both the shear condition derivation and the bend condition derivation. I went through and rederived all of the formulas without this assumption, and obtained the formulas given in the appendix. After all of this work, there was still a problem, which proved to be the error in Macri's versions of equations (54) and (57).

After finding this error in Macri's work, I returned to the unit magnitude issue. I compared the results of my simulation using the formulas given in the body of the paper (using Baraff & Witkin's approximation) to the results using the formulas given in the appendix (without the approximation). In my test case, I found that the approximation gave more pleasing results! In total, stretch energy was halved, and shear and bend energy dropped by a small amount. The formulas for both stretch and bend energy are unchanged between the two cases, so this is a valid comparison of the energy in the cloth. The formula for shear energy is slightly changed by removing the unit magnitude approximation, however, so comparisons of energy aren't very meaningful.

To be honest, I was surprised that such a drastic change to the formulas had a relatively small effect on the final result. I suspect that there may be a lot of room for optimisation of Baraff & Witkin's system, achieving comparable results at less computational cost.

In the interests of presenting a direct implementation of Baraff & Witkin's paper, and on the basis of the better results, I have therefore presented the formulas with the constant magnitude approximation.

5.4 Parameters

For parameters, I adopted values that gave good results. I looked at Alias $\vert$ Wavefront's Maya cloth software, but I suspect that they scale their parameters internally by an unknown amount. For damping forces, I multiplied the associated regular force's constant by the damping constant (e.g. $k_{stretch} k_{damp}$ was used to damp stretch forces).

Table 1: Parameters used for simulation.

Name	Symbol	Value
Stretch force strength	$k_{stretch}$	$5.0 \times 10^3$
Shear force strength	$k_{shear}$	$0.5 \times 10^3$
Bend force strength ($u$ )	$k_{bend_u}$	$0.01 \times 10^{-3}$
Bend force strength ($v$ )	$k_{bend_v}$	$0.01 \times 10^{-3}$
Damping force strength	$k_{damp}$	0.2
Rest stretch ($u$ )	$b_u$	1
Rest stretch ($v$ )	$b_v$	1
Time step	$h$	0.02
Density	$\rho$	0.1
Gravity	$g$	9.81

6 Results

To test out the simulator, I used a variety of cloth sizes and constraints. These can be seen at the Freecloth website http://davidpritchard.org/freecloth, in the screenshots section.

The simulator's performance is reasonable, with a $66\times66$ cloth typically requiring about thirty to forty minutes of computation per second of animation on a 3 GHz Pentium IV system using a step size of $0.02$ . Smaller cloth sizes are much faster, and there is still a fair bit of room for optimisation, especially in terms of adaptive timesteps.

My implementation does exhibit some instability, but almost always works with a step size of 0.01 seconds; sometimes, a smaller stepsize is needed.

7 Conclusion

Baraff & Witkin's paper does indeed provide a solid framework for cloth simulation. It appears to allow relatively large timesteps, and with suitable optimisation should run quite quickly.

Hopefully, this paper should clear up some of the confusions of both Macri's paper, and the original paper. Implementing Baraff & Witkin's paper should be straightforward using the equations shown here.

8 Acknowledgement

This work was supported in part by a scholarship from the Natural Sciences and Engineering Research Council of Canada, and by the British Columbia Advanced Systems Institute.

Thanks to Sebastian Lipponer for pointing out a correction to the bend condition formulas.

Bibliography

BW98

D. Baraff and A. Witkin.
Large steps in cloth simulation.
In ACM SIGGRAPH Conference Proceedings, pages 43-54, 1998.

Mac00

D. Macri.
Real-time cloth.
In Game Developers Conference Proceedings, 2000.

A. Appendix

The equations given in this appendix are replacements for certain equations given in the body of the paper. By using these, several of Baraff & Witkin's approximations are avoided. However, my tests show that the results are less pleasing.

A..1 Per-Triangle Quantities

The first and second derivatives of the normalised vectors ${\bf {\hat w}}_u$ and ${\bf {\hat w}}_v$ are

$$\frac{\partial {\bf {\hat w}}_u}{\partial {\bf x}_{m_s}} = \frac{1}{\begin{Vmatrix}{\bf w}_u\end{Vmatrix}} \frac{\partial ... ...al {\bf x}_{m_x}} \left( {\bf I}_s- {\bf {\hat w}}_{u_s}{\bf {\hat w}}_u\right)$$ (A.1)

$$\frac{\partial {\bf {\hat w}}_v}{\partial {\bf x}_{m_s}} = \frac{1}{\begin{Vmatrix}{\bf w}_v\end{Vmatrix}} \frac{\partial ... ...al {\bf x}_{m_x}} \left( {\bf I}_s- {\bf {\hat w}}_{v_s}{\bf {\hat w}}_v\right)$$ (A.2)

$$\frac{\partial^2 {\bf {\hat w}}_u}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = -\frac{1}{\begin{Vmatrix}{\bf w}_u\end{Vmatrix}} \frac{\partial {... ...\partial {\bf {\hat w}}_{u_s}}{\partial {\bf x}_{n_t}} {\bf {\hat w}}_u \right)$$

$$- \frac{\partial {\bf w}_{u_x}}{\partial {\bf x}_{m_x}} \frac{\pa... ...\bf {\hat w}}_u) {\bf {\hat w}}_{u_t}}{\begin{Vmatrix}{\bf w}_u\end{Vmatrix}^2}$$ (A.3)

$$\frac{\partial^2 {\bf {\hat w}}_v}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = -\frac{1}{\begin{Vmatrix}{\bf w}_v\end{Vmatrix}} \frac{\partial {... ...\partial {\bf {\hat w}}_{v_s}}{\partial {\bf x}_{n_t}} {\bf {\hat w}}_v \right)$$

$$- \frac{\partial {\bf w}_{v_x}}{\partial {\bf x}_{m_x}} \frac{\pa... ...\bf {\hat w}}_v) {\bf {\hat w}}_{v_t}}{\begin{Vmatrix}{\bf w}_v\end{Vmatrix}^2}$$ (A.4)

A..2 Shear Condition

If we remove the assumption that the ${\bf {\hat w}}_u$ and ${\bf {\hat w}}_v$ vectors are unstretched, then equation (19) should be replaced with

$$C = \alpha {\bf {\hat w}}_u\cdot {\bf {\hat w}}_v$$ (A.5)

$$\frac{\partial C}{\partial {\bf x}_{m_s}} = \alpha \left( \frac{\partial {\bf {\hat w}}_u}{\partial {\bf x}_{... ...hat w}}_u\cdot \frac{\partial {\bf {\hat w}}_v}{\partial {\bf x}_{m_s}} \right)$$ (A.6)

$$\frac{\partial^2 C}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} = \alpha \bigg( \frac{\partial^2 {\bf {\hat w}}_u}{\partial {\bf x}... ...l {\bf x}_{m_s}} \cdot \frac{\partial {\bf {\hat w}}_v}{\partial {\bf x}_{n_t}}$$

$$+ \frac{\partial {\bf {\hat w}}_u}{\partial {\bf x}_{n_t}} \cdot ... ...tial^2 {\bf {\hat w}}_u}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}} \biggr)$$ (A.7)

A..3 Bend Condition

If we remove the assumption that the normal and edge vectors have constant length, we must replace equations (40) and (43) with those shown below.

$$\frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{m_s}} = \left({\bf I}- {\bf {\hat n}}^A({\bf {\hat n}}^A)^T\right) \frac{ S_s( {\bf q}^A_m ) }{ \begin{Vmatrix}{\bf n}^A\end{Vmatrix} }$$ (A.8)

$$\frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{m_s}} = \left({\bf I}- {\bf {\hat n}}^B({\bf {\hat n}}^B)^T\right) \frac{ S_s( {\bf q}^B_m ) }{ \begin{Vmatrix}{\bf n}^B\end{Vmatrix} }$$ (A.9)

$$\frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{m_s}} = ({\bf I}- {\bf {\hat e}}{\bf {\hat e}}^T) \frac{ q^e_m {\bf I}_s}{ \begin{Vmatrix}{\bf e}\end{Vmatrix} }$$ (A.10)

$$\frac{\partial^2 {\bf {\hat n}}^A}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$= - \left( \frac{\partial {\bf {\hat n}}^A}{\partial {\bf x}_{n_t}}... ...right)^T \right) \frac{S_s({\bf q}^A_m)}{\begin{Vmatrix}{\bf n}^A\end{Vmatrix}}$$

$$+ \frac{{\bf I}- {\bf {\hat n}}^A({\bf {\hat n}}^A)^T}{\begin{Vma... ...{Vmatrix} S_s\left( \frac{\partial {\bf q}^A_m}{\partial {\bf x}_{n_t}} \right)$$

$$\quad\quad\quad\quad\quad\quad\quad - \left( S_t({\bf q}^A_n) \cdot {\bf {\hat n}}^A \right) S_s({\bf q}^A_m) \bigg)$$

$$\frac{\partial^2 {\bf {\hat n}}^B}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$= - \left( \frac{\partial {\bf {\hat n}}^B}{\partial {\bf x}_{n_t}}... ...right)^T \right) \frac{S_s({\bf q}^B_m)}{\begin{Vmatrix}{\bf n}^B\end{Vmatrix}}$$

$$+ \frac{{\bf I}- {\bf {\hat n}}^B({\bf {\hat n}}^B)^T}{\begin{Vma... ...{Vmatrix} S_s\left( \frac{\partial {\bf q}^B_m}{\partial {\bf x}_{n_t}} \right)$$

$$\quad\quad\quad\quad\quad\quad\quad - \left( S_t({\bf q}^B_n) \cdot {\bf {\hat n}}^B \right) S_s({\bf q}^B_m) \bigg)$$

$$\frac{\partial^2 {\bf {\hat e}}}{\partial {\bf x}_{m_s} \partial {\bf x}_{n_t}}$$

$$= - \left( \frac{\partial {\bf {\hat e}}}{\partial {\bf x}_{n_t}} {... ...{n_t}} \right)^T \right) \frac{S_s(q^e_m)}{\begin{Vmatrix}{\bf e}\end{Vmatrix}}$$

$$- \frac{{\bf I}- {\bf {\hat e}}{\bf {\hat e}}^T}{\begin{Vmatrix}{\bf e}\end{Vmatrix}^2} q^e_m q^e_n {\bf {\hat e}}_t {\bf I}_s$$ (A.11)