Cloth Simulation and Collision Detection using Geometry Images Nico Zink∗ University of Johannesburg

Abstract The simulation and animation of cloth has attracted considerable research interest by the computer graphics community. Cloth that behaves realistically is already expected in animated films, and real-time applications are certain to follow. A common challenge faced when simulating the complex behaviour of cloth, especially at interactive frame rates, is maintaining an acceptable level of realism while keeping computation time to a minimum. A common method of increasing the efficiency is a decrease in the number of nodes controlling the cloth movement, sacrificing details that could only be obtained using a dense discretization of the cloth. A simple and efficient method to simulate cloth is the mass-spring system which utilises a regular grid of vertices, representing discrete points along the cloth’s surface. The structure of geometry images is similar, which makes them an ideal choice for representing arbitrary surface meshes in a cloth simulator whilst retaining the efficiency of a mass-spring system. In this paper we present a novel method to apply geometry images to cloth simulation in order to obtain cloth motion for surface meshes of arbitrary genus, while retaining the simplicity of a mass-spring model. We also adapt an implicit/explicit integration scheme, utilising the regular structure of geometry images, to improve performance. Additionally, the cloth is able to drape over other objects, also represented as geometry images. Our method is efficient enough to allow for fairly dense cloth meshes to be simulated in real-time. CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modelling—Physically based modelling; I.3.7 [Computer Graphics]: Three Dimensional Graphics and Realism—Virtual Reality Keywords: cloth simulation, collision detection, geometry images, real-time rendering

1

Introduction

The simulation and animation of cloth has interested the computer graphics community for the last 20 years since this is a field that is both very challenging and applicable to many problems. Extensive work has been done and many groups and individuals researching the subject have made a great deal of progress. The most common demand for simulated cloth has come from production animation, since an animator is able to apply visually attractive and plausible clothing to a character without having to painstakingly animate it by hand. Virtual cloth is also useful to garment and fashion designers who can see and interact with the dynamics of their designs before manufacture. Another field benefiting from the simulation of cloth is the e-commerce field, where plans exist to ∗ e-mail:

[email protected]

† e-mail:[email protected]

