Highlights

• Arbitrary model loading and handling, including both convex and non-convex triangle models
• Support for multiple bodies
• Implementation and comparison of incremental and linear axis-aligned bounding box update for use in collision detection between objects and surrounding world cube
• Four types of body-to-body collision detection algorithms of varying degrees of accuracy
• Implementation and comparison of five different ODE solvers
• Simplified rotational dynamics
• Full-featured interactive graphical simulation
• Fully configurable simulation script files
• Compiles and runs on both SGIs and PCs without modification using OpenGL and GLUT
• All routines (except for basic OpenGL viewer) were written specifically for this project; no collision detection or ODE solver libraries were used (i.e.: I-COLLIDE, RAPID, and Numerical Recipes code)
• Implemented a full particle simulation for further analysis of ODE solvers and the force approximations.

• Demonstration and Screen Shots of interactive simulation running on an SGI Infinite Reality.

Source Files

• Source Code and PC Executable
• Input Script Files (bodies.dat)
• Convex Triangle Models (*.tri)

Project Highlight Descriptions

Model Handling

• Loads arbitrary "triangle soup" models with per-triangle information
• Reads in models' output from a public domain modeler called AC3D
• Converts models to shared-vertex format with vertex neighbor information for use in the incremental bounding box update routine

Collision Detection with Surrounding World Cube

• This system supports both linear O(n) and incremental O(1) axis-aligned bounding box update. The problem of finding closest features has been reduced to finding the 6 extreme points that form a tight-fitting axis-aligned bounding box. The AABB approach greatly simplifies collision detection by reducing the overlap test to a simple min-max test. This technique seems to be more general since it can be used in more sophisticated sweep-and-prune collision detection systems (like I-COLLIDE). See previous assignment for details on the incremental algorithm.
• The linear and incremental approaches were nearly equal in performance until the size of the models became quite large. The incremental approach performed noticeably faster on simulations containing 2000 triangles each.
• Non-conves objects are supported in this simulation by using the brute-force linear AABB update; an incremental approach requires a convex decomposition of the model, making it prohibitively complex in theory and practice.

Body-to-Body Collision Detection

• Using a simple, 3-level hierarchical approach, the system provides 3 levels of accuracy in body-to-body collision detection to trade accuracy for speed.
• The 3-level hierarchical approach proceeds as follows:
• Test 1: Pairwise O(n²) fast axis-aligned object bounding box overlap test (n being the number of objects, typically quite small, 1 to 40 in this simulation). For the coarsest and fastest level of collision detection, these are the only tests.
• Test 2: For objects overlapping in Test 1, we then perform pairwise testing of the triangle AABBs. The second level of accuracy stops here.
• Test 3: Pairs of triangles passing Test 2 are then tested for exact overlap using edge-triangle intersection tests.
• This system has the advantage of being conceptually simple, easy to implement, easily configurable for performace tweeking (oder of operations, additional levels of hierarchy, etc.), quite accurate, and robust.
• Numerical issues do arise concerning intersecting nearly parallel edges and plains. This did not seem to be a problem in this simulation. Anyway, using a more robust technique based on Stefan Gottschalk 's OBBtree separating axis theorem could remove this problem.

Forces Used in the Simulation

• Gravity: always directed downward with respect to the viewer.
• Viscous Drag: force acting opposite the direction an object is travelling. This force is approximated by negative linear scaling term on the object's velocity. Some variations use a velocity squared term in determining the drag force; this was tried and did not yield any significant results to justify the extra computational expense.
• Restitution: dampening forces due to elastic collisions. These forces are required to make the system more stable and realistic; purely inelastic collisions appear too "bouncy". This dampening force is approximated by a linear scaling term on the normal component of the reflected velocity.
• Lateral Friction: this force dampens the lateral components of the reflected velocity. This is approximated by a simple linear scaling factor on the lateral components of the reflected velocity. This force should normally be proportional to the normal component of the incident velocity.

Differential Equation Solvers

• Euler: The step, or increment, to the next state vector is equal to the derivative of the current step vector. This is a first-order method so global truncation error is proportional to O(h), where h is the step size.
• Two-Step Averaged Euler: The step for the current state vector is equal to its derivative. The next state vector is computed using an Euler step, and the step for this state vector is also computed. The average of these two steps is used as the step for the original state vector. Global truncation error is proportional to O(h²).
• Backwards Euler: Two Euler steps are performed and the second step is used as a step for the original state vector. Global truncation error is proportional to O(h²).
• Midpoint: One half of an Euler step is performed, the step is computed at this midpoint, and used as a step from the original state vector. Global truncation error is proportional to O(h²).
• 4th-Order Runge-Kutta:
• The step vector is computed as follows:
• Step1 is computed as a normal Euler step
• Step2 is computed as the derivative of the state vector, computed by travelling half of step 1 from the original state vector
• Step3 is computed as the derivative of the state vector, computed by travelling half of step 2 from the original state vector.
• Step4 is computed as the derivative of the state vector, computed by travelling a full step 3 from the original state vector.
• The resulting step is computed as a linear combination of steps 1-4 as follows: one-sixth of step1 + two-thirds of step2 + two-thirds of step3 + one-sixth of step4.
• The global truncation error is proportional to O(h⁴).
• The basic Euler method proved to be unstable without the presence of sufficient damping forces; objects seemed to "jitter" when in resting contact.
• The remaining four solvers behaved very similarly in terms of stability and accuracy.
• The 4th-Order Runge-Kutta method proved to be unjustifiably computationally expensive for use in a particle system where the number of bodies tends to be much higher. Nevertheless, there was no noticeable difference between this method and the O(h²) methods.
• The best choice in general seems to be the Backwards Euler method, since it is second only to the Euler method in efficiency, and proved to be quite stable and accurate in our simulation.

