PyCG 6: Cloth Simulation

17 Aug 2019

The computer graphics lecture at Berkeley constantly uses awesome cloth simulation videos (e.g. this one) to demonstrate the beauty of computer graphics, and I was truly absorbed by it. As it turns out, the mathematics behind it is not very difficult to understand, and it can be fairly easy to build a simple simulation.

The mass-spring model

The nature of cloth behavior can be approximated by the mass-spring model¹ at a low cost, albeit the existence of more realistic and complicated methods such as Finite Element Analysis². The cloth is described as an \(m\times n\) grid of point mass, and the point masses are connected by springs. To simulate cloth more accurately, for a given point mass \((i, j)\), there are three types of springs:

1. structural spring: the links between the point and \((i + 1, j)\), \((i - 1, j)\), \((i, j + 1)\), \((i, j - 1)\).
2. shear spring: the links between the point and \((i + 1, j + 1)\), \((i + 1, j - 1)\), \((i - 1, j + 1)\), \((i - 1, j - 1)\).
3. bend spring: the links between the point and \((i + 2, j)\), \((i - 2, j)\), \((i, j + 2)\), \((i, j - 2)\).

With symbols described in the table below, the forces considered in this model are listed as follows:

Symbol	Description
\(\bm{P}_{i, j}\)	The spatial position on point \((i, j)\)
\(\bm{v}_{i, j}\)	The velocity of point \((i, j)\)
\(\bm{a}_{i, j}\)	The acceleration of point \((i, j)\)
\(\bm{n}_{i, j}\)	The surface normal at point \((i, j)\)
\(\mu\)	The mass of each point

1. Internal force:

   Denote the position of point \((i, j)\) by \(\bm{P}_{i, j}\), the internal forces caused my the spring grid for \((i, j)\) is given by

   \[\begin{align*} \bm{F}_{\mathrm{int}}(i, j) = \sum_{(k, l) \in \mathcal{R}_{i, j}} K \left( \bm{l}_{i, j, k, l} - l^0_{i,j, k, l} \frac{\bm{l}_{i,j, k, l}}{\lVert \bm{l}_{i,j, k, l} \rVert} \right), \end{align*}\]

   where \(\mathcal{R}_{i, j}\) is the set of all existing links of point \((i, j)\); \(K\) is stiffness; \(\bm{l}_{i, j, k, l} = \overrightarrow{\bm{P}_{i, j}\bm{P}_{k, l}}\), and \(l^0_{i,j, k, l}\) is the natural length of spring \((i, j) \rightarrow (k, l)\).

2. Viscous damping:

   To simulate the dissipation of mechanical energy in the model, the viscous damping for point \((i, j)\) is given by

   \[\begin{align*} \bm{F}_{\mathrm{dis}}(i, j) = -C_{\mathrm{dis}} \bm{v}_{i, j}, \end{align*}\]

   where \(C_{\mathrm{dis}}\) is the damping coefficient.

3. Gravity:

   The force of gravity for point \((i, j)\) is given by

   \[\begin{align*} \bm{F}_{\mathrm{gr}}(i, j) = \mu g, \end{align*}\]

   where \(g\) is the gravity constant.

4. Wind:

   Assume the wind has a constant speed of \(\bm{v}_{\mathrm{wind}}\), the force of wind for point \((i, j)\) is given by

   \[\begin{align*} \bm{F}_{\mathrm{wind}}(i, j) = C_{\mathrm{wind}} \left[\bm{n}_{i, j} \cdot (\bm{v}_{\mathrm{wind}} - \bm{v}_{i, j}) \right] \bm{n}_{i, j}, \end{align*}\]

   where \(C_{\mathrm{wind}}\) is wind strength constant, and “\(\cdot\)” denotes dot product.

The force that the point \((i, j)\) receives (denoted by \(\bm{F}_{i, j}\)) is the sum of all forces listed above.

Solving for the next state

Given the state of the entire system at time \(t\), which contains \(\bm{a}_{i, j}(t)\), \(\bm{v}_{i, j}(t)\), \(\bm{P}_{i, j}(t)\), we would like to compute the next state after a small time increment \(\Delta t\). By using a simple Euler method, it can be seen that

\[\begin{equation*} \begin{cases} \bm{a}_{i, j}(t + \Delta t) = \bm{F}_{i, j}(t) / \mu\\ \bm{v}_{i, j}(t + \Delta t) = \bm{v}_{i, j}(t) + \Delta t \bm{a}_{i, j}(t + \Delta t)\\ \bm{P}_{i, j}(t + \Delta t) = \bm{P}_{i, j}(t) + \Delta t \bm{v}_{i, j}(t + \Delta t) \end{cases}. \end{equation*}\]

It is very straightforward and computationally cheap. However, one has to care about the choice of \(\Delta t\). The linear system above is divergent when \(\Delta t\) is greater than the natural period of the system, which is given by

\[\begin{align*} T_0 \approx \pi \sqrt{\frac{\mu}{K}}. \end{align*}\]

That is why the demo gif of this article needs to be accelerated five times, because I need to select a small \(\Delta t\) in order to keep the system stable. But doing this will increase the duration of the animation.

Implementing these in Python

Again, when the mathematical basis of this problem is all set, implementing the visualization can be trivial. However, because Python is a language of inferior performance, a good implementation needs to avoid using loops and try to use numpy intensively. I try to do this throughout my code, and the only part that still needs Python loop is the computation of surface normals, for the process needs to find the surface normal of a point, which is the averaged on a variant amount of surface normals acquired in triangulation. Nonetheless, the performance of the program is still not very ideal. When using a grid larger than \(20 \times 20\), the frame rate can drop to less than 10 fps. To make a more appealing animation, I only apply a small regional wind on the bottom-right region of the cloth.

This is actually simplified version of an assignment of Berkeley CS184 in spring 2018. I think it is impractical to implement collision management in Python, because it is too computationally intensive for the language. But I think the simple Python solution to this problem is still very beautiful - with only around 250 lines of code, we can have a (slow) solution to this problem. And OpenGL is the hardcore way of visualizing the result, it could have been done with matplotlib with several lines of code. Therefore, I still think that for a project like this which is seemingly not suitable for Python, it is always OK to start with Python. After all, if not for the fancy numpy code that I tried to write and debug for 2-3 hours, the entire solution could have been finished within 1-2 hours. It is generally awesome to have a prototype solution for a problem in such a short time.

Demo

The source code of the demo can be found on GitHub.

Basic usages:

• Use mouse and “wasd” to look/walk around
• Press Esc to exit
• Press “p” to take screenshots
• Press “o” to toggle fill/wireframe mode