Alexandre Hardy† University of Johannesburg implement a system where a user can upload a model representing his/her own body and test the particular look of various clothing [Protopsaltou et al. 2002]. Cloth simulation is also applicable to the automation of textile production lines [Eischen et al. 1996]. In order to meet these demands, the computer graphics community started researching the simulation of cloth, and have identified two common bottlenecks for the performance of creating virtual cloth: the numerical simulation of the cloth’s properties and collision detection with the objects in the scene and with the cloth itself. The numerical simulation of cloth concerns capturing the properties describing cloth in a mathematical model, and deriving the shape of a cloth object in a manner that resembles real-life cloth. This is very difficult because the behaviour of cloth is very complex and highly non-linear. Adding to this problem is the fact that cloth hardly stretches under normal conditions, giving rise to very stiff equations. This generally requires more computational time to ensure convergence of the solution, as most numerical solvers cannot deal with these conditions efficiently. Collision detection is a bottleneck in cloth simulation since testing the cloth for collision with all of the objects present in the scene requires a lot of time. Also, for the cloth to be able to look realistic, self-collision detection must be performed. This means that, in order for the cloth to interact with itself (e.g. fall onto itself), the surface of the cloth must be tested for collision with the rest of its surface. This has proven to be quite challenging, with collision and self-collision detection taking up huge slices of the time required to calculate each frame of a simulation. Visually attractive and plausible animation of cloth has been increasing in priority and is now a requirement for most characters produced by the animation and special effects industry. These requirements will almost certainly extend to virtual reality applications such as games in the future as techniques improve and increasingly powerful computational resources become available. The computational costs are still high, and various problems continue to make simulating virtual cloth a difficult task. There are many factors that make the simulation of cloth with the consideration of all these criteria difficult. This means that generally some kind of compromise must be made. This is especially true in highly constrained situations: such as when demands for real-time performance have to be met, or when collisions and self-collisions in complex environments create complicated cloth movement. Many advances in the simulation of cloth provide realistic motion of cloth, but the amount of computational time limits many applications to off-line rendering. The problem of real-time cloth simulation remains a challenge, and many solutions suffer from instability. Implicit integration [Baraff and Witkin 1998; Desbrun et al. 1999] resolves the but methods generally require ´ ` instability, a running time of O n2 where n is the number of particles present in the system, limiting the scalability of the system. A common method used to reduce computational time is to decrease the resolution of the cloth, sacrificing valuable detail that can only be obtained using a dense discretization of the cloth. For real-time applications, the mass-spring system proves to be the most popular cloth model. The cloth’s particles are arranged in a regular manner with simple linear springs generating forces to determine the displacement of each particle. The structure of

Geometry images were introduced by Gu et al. [2002] to store surface geometry in a completely regular structure. A geometry image is a mesh stored in an ordinary image file consisting of an n × n array of x, y and z values. Surface signals, such as normal maps, textures, etc., are stored in similar image files. Geometry images can be created by cutting a mesh such that the topology resembles that of a disk, and then parameterising the surface of this disk onto the domain of a rectangle, which is sampled at regular intervals with the positions of the samples stored in the geometry image. (a)

(b)

Figure 1: Torus model: (a) Geometry image; (b) Normal mapped reconstruction.

geometry images is similar, which makes them an ideal choice for representing arbitrary surface meshes in a cloth simulator while retaining the efficiency of a mass-spring system. In this paper we present a new method of simulating cloth with geometry images. Geometry images present a way of representing a triangulated mesh as a regular 2-dimensional image file, as shown in figure 1. Our cloth simulator uses a multilevel approach based on an implicit/explicit mass-spring model to efficiently handle fairly dense meshes stored as geometry images. The rest of this paper is organised as follows: In section 2, a brief summary of the previous work in the areas of cloth modelling, geometry images and collision detection is given. Section 3 describes the process of creating geometry images. The cloth model is presented in section 4 and section 5 describes the collision detection. The results and conclusions are presented in sections 6 and 7.

2 2.1

Related Work Geometry Images

Surface geometry is often stored as an irregular and unsorted series of triangles, also referred to as “polygon soup”. This structure is the most flexible, since it allows for any mesh of arbitrary size, shape and connectivity to be stored. Due to the unordered nature of the information, however, mesh traversal is complicated. As a result, many operations cannot be accomplished efficiently. This representation generally requires some kind of processing of the information to establish the connectivity of the mesh. The information is arranged into various interconnected data structures. The surface of the mesh can be explored by performing operations on neighbouring nodes in the data structure used to store the mesh. The data structure needs to be highly dynamic to be able to store a generic mesh arranged in any order. Since meshes often contain a very large number of triangles, traversal of the data structure needs to be efficient. The surface can be re-meshed into a structure that contains semiregular connectivity. This has been done by separating the mesh into several charts and packing them into an atlas. The individual charts are regular in the sense that neighbouring nodes in 3D space are also neighbouring nodes in storage and connectivity within the charts is standard. This property greatly improves operations like texture mapping. Textures can be separated into charts equivalent to the corresponding charts of the mesh. Although the individual charts are all structured regularly, the overall geometry still needs to be constructed by re-assembling the separate charts.

The triangle mesh is thus converted into a completely regular structure that can be stored in a bitmap type file. This bitmap can then be sent directly to the graphics hardware, thereby simplifying the rendering pipeline by avoiding expensive pointer structures. Since vertices are stored in a two-dimensional array-like structure with neighbouring points next to each other, many operations can be performed more efficiently. The finding of adjacent vertices, faces and the like can be performed by inspection without complex data structures or methods of graph traversal. Normal and texture mapping are also simplified since texture coordinates are implicit. Textures and normal maps are stored in similar files to the geometry images so mapping onto the geometry is straight forward. Image manipulation, such as compression, can be applied to the geometry image. Naive level-of-detail can be performed by simply mipmapping the geometry image. The process of creating geometry images is presented in detail in section 3. The rest of this section presents an overview of the applications that geometry images have been applied to. The parameterisation of a surface mesh onto a rectangular domain is difficult in some cases and impossible in others. Multi-chart geometry images [Sander et al. 2003] split these surfaces into charts, which are sampled and fit onto a rectangular bitmap. A zippering algorithm is applied to connect the charts to form the representation of the mesh. GI-COLLIDE [Beneˇs and Villanueva 2005] is a method of performing collision detection with (2n + 1) × (2n + 1) geometry images. Spheres are used as a bounding volume hierarchy. A sphere tree is built from the bottom up by placing spheres around the area between every quadruple of vertices, every quadruple of spheres is bounded by a larger sphere, and so forth, until the whole geometry image is bounded by a single sphere. Two geometry images can be tested for collision by alternately recursing the two sphere trees and testing the spheres for intersection. When the lowest level is reached the exact collision test is performed.

2.2

Cloth Simulation

Since most of the demand for virtual cloth came from production animation, the prominent design criteria for many researchers has been the visually plausible and attractive animation of cloth, while maintaining an acceptable level of realism. Many proposed cloth models have achieved complex and realistic results. Interest in real-time cloth has been increasing, where efficiency is the biggest priority. Mass-spring models are the most efficient as well as the simplest of the cloth models. This method is one of the most popular techniques for simulating cloth, especially when interactive frame rates are required. Provot [1995] presented a mass-spring model where cloth is discretized by using a regular grid of nodes. Stretch springs connect every node to its adjacent neighbours, shear springs connect each node to its diagonal neighbours and bend springs connect each

particle to its second neighbours. The forces fi and fj generated by a spring connecting particles i and j are equal to « „ li,j fi = −ki,j li,j − r0 kli,j k fj = −fi where ki,j is the stiffness constant of the spring, r0 is the predefined rest length, and li,j = xj − xi , where xi and xj are the displacements of the nodes i and j. The forces are integrated using an explicit Euler method, and the resulting stiffness problems are dealt with using an inverse dynamics procedure: whenever a spring elongates above a certain threshold from its rest length, the spring is artificially shortened by placing the nodes closer together. Unfortunately, the explicit Euler scheme is not very accurate. The accuracy decreases even more in cases of large values of the time step size ∆t. As a result, explicit integration diverges in stiff systems with large time steps. This is due to excessive forces overcompensating for the position of particles due to them being kept constant over too long a time step. The stiffness problem can be addressed by turning to implicit integration. Implicit integration is able to remain unconditionally stable independent of the time step used. The disadvantage is that a large linear system needs to be solved. Baraff and Witkin [1998] use the sparse linear equation „ « „ « ` ´ ∂f ∂f ∂f M − ∆t − ∆t2 ∆v = f xtn , vtn + ∆t vtn ∂v ∂x ∂x ∂f ∂f where the Jacobian matrices ∂x and ∂v are ` evaluated ´ based on tn tn the current time step, i.e. x , v and f xtn , vtn . At every time step, this equation is efficiently solved for ∆v using their iterative modified preconditioned conjugate gradient algorithm. This method enables the simulation of stiff springs, even when using large time steps. Unfortunately, whilst computational time can be low enough ` ´ for interactive frame rates, the size of the linear system is O n2 where n is the number of particles present in the system. This means that an increase in the number of particles drastically reduces efficiency.

