Geometric Foundations for Graphics and Visualization

Physics-Based Modeling, Simulation, and Computing

CSE 530 + CSE631

Final Project

by Evan Closson

5/19/2003

- This in word: collision_detection_and_response.doc
- Presentation: collision_detection_and_response.ppt
- Program: program.zip

Collision detection and the responses after collision are important elements in the simulation of virtual environments. When working with large polygonal models in a scene where a "real-time" simulation is paramount, it is important to quickly determine collisions and their responses.

This final report discusses different techniques used to speed up rigid body collision detection with bounding volume hierarchies and the details involved in calculating accurate responses to the collisions. The theoretical and research details are discussed as well as the implementation details encountered while working with the project.

Collision detection and response are crucial, elements when creating a realistic simulated environment. When working with real-time simulation constraints, the speed with which these computationally expensive operations work must be greatly improved in order to retain the correct time and still keep the realist results of the collision.

This paper focuses on speeding up collision detection and providing an accurate response for rigid bodies by using bounding volume hierarchies. There has been much work done with bounding volume hierarchies and collision response.

Many efficient bounding volume hierarchies have been created with spheres [HUBBARD96], AABBs [VAN DEN BERGEN97], k-dops also known as discrete orientation polytopes [KLOSOWSKI98], and OBB trees [GOTTSCHALK96]. Each of these techniques has a tradeoff between tight fitting bounding volumes which require less intersection tests, or looser fitting bounding volumes like spheres that require more bounding volumes to approximate a surface but require very inexpensive intersection tests.

Much research has also been done in numerical computing techniques geared for rigid and flexible body dynamics in [BARAFF89,92,94].

Collision detection is a well researched area, but it is very fundamental in many applications and having a thorough understanding of the techniques and methods used as well as the problems involved in collision detection is crucial for any simulated environment. The focus of this paper is on collision detection using bounding volume hierarchies the general algorithms for bounding volume hierarchies are explained, and two types of bounding volume structures are discussed with an emphasis on sphere trees.

A collision volume hierarchy is a way to approximate a model with an efficient data structure that can be used for making queries of the type: do these two bodies A and B collide and at what point. The basic design that bounding volume hierarchies follow is to approximate a model at different levels of detail. Each level of detail is a different depth in a search tree. The top level - or root level - has the lowest approximation and is used to quickly determine that two bodies are not in contact. If there is a contact the tree is traversed and the points of contact are returned.

Hierarchies of axis aligned bounding boxes can be used to create collision detection trees. These bounding volume structures make use of binary tree's to decompose a model. Axis aligned bounding hierarchies can also handle deformable models of most types well [VAN DEN BERGEN97].

Instead of using the AABBs to approximate the object it is actually used instead to quickly find the primitives of the object. The tree is created by first creating a bounding box around the entire object. Then a partitioning plane is chosen with the emphasis of keeping the bonding boxes as `cube like' as possible. This is done because the more `cube like' the bounding boxes the less overlap they will contain. The underlying primitives - usually triangles - are split into two sets by the partitioning plane. The portioning step creates two new children for the parent node each one with its respective new set of primitives. Now both of the new children place a tight fitting AABB around the primitives it is associated with. Each child now repeats this partitioning step repeatedly until every node contains only one primitive.

For AABBs the size of the tree structure is dependent upon the geometry of the model. If the model has n primitives then there will be n leaf nodes. The benefit of this however is that it makes it easier to rebuild parts of the tree when the underlying geometry gets deformed.

The AABB hierarchy is then traversed and checks for collisions by finding overlap between two AABBs under consideration. If there is overlap their children are traversed. When a two leaf nodes are considered their primitives are checked against each-other the result of the test will be whether or not an intersection has occurred and at what point.

Sphere trees are very efficient bounding volume hierarchies. What they lack in their ability to easily conform to arbitrary objects they make up for in their simplistic collision calculations. Presented here are two methods for creating bounding sphere hierarchies, one that makes use of octrees, and another that uses a medial axis approximation.

Creating sphere trees with the use of octrees is probably the simplest creation method.

