Narrow Phase: Sweeping a Sphere Against a Polygon

Before sweeping a sphere against a polygon, you should determine if the sphere is already intersecting the polygon. I will discuss the case where the sphere is not initially intersecting the polygon.

We will use a 2D representation of the hypothetical setup, since it's easier to draw 2D pictures.

In this setup, there is a sphere with radius r located at point c, a relative velocity vector v that represents the difference of velocity between the sphere and the polygon. This relative velocity can be considered as the velocity of the sphere from the point of view of the polygon. The velocity is given in units that represent how much the sphere will move over the course of a given frame. Thus, at the end of the frame, the center of the sphere should be positioned at the tip of the velocity vector. The polygon lies entirely within some plane, and we are projecting the setup in such a way that the plane occupies a single line. The polygon itself is the bold section of the line. The normal n of the plane is pointing out of the front side of the polygon.


Before determining collision with the polygon, we will first determine collision with the plane. The first step is to find the point on the sphere that will first come into contact with the plane. This can be done scaling the normal of the plane by the radius of the sphere, and subtracting this from the center of the sphere.

We will sweep the point p1 forward. The sweep is accomplished by finding a line that passes through p1, and is parallel to v. We then find the point where the line intersects the plane. The point on the plane is the point where the sphere will first make contact with the plane.



The point p2 can be calculated using the following expression

Here d is the plane parameter, and can be found by dotting the plane normal with any point that resides on the plane, such as one of the vertices of the polygon.

Before calculating p2, check if s is less than 0, or greater than 1. If it is then the collision is not going to happen during this frame, so you can return a false on your collision routine.

If s is between 0 and 1, then go ahead and calculate p2 and continue with the operation.

The next step is to find a point inside the polygon that is nearest to p2. The method for doing this is dependant on the polygon. If p2 already resides within the polygon, then the nearest point to p2 is just p2. Otherwise, the nearest point lies on the boundary of the polygon.

NOTE: Finding the nearest point on the polygon to p2 does not actually work in 3D, since it is possible to have components of velocity that are tangent to the polygon plane. To correct for this, you should try to find the point on the polygon that is closest to the ray that is cast from p1. This is done by finding which feature of the polygon (vertex, edge, or face) is closest to the line, and then finding the closest point on this feature to the line.

The point p3 is the point on the polygon that will first be touched by the sphere.

Now we will sweep the point p3 backward, and see if it touches the sphere. This sweep is performed by finding a line that passes throw p3 and is parallel with v. This line will intersect the sphere in 0, 1, or 2 places.


The formula for this sweep involves a square root. If the stuff we are trying to take the square root of is negative, then there is no intersection. So, we first check to see if this stuff is negative.

If this is negative, then there is no collision and we can return false. Otherwise, we will continue to calculate the fraction of the time step t at which the collision takes place

This quantity should not be less than 0. It may however be greater than 1. If it is, then there will be a collision, but not during this time interval, so return false. If however, t is less than 1, you can proceed to determine the point p4 that represents the point that the sphere will first come into contact with the polygon.

At this point, we can move the sphere forward so that it just touches the polygon. The task now is to determine how to resolve this collision, by adjusting the velocity of the sphere. Once we have determined the new, resolved, velocity we continue to sweep the sphere for the remainder of the time step - checking for collisions with other polygons along the way.

Collision Overview

Collision detection is one of the very key components of a physics engine. However, collision detection on its own is a much more broad topic. It is important to know what kinds of information your collision detection code should be providing, and the different means of accessing this information.

Generally, the primary thing we want the collision system to do is notify us if two objects are intersecting. As a secondary part of this computation, we often would like to know how to keep the objects from intersecting.

There are a couple of methods I have used, when writing collision systems, to deliver this information.

The first is a direct query method, which has the possible functional prototype

bool CollisionSystem::CollideObjects(const PhysObj& A, const PhysObj& B, ColResData* resolve = NULL);

Such a function would return true if the two objects were intersecting. We may optionally pass in a data structure to retrieve information about how to resolve the collision.

The second method is a bit more automatic. All collision information is passed to the user through callback functions. The collision system accumulates all motion of all objects, and automatically determines which objects need collision processing. If a collision event occurs, a callback function is called.

typedef void (*CollisionHandler)(PhysObj&, PhysObj&, ColResData*);
void collisionHandler(PhysObj& A, PhysObj& B, ColResData* resolve);

The physics objects are not passed in as const because the collision handler may change the state of these objects in order to resolve the collision.

We may want to have several different collision event handlers. For instance, we may want to handle an event which corresponds to objects touching, or an event which corresponds to the moment that they cease to touch. There are several different objects in the world, and we may want a different handler for different object types. Therefore, when we proved the collision system with collision handlers, we must supply the object types, and event types for which the collision handler applies