Desbrun et al. [1999] present a method that approximates the implicit integration found in the Baraff and Witkin [1998] cloth model. The precomputed inverse of the Hessian matrix „ «−1 ∆t W = I− H m is used as a force filter at every time step to smooth out any irregular cloth movement. This method is far more efficient than the original model, but it still needs an n × n matrix which is used to filter the forces, so the algorithm does not scale well.

2.3

Collision Detection

Collision detection with cloth is an expensive task. The cloth must be tested against a wide variety of objects. It is desirable, though very expensive, to have the cloth also perform self-collision detection. To test for self collision each part of the cloth must be tested against each other part of the cloth. Traditional rigid body methods for collision detection are too slow to make this approach feasible with high frame rates. This section presents a brief survey of existing collision detection methods. Collision detection for rigid bodies has been well researched. Methods generally rely on a hierarchy of various bounding volumes to cull the number of low level intersection tests that need to be

performed. Common methods include I-COLLIDE [Cohen et al. 1995], V-COLLIDE [Hudson et al. 1997], RAPID [Gottschalk et al. 1996], k-DOPs [Klosowski et al. 1998] and many more. A survey of rigid collision detection methods as well as a breakdown of certain properties of collision detection algorithms was presented by Lin and Gottschalk [1998]. The bounding volumes used in the collision testing are usually created in a pre-processing phase. They are optimised to cover as little area as possible and for inexpensive intersection tests against one another. After these bounding volumes have been precomputed they can be efficiently accessed at run time. This is generally not efficient for deformable models. Collision detection with deformable models is more problematic since changes in the underlying positioning would force these bounding volumes to have to be updated or re-calculated at each time step. Also, deformable models may collide with themselves, creating a need for self-collision detection. Some algorithms that move away from hierarchical bounding volumes are being explored, such as when spacial subdivision [Teschner et al. 2003] is being used. Spacial subdivision divides the area surrounding the models into smaller areas. A hashing function is used to quickly assign objects into groups. This method is often fast enough to support self-collision detection. Lately several collision detection techniques based on imagespace techniques have been proposed. These methods work with projections of the meshes to decide which parts of the mesh overlap. These methods do not use any pre-computed data and can therefore easily be applied to deformable models. Another advantage is the fact that they can be implemented on programmable graphics hardware. An example of this is Quick-CULLIDE [Manocha 2005]. Based on CULLIDE [Chen et al. 2004], which tested for collisions between multiple deformable and breakable objects, Quick-CULLIDE tests for self-collision as well. CULLIDE works by constructing a potentially colliding set (PCS) of all objects that could possibly collide. The scene to be tested for collision is sent to the graphics hardware. From there it is rasterized from two sides: first the front and then the back. The results are checked to see which objects are occluded by other objects. Objects that are fully visible are removed from the PCS. Next all sub-objects of all the objects in the PCS are tested for visibility. Then exact collision detection is applied to the remaining objects left in the PCS. These tests are very fast and leave a relatively small PCS. Also, another benefit is that no pre-computation on the model is needed and the method is even feasible to use with “polygon soup”. This method also works with deformable objects. There is no support for self collision though. This makes collision detection fast enough for real time applications and removes the burden from the CPU. Teschner et al. [2005] present a survey of techniques for collision testing with deformable models.

3

Geometry Images

We will describe the process of creating a geometry image from a given triangle mesh based on the method presented by Gu et al. [2002]. Geometry images are created by parameterising the mesh onto a plane, sampling the surface of the plane at regular intervals, and finally storing the positions of these samples in a regular image file. In order for the planar parametrisation to be possible, the surface topology of the mesh must resemble that of a disk. Any mesh, with the exception of non-manifolds, can be

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 2: Mannequin model: (a) Original model with the basis of the cut path drawn in white; (b) Original Floater parametrisation; (c) Improved cut path; (d) Improved Floater parametrisation; (e) Stretch-minimising parametrisation; (f) Geometry image; (g) Normal map; (h) Normal mapped reconstruction.