First create a bounding box enclosing the entire model. Split this box with three evenly divided planes to create eight child boxes within the top box. Repeat this on each child until the desired depth is reached.

Once the octree is created we need to traverse the octree and create a sphere at each node that encompasses all of the geometry that is contained inside that bounding box. There are different methods to do this and they generally depend upon the type of sampling that is done for the model. One simple approach is to sample the model uniformly and create points that are inside of it. There are many other methods some that don't involve sampling but instead make use of collision detection on the geometry to find the sections of the node that represent the model.

The medial axis is defined as the Voronoi diagram of the edges of a polygon, which are interior to that polygon. Note: the medial axis can actually be computed in linear time with respect to the number of edges of a polygon [BERG00] it would be interesting to see Hubbards algorithm [HUBBARD94, 96] updated to make use of this. Here however an approximation is used based upon the Voronoi diagrams of the vertices.

Using the medial axis approximation method is no trivial manner and for detailed information one should see [HUBBARD94, 96]. A brief overview is given here.

The idea is to fit spheres to the medial axis - or skeleton - of a model. These spheres will more tightly fit the structure because their placement depends upon the geometric structure of the model.

To create the sphere tree, first a Voronoi diagram of the model is computed, this is then used to approximate the medial axis. The lowest level of spheres is then created using the interior Voronoi vertices. Finally the sphere tree is created from these lowest level vertices by a root node and then merging sections of the lowest level of spheres together to obtain the desired number of children for the root node, along the way the rest of the tree is build up. The spheres are merged through adjacency info that they have access to through the Voronoi diagram.

One of the benefits for using sphere trees are the very simple and efficient calculations involved when checking for collisions between two spheres. Figure 2.1.2.3-1 demonstrates the calculations.

These calculations are simple enough and can be sped up even more by not taking the square root of the distance. Figure 2.1.2.3-2 shows code checking for collisions.

This code checks to see if distance between the two centers of the spheres is less then the sum of the two radii. If so a collision has occurred.

Since the motion of the bodies when running a simulation is estimated by minute time steps we are not dealing with continuous time, but time that is sampled. This sampled time can make the collision be detected too late, as in after a collision has already occurred.

To get around this sampled time problem different methods can be used. First one of the obvious things to try is to decrease the time step used in the simulation. This will increase the accuracy of the simulation causing the exact time of collision to be more accurate. There is a large computational tradeoff here though. The more time steps there are the more computing power is need and that power might not be available.

Another method would be to estimate the point of collision to some degree of accuracy specified by the user using a binary search type approach. If at simulation time ti + 1 a collision is detected we move the simulation time back to ti +.5 and test again. If a collision is again detected try between ti and ti +.5 else try between ti +.5 and ti + 1. At the limit this will find the exact point of collision between ti and ti + 1, so the more we subdivide the more accurate the result. This method is better then the first presented since the extra time steps only occur when a collision occurs, which may not be often.

There are other methods that involve predicting where bodies will collide ahead of time; these are presented in [BARAFF89].

Detecting when a collision occurs is not all that needs to be done, the point of contact must be computed in order to pass this information onto the collision response system.

For sphere trees when leaf-leaf spheres are detected in a collision this becomes a collision point. There can be many collision points at a time but to simplify the collision response system it is easier to have only one point of collision. To do this, all the colliding spheres are averaged to estimate the collision point.

Once the collision points for both colliding bodies are known the line of action - normal to each surface since we are working with spheres - between them must be computed. This gives the direction along which the collision has occurred. It can be found by simply finding the normalized vector from the distance between the points of contact on both A and B.

So now we pass the collision normal, the collision point on both A and B in their local space, and the relative velocity between the two points onto the collision response system.

The tree hierarchy that I used made use of spheres. I felt that the simplicity of the collision tests would greatly help cut down on the amount of time checking for collisions even though it has a generous boundary for most shapes.