bool CollisionSystem::AddHandler(CollisionHandler handler, int AType, int BType, ColEventType event);

This function passes in a collision handler, and an int for the types of objects A and B. An int is used so that the user may map these types to any meaning desired. If the function returns false, it was not able to map the collision handler, probably because there is already a handler mapped to the given types and event. The ColEventType is an enum that may contains events like

enum ColEventType
{
ColEvt_OnNear,
ColEvt_Near,
ColEvt_OnIntersect,
ColEvt_Intersect,
ColEvt_Resolve,
ColEvt_OnSeperation,
ColEvt_OnFar
};

We may want to add more collision events, but these cover many of the bases. We want to know if two objects have just gotten near to each other, if they are near to each other, if they have just begun to intersect, if they are intersecting, if we want to resolve the collision, if they have just separated, and if they have just moved away from each other.

The way that we actually process collisions is broken up into three distinct portions

1. Broadphase - rough estimation if objects are even close to each other
2. Midphase - determines which parts of a complex object might be colliding
3. Narrowphase - performs collision on convex primitives, and calculates resolution information

In the following posts, I will discuss each of these topics.

Quaternion Tricks

Trick #1 : Optimizing the Quaternion Rotation Formula

The quaternion rotation formula is given as

(0, v') = R(0, v)R^†

This formula can be expressed in terms of standard 3D vector products such as scaling, dot products, and cross products. This equation involves two quaternion products, each of which requires 16 multiplies and 12 additions. Thus there are 56 operations in all.

We can express the vector cross product in terms of the quaternion product. We know that if you switch the factors in a cross product, the product switches sign. Knowing this, you might accept without proof the following formula.

q_V × p_V = (1/2) (qp - pq)

Indeed if you flip the order of the arguments, the product will change sign.

We will use this relationship, in a rearranged form, on the first two factors in the quaternion rotation formula.

R(0, v) = (0, 2R_V × v) + (0, v)R

Inserting the expression on the right hand side into the quaternion rotation formula yields

(0, v') = (0, 2R_V × v)R^† + (0, v)RR^†

When examining the first term on the right hand side, it can be shown that the scalar portion is always zero. For the second term, we can use the fact that the rotation quaternion R obeys the identity RR^† = 1. Since all 3 terms have zero scalar component, we can now rewrite this expression in terms of vector operations.

v' = v + 2[ R_S( R_V × v ) + R_V × ( R_V × v ) ]

This expression employs 18 multiplications and 12 adds, totalling 30 unique operations in all. This optimisation basically cuts the expense of the quaternion rotation formula in half.

Here is a function that employs this omptimisation

Vector rotateVector(const Vector& v, const Quat& q)
{
     Vector result;
     float x1 = q.y*v.z - q.z*v.y;
     float y1 = q.z*v.x - q.x*v.z;
     float z1 = q.x*v.y - q.y*v.x;

     float x2 = q.w*x1 + q.y*z1 - q.z*y1;
     float y2 = q.w*y1 + q.z*x1 - q.x*z1;
     float z2 = q.w*z1 + q.x*y1 - q.y*x1;

     result.x = v.x + 2.0f*x2;
     result.y = v.y + 2.0f*y2;
     result.z = v.z + 2.0f*z2;

     return result;
}

Let's express the scalar and vector portions of the quaternion in terms of the axis and angle of the rotation

R_S = cos( θ/2 )

R_V = n sin( θ/2 )

Pluging these values into the optimized quaternion rotation formula, and using a little bit of trigonometric magic (double angle identity), the formula reduces to

v' = v + [ sin(θ) (n × v) + ( 1 - cos(θ) ) n × (n × v) ]

Using the vector triple product identity, and the fact the n is normalized, we can reduce it a bit further.

v' = cos(θ)v + sin(θ)(n × v) + (1 - cos(θ)) (n•v) n

This operation has 19 mults and 12 adds, totalling 31 operations. We didn't do anything to optimize the formula, but we were able to transform it into a form involving axis and angle.

Trick #2: Determining a Rotation From an Initial and Final Vector

One of the questions which I am asked most frequently is how to find a rotation which takes some initial vector into some final vector. There are several variants of this problem, including "I am pointing my gun at Bob, and I would like to point my gun at Jim." I have seen this problem solved in many different ways, many times involving inverse trig functions. We are going to solve this problem with quaternions - and it's really easy!

We are going to call the initial vector A, and the final vector B. Both A and B are normalized 3D vectors.

We are going to multiply these two 3D vectors together, and if you remember, this is accomplished with a quaternion multiplication. To reitterate, our input vectors are just standard 3D vectors, but the output of the product is a quaternion, which represents the product of the inputs. Since the ordering matters, we are going to place A on the right, and B on the left.

BA = -A•B - A × B

We can relate the dot and cross products to functions involving the angle between A and B. If A and B are normalized, then their product has the following form

BA = - cos(θ) - n sin(θ)

The angle θ is the angle between the initial and final vector, and the vector n is a normalized vector that is perpendicular to the input vectors. This is the exact form of a rotation quaternion. Don't worry about the minus sign, since we can flip the sign without changing the rotation we are representing. In fact, this almost represents the quaternion that we are looking for. We need a quaternion with a half angle instead of a full angle.

Instead of multiplying the initial vector A by the final vector B, let's instead find a vector that represents a halfway rotation between A and B. This is very easy to do - just find the average of A and B!

H = (1/2) (A + B)

Now, this H is not normalized, and it needs to be. If we are going to normalize it anyway, we don't need to worry about the factor of 1/2. The normalized version of H is

H = (A + B) / A + B

Knowing that the angle between A and H is half of the angle between A and B, we can use our previous result to show that the rotation which takes A into B is given by

R = HA

Here is a function that will give you the rotation quaternion that will rotate some initial vector into some final vector

Quat getRotationQuat(const Vector& from, const Vector& to)
{
     Quat result;

     Vector H = VecAdd(from, to);
     H = VecNormalize(H);

     result.w = VecDot(from, H);
     result.x = from.y*H.z - from.z*H.y;
     result.y = from.z*H.x - from.x*H.z;
     result.z = from.x*H.y - from.y*H.x;
     return result;
}

It doesn't get much easier than that folks.

Trick #3: Tangent Space Compression / Quaternion Maps

Tangent space is a set of 3 vectors - normal, tangent, and bi-normal - which are used in lighting calculations such as bump mapping, or anisotropic shading techniques, like fur or hair shaders.

The tangent space is a representation of the surface of the mesh. The tangent space vectors are stored per-vertex, and interpolated across the polygon for per-pixel operations. If you use floats to represent the components of the tangent space vectors, then you are dedicating 3 x 3 x 4 bytes = 36 bytes per vertex to the storage of tangent space.

The tangent space vectors form an ortho-normal frame. Consider the rotation which will transform the ortho-normal frame of tangent space into the standard frame in world space. If this rotation is expressed as a matrix, then the 3 columns of this matrix represent the 3 tangent space vectors.

If you can express a rotation as a matrix, then you can also express it as a quaternion. This immediately reduces the tangent space representation to 4 components = 16 bytes. When decompressing the vertex, use the quat->matrix conversion, and pull out the 3 columns of the resulting matrix as your tangent space vectors. The quat->matrix conversion costs 13 multiplies and 15 subtractions. It may be worthwhile to determine if it is possible to optimize this conversion for shader operations.

If you want to save more space, you can go even further, if you have the time - computationally speaking. Since rotation quaternions are normalized, there are only 3 degrees of freedom instead of 4. Thus you only really need to store 3 components of the quaternion, and the 4th can be calculated using the formula

w = sqrt(1 - x² - y² - z²)

This brings the tangent space size down to 12 bytes.

Finally, if you are super hard pressed for space, you can take advantage of the fact that the components of a quaternion are within the range [-1,1], which may not require 32 bits to have a sufficient representation. You may possibly use 16 bit floats to represent these components. Doing so will reduce the size further - 6 bytes.

With a 6-byte tangent space representation, it seems conceivable to employ quaternion texture maps, which store the tangent space directly in the texture, rather than using bump map data to perturb the interpolated tangent vectors. You could then apply a texture compression, like 3DC which gets 2:1 compression. Now tangent space only consumes 3-bytes of texture, and 0-bytes on the vertex. If you were previously using bump maps, and you replace them with quaternion maps, then you aren't really increasing the amount of texture memory used, and you have successfully removed all 36 bytes of tangent space data from the vertex.

Of course, doing quaternion map lookups in the pixel shader would also mean doing the quat->matrix conversion in the pixel shader - or, alternatively recasting the pixel shader lighting equation in terms of quaternion algebra, rather than matrix algebra - whichever ends up being cheaper. To gain perspective on the cost of this conversion, you can perform 5 quat->matrix conversions in roughly the same amount of time as a single matrix multiply. A single quat->matrix conversion is about as costly as transforming 3 vectors with a matrix.

The quaternion map would take the place of bump maps altogether, and would provide the added feature of having fine tuned control over tangent space on a per-pixel level. For anisotropic shaders, like hair, fur, and velvet shaders, you could introduce swirls and wiggles in the middle of a polygon.

There are two rather large problems facing the idea of a quaternion map: 1) a lot of people get freaked out when you use the word quaternion. 2) I can't off the top of my head conceive of a way for an artist or content creator to generate such a map.

At any rate, it seems like a very intriguing idea.