unwrapped when cut open along an appropriate path. For any mesh where the topology does not resemble that of a disk, the algorithm must find a suitable cut along a number of edge paths. These edge paths, or cuts, allow the mesh to be opened up and thus be parametrised onto a rectangle. Since the resulting geometry will be stored in an image format, we chose to scale the object to fit inside the unit cube. This is done in order to conform to most image file standards, that do not allow negative values as well as placing a maximum bound on the colour signal size. The goal of the first step in creating geometry images is to form an appropriate cut path along the edges of the mesh in such a manner that the topology of the mesh resembles that of a disk, enabling the mesh to be parametrised onto a square. This cut path is found by removing a disk-like section from the mesh, with the remaining edges forming the cut path. The first step in cutting a mesh is to detect all boundary edges. These edges must remain fixed throughout the procedure and must be a subset of the final cut path. Then a random seed triangle is chosen and removed from the mesh. The algorithm now proceeds in two phases. In the first phase edges adjacent to only one face, but not on the mesh boundary, are chosen one at a time for removal. Both the chosen edge and the connected face are removed from the mesh. Removals are ordered according to the geodesic distance [Surazhsky et al. 2005] between the edge and the original seed triangle. When this phase terminates, all of the faces will have been removed. In the second phase all unnecessary edges are removed. We repeatedly identify a vertex adjacent to only one edge. Both the vertex and the edge are removed. When this phase terminates all that is left is one or more connected loops. All of the remaining

edges form the basis of the cut. All that remains is to shorten the cut path where possible. Special care must be taken when cutting genus zero surfaces. If this is the case, all of the edges in the mesh will have been removed, leaving only a single vertex behind. Two adjacent edges are then added back to the vertex, so that they form the cut. Since this vertex was found arbitrarily, a better cut may be obtained by choosing another vertex. This is repeated until the cut no longer improves. The cut can now be optimised by iteratively augmenting the cut. It is desirable to have the cut pass through various areas that result in unnecessary geometric stretch [Sander et al. 2001] when the current cut is being used. Floater’s [1997] shape-preserving parametrisation is used to perform the parametrisation. The triangle with the most geometric stretch is located, and one of its vertices is connected to the cut path via the shortest path. An example of this applied to a mannequin model is shown in figure 2(a-d). This process is repeated until the geometric stretch under the Floater parametrisation no longer improves. This will yield the final cut path. The first step in obtaining a planar parametrisation is parameterising the boundaries. The cut consists of various cut nodes, and cut paths [Gu et al. 2002] connecting them. This is done by allocating an amount of the boundary of the rectangle to each cut path proportional to its length in the mesh. In order for the border of the geometry image to be reconstructed properly, the cut nodes are placed at the sample points of the geometry image, while all adjacent cut paths are allocated the same amount of space along the boundary. Each edge placed along one of the four corners is split into two new edges, while the same splitting is applied to its boundary mate. Now the interior of the mesh is parametrised. A simplified representation of the mesh is created and optimally placed within the rectangle by brute force. The simplified mesh is refined and

the new vertices are inserted while shifting neighbouring vertices to minimise stretch [Yoshizawa et al. 2004]. The parametrisation maps the mesh onto a rectangle. Within this rectangle, nodes along the boundary all correspond to nodes on the cut path, with the rest of the nodes on the inside. Each vertex on the rectangle corresponds to a vertex on the cut mesh. Each vertex on the rectangle is given the vertex’s position in the original mesh as a colour. The geometry image can now be sampled simply by rendering the rectangle and writing the result into an image file. Normal maps are created in a similar manner. The results are shown for a mannequin model in figure 2. The process of rendering the geometry image to a file often causes cracks to appear along the boundary. In order to be able to reconstruct the topology of the original geometry, the borders will have to be matched. We let the set b be all boundary nodes, consisting of b = {b1 ,b2 , . . . ,bn } where each bi consists of the corresponding boundary mates bi = {p1 ,p2 , . . . ,pni } that need to be matched during the reconstruction, where ni is the cut valence of bi .

4

Cloth Simulation

Since our primary concern is real-time performance, we have chosen to use a mass-spring model [Provot 1995]. The surface of the cloth is represented as a rectangular array of particles. Connectivity is maintained using three types of springs: • Stretch springs connect vertical and horizontal neighbours • Shear springs connect diagonal neighbours • Bend springs connect second vertical and horizontal neighbours. These springs generate linear forces on pairs of particles based on the spring’s elongation or compression relative to a predefined rest length. Mass-spring systems have important drawbacks. By dropping non-linear forces in favour of linear ones, springs are often given excessive permission to stretch, resulting in visibly elastic cloth. Another problem is the locality of the applied forces, as nodes only interact with their immediate neighbours. The locality becomes even worse with the use of a denser discretization of the cloth mesh. Despite these shortcomings, mass-spring systems are able to simply and efficiently create visually pleasing cloth, and have remained a popular method. An important part in the simulation of cloth is the numerical solution of the underlying mathematical model. The forces acting on each node are gathered by inspecting each spring, and are used to update the positions of the nodes. Integration is used to track the position of each cloth node by making their positions a function of time, which is successively evaluated at regular time intervals. Based on Newton’s second law of dynamics (f = m×a), the forces acting on the cloth give rise to the second-order ordinary differential equation x ¨ = M−1 f (x, x) ˙ where M is a diagonal 3n × 3n matrix with mii being equal to the mass of particle i. The force f consists of internal and external forces, i.e. f = fint + fext . The displacement x and its derivatives, the velocity v = x˙ and the acceleration a = x ¨, are vectors of size 3n and contain the displacement, velocity and acceleration of each node. The