In order to create the sphere tree's I used two methods. First an auxiliary program was created to make simple sphere tree hierarchies. These were not very tight fitting approximations of the object but did the job of being able to get the hierarchy in place and tested. The second method was making use of [BRADSHAW03]'s program that can create sphere trees. I wrote converters to go back and forth between the file types I used and the file types used in that program. This proved very beneficial for I had the use of many different types of methods for the creation of sphere tree's that could then be use to test out how efficient different hierarchical sphere structures were. The method for creating sphere trees that I felt performed the most efficiently during collision detection was the medial axis method [HUBBARD94, 96].

The results that were obtained for the collision detection were adequate. I learned that the problem of collision detection is a very hard one to get looking `just right' when trying to keep speed optimal. It is for this reason that penetration and computing the exact time of collision is not as accurate as I would like them to be. Initially when a collision was detected the last time step and the one in which a collision was detected would be cut in half multiple times in a 'binary search' type of approach to find a more accurate time of collision. The results of this were slower then expected so this was eventually removed - a faster tree hierarchy was needed. Asides from this, the collision detection looks ok.

Another problem that was also encountered involved the collision detection results that get passed onto the collision response system. Many collision points - each leaf sphere - can be found to be in contact for a collision. The collision response method I used could use only one point of contact. If many were recorded the impulse force calculation is dependent upon the order in which the collisions were detected. To approximate the collision point found, all the colliding spheres were averaged together to create one collision point approximation, this had rather good results. There are ways to more accurately account for multiple collision points during impulse calculations [BARAFF89] but the results I received by averaging the colliding points were rather good.

Pictures of the sphere trees and their use during collision detection can be seen in Figures 1 and 2 at the end of this paper.

So much work has been done in collision detection hierarchies that it has become a very saturated field, however with the steadily increasing availability of fast consumer grade graphics hardware more research is being done in the area of hardware assisted collision detection. Using graphics hardware collision detection can be computed much faster on the graphics hardware then on regular software and this is one direction of research that collision detection is starting to take. Also more work on simulating deformable model collision detection with updatable bounding volume hierarchies is being done and although deformable model collision detection has also been researched a lot there is still room for improvement.

Simply put, what happens between two colliding body's is a collision response. There are different ways to handle collision response, this project used impulse based collision response, but there are other types of collision response such as penalty methods that inserts a repellent force between bodies so that as they collide the force pushes them apart again.

To simulate an accurate collision response one must first understand the mass properties of a body, the equations of motion and how to compute them, and how to compute the impulse force from collision information determined during collision detection.

Mass, center of mass, and moment of inertia make up the mass properties of a body. These mass properties are intimately related to the study of a body's motion. The responses that a body has to the forces that are exerted upon it are dependent upon its mass properties for they determine a body's resistance to motion.

Simply put, mass can be defined as the total number of elemental particles in a body. The mass can be found by taking the integral over the volume of the body while taking into consideration the density throughout the body as well:

M = ? ?dV

*Equation 5.1-1:* The integral of the volume to calculate the mass of a body where ? is the density and V is the volume.

Using the integral over the volume is not needed in order to find the mass of a body. The body can be broken up into masspoints, and these points can be summed to compute the total mass:

M = ä mi

*Equation 5.1-2:* The calculation of total mass given i finite mass points.

Usually a masspoint is just a vertex in the mesh of a model with an extra mass value associated with it.

To relate mass to the motion of a body, mass can be thought of as the resistance that a body has to linear motion. Looking at the formula for acceleration in Equation 5.1-3 it becomes obvious that the greater the mass of a body the larger a force that is needed to accelerated the body faster.

a = F/m

*Equation 5.1-3:* The formulation for acceleration where a is the acceleration of a body F is the force acting upon the body and m is the mass of the body.

The center of mass for a body is a point about which the body's mass is evenly distributed. The center of mass has some interesting properties of note in body dynamics: if a force acts directly upon the center of mass the body will have no rotational movement and if a force acts off-center of the center of mass then there is a rotation of the body about its center of mass. Center of mass is formally defined over a continuous body in Equation 5.1-4.

COM = ? xyz dm

*Equation 5.1-4:* This is the equation to find the center of mass over a continuous body. Here xyz is a position vector and so is COM.

As before when calculating the mass of a body, the center of mass can also be approximated using masspoints:

COM = (ä pimi)/M

*Equation 5.1-5:* The calculation of the center of mass for a body where pi is the location of masspoint mi and M is the total mass of the body.

Inertia is defined as the resistance to the change in motion of an object. Mass is also known as inertial mass. The moment of inertia is the inertia for a point mass about a given axis. So since bodies rotate about their center of mass the moment of inertia defines the inertia mass for the points that don't lie on the center of mass. The moment of inertia over a continuous body is formally defined in Equation 5.1-6.

I = ? r^2 dm

*Equation 5.1-6:* This is the equation to find the moment of inertia for a continuous body where r is the radius from a given axis of rotation.

When working with a body in 3D the mass may not be symmetric. If this is the case the inertia needs to be represented with a second rank tensor to account for the rotational velocity's not always being parallel. This is called the inertia tensor. The inertia tensor shows the relationship of mass over the body with respect to the center of mass. This tensor which is represented by a 3x3 matrix consists of the moment of inertia values along the diagonal and the products of inertia making up the rest of the elements of the matrix.

Calculating the inertia tensor for a 3D body with a simple known shape - cylinder, sphere, or cube - can be done using simple formulas for inertia found in any introductory college level physics text book - the products of inertia for these tensors will be zero. However, to approximate the inertia tensor for an arbitrary shape we can once again break down the arbitrary shape into many smaller masspoints and assign each mass point its own local inertia. The local inertia for a masspoint can be calculated by assuming it is a simple shape and then using the formula for calculating that shape's moment of inertia.

Once the local inertias for the masspoints are known we can use the following Equation 5.1-7 to estimate the inertia tensor:

Like mass, the moment of inertia can be thought of as a resistance to motion. Unlike mass it is not a resistance to linear motion, but to a body's rotational motion.

A rigid body has two types of motion: linear motion and angular motion. A force applied to a body will make it move. A torque applied to a body will make it rotate.

F = ma

ç = r x F

ç = Ià

*Equation 5.2:* The equations for force and torque. Here r x F is the cross product with the force and r which is the distance to the axis of rotation. The can also be defined in a manner similar to the force by the moment of inertia times the angular acceleration.

ç = r x F

ç = Ià

Linear motion is a very simple concept and for a rigid body can be thought of as just the motion of the center of mass for the body. The equations for linear motion are nothing that anyone hasn't seen before and are shown below in Equation 5.2.1.

a = F/m

v = ds/dt

s = ? v dt

*Equation 5.2.1:* The equations for linear motion. The acceleration is force over mass, the velocity is the distance over the time, and the displacement is s.

v = ds/dt

s = ? v dt

Also it must be noted that in the actual simulation the position or displacement is updated from the velocity in very small time steps.

Angular motion describes the motion of the body about its center of mass. Instead of being governed by the force and mass it is controlled by the torque - derived from the force - and moment of inertia. Figure 5.2.2 shows the equations for rotational motion.

à = I^-1ç

? = dê/dt

ê = ? ? dt

*Equation 5.2.2:* The equations for angular motion. Here I^-1 is the inverse inertia tensor and ê is the orientation or angular displacement.

? = dê/dt

ê = ? ? dt

These formulas are very similar to the formulas for linear motion. Like linear motions position the current orientation of the body is updated from the angular velocity ?. The orientation is often stored as either a rotation matrix or a quaternion and is updated with time steps.

Real-time simulations are limited by the amount of time they have to solve the ordinary differential equations used in motion. This requires that quick approximation schemes for integration are needed in order estimate the solution without much computational overhead. When working with real-time simulations it is important to note that we are dealing with initial value problems, these problems can be solved incrementally over time. One of these approximation schemes - also the simplest and rather error prone - that can be used to solve the equations of motion is Euler's Method.

A Taylor series expansion can be used to approximate an initial value problem. Euler's Method is the Taylor series expansion of order one.

f(x) = ä(fn(a)/n!)(x - a)^n

*Equation 5.3-1:* Equation of the Taylor series where f(x) is centered at a, n is the order and fn(a) is the derivative of f to the n'th order.

To give a more general explanation of what Euler's Method does it helps to look at a simple example.

v1 = v0 + a*t

*Equation 5.3-2:* Example of integration using Euler's method for updating the velocity. Here v1 is the new velocity v0 is the old velocity a is the acceleration and t is the time-step.

This equation shows the Euler's Method solution for the initial value problem for when the velocity v0 is known and we want to find out the velocity for v1. Acceleration is velocity over time - or the derivative of the velocity with respect to time - and here the new velocity is approximated with a finite time-step t. Expanding the Taylor series to solve the velocity as a function of time shows us how to approximate these initial value problems.

v(ti + dt) = v(ti) + dt v'(ti) + (dt^2/2!)v''(ti) + (dt^3/3!)v'''(ti) + . + (dt^n/n!)vn(ti) + .

*Equation 5.3-3:* The expanded Taylor series solving for velocity as a function of time.

When using Euler's Method which has order one, it is easily seen how the Taylor series expansion in Equation 5.3-4 is related to the Equation 5.3-2.

v(ti + dt) = v(ti) + dt v'(ti)

*Equation 5.3-4:* The Taylor series of order one, where v'(t) is acceleration.

Euler's Method for integration is an approximation of the actual solution. This approximation introduces errors and the errors in a simulation can build up over time to cause a divergence from the actual solution. There are many ways to improve these approximations so that they are more accurate. Taking a look at Euler's Method it is obvious that if we increase the time-step that is used the solution will have less errors. However this is computationally expensive when working with real-time constraints. There are other methods for integrating that are more accurate then Euler's Method and some are briefly discussed below.

Looking at the example of finding the velocity in Equation 5.3-3 and 5.3-4 one can see that the Taylor expansion of Euler's Method is only order one. Increasing the order of the Taylor series is one way to increase the stability of the solution. This method is often used however the problem of calculating the higher order derivatives can be time consuming.

Runge-Kutta methods, which work by approximating the Taylor series methods using linear combinations to approximate their derivatives, can be used to more accurately integrate a solution. The more common Runge-Kutta method is the fourth order method which approximates the Taylor series up to order four.

There are many integration methods that can be used and references to the topic are included at the end of this paper.

Impulse is the change in momentum. The equations for momentum and impulse can be seen in Equation 5-4.

p = mv

J = ?p = Favg * ?t

*Equation 5.4:* The equations for momentum and impulse where Favg is the average force exerted over a period of time.

J = ?p = Favg * ?t

Impulses are used to determine the collision response between two bodies A and B and are calculated with the information provided by the collision detection scheme. The result of the collision detection routine that is used in this project returns one point of contact between the two objects that are colliding. In the case of multiple points colliding, the collision detection routine averaged them to find an approximate point of contact. Once the point of contact is known the normal is calculated along with the relative velocity of the colliding points.

First the linear impulse is explained but this is only suitable for a linear system. The impulse that is needed when rotations are taken into consideration is explained afterwards.

When working with arbitrarily shaped rigid bodies linear impulse alone is not enough, however point masses can be modeled with linear impulse alone and to set the stages for impulse that takes into consideration rotations it helps to first take a look at the simple case of linear impulse.

To calculate a new linear velocity based upon an impulse force we can use the following definition of the change in velocity:

?v = J/m

vi = vi-1 + J/m

*Equation 5.4.1-1:* Calculating the change in velocity from the impulse.

vi = vi-1 + J/m

Now to solve for J it is first broken up into two parts consisting of a magnitude and a direction - the normal vector between two colliding bodies. Since are taking into consideration two bodies that are colliding, once the impulse is calculated one body get the calculated impulse and the other gets its reverse.

J = jn

*Equation 5.4.1-2:* Calculating the impulse based upon the magnitude and the normal direction.

To solve for j we must look at the equation for a frictionless collision, or reflection shown in Equation 5.4.1-3. Here we have the relative velocity between two objects A and B. The new velocity after a collision will just be the reflection of the initial relative velocity.

vreli = -evreli-1

vrel = n (va - vb)

*Equation 5.4.1-3:* Calculating the new relative velocity for a colliding object which is just a simple reflection of the previous relative velocity scaled by the coefficient of restitution e. The relative velocity is computed by taking the difference between the velocity of the colliding point on object a with the velocity of the colliding point on b and reflecting it along the normal of collision.

vrel = n (va - vb)

Here, e is the coefficient of restitution and is basically a factor estimating the softness of the bodies in collision. It is a value between 0 and 1 and shows the kinetic energy that is lost during the collision due to deformation. A value of 1 means no energy is lost a value of 0 means all energy is lost.

Looking at Equation 5.4.1-3 we can see that if we somehow can define vreli in terms of the previous known velocities of objects A and B at i-1 and of the magnitude j then since vreli-1 is known we can compute j.

We know from Equation 5.4.1-1 and Equation 5.4.1-2 that we can write the new velocity of the colliding point as:

vi = vi-1 + jn/m

*Equation 5.4.1-4:* Calculating the new velocity based upon old velocity and the impulse.

If we do this for the new velocities for both colliding objects A and B we can plug this back into the formula for relative velocity in Equation 5.4.1-3:

vreli = n[(vai-1 + jn/ma) - (vbi-1 - jn/mb)]

*Equation 5.4.1-5:* Calculates the new relative velocity based upon j and the old velocity of object A and B. Here - jn is used since we apply an opposite impulse to object B.

Finally j can now be solved for by using the Equation 5.4.1-3 for the new relative velocity.

-evreli-1 = n[(vai-1 + jn/ma) - (vbi-1 - jn/mb)]

-evreli-1 = n(vai-1 - vbi-1) + j(n/ma + n/mb)

-evreli-1 = vreli-1 + j(n/ma + n/mb)

-evreli-1 - vreli-1 = j(n/ma + n/mb)

-(e + 1)vreli-1 = j(n/ma + n/mb)

j = [-(e + 1)vreli-1]/ (nn * (1/ma + 1/mb))

*Equation 5.4.1-6:* Solving for j.

-evreli-1 = n(vai-1 - vbi-1) + j(n/ma + n/mb)

-evreli-1 = vreli-1 + j(n/ma + n/mb)

-evreli-1 - vreli-1 = j(n/ma + n/mb)

-(e + 1)vreli-1 = j(n/ma + n/mb)

j = [-(e + 1)vreli-1]/ (nn * (1/ma + 1/mb))

Now that j is known we can compute J using equation 5.4.1-2 and then plug the impulse into the Equation 5.4.1-1 to obtain the new linear velocity of the body.

Calculating the impulse while taking into consideration the angular movement involves a similar process to calculating the impulse when only dealing with linear motion. Looking at Equation 5.4.1-1 we can see how to calculate the change in linear velocity from the impulse, in Equation 5.4.2-1 the calculation for the change in angular velocity is shown.

?? = çimp * I^-1

?i = ?i-1 + çimp * I^-1

çimp = r x J

*Equation 5.4.2-1:* Calculating change in angular velocity where çimp is the impulsive torque and is defined as the corss-product of r - the distance from the point of contact and the axis of rotation - and J.

?i = ?i-1 + çimp * I^-1

çimp = r x J

Now J can be broken down the same way here as in Equation 5.4.1-2. The difference is that when we calculate the velocity at the points of collision we must account for the rotational motion. Looking at Equation 5.4.1-3 va and vb, which were the velocities of object A and B now have to take into account the angular velocity of A and B as well. This is shown in Equation 5.4.2-2 now va and vb have been replaced with pa and pb which is the velocity at the point of collision for A and B respectively.

vrel = n (pa - pb)

p = v + ? x r

*Equation 5.4.2-2:* The relative velocity is computed by taking the difference between the velocity of the colliding point on object a with the velocity of the colliding point on b and reflecting it along the normal of collision, the velocity of the colliding point also takes into account the angular velocity. Here r is the distance from the colliding point to the axis of rotation.

p = v + ? x r

The same steps are now taken to solve for j, only not instead of substituting for va and vb, we substitute for pa and pb using the following formulas:

vi = vi-1 + jn/m

?i= ?i-1 + [r x (jn/m)] * I^-1

pi = vi + ?i x r

*Equation 5.4.2-3:* Substituting for J and çimp so that pi can be re-written in known terms so that j can be computed.

?i= ?i-1 + [r x (jn/m)] * I^-1

pi = vi + ?i x r

Now that we have substituted the impulse and impulsive torque into the equation we can solve for j yielding:

t1 = 1/ma + 1/mb

t2 = n [(Ia^-1 * (ra x n)) x ra]

t3 = n [(Ib^-1 * (rb x n)) x rb]

j = [-(e + 1)vreli-1]/ (t1 + t2 + t3)

*Equation 5.4.2-4:* This equation shows the solution for j broken down into three components.

t2 = n [(Ia^-1 * (ra x n)) x ra]

t3 = n [(Ib^-1 * (rb x n)) x rb]

j = [-(e + 1)vreli-1]/ (t1 + t2 + t3)

Finally now that j is known we can solve for the new linear and angular velocity using Equations 5.4.1-1 and 5.4.2-1.

When two bodies collide; while there surfaces are in contact, friction occurs. The friction can be difficult to calculate and may often be ignored with decent results, however if friction is needed to gain a more realistic response it can be estimated by taking into account a tangential force along the collision normal.

t = (n x v) x n

t = t/|t|

*Equation 5.5-1:* Calculating the tangent from the normal and relative velocity at the point of contact. X is the cross-product.

t = t/|t|

Using the tangent vector, the change in linear and angular velocity during collision can now be calculated while taking into consideration a frictional force æ.

vi+1 = vi + (jn + (æj)t)/m

?i+1 = ?i + [l x (jn + (æj)t)]/I

*Equation 5.5-2:* Calculating the new linear and angular velocity while taking into account a frictional force æ. Here l is the location of collision for the body with respect to its center of mass.

?i+1 = ?i + [l x (jn + (æj)t)]/I

The results I had for the collision response portion of this project were greatly dependent upon the accuracy of my collision detection scheme. Throughout the project as the collision detection was repeatedly and steadily improved the collision response results greatly improved. At present since the collision detection routines are not catching penetration, some of the impulse calculations that occur because of repeated collisions can give unrealistic results to the collision response. Fixing this problem depends upon fixing the penetration problems for my collision detection scheme.

I used the simple Euler's Method for integration the equations of motion, this also led to an increase in instability. To help keep the system stable I introduced energy loss damping to prevent the angular motion from getting out of hand. I also introduced some energy loss in the impulse calculations to simulate the minute loss of energy due to surface deformation during contact; this also helped to keep the system more stable.

Since the results can only truly be viewed by seeing the animations of the collision response simulation the second best approach is to show a couple of pictures showing the steps of the response. These pictures in Figure 3 can be seen at the end of this paper.

Possible future work that I would like to do involving Collision response would be to work on accurate collision response for deformable models with energy loss and to speed up the integration of the equations of motion for models using vertex shaders.

Deformable models are much more complicated computationally from rigid body models since many of the constant value mass properties change over time. This makes the computations much more complex and harder to implement in a real-time virtual environment.

Vertex shaders replace the fixed-functional transformation pipeline with programs that can be run on a per vertex basis on programmable graphics hardware. The graphics hardware could then be used to quickly transform the vector and matrix values for motion by pushing them through the vertex shaders.

While working on the project Efficient and Accurate Collision Detection Hierarchies with Collision Response I learned a great deal about the intricacies of collision detection and the dependence of collision response on the accuracy of the collision detection. I gained a great respect for the difficulties in accurately detecting collisions robustly and in an efficient manner.

The problems inherent in numerical computation are of great importance when working with collision response, throughout this project I learned a great deal about stability issues that arise through errors in numerical computation and I gained a glimpse into the field of numerical computation as a whole. I would defiantly like to learn more about this topic and work more with the computation of numerical solutions for virtual environment simulations.

Overall, collision detection and collision response are very well researched areas, specifically when talking about rigid body motion. When working with deformable models there is a lot of work that has also been done but there is still a lot of room for improvements and work yet to be done.

[HECHT96]: "Physics Calculus," Eugene Hecht. Brooks/Cole Publishing Company, 1996.

[BOURG02]: "Physics for Game Developers," David M. Bourg. O'Reilly & Associates, Inc, 2002.

[HUBBARD94]: "Collision Detection for Interactive Graphics Applications," Philip M. Hubbard. Dept. of Comp. Sci. Brown University. Ph.D. Thesis, 1994.

[HUBBARD96]: "Approximating Polyhedra with Spheres for Time-Critical Collision Detection," Philip M. Hubbard. ACM Transactions on Graphics, Vol. 15, No 3, July 1996, Pages 179-210.

[BERG00]: "Computational Geometry: Algorithms and Applications," Mark de Berg, Marc van Kreveld, Mark Overmars, Otfried Schwarzkopf. Springer-Verlag, 2nd rev. ed. 2000.

[O'ROURKE98]: "Computational Geometry in C (2nd Ed.)," Joseph O'Rourke. Cambridge University Press, September 1998.

[LANDER00]: "Collision Response: Bouncy, Trouncy, Fun," Jeff Lander. Gamasutra, 2000.

[BOBIC00]: "Advanced Collision Detection Techniques," Nick Bobic. Gamasutra, 2000.

[KLOSOWSKI98]: ``Efficient Collision Detection Using Bounding Volume Hierarchies of k-DOPs'', J.T. Klosowski, M. Held, J.S.B. Mitchell, H. Sowizral, K. Zikan, IEEE Trans. Visualizat. Comput. Graph. 4(1):21-36, Jan 1998.

[BARAFF89]: "Analytical Methods for Dynamic Simulation of Non-Penetrating Rigid Bodies," David Baraff. Computer Graphics, Volume 23 Number 3, July 1989, SIGGRAPH 1989.

[BARAFF94]: "Fast contact force computation for nonpenetrating rigid bodies ," D. Baraff. Computer Graphics Proceedings, Annual Conference Series: 23-34, 1994.

[BARAFF92]: "Dynamic simulation of non-penetrating flexible bodies," D. Baraff and A. Witkin. Computer Graphics 26(2): 303-308, 1992.

[BARAFF97]: "Physically Based Modeling: Principles and Practice," Andrew Witkin and David Baraff. Online Course Notes for SIGRAPH 97, http://www-2.cs.cmu.edu/~baraff/sigcourse/index.html.

[BRADSHAW03]: "Sphere-Tree Construction Toolkit," Gareth Bradshaw http://isg.cs.tcd.ie/spheretree/, 2003.

[BRADSHAW02]: "Sphere-Tree Construction Using Dynamic Medial Axis Approximation," Bradshaw, G. and O'Sullivan, C. ACM SIGGRAPH Symposium on Computer Animation SCA 2002.

[VAN DEN BERGEN97]: "Efficient Collision Detection of Deformable Models using AABB Trees," Gino van den Bergen, Journal of Graphics Tools, 2(4):1-14, 1997.

[GOTTSCHALK96]: "OBB-Tree: A Hierarchical Structure for Rapid Interference Detection," S. Gottschalk, M. Lin, and D. Manocha. Proc. ACM SIGGRAPH, 1996.

[KAISER00]: "3D Collision Detection," Kevin Kaiser. Game Programming Gems, Charles River Media, 2000. Pgs 391-402.

[OLSEN00]: "Interpolation Methods," John Olsen. Game Programming Gems, Charles River Media, 2000. Pgs 141-149.

[GOMEZ00a]: "Integrating the Equations of Rigid Body Motion," Miguel Gomez. Game Programming Gems, Charles River Media, 2000. Pgs 150-160.

[GOMEZ00b]: "Using Implicit Euler Integration for Numerical Stability," Miguel Gomez. Game Programming Gems, Charles River Media, 2000. Pgs 177-181.