Translational Collision Response with Walls of Axis-Aligned Surrounding World Cube

Using the current axis-aligned bounding box that surrounds the object, we simply handle each component (or dimension) separately as a one dimensional interval overlap problem. If the box extends beyond the world cube in a particular dimension, we push it back to the wall only along that axis and flip the velocity; restitution and frictional forces are applied here as a dampening on the normal and lateral components of the flipped velocity respectively.

Translational Collision Response Between Bodies

Because our collision detection system only returns contact pairs and rotational dynamics was not implemented for body-to-body collision (because of time), we simply query the collision system for overlap only. We then create a contact plane using the vector formed between the object's centroids as the normal to the plane. All exchange and reflection of normals occurred about this plane. Perfectly inelastic collision responses were performed; so, if two objects are moving towards each other during contact, they exchange velocities (the objects are assumed to have equal mass).

Simplified Rotational Dynamics

Each body in our simulation has the following:

• Position in world space (Pos): 3D vector
• Linear Velocity (Vel): 3D vector
• Angular Velocity (w): 3D vector that forms the axis about which to rotation. the magnitude is the rotation rate in revolutions / time.
• A set of orientation axes (X,Y,Z): normalized and orthogonal. A 3x3 with the axes as columns forms the orientation matrix to properly orient an object in world space.

For each frame of animation, the position and velocity is updated in the normal fashion completely independent of rotational effects. The orientation axes of the object are updated as follows:

X += (w x X)
Y += (w x Y)
Z += (w x Z)

Where the cross product between the angular velocity and an axis is the velocity for the vector. This is an Euler step for the orientation. to take collision into account we need only properly modify the angular velocity (w).

We have effectively factored out angular momentum and the inertia tensor by assuming the mass is equally distributed about the centroid (a uniform density sphere). The inertia tensor then becomes the identity matrix, and the angular momentum, which is the inertia tensor times the angular velocity, just becomes the angular velocity. The benefits of this will become apparent below.

In order to model the change in the angular velocity, we need the following:

• vertex point of contact (Pi) in the object space with the centroid of the object at the origin. The average position of all vertices should be the centroid (assuming equal masses for each vertex). It is important to note, that Pi is in the object space, but it must still be properly oriented. We can obtain this by simply applying the orientation rotation only to the point of contact.
• velocity (Vi) of the point of contact
• resulting velocity (Rvi) of the point of contact after reflection off of contact plane
• instantaneous impulse "force" (I) exerted on the object at the point of contact that will cause the change in rotation
• torque (T) that is the change in angular momentum, but in our case it is the change in angular velocity directly since the inertia tensor is assumed to be the identity matrix.

Now given a point of contact Pi (which is just a vertex of the object), we can compute the velocity of the particle (Vi) as follows:

Vi = CrossProd(Angular Velocity, Point of contact) + linear velocity of the body
Vi = (w x Pi) + Vel

The resulting velocity (Rvi) is then computed by simply treating the point of contact like a particle and reflecting its velocity about the normal of the contact plane. Again, normal and lateral components of the reflection can be dampened to account for restitution and lateral friction.

Normal Component of Rvi = -DOT(Vi,N) * N
Lateral Component of Rvi = Vi + Normal Component of Rvi
Resulting Rvi = Lateral Component + Normal Component
Resulting Rvi w/ dampening = (Lateral Component * kFriction) + (Normal Component * kRestitution)

The impulse (I) is simply the dampened normal component of the resulting velocity (Rvi) which constututes an instantaneous force applied perpendicularly to the contact plane into the contacting body.

I = Normal Component of Rvi * kRestitution

The torque (T) is simply computed as follows:

T = CrossProd(Point of Contact, Impulse)
T = Pi x I

The torque IS the desired change in the angular velocity (w), so the angular velocity is immediately changed with an Euler step as follows:

New Angular Velocity = Old Angular Velocity + Torque
w += Torque

Existing Problems (remaining bugs in implementation)

• Collision response between the bodies is sometimes incorrect. Since we exchange velocities upon colliding contact, sometimes objects will interpenetrate when they are moving slowly or in the same direction.
• Rotational dynamics is wrong despite the correctness of the simplification derivation. I have carefully drawn all of the vectors involved and they seem to be correct, but the angular velocity does not properly dampen and come to a rest.
• Collision detection system could use more levels of hierarchy to make it more efficient.

Summary

I have thoroughly completed the three minimum guidelines to this project, and have sufficiently justified majority of all three of the extra credit requirements. Despite the presence of bugs in both the N-body collision response and the rotational dynamics, both of these requirements were thoroughly covered. I feel that I have more than compensated for these shortcomings with a complete implementation of a multi-level of accuracy collision detection system and with a compelling interactive fully-configurable simulation.