relationship between the displacement x, velocity v and force f is formulated as the coupled first-order ordinary differential equation ! ! v d x . = M−1 f (x, v) dt v At time tn , given values for the displacement xtn and velocity vtn , the new position xtn +∆t and velocity vtn +∆t at time tn + ∆t need to be calculated. Most of the early work on cloth simulation used explicit time integration to advance the system. The explicit Euler method is the simplest of the integration schemes, approximating a function x at time tn + ∆t as xtn +∆t = xtn + ∆tx˙ tn . Values for xtn +∆t and vtn +∆t are thus calculated as ! ! ! xtn +∆t xtn vtn = + ∆t vtn +∆t vtn M−1 f (xtn , vtn ) Unfortunately, the explicit Euler scheme is not very accurate. The accuracy decreases even more in cases of large values of ∆t. Explicit integration methods suffer from instability, making the size of the time step dependant on the stiffness of the springs. This is due to excessive overcompensation on the positions of the particles due to the forces being kept constant over too long a time period. This means that the system moves forward in time blindly, as there is no guarantee that the force model ends up being in a consistent state at the end of the time step. As a result, explicit integration diverges in stiff systems with large time steps. This problem has been addressed by turning to implicit integration [Baraff and Witkin 1998]. Implicit integration is formulated as ! ! ! xtn +∆t xtn vtn +∆t = + ∆t vtn +∆t vtn M−1 f (xtn +∆t , vtn +∆t ) By approximating the forces f at time tn + ∆t, Baraff and Witkin [1998] derive the sparse linear equation „ « „ « ` ´ ∂f ∂f ∂f M − ∆t − ∆t2 ∆v = f xtn , vtn + ∆t vtn ∂v ∂x ∂x ∂f ∂f where the Jacobian matrices ∂x and ∂v based on the ` are evaluated ´ tn tn current time step, i.e. x , v and f xtn , vtn .

At every time step, this equation is efficiently solved for ∆v using their iterative modified preconditioned conjugate gradient algorithm. Given a value for ∆v, the new velocity vtn +∆t and position xtn +∆t are simply solved as vtn +∆t = vtn + ∆v xtn +∆t = vtn + ∆tvtn +∆t . What this means is that unlike explicit integration, the approximated forces at the end of the time step are used to update the velocity. Although forces are still assumed constant within the time step, this ensures a dramatic increase in stability. Implicit integration is able to remain unconditionally stable independent of the time step and spring stiffness used. This method enables the simulation of stiff springs, even when large time steps are being used. The disadvantage is that a large but sparse linear system needs `to be´ solved at every time step. The size of the linear system is O n2

(a)

(b)

Figure 3: (a) Wl propagates forces locally; (b) Wg propagates forces globally. where n is the number of particles present in the system. While the method provides interactive frame rates at low resolutions of the cloth, an increase in the number of particles drastically reduces efficiency. Desbrun et al. [1999] present a method that approximates the implicit integration found in the Baraff and Witkin [1998] cloth model. The precomputed inverse of the Hessian matrix «−1 „ ∆t2 H W = In − m is used to filter the forces at every time step. The n × n matrix H is computed per connected node pair i and j as ( kij if i 6= j P Hij = − i6=j kij if i = j. The matrix W is symmetric, and the sum of the entries in each row and column are all equal to 1. This matrix is used to propagate the individual forces through the entire cloth. For each particle i, the entries contained in row i are used to compute a weighted average of all particles, i.e. Particle j contributes Wij of its force to the total force of particle i. This is expressed as fi f iltered = Wi,1 ×f1 +Wi,2 ×f2 +. . .+Wi,n−1 ×fn−1 +Wi,n ×fn . W is used to distribute forces, pointing all forces in a similar direction, minimising local differences and improving stability by pointing particle velocities in similar directions. This method is far more efficient than the original model proposed by Baraff and Witkin [1998]. Even though H is sparse, W tends to ` be´ a dense matrix. Therefore, the running time tends to be O n2 , which creates problems for larger resolutions of the cloth. Our method is based on the implicit/explicit integration presented by Desbrun et al. [1999]. We use a multilevel approach to distribute the forces through the cloth more efficiently, in order to improve the performance for large cloth meshes. The cloth is stored as a geometry image of size n × n, giving us a total of n2 particles. Since we are now working with a grid, we let fi,j denote the force of the particle located in the ith row and the jth column of the cloth. Also, to simplify the following operations, indexing will start at 0 and end at n − 1. After all forces f have been computed, a small version of W, which we will refer to as Wl , is used to gather forces locally. This matrix is computed as „ «−1 ∆t l l , H W = In2 − l m

where Hl is computed based on a standard nl × nl mass-spring model, giving an n2l × n2l matrix, with nl defined by the user. The matrix Wl is used to gather the forces f into f l . For every particle, l fi,j is calculated by placing this nl × nl mass-spring model over the cloth, centred over particle i,j. An example of such a placement is shown in figure 3(a). It is required that nl must be odd to ensure that the mass-spring model has a centre that can be used to gather forces for a given particle. If the centre particle cannot be used due to there being an overlap with the border of the cloth, the matrix position is adjusted by aligning the side of the mass-spring model with the border, and using the relevant row of the matrix for the l calculation of fi,j . The mass-spring model is placed over the cloth min with its corner placed at xmin and yi,j , calculated as i,j 8 > if i < (nl − 1)/2 <0 xmin = (n − 1) − (nl − 1) if i > (n − 1) − nl2−1 i,j > :i − nl −1 otherwise 2 8 > if j < (nl − 1)/2 <0 min yi,j = (n − 1) − (nl − 1) if j > (n − 1) − nl2−1 > :j − nl −1 otherwise. 2 After the placement has been performed, the matrix Wl is used to gather neighbouring forces for the current particle. The force calculation for particle i,j is done as nl −1 nl −1 l fi,j =

