Real-Time Simulation of Rigid Body Dynamics

Final Project: CSE 557

Charles Gordon



Overview

A great deal of work has been done on creating realistic simulations of rigid body dynamics at frame rates. This project attempts to combine the results of a number of different, state of the art algorithms for modelling realistic physics and collisions. The goal was to create realistic approximations to real-world physics while maintaining at least the 30 frames/second needed for animation purposes.

References

Chris Hecker has done a series of articles for Game Developer magazine and put up a wonderful web page detailing his results and the library-full of references he used.

One of those references is to David Baraff's set of papers on physicaly based modelling, which were used as a tutorial for SIGGRAPH '95. These include a simple algorithm for collision detection and response which turned out to be monstrously difficult to implement (see my notes below).

The collision detection algorithm used in my program is the same one used for Q-COLLIDE. This algorithm combines some interesting results from computation geometry with some (fairly difficult to understand) algorithms for finding closest points on arbitrary convex polymeshes. Kelvin Leung's master's thesis on QCOLLIDE can be found here.

Details

The most common cause of aliasing problems in physics modelling is the numerical integrator. The easiest way to approximate the closed form solutions of classical physics is to use an Euler integration step, which is really just a first order Taylor polynomial. In order to increase the accuracy of the integrator and to also increase the size of the time step allowable, I used an adaptive midpoint approximation, which is a second order Taylor approximation that detects when the error is beyond some tolerance and subdivides the time step (this is a surprisingly simple operation). Where the Euler integrator required time steps on the order of hundredths of a second to be anywhere near accurate, the midpoint evaluator uses time steps right around the 1/30th of a second necessary for real-time animation.

The program simulates both linear and angular dynamics. Quaternions are used to specify the orientation of the bodies in space (and also because they are much easier to handle with the integrator). Initial conditions can be set and the bodies follow all of the normal constraints of gravity. Special objects I called "attractors" can be placed in scenes to act like small planetary bodies, giving out an attractive force that is inversely proportional to the distance squared and directly proportional to the mass of the object. These can be set to any multiple of the gravitational constant so that they either act as strong gravity wells or, if the scalar is negative, as strongly repulsive forces.

Collision detection turned out to be an immensely more difficult problem than I had anticipated. This is especially true of trying to make it accurate in real-time, as the calculations involved can be very time consuming. In order to reduce the number of calculations there must be a simple, high granularity step using bounding boxes or spheres. When the bounding boxes of two objects collide there is another step to determine if the objects themselves are actually in contact. This is followed by a determination of exactly which points on each object are in contact and then the calculation of an impulse force to be applied there. The algorithms I used for each of these steps is summarized below.

Bounding Boxes: The first step in the CD algorithm is to put axis-aligned bounding boxes around all of the objects in the scene. To avoid the overhead of determining the maximum extent of each object in three-space, each object also had an object-aligned bounding box that underwent the same rotations and translations as the object. The axis aligned bounding box was determined from the eight points of this object aligned box. To determine collisions, the axis-aligned boxes were intially sorted along each axis. At every time step the coordinate lists were traversed and resorted using a simple insertion sort (so that only objects that had actually changed position were considered). Any objects that had overlapping boxes were added to a list of objects to be checked in the next phase of collision detection. The bounding box algorithm required a small overhead at every time step to keep the object-aligned boxes in the correct orientation, followed by an expected linear time to re-sort the axis aligned boxes.

Separating Vectors: The separating plane algorithm used by Baraff can take up to O(mn) time to execute, where m and n are the number of vertices of the colliding polymeshes. A better algorithm using a separating vector is used by the QCOLLIDE library and takes time O(m+n). By caching the vector (or plane) between time steps this can be reduced to nearly constant time. In this program a separating vector is found between polymeshes everytime their bounding boxes overlap. This vector can be found very quickly by using incremental updates to its position based on the results of a set of dot products at the last position. It is detailed in the thesis mentioned above. Unfortunately, due to time constraints, I was unable to implement the caching algorithm, so this step always takes O(m+n) time every time bounding boxes are found to intersect.

Closest Features: This is by far the most difficult part of collision detection and it took the bulk of my time (this is an immense understatement, I have over 30 hours into this). The first approach I tried was an O(mn) algorithm that simply compared the distance between every pair of features on both polymeshes until it found the closest pair. This approach is neither effective nor fast. It could only find that, say, a point intersected a face, but not where on the face they intersected. This led to very unrealistic looking collisions between objects with large faces (such as squares). Although it worked well for spheres with a large number of small faces. It also ran very slowly, bringing the simulation to its knees. The next approach I tried used the results from the previous step along with an internal graph of neighboring vertices for each vertex in the polymesh. This allowed a much faster lookup of the colliding features, but did not fix the problem of not knowing exactly which points were colliding. It was at this point that I came across a description of the Gilber-Johnson-Keerthi algorithm for collision detection. It is explained in the master's thesis mentioned above as well as in:

Bergen, Gino Van Den. A Fast and Robust GJK Implementation for Collision Detection of Convex Objects. Department of Mathematics and Computer Science, Eindhoven University of Technology.

This algorithm does an amazing job of finding the closest points between any two arbitray polymeshes in expected linear time (using temporal coherence and caching). The basic idea is to do a pairwise subtraction of all the points in both polymeshes and then find the points that are closest to the origin. Rather than actually performing the pairwise subtraction, they devised a method for incrementally descending to the correct pair of points starting from the convex hull of any four points on the surface. The Johnson algorithm, which is used as a subalgorithm, finds the closest point from these four points to the origin. This point is used to find a new point on the original surface and the process is repeated. Empirical results in the paper (and repeated in my program) show that the algorithm generally terminates after 4 to 6 of these iterations. With some selected use of caching this can be reduced even further.

The results of using this algorithm in my project leave a lot to be desired. Besides the mathematical difficulty of the method, very little was said in either of the papers I found about parameter tuning, which is the most difficult part to making it work correctly. The GJK algorithm includes a number of small error thresholds that must be tuned very precisely to insure that the algorithm terminates correctly. After 15 hours of watching it go into infinite loops I can assure you that this is no triviality. I am fairly confident that my current implementation is correct, but it does not interface with the collision response algorithm very well, which means that objects will occasionally get "glued together" because the physics engine can't generate the correct force to get them unstuck. This is a fairly difficult problem, but not unsolvable. In the meantime I have fudged a small repulsive force to guarantee that they do not stick. This makes the result less than realistic, but it was all I had time for.