Hello and welcome back to my blog!

This article is a kind of companion article to Physics engines for dummies and talks about the act of actually detecting a collision between two shapes.

This article assumes the reader has a basic grasp of maths and geometry.

**Definitions**

I will define matrices in bold uppercase, vectors in uppercase and scalars in lowercase. The dot product will be shown as . and cross product as x.

||A|| means the magnitude of A.

|A| means the absolute value of A.

When I want to express individual components of a vector, I will write A_x to mean the x component of A.

I will use * to denote multiplication, but also note that sV might be used to show the scalar s multiplied by the vector V.

## Collision Detection

Briefly, collision detection is the process of detecting when two shapes are about to collide or have already collided. There is a distinction there because there are chiefly two different types of collision detection methods:

**Discrete**

This type will only tell you whether two shapes are intersecting at the current time.

At frame 1, the objects are not colliding, but by frame 2 they are.

**Continuous**

This type tries to predict the exact point and time of collision between the two objects.

At frame 1, the objects are *predicted* to collide before frame 2 arrives. The system estimates the point and time of the collision so that when frame 2 does arrive the objects are not intersecting.

Obviously the continuous method seems the better choice, since we don’t want our objects to be intersecting on screen – however before we can understand the way this method works, we must first explore the other method.

## Discrete collision detection

**Penetration distance**

In discrete collision detection, it’s very important to know exactly how much one shape has penetrated into another and what the shortest direction to push the two apart again is.

This is important because we will need to be able to *resolve *the collision and push the two shapes apart so they are just touching.

In Figure 1, the penetration distance is *d* and the direction to push is *n*.

Let’s start having a look at some common collision detection cases.

**Circle vs circle**

Two circles *A *and *B *are defined as having a centre and a radius. They collide when the distance between the two centres is less than the sum of their radii:

||A-B|| < (r1+r2)

Of course you can speed this up by not having the compute the magnitude ||A-B|| which involves a sqrt by squaring both sides of the equation:

||A-B||² < (r1+r2)²

So to calculate a measure of the penetration we can say:

p = ||A-B|| – (r1+r2)

You can see that this is a simple rearrangement of the collision test inequality

||A-B|| < (r1+r2)

=

||A-B|| – (r1+r2) < 0

All I have done is to subtract (r1+r2) from both sides of the inequality which is handily also the measure of penetration between the two circles. So, when

p < 0

The two spheres are penetrating by distance p. We would also like the *penetration vector *so that we can correct the penetration once we discover it. This is the vector that moves both circles to the point where they just touch, correcting the penetration. Importantly it is not only just a vector that does this, it is the *only *vector which corrects the penetration by moving the minimum amount. This is important because we only want to correct the error, not introduce more by moving too much when we correct, or too little.

N = (A-B) / ||A-B||

P = N*p

Here we have calculated the normalised vector N between the two centres and the penetration vector P by multiplying our unit direction by the penetration distance.

At this point you may have noticed that I haven’t tried to optimise this part of the calculation by using squared lengths etc. This is generally ok because we can still do the initial intersection test using squared maths and then once we have detected an intersection, we can go for the heavier maths since it will happen far less often.

You may also have noticed that during our calculations we have calculated the collision normal N.

**AABB vs AABB**

Similarly as useful as the mighty circle, the Axis Aligned Bounding Box, or AABB is defined by its corner points and represents a un-rotatable rectangle in 2d. These are generally used to bound more complex objects so that a simple overlap check can be done to determine broad overlap, and then a more expensive fine grained check can be performed later.

*Figure 2* shows complex object *A *bounded by its AABB.

For two AABBs *A *and *B *defined by their centres and half-extents (half width, half height), overlap can be determined quite simply:

D = |centreB-centreA| – (halfExtentsA+halfExtentsB)

binary overlap = D_x < 0 && D_y < 0;

In the example below the half extents are shown in red.

Fun AABB facts:

To build an AABB from a complex object, just loop through all the vertices keeping a track of the minimum and maximum coordinate in each axis:

/// <summary> /// /// </summary> static public function BuildAABB(points:Vector.<Vector2>):AABB { var min:Vector2 = new Vector2(Number.MAX_VALUE, Number.MAX_VALUE); var max:Vector2 = min.m_Neg; for each(var p:Vector2 in points) { min = p.Min(min); max = p.Max(max); } return new AABB ( min.Add(max).MulScalar(0.5), max.Sub(min).MulScalar(0.5) ); } |