X X

Wwl i,j ,wk,l fk+xmin ,l+ymin i,j

i,j

k=0 l=0

` ´ ` ´ min where wi,j = i − xmin nl + j − yi,j and wk,l = knl + l. i,j We precompute another version of W, called Wg , to globally distribute the gathered forces from f l into f g . Wg is calculated similarly to Wl with „ «−1 ∆t g Wg = In2g − H , m and is of size n2g × n2g . The matrix Hg is computed based on the mass-spring system that is found by sampling the geometry image at a lower resolution of ng × ng , as shown in figure 3(b). It is required that n − 1 must be a multiple of ng − 1, so that the global nodes are spaced evenly, ensuring that: • Forces are distributed evenly • Borders remain matched. The construction of Wg must take borders into account, so rows and columns are merged for adjacent boundary nodes during the construction of Hg . For each set of boundary mates bi = {p1 ,p2 , . . . ,pni } in b, the relevant rows of Hg are merged. Wg is used to gather forces of every (n−1) th particle along the cloth’s (ng −1) g rows and columns. The force fi,j is calculated as g fi,j =

ng ng X X k=1 l=1

Wwg i,j ,wk,l fkl

n−1 ,l n−1 ng −1 ng −1

“ ” ng −1 ng −1 where wi,j = i n−1 − 1 nl +j n−1 and wk,l = (k − 1) nl +l. The rest of the global forces are then linearly interpolated to obtain forces for each particle, and are added to the gathered forces in f l to obtain f t , which is used to perform the integration. Since the matrices involved are all smaller than W, our approach scales much more efficiently than the original method.

As was noted by Desbrun et al. [1999], such an approximation to implicit integration does not preserve angular momentum. They provide a simple and efficient post-correction method to preserve this quantity, which we also follow. The global torque is calculated as n X δT = f tn × xi , i=1

which is used to correct the position of each particle i by calculating the term fc = (xG − xi ) × δT∆t, 2

and added to i’s position. which is integrated as fc ∆t m Since we are using linear spring forces, these often allow for large spring displacements when large forces are present. We use the method first used by Provot [1995] to correct for this. We iteratively inspect each spring. Whenever a spring exceeds its rest length r by a predefined threshold , both ends are pulled back together so that the elongation rate is exactly . Our method is presented as pseudocode in figure 4.

5

Collision Detection

Any object present in the scene will be represented as a geometry image. Each of these objects is located inside a unit cube, bound by the values (0,1) on all coordinate axes. We place a√sphere around each object with centre (0.5,0.5,0.5) and radius 22 so that the centre of the sphere corresponds to the centre of the unit cube, and the unit cube is bound by the sphere. Any cloth particle located inside this sphere will be tested for collision. In a pre-processing phase all vertex coordinates are converted into polar coordinates relative to the centre of the sphere. The two angular components of the polar coordinates serve as an index into an n×n array, which will be used for proximity queries. According to their polar coordinates, each vertex of the geometry image is inserted into the array. In the case that two or more entries are placed into the same position in the array, they are sorted according to their distance from the centre of the sphere. Additionally, the interior of each triangle is rasterized into the array so that no gaps appear. The collision queries are done per vertex of the cloth’s surface. After the vertex to be checked has been transformed into objectspace, if the vertex is located outside the bounding sphere, there cannot be a collision. As was done in the pre-processing phase, the polar coordinate of each vertex of the cloth is calculated and used to index a point in the array. We find the closest node on the surface of the object xs and its normal signal ns by inspecting the distance part of the polar coordinates of the object and the cloth. If the sign of ns · (xs − xi,j ) is larger than −, the vertex is on the wrong side of the object, and must be pushed back to the other side of the object. The position update of the cloth is formulated as tn +∆t xi,j = xti,jn +∆t + (ns · (xs − xi,j ) + ) ns .

By adding  to the collision response term, the cloth is pulled back past the object, minimising the occurence of sharp parts of the object passing through the cloth. In our collision detection scheme, cloth self-intersection is not dealt with.

6

Results

We implemented our cloth model on a Packard Bell Easynote R9502 laptop with a Intel Centrino 2.0 Ghz processor, 1 gigabyte of memory and a Nvidia 6200 Go graphics card. Our simulations

“ ”−1 Precompute Wl = In2 − ∆t Hl m “ l ”−1 g Precompute Wg = In2g − ∆t H m At each time step dt For every particle i, j For every particle k, l connected to i, j d = xk,l − xi,j d˙ = vk,l − vi,j “ ” d ˙ fi,j = fi,j + ki,j,k,l (kdk − ri,j,k,l ) kdk + d∆t For every particle i, j Pnl −1 Pnl −1 l l fi,j = k=0 l=0 Wwi,j ,wk,l fk+xmin ,l+ymin For every nn−1 nth particle i, j −1 Pgng −1 Png −1 g g l = k=0 fi,j W wi,j ,wk,l fk n−1 ,l n−1 l=0 ng −1

ng −1

Interpolate f g For every particle i, j g l t = fi,j + fi,j fi,j For every particle i, j xG = xG + xti,jn xG = xnG2 For every particle i, j t × xti,jn ∂T = ∂T + fi,j For every particle i, j` ´ tn +∆t tn t ∆t vi,j = vi,j + fi,j m tn tn +∆t tn xi,j = xi,j + vi,j ∆t For every ` particle i,´j fc = xG − xti,jn × ∂T∆t 2 xti,jn +∆t = xti,jn +∆t + fc ∆t m Repeat maxiter times For every particle i, j For every particle k, l connected to i, j n +∆t d = xtk,l − xti,jn +∆t If kdk > ri,j,k,l (1 + ) c = kdk − (ri,j,k,l (1 + )) If kdk < ri,j,k,l (1 − ) c = kdk − (ri,j,k,l (1 − )) xti,jn +∆t = xti,jn +∆t + cd n +∆t n +∆t = xtk,l + cd xtk,l tn +∆t tn +∆t cd vi,j = vi,j + ∆t tn +∆t tn +∆t cd vk,l = vk,l + ∆t For every particle i, j Find nearest point xs with normal ns If ns · (xs − xi,j ) > − xti,jn +∆t = xti,jn +∆t + (ns · (xs − xi,j ) + ) ns xtn = xtn +∆t tn = tn + ∆t Figure 4: Algorithm Pseudo-code

Cloth Size 11 × 11 21 × 21

Desbrun et al. [1999] n Precompute Integration W 121 1.591 0.00414493 441 366.508 0.0483582

31 × 31

961

12609.9

0.182559

41 × 41

1681

−

0.24063

51 × 51

2601

−

0.58068

61 × 61

3721

−

1.0242

71 × 71

5041

−

1.93297

81 × 81

6561

−

50.5037

Our method ng nl Precompute Wg Wl 6 5 0.026 0.008 6 9 0.026 0.379 11 5 1.591 0.008 7 11 0.071 1.588 11 7 1.602 0.069 16 5 33.58 0.008 9 11 0.379 1.576 11 9 1.601 0.383 21 5 364.699 0.009 11 11 1.593 1.578 26 5 2166.05 0.008 16 9 54.266 0.631 21 7 366.654 0.07 31 5 12609.9 0.009 15 11 18.518 1.611 36 5 − 0.01 17 11 56.439 1.59 21 9 367.489 0.384 41 5 − 0.01

Integration 0.00333784 0.0146533 0.0154459 0.0367027 0.0355135 0.0405753 0.0676944 0.0661781 0.0994028 0.111589 0.191877 0.172863 0.221257 0.345918 0.251877 0.578284 0.3575 0.409644 0.932329

Figure 5: Algorithm runtime comparison in seconds.

consist of cloth free-falling under the influence of gravity, until finally draping over an object present in the scene. Our implementation achieves interactive frame rates for cloth meshes containing more than 6000 particles and yields visually attractive and interesting results. A performance comparison between our method and that of Desbrun et al. [1999] is presented in figure 5. Notable is the difference in the amount of time taken to precompute the matrices, as the method proposed by Desbrun et al. [1999] requires over 3 hours (12609.9 seconds) to calculate the 961 × 961 matrix needed for a 31 × 31 cloth mesh. We did not attempt to calculate any matrices larger than that, hence the entries are left blank. Instead, we filled large matrices with arbitrary values, so that the integration time taken per frame could still be measured. The integration time was measured as the average computational time needed per frame to simulate a piece of cloth free-falling under gravity for 60 frames of animation. From this table it can be seen that the performance of our method is a vast improvement, both in terms of efficiency, and in terms of the scalability as the number of particles increases. The main drawback with our model is that extra parameters are needed that are not physically motivated, such as the values of nl and ng , as well as the stiffness values used when precomputing the matrices. As a result, the fine-tuning of the simulation parameters needed by our model is less intuitive when compared to that of the original method [Desbrun et al. 1999]. Some screenshots from our implementation of cloth draping over various objects are shown in figure 6.

7

Conclusions

In this paper we improved the performance of an existing cloth simulation method. Utilising the regular structure of geometry images, we used a multilevel scheme to propagate forces through the cloth. By changing the way in which forces are distributed, we achieved a significant performance gain over the existing method. Since the size of the matrices used is far smaller, both the required precomputation time as well as the integration time per frame are reduced. Even though the fine-tuning of the simulation parameters becomes less intuitive, our method produces attractive draping behaviour of dense cloth meshes at interactive frame rates.