To build an AABB from two separate AABBs, do the same thing, but form the min and max points of each AABB and use them as the vertices:

static public function CombineAABBs(a:AABB, b:AABB):AABB { const minA:Vector2 = a.m_centre.Sub( a.m_halfExtents ); const maxA:Vector2 = a.m_centre.Add( a.m_halfExtents ); const minB:Vector2 = b.m_centre.Sub( b.m_halfExtents ); const maxB:Vector2 = b.m_centre.Add( b.m_halfExtents ); const min:Vector2 = minA.Min(minB); const max:Vector2 = maxA.Max(maxB); const centre:Vector2 = min.Add(max).MulScalar(0.5); const halfExtents:Vector2 = max.Sub(min).MulScalar(0.5); return new AABB( centre, halfExtents ); } |

**Circle vs AABB**

Not terribly useful in the real world, but very useful as tool for understanding, the circle vs AABB test is also reasonably simple – the trick is to convert this problem into the simpler problem of circle vs circle. For an AABB A and a circle B:

The first step is to get the centre of the circle and project it on to the boundary of the AABB. Fortunately this is as easy as forming the vector *V *between AABB and circle and simply clamping it against the half extents of the AABB. This yields point B_{c}. In *Figure 2*, the clamped x, y region of V is shown in red.

B_{c} is actually the closest point on the AABB to the circle, which means the distance from the AABB to the circle is simply the length of the vector from B_{c} to B minus the radius of the circle.

D = B – B_{c}

distance = ||D|| – radius

**Circle vs OBB**

Circle vs Oriented Bounding Box or OBB is much more useful in the real world, but more tricky because the box now has an orientation to worry about. But, once again we can convert this problem into a simpler one, that of circle vs AABB by simply transforming the circle into the space of the OBB first.