References BARAFF , D., AND W ITKIN , A. 1998. Large Steps in Cloth Simulation. In Proceedings of ACM SIGGRAPH 98, ACM Press / ACM SIGGRAPH, New York, M. Cohen, Ed., Computer Graphics Proceedings, Annual Conference Series, ACM, 43–54. B ENE Sˇ , B., AND V ILLANUEVA , N. G. 2005. GI-COLLIDE: Collision Detection with Geometry Images. In Proceedings of the 21st Spring Conference on Computer Graphics, ACM Press, New York, NY, USA, 95–102,218. C HEN , W., WAN , H., Z HANG , H., BAO , H., AND P ENG , Q. 2004. Interactive Collision Detection for Complex and Deformable Models Using Programmable Graphics Hardware. In Proceedings of the ACM Symposium on Virtual Reality Software and Technology 2004, ACM Press, New York, NY, USA, 10–15. C OHEN , J. D., L IN , M. C., M ANOCHA , D., AND P ONAMGI , M. 1995. I-COLLIDE: An Interactive and Exact Collision Detection System for Large-Scale Environments. In Proceedings of the 1995 Symposium on Interactive 3D Graphics, ACM Press, New York, NY, USA, 189–196. ¨ D ESBRUN , M., S CHR ODER , P., AND BARR , A. 1999. Interactive Animation of Structured Deformable Objects. In Proceedings of Graphics Interface 1999, Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1–8. E ISCHEN , J. W., D ENG , S., AND C LAPP, T. G. 1996. FiniteElement Modeling and Control of Flexible Fabric Parts. IEEE Computer Graphics and Applications 16, 5, 71–80. F LOATER , M. S. 1997. Parametrization and Smooth Approximation of Surface Triangulations. Computer Aided Geometric Design 14, 3, 231–250. G OTTSCHALK , S., L IN , M. C., AND M ANOCHA , D. 1996. OBBTree: A Hierarchical Structure for Rapid Interference Detection. In Proceedings of ACM SIGGRAPH 96, ACM Press / ACM SIGGRAPH, New York, H. Rushmeier, Ed., Computer

Figure 6: Snapshots taken from our cloth simulator.

Graphics Proceedings, Annual Conference Series, ACM, 171– 180.

Proceedings of Graphics Interface 1995, Canadian ComputerHuman Communications Society, 147–154.

G U , X., G ORTLER , S. J., AND H OPPE , H. 2002. Geometry Images. ACM Transactions on Graphics 21, 3, 355–361.

S ANDER , P. V., S NYDER , J., G ORTLER , S. J., AND H OPPE , H. 2001. Texture Mapping Progressive Meshes. In Proceedings of ACM SIGGRAPH 2001, ACM Press / ACM SIGGRAPH, New York, E. Fiume, Ed., Computer Graphics Proceedings, Annual Conference Series, ACM, 409–416.

H UDSON , T. C., L IN , M. C., C OHEN , J., G OTTSCHALK , S., AND M ANOCHA , D. 1997. V-COLLIDE: Accelerated Collision Detection for VRML. In Proceedings of the Second Symposium on Virtual Reality Modeling Language, ACM Press, New York, NY, USA, 117–124. K LOSOWSKI , J. T., H ELD , M., M ITCHELL , J. S. B., S OWIZRAL , H., AND Z IKAN , K. 1998. Efficient Collision Detection Using Bounding Volume Hierarchies of k-DOPs. IEEE Transactions on Visualization and Computer Graphics 4, 1, 21–36. L IN , M., AND G OTTSCHALK , S. 1998. Collision Detection between Geometric Models: A Survey. In Proceedings of the 8th IMA Mathematics of Surfaces Conference, Information Geometers, Ltd., Winchester, UK, 37–56. M ANOCHA , D. 2005. Quick-CULLIDE: Fast Inter- and Intra-Object Collision Culling Using Graphics Hardware. In Proceedings of the IEEE Virtual Reality Conference 2005, IEEE Computer Society, Washington, DC, USA, 59–66, 319. P ROTOPSALTOU , D., L UIBLE , C., A REVALO , M., AND M AGNENAT-T HALMANN , N. 2002. A Body and Garment Creation Method for an Internet Based Virtual Fitting Room. In Proceedings of Computer Graphics International 2002, Springer-Verlag, 105–122. P ROVOT, X. 1995. Deformation Constraints in a MassSpring Model to Describe Rigid Cloth Behavior. In

S ANDER , P. V., W OOD , Z. J., G ORTLER , S. J., S NYDER , J., AND H OPPE , H. 2003. Multi-Chart Geometry Images. In Proceedings of the 2003 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing, Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 146–155. S URAZHSKY, V., S URAZHSKY, T., K IRSANOV, D., G ORTLER , S. J., AND H OPPE , H. 2005. Fast Exact and Approximate Geodesics on Meshes. ACM Transactions on Graphics 24, 3, 553–560. T ESCHNER , M., H EIDELBERGER , B., M UELLER , M., P OMERANETS , D., AND G ROSS , M. 2003. Optimized Spatial Hashing for Collision Detection of Deformable Objects. In Proceedings of Vision, Modeling, and Visualization 2003, 47–54. T ESCHNER , M., K IMMERLE , S., H EIDELBERGER , B., Z ACHMANN , G., R AGHUPATHI , L., F UHRMANN , A., C ANI , M. P., FAURE , F., M AGNENAT-T HALMANN , N., S TRASSER , W., AND VOLINO , P. 2005. Collision Detection for Deformable Objects. Computer Graphics Forum 24, 1, 61–81. YOSHIZAWA , S., B ELYAEV, A., AND S EIDEL , H.-P. 2004. A Fast and Simple Stretch-Minimizing Mesh Parameterization. In Proceedings of Shape Modeling International 2004, IEEE Computer Society, Washington, DC, USA, 200–208.