public class OBB { /// <summary> /// /// </summary> public override function Collide( vis:Visual ):Contact { if ( vis is Circle ) { const c:Circle = Circle( vis ); // transform circle into OBB space const v:Vector2 = m_matrix.TransformIntoSpaceOf( c.m_Pos ); const clamped:Vector2 = v.Clamp( m_halfExtents.m_Neg, m_halfExtents ); // form point on OBB const point:Vector2 = m_matrix.TransformBy( clamped ); // get distance to circle const d:Vector2 = c.m_Pos.Sub( point ); const dist:Number = d.m_Len - c.m_Radius; return new Contact( d.m_Unit, dist, point ); } else { throw new NotImplementedException( ); } } } |

That’s the code snippet.

**OBB vs OBB**

Now we’re starting to get into the interesting stuff. Oriented Bounding Box vs Oriented Bounding Box is actually pretty much the same complexity as convex polygon vs convex polygon. A lot of people advocate using the Separating Axis Test in order to solve this problem, but personally I’ve never been able to visualise it and generally, if I can’t visualise something, I’m not going to use it as a solution.

A concept I have been able to visualise and the one I’m going to demonstrate here is that of the Minkowski Difference. Once again the idea behind using this technique is that it enables us to covert this difficult problem into a more simple one.

The Minkowski Difference consists of the convex hull of every point on the boundary of one shape *subtracted* from every point on the boundary of the other shape. Its easier to imagine this as shrinking one of the shapes down until it becomes a point, and expanding the other shape by the shape of the first. *Configuration space* is the space that the Minkowski difference resides within – the set of vector differences of two objects.

*Figure 4* shows objects *A *and *B*, and then also shows object *B *shrunk down to a point, while *A *grows by the shape of *B*. Of course, because its a subtraction of one from the other, the position of the resulting shape changes, and it’s the relationship with the origin (0,0) which then becomes very interesting.

It turns out that when the Minkowski Difference (written *MD *from now on) of two shapes contains the origin, the two shapes are intersecting. Furthermore, the distance from the origin to the MD is actually the distance between the two shapes, and the vector between the two gives you the *Minimum Translational Distance* which is the globally shortest distance you can move either shape from penetration until they both just touch.

In the above demo, you can manipulate the OBBs *A *and *B *by moving them, or you can drag the handles to rotate them. Note the MD *B-A* shown and how it reacts as the objects are moved. The origin is shown as the black dot in the centre – note that when *A *and *B *intersect, the MD contains the origin. The *Minimum Translational Distance* between the two shapes is shown as the blue line.

The observant among you will also notice the coloured edges in the MD – these represent the shape they originated from; red edges come from the red shape *A *and green from *B*.

So, how do you build this MD for two OBBs, then? Each edge in the MD is made up of a vertex from one shape subtracted from an edge from the other. In order to know which vertices we subtract from which edges, we need to define the concept of a *supporting vertex*.

A *supporting vertex* is simply a vertex which is “most in the direction of” a given direction vector. Mathematically, this can be found as the vertex which has the greatest dot product with a given direction vector.

Once we have this we can simply loop through all the edges of *A*, finding the supporting vertex in B for each edge by using reversed edge normal. The edge normal is reversed because we want to find the vertex “most opposite” this edge. As we do so, we project the origin onto each edge, and keep track of the distance – we’re looking for the least negative distance, which will be the *Minimum Translational Distance*.

Once we’ve gone through edges of *A*, we do the same thing with *B*, this time using supporting vertices from A, again tracking the least negative distance. If we’re only concerned about penetration we can terminate our search as soon as we find a positive distance – since this represents a separating axis.

var p0:Vector2 = new Vector2( ); var p1:Vector2 = new Vector2( ); var leastPenetratingDist:Number = -Number.MAX_VALUE; var leastPositiveDist:Number = Number.MAX_VALUE; // face of A, vertices of B for ( var i:int = 0; i<A.m_numPoints; i++ ) { // get world space normal var wsN:Vector2 = A.GetWorldSpaceNormal( i ); // world space edge var wsV0:Vector2 = A.GetWorldSpacePoint( i ); var wsV1:Vector2 = A.GetWorldSpacePoint( ( i+1 )%A.m_numPoints ); // get supporting vertices of B, most opposite face normal var s:Vector.<SupportVertex> = B.GetSupportVertices(wsN.m_Neg); for (var j:int=0; j<s.length; j++) { // form point on plane of minkowski face var mfp0:Vector2 = s[j].m_v.Sub(wsV0); var mfp1:Vector2 = s[j].m_v.Sub(wsV1); var faceDist:Number = mfp0.Dot( wsN ); // project onto minkowski edge var p:Vector2 = ProjectPointOntoEdge( new Vector2( ), mfp0, mfp1 ); // get distance var dist:Number = p.m_Len*Scalar.Sign(faceDist); // track negative if ( dist>leastPenetratingDist ) { p0 = p; leastPenetratingDist = dist; } // track positive if ( dist > 0 && dist<leastPositiveDist ) { p1 = p; leastPositiveDist = dist; } } } // face of B, vertices of A for ( i = 0; i<B.m_numPoints; i++ ) { // get world space normal wsN = B.GetWorldSpaceNormal( i ); // world space edge wsV0 = B.GetWorldSpacePoint( i ); wsV1 = B.GetWorldSpacePoint( ( i+1 )%B.m_numPoints ); // get supporting vertices of A, most opposite face normal s = A.GetSupportVertices(wsN.m_Neg); for (j=0; j<s.length; j++) { // form point on plane of minkowski face mfp0 = wsV0.Sub(s[j].m_v); mfp1 = wsV1.Sub(s[j].m_v); faceDist = -mfp0.Dot( wsN ); // project onto minkowski edge p = ProjectPointOntoEdge( new Vector2( ), mfp0, mfp1 ); // get distance dist = p.m_Len*Scalar.Sign(faceDist); // track negative if ( dist>leastPenetratingDist ) { p0 = p; leastPenetratingDist = dist; } // track positive if ( dist > 0 && dist<leastPositiveDist ) { p1 = p; leastPositiveDist = dist; } } } if ( leastPenetratingDist<0 ) { // penetration by leastPenetratingDist } else { // separated by leastPositiveDist } |

The above snippet also tracks positive distance so can be simplified if you only want to track penetration. One caveat to watch out for is the case where the two shapes are exactly aligned – as they are in the demo above. In this case, you have to deal with the possibility of there being two supporting vertices for an edge instead of one, and it’s important to re-project onto the edge from the origin to make sure you don’t accidentally pick the wrong edge-vertex pair. The above code deals with this case.

**Poly vs Poly**

The really good thing about this technique is that it applies directly to arbitrary convex polygons as well.

Again, manipulate the live demo by dragging the shapes and their handles…

So, we can now calculate the penetration distance for two arbitrary convex polygons.

## Continuous

So, finally we get around to talking about continuous collision detection. Now we start to be able to reap the rewards of all this *Minkowski Difference* nonsense

**Continuous Linear**

It turns out that if all we want to do is a linear sweep of our shapes to find the first time of collision ignoring rotation, the only thing we have to do is to ray-cast from the origin towards the MD using the relative velocity of the two objects as the ray segment end!

In *Figure 5*, the first time of collision of *B *with *A *can be found by shooting a ray from the origin towards the MD of *B-A*.

For polygons this is the same amount of work as calculating the penetration distance – you fire the ray at the edges of the MD, which you build as usual. You can even avoid testing edges if the ray is back-facing with respect to the edge normal.

In the above demo, A is being swept against B along A’s velocity; the pink A’ shows where A hits B. The MD shows what this looks like in Configuration Space – don’t forget to rotate/move A and B to see what the result looks like. Notice how the ray at the origin fired against the MD corresponds exactly to when A hits B?

One caveat is that A must not intersect B to start with, or the test will fail. Gino Van Den Bergen has generalised this approach to work with the GJK algorithm, this is well worth taking a look at if you want to cast against non-polygonal shapes, like circles or capsules.

The basic idea is this:

You provide a function which can return the closest distance between two shapes and the normal at that point. These two return values are also the distance from the MD to the origin and the direction that distance is in. In this example, the function used is returning distance and normal between two circles – shown in the diagram is the MD of two circles.

In *Figure 6*, two iterations of Gino Van Den Bergen’s algorithm are shown, I_{1} and I_{2}. The black point represents the current estimate, and starts at the origin in I_{1} and moves along the velocity shown as the arrow. The green point represents part of the information returned by the function – the closest point on the MD, which is at the distance returned, also the normal at that point has been found as well (again, returned by the function).

The algorithm intersects the velocity ray against the returned normal and produces the red point shown in the diagram. Then it moves on to the next iteration. The red point becomes the new estimate and the function is queried again for the distance to the MD and the normal. The information returned results in the new red intersection point in I_{2} and at this point the diagram becomes less useful because we’re so close to the visible tolerance. In practice the algorithm will terminate when the distance returned from the function is less than some predefined tolerance.

I’m describing this algorithm because it will be useful when we think about intersection including angular velocity.

**Continuous Angular/Linear**

Things start to get a little more tricky when we want to calculate the exact hit time where there is rotation present. To solve this problem, we must turn to the work of Brian Mirtich, Gino Van Den Bergen and Erwin Coumans.

Brian Mitch’s thesis presents some of the fundamental building blocks used in solving this problem:

In *Figure 7*, (lifted directly from his thesis) he shows two bodies, with their closest points shown and the two actual contact points *x1 *and *x2 *for this collision. Given the fact that these bodies are convex, he shows that the distance between the closest points represents the first possible time of collision between the two shapes under **no rotation**. Further, he shows that, given the point with the largest distance from the centre for each body, it’s possible to bound the earliest time of collision between two rotating shapes.

You can see why by looking at *Figure 8* – on the left, the thin rectangle rotates and extreme point p moves to p’ following the rotation. The length of the arc travelled by p can be found by simply using the circular arc length formula from geometry – L=theta*r as shown on the right in *Figure 8*. Therefore if we know the maximum distance from the origin of the shape, and we know the angular velocity, we can find an *upper bound* for the linear component of the angular rotation of a shape, and therefore we can simply add this linear distance value on to the bound that we had for non-rotating shapes.

Erwin Coumans took this idea and using a method similar to the linear cast of Gino Van Den Bergen, produced a new algorithm which determines the time of impact between translating/rotating bodies:

The live demo above shows A and A’, the position of A after some linear and angular motion, 90 degrees to be precise. The various iterations of the MD are shown on the right. Don’t forget to move the shapes around!

The thing to note is that the **shape **of the MD changes under rotation of either body, unlike with linear only motion where it just translates.

Each of the MDs shown represent one iteration of the algorithm. The algorithm itself is exactly the same I described using Figure 6 above, except that there is an extra distance term added which means that instead of the ray intersecting the plane and the estimate getting moved directly there, it only approaches it instead – taking care not to advance more than the extra distance terms allows each iteration. This way, the correct result is approached each iteration, or the algorithm terminates with a miss.

The actual code behind the algorithm is quite small:

public function AngularCast( B:Polygon, velA:Vector2, velB:Vector2, avA:Number, avB:Number ):LinearCastResult { const A:Polygon = this; const origin:Vector2 = new Vector2( ); var t:Number = 0; // relative linear velocity const ray:Vector2 = velB.Sub( velA ); var mA:Matrix23 = new Matrix23( A.m_Pos.Add( velA.MulScalar( t ) ), A.m_angle+avA*t ); var mB:Matrix23 = new Matrix23( B.m_Pos.Add( velB.MulScalar( t ) ), B.m_angle+avB*t ); var contact:Contact = GetDistance( B, mA, mB, parent, screenCentre ); while ( contact.m_dist>kAngularTollerance ) { // intersect velocity against normal const rayN:Number = ray.Dot( contact.m_normal ); const dRel:Number = rayN + A.m_r*Math.abs(avA) + B.m_r*Math.abs(avB); // compute conservative advancement t += contact.m_dist/dRel; if ( t<0 || t>1 ) { // never hit return new LinearCastResult( false, t, new Vector2( ) ); } // interpolate mA = new Matrix23( A.m_Pos.Add( velA.MulScalar( t ) ), A.m_angle+avA*t ); mB = new Matrix23( B.m_Pos.Add( velB.MulScalar( t ) ), B.m_angle+avB*t ); // get distance again contact = GetDistance( B, mA, mB, parent, screenCentre ); } return new LinearCastResult( true, t, contact.m_normal ); } |

As you play around with the demo, you’ll notice it sometimes takes this algorithm a while to converge – in production code you would most likely limit the number of allowed iterations. Something like having it terminate if its not making ‘enough’ progress in each iteration might be a good idea.

## Conclusion

In this article, we have explored the various forms of collision detection and the algorithms behind them, from the simplest possible case of discrete circle vs circle to the complex case of continuous translating/rotating polygons. I hope this article has been intelligible and will inspire people to start diving in and exploring these techniques in more detail.

## Source code

As always, you can buy the source-code accompanying this article if you wish – it is written in Actionscript 3.0 but the techniques are applicable to all OO languages. All purchases are very much appreciated and ensure that I can afford to make more like this one

You get code for every single demo on this page, containing useful functions for computing the distance between arbitrary convex polygons, penetrating or not, and for computing the exact time of contact between rotating/translating polygons and best of all, the licence allows you to use this code for commercial and non-commercial applications alike

USD 5.99

Subscribers can access the source here

Until next time, have fun!

Cheers, Paul.

Hi Paul,

A happy new year to you.

Thank you as ever for a very interesting and inspiring read; I used some of your ideas (along with many others) to write an (almost) complete 2d physics engine ..it’s taken a long time!

I have just revisted some of my code which is based on some of the ideas in this article, because I found a problem, and I believe that I can present to you a (rare) edge case where the supported vertices approach detailed here can fail to resolve the correct penetration distance.

I have a slightly tilted falling box whose motion bounds eventually intersect the bounds of a static triangle; the lower box edge is wrongly selected as the edge feature. Because i’m using motion bounds, one of the box corners is slightly higher than the top point of the triangle, and just misses it in x- (the top point of the triangle is a supporting vertex because it obviously most opposes the bottom edge). the projected point is clamped to the lower edge as normal- and it’s distance is certainly the smallest- but it is large enough combined with the velocity in my speculative-contact-impulse calulation to determine that a corrective impulse is to be applied, even though the polys completely miss each other! this code has been given *a lot* of testing, and appears to be the only outstanding problem.

Im aware that even though i’ve had a go, it’s very difficult to explain this in words. I can send you small images of the complete setup if you wish.

I wondered if perhaps you had time to pass your eye over this problem with me, out of interest sake? I like this appoach, it has elegance- and I really don’t want to revert to something like SAT unless i can help it..

cheers,

Jon

Hi Paul,

I’ve been debugging this thing- and realised i missed a wrong negative sign, which fixes my issue.

thanks again for an excellent article.

cheers,

Jon

Hi Jon,

Glad you got your problem sorted, all the best with it!

Cheers, Paul.

Hi Paul

I am currently developing a haptic stage, I import models in 3DS, I have a problem of collision detection between objets.j I need help.

thank you

Hey Paul!

Got a quick question, its about the angular velocities. I am trying to understand how the rotational information gets tacked onto linear velocity information. I’ve been examining it for a while, and have made little to know break-through’s

Hi there,

I recommend taking a look at this article by Chris Hecker which explains it quite well: http://chrishecker.com/images/c/c2/Gdmphys2.pdf

Cheers, Paul.