## Game Physics Tutorial Homework

This is Part II of our three-part series on video game physics. For the rest of this series, see:

Part I: An Introduction to Rigid Body Dynamics
Part III: Constrained Rigid Body Simulation

In Part I of this series, we explored rigid bodies and their motions. In that discussion, however, objects did not interact with each other. Without some additional work, the simulated rigid bodies can go right through each other, or “interpenetrate”, which is undesirable in the majority of cases.

In order to more realistically simulate the behavior of solid objects, we have to check if they collide with each other every time they move, and if they do, we have to do something about it, such as applying forces that change their velocities, so that they will move in the opposite direction. This is where understanding collision physics is particularly important for game developers.

In Part II, we will cover the collision detection step, which consists of finding pairs of bodies that are colliding among a possibly large number of bodies scattered around a 2D or 3D world. In the next, and final, installment, we’ll talk more about “solving” these collisions to eliminate interpenetrations.

For a review of the linear algebra concepts referred to in this article, you can refer to the linear algebra crash course in Part I.

## Collision Physics in Video Games

In the context of rigid body simulations, a collision happens when the shapes of two rigid bodies are intersecting, or when the distance between these shapes falls below a small tolerance.

If we have n bodies in our simulation, the computational complexity of detecting collisions with pairwise tests is O(n2), a number that makes computer scientists cringe. The number of pairwise tests increases quadratically with the number of bodies, and determining if two shapes, in arbitrary positions and orientations, are colliding is already not cheap. In order to optimize the collision detection process, we generally split it in two phases: broad phase and narrow phase.

The broad phase is responsible for finding pairs of rigid bodies that are potentially colliding, and excluding pairs that are certainly not colliding, from amongst the whole set of bodies that are in the simulation. This step must be able to scale really well with the number of rigid bodies to make sure we stay well under the O(n2) time complexity. To achieve that, this phase generally uses space partitioning coupled with bounding volumes that establish an upper bound for collision.

The narrow phase operates on the pairs of rigid bodies found by the broad phase that might be colliding. It is a refinement step where we determine if the collisions are actually happening, and for each collision that is found, we compute the contact points that will be used to solve the collisions later.

In the next sections we’ll talk about some algorithms that can be used in the broad phase and narrow phase.

In the broad phase of collision physics for video games we need to identify which pairs of rigid bodies might be colliding. These bodies might have complex shapes like polygons and polyhedrons, and what we can do to accelerate this is to use a simpler shape to encompass the object. If these bounding volumes do not intersect, then the actual shapes also do not intersect. If they intersect, then the actual shapes might intersect.

Some popular types of bounding volumes are oriented bounding boxes (OBB), circles in 2D, and spheres in 3D. Let’s look at one of the simplest bounding volumes: axis-aligned bounding boxes (AABB).

AABBs get a lot of love from physics programmers because they are simple and offer good tradeoffs. A 2-dimensional AABB may be represented by a struct like this in the C language:

The field represents the location of the lower left corner of the box and the field represents the top right corner.

To test if two AABBs intersect, we only have to find out if their projections intersect on all of the coordinate axes:

This code has the same logic of the function from the Box2D engine (version 2.3.0). It calculates the difference between the and of both AABBs, in both axes, in both orders. If any of these values is greater than zero, the AABBs don’t intersect.

Even though an AABB overlap test is cheap, it won’t help much if we still do pairwise tests for every possible pair since we still have the undesirable O(n2) time complexity. To minimize the number of AABB overlap tests we can use some kind of space partitioning, which which works on the same principles as database indices that speed up queries. (Geographical databases, such as PostGIS, actually use similar data structures and algorithms for their spatial indexes.) In this case, though, the AABBs will be moving around constantly, so generally, we must recreate indices after every step of the simulation.

There are plenty of space partitioning algorithms and data structures that can be used for this, such as uniform grids, quadtrees in 2D, octrees in 3D, and spatial hashing. Let us take a closer look at two popular spatial partitioning algorithms: sort and sweep, and bounding volume hierarchies (BVH).

### The Sort and Sweep Algorithm

The sort and sweep method (alternatively known as sweep and prune) is one of the favorite choices among physics programmers for use in rigid body simulation. The Bullet Physics engine, for example, has an implementation of this method in the class.

The projection of one AABB onto a single coordinate axis is essentially an interval [b, e] (that is, beginning and end). In our simulation, we’ll have many rigid bodies, and thus, many AABBs, and that means many intervals. We want to find out which intervals are intersecting.

In the sort and sweep algorithm, we insert all b and e values in a single list and sort it ascending by their scalar values. Then we sweep or traverse the list. Whenever a b value is encountered, its corresponding interval is stored in a separate list of active intervals, and whenever an e value is encountered, its corresponding interval is removed from the list of active intervals. At any moment, all the active intervals are intersecting. (Check out David Baraff’s Ph. D Thesis, p. 52 for details. I suggest using this online tool to view the postscript file.) The list of intervals can be reused on each step of the simulation, where we can efficiently re-sort this list using insertion sort, which is good at sorting nearly-sorted lists.

In two and three dimensions, running the sort and sweep, as described above, over a single coordinate axis will reduce the number of direct AABB intersection tests that must be performed, but the payoff may be better over one coordinate axis than another. Therefore, more sophisticated variations of the sort and sweep algorithm are implemented. In his book Real-Time Collision Detection (page 336), Christer Ericson presents an efficient variation where he stores all AABBs in a single array, and for each run of the sort and sweep, one coordinate axis is chosen and the array is sorted by the value of the AABBs in the chosen axis, using quicksort. Then, the array is traversed and AABB overlap tests are performed. To determine the next axis that should be used for sorting, the variance of the center of the AABBs is computed, and the axis with greater variance is chosen for the next step.

### Dynamic Bounding Volume Trees

Another useful spatial partitioning method is the dynamic bounding volume tree, also known as Dbvt. This is a type of bounding volume hierarchy.

The Dbvt is a binary tree in which each node has an AABB that bounds all the AABBs of its children. The AABBs of the rigid bodies themselves are located in the leaf nodes. Typically, a Dbvt is “queried” by giving the AABB for which we would like to detect intersections. This operation is efficient because the children of nodes that do not intersect the queried AABB do not need to be tested for overlap. As such, an AABB collision query starts from the root, and continues recursively through the tree only for AABB nodes that intersect with the queried AABB. The tree can be balanced through tree rotations, as in an AVL tree.

Box2D has a sophisticated implementation of Dbvt in the class. The class is responsible for performing the broad phase, and it uses an instance of to perform AABB queries.

## Narrow Phase

After the broad phase of video game collision physics, we have a set of pairs of rigid bodies that are potentially colliding. Thus, for each pair, given the shape, position and orientation of both bodies, we need to find out if they are, in fact, colliding; if they are intersecting or if their distance falls under a small tolerance value. We also need to know what points of contact are between the colliding bodies, since this is needed to resolve the collisions later.

### Convex and Concave Shapes

As a video game physics general rule, it is not trivial to determine if two arbitrary shapes are intersecting, or to compute the distance between them. However, one property that is of critical importance in determining just how hard it is, is the convexity of the shape. Shapes can be either convex or concave and concave shapes are harder to work with, so we need some strategies to deal with them.

In a convex shape, a line segment between any two points within the shape always falls completely inside the shape. However for a concave (or “non-convex”) shape, the same is not true for all possible line segments connecting points in the shape. If you can find at least one line segment that falls outside of the shape at all, then the shape is concave.

Computationally, it is desirable that all shapes are convex in a simulation, since we have a lot of powerful distance and intersection test algorithms that work with convex shapes. Not all objects will be convex though, and usually we work around them in two ways: convex hull and convex decomposition.

The convex hull of a shape is the smallest convex shape that fully contains it. For a concave polygon in two dimensions, it would be like hammering a nail on each vertex and wrapping a rubber band around all nails. To calculate the convex hull for a polygon or polyhedron, or more generally, for a set of points, a good algorithm to use is the quickhull algorithm, which has an average time complexity of O(n log n).

Obviously, if we use a convex hull to represent a concave object, it will lose its concave properties. Undesirable behavior, such as “ghost” collisions may become apparent, since the object will still have a concave graphical representation. For example, a car usually has a concave shape, and if we use a convex hull to represent it physically and then put a box on it, the box might appear to be floating in the space above.

In general, convex hulls are often good enough to represent concave objects, either because the unrealistic collisions are not very noticeable, or their concave properties are not essential for whatever is being simulated. In some cases, though, we might want to have the concave object behave like a concave shape physically. For example, if we use a convex hull to represent a bowl physically, we won’t be able to put anything inside of it. Objects will just float on top of it. In this case, we can use a convex decomposition of the concave shape.

Convex decomposition algorithms receive a concave shape and return a set of convex shapes whose union is identical to the original concave shape. Some concave shapes can only be represented by a large number of convex shapes, and these might become prohibitively costly to compute and use. However, an approximation is often good enough, and so, algorithms such as V-HACD produce a small set of convex polyhedrons out of a concave polyhedron.

In many collisons physics cases, though, the convex decomposition can be made by hand, by an artist. This eliminates the need to tax performance with decomposition algorithms. Since performance is one of the most important aspects in real-time simulations, it’s generally a good idea to create very simple physical representations for complex graphic objects. The image below shows one possible convex decomposition of a 2D car using nine convex shapes.

### Testing for Intersections - The Separating Axis Theorem

The separating axis theorem (SAT) states that two convex shapes are not intersecting if and only if there exists at least one axis where the orthogonal projections of the shapes on this axis do not intersect.

It’s usually more visually intuitive to find a line in 2D or a plane in 3D that separates the two shapes, though, which is effectively the same principle. A vector orthogonal to the line in 2D, or the normal vector of the plane in 3D, can be used as the “separating axis”.

Game physics engines have a number of different classes of shapes, such as circles (spheres in 3D), edges (a single line segment), and convex polygons (polyhedrons in 3D). For each pair of shape type, they have a specific collision detection algorithm. The simplest of them is probably the circle-circle algorithm:

Even though the SAT applies to circles, it’s much simpler to just check if the distance between their centers is smaller than the sum of their radii. For that reason, the SAT is used in the collision detection algorithm for specific pairs of shape classes, such as convex polygon against convex polygon (or polyhedrons in 3D).

For any pair of shapes, there are an infinite number of axes we can test for separation. Thus, determining which axis to test first is crucial for an efficient SAT implementation. Fortunately, when testing if a pair of convex polygons collide, we can use the edge normals as potential separating axes. The normal vector n of an edge is perpendicular to the edge vector, and points outside the polygon. For each edge of each polygon, we just need to find out if all the vertices of the other polygon are in front of the edge.

If any test passes – that is, if, for any edge, all vertices of the other polygon are in front of it – then the polygons do not intersect. Linear algebra provides an easy formula for this test: given an edge on the first shape with vertices a and b and a vertex v on the other shape, if (v - a) · n is greater than zero, then the vertex is in front of the edge.

For polyhedrons, we can use the face normals and also the cross product of all edge combinations from each shape. That sounds like a lot of things to test; however, to speed things up, we can cache the last separating axis we used and try using it again in the next steps of the simulation. If the cached separating axis does not separate the shapes anymore, we can search for a new axis starting from faces or edges that are adjacent to the previous face or edge. That works because the bodies often don’t move or rotate much between steps.

Box2D uses SAT to test if two convex polygons are intersecting in its polygon-polygon collision detection algorithm in b2CollidePolygon.cpp.

### Computing Distance - The Gilbert-Johnson-Keerthi Algorithm

In many collisions physics cases, we want to consider objects to be colliding not only if they are actually intersecting, but also if they are very close to each other, which requires us to know the distance between them. The Gilbert-Johnson-Keerthi (GJK) algorithm computes the distance between two convex shapes and also their closest points. It is an elegant algorithm that works with an implicit representation of convex shapes through support functions, Minkowski sums, and simplexes, as explained below.

Support Functions

A support functionsA(d) returns a point on the boundary of the shape A that has the highest projection on the vector d. Mathematically, it has the highest dot product with d. This is called a support point, and this operation is also known as support mapping. Geometrically, this point is the farthest point on the shape in the direction of d.

Finding a support point on a polygon is relatively easy. For a support point for vector d, you just have to loop through its vertices and find the one which has the highest dot product with d, like this:

However, the real power of a support function is that makes it easy to work with shapes such as cones, cylinders, and ellipses, among others. It is rather difficult to compute the distance between such shapes directly, and without an algorithm like GJK you would usually have to discretize them into a polygon or polyhedron to make things simpler. However, that might lead to further problems because the surface of a polyhedron is not as smooth as the surface of, say, a sphere, unless the polyhedron is very detailed, which can lead to poor performance during collision detection. Other undesirable side effects might show up as well; for example, a “polyhedral” sphere might not roll smoothly.

To get a support point for a sphere centered on the origin, we just have to normalize d (that is, compute d / ||d||, which is a vector with length 1 (unit length) that still points in the same direction of d) and then we multiply it by the sphere radius.

Check Gino van den Bergen’s paper to find more examples of support functions for cylinders, and cones, among other shapes.

Our objects will, of course, be displaced and rotated from the origin in the simulation space, so we need to be able to compute support points for a transformed shape. We use an affine transformationT(x) = Rx + c to displace and rotate our objects, where c is the displacement vector and R is the rotation matrix. This transformation first rotates the object about the origin, and then translates it. The support function for a transformed shape A is:

Minkowski Sums

The Minkowski sum of two shapes A and B is defined as:

That means we compute the sum for all points contained in A and B. The result is like inflatingA with B.

Similarly, we define the Minkowski difference as:

which we can also write as the Minkowski sum of A with -B:

For convex A and B, A⊕B is also convex.

One useful property of the Minkowski difference is that if it contains the origin of the space, the shapes intersect, as can be seen in the previous image. Why is that? Because if two shapes intersect, they have at least one point in common, which lie in the same location in space, and their difference is the zero vector, which is actually the origin.

Another nice feature of the Minkowski difference is that if it doesn’t contain the origin, the minimum distance between the origin and the Minkowski difference is the distance between the shapes.

The distance between two shapes is defined as:

In other words, the distance between A and B is the length of the shortest vector that goes from A to B. This vector is contained in A⊖B and it is the one with the smallest length, which consequently is the one closest to the origin.

It is generally not simple to explicitly build the Minkowski sum of two shapes. Fortunately, we can use support mapping here as well, since:

Simplexes

The GJK algorithm iteratively searches for the point on the Minkowski difference closest to the origin. It does so by building a series of simplexes that are closer to the origin in each iteration. A simplex – or more specifically, a k-simplex for an integer k – is the convex hull of k + 1affinely independent points in a k-dimensional space. That is, if for two points, they must not coincide, for three points they additionally must not lie on the same line, and if we have four points they also must not lie on the same plane. Hence, the 0-simplex is a point, the 1-simplex is a line segment, the 2-simplex is a triangle and the 3-simplex is a tetrahedron. If we remove a point from a simplex we decrement its dimensionality by one, and if we add a point to a simplex we increment its dimensionality by one.

GJK in Action

Let’s put this all together to see how GJK works. To determine the distance between two shapes A and B, we start by taking their Minkowski difference A⊖B. We are searching for the closest point to the origin on the resulting shape, since the distance to this point is the distance between the original two shapes. We choose some point v in A⊖B, which will be our distance approximation. We also define an empty point set W, which will contain the points in the current test simplex.

Then we enter a loop. We start by getting the support point w = sA⊖B(-v), the point on A⊖B whose projection onto v is closest to the origin.

If ||w|| is not much different than ||v|| and the angle between them didn’t change much (according to some predefined tolerances), it means the algorithm has converged and we can return ||v|| as the distance.

Otherwise, we add w to W. If the convex hull of W (that is, the simplex) contains the origin, the shapes intersect, and this also means we are done. Otherwise, we find the point in the simplex that is closest to the origin and then we reset v to be this new closest approximation. Finally, we remove whatever points in W that do not contribute to the closest point computation. (For example, if we have a triangle, and the closest point to the origin lies in one of its edges, we can remove the point from W that is not a vertex of this edge.) Then we repeat these same steps until the algorithm converges.

The next image shows an example of four iterations of the GJK algorithm. The blue object represents the Minkowski difference A⊖B and the green vector is v. You can see here how v hones in on the closest point to the origin.

For a detailed and in-depth explanation of the GJK algorithm, check out the paper A Fast and Robust GJK Implementation for Collision Detection of Convex Objects, by Gino van den Bergen. The blog for the dyn4j physics engine also has a great post on GJK.

Box2D has an implementation of the GJK algorithm in b2Distance.cpp, in the function. It only uses GJK during time of impact computation in its algorithm for continuous collision detection (a topic we will discuss further down).

The Chimpunk physics engine uses GJK for all collision detection, and its implementation is in cpCollision.c, in the function. If the GJK algorithm reports intersection, it still needs to know what the contact points are, along with the penetration depth. To do that, it uses the Expanding Polytope Algorithm, which we shall explore next.

### Determining Penetration Depth - The Expanding Polytope Algorithm

As stated above, if the shapes A and B are intersecting, GJK will generate a simplex W that contains the origin, inside the Minkowski difference A⊖B. At this stage, we only know that the shapes intersect, but in the design of many collision detection systems, it is desirable to be able to compute how much intersection we have, and what points we can use as the points of contact, so that we handle the collision in a realistic way. The Expanding Polytope Algorithm (EPA) allows us to obtain that information, starting where GJK left off: with a simplex that contains the origin.

The penetration depth is the length of the minimum translation vector (MTV), which is the smallest vector along which we can translate an intersecting shape to separate it from the other shape.

When two shapes are intersecting, their Minkowski difference contains the origin, and the point on the boundary of the Minkowski difference that is closest to the origin is the MTV. The EPA algorithm finds that point by expanding the simplex that GJK gave us into a polygon; successively subdividing it’s edges with new vertices.

First, we find the edge of the simplex closest to the origin, and compute the support point in the Minkowski difference in a direction that is normal to the edge (i.e. perpendicular to it and pointing outside the polygon). Then we split this edge by adding this support point to it. We repeat these steps until the length and direction of the vector doesn’t change much. Once the algorithm converges, we have the MTV and the penetration depth.

Using GJK in combination with EPA, we get a detailed description of the collision, no matter if the objects are already intersecting, or just close enough to be considered a collision.

The EPA algorithm is described in the paper Proximity Queries and Penetration Depth Computation on 3D Game Objects, also written by Gino van den Bergen. The dyn4j blog also has a post about EPA.

## Continuous Collision Detection

The video game physics techniques presented so far perform collision detection for a static snapshot of the simulation. This is called discrete collision detection, and it ignores what happens between the previous and current steps. For this reason, some collisions might not be detected, especially for fast moving objects. This issue is known as tunneling.

Continuous collision detection techniques attempt to find when two bodies collided between the previous and the current step of the simulation. They compute the time of impact, so we can then go back in time and process the collision at that point.

The time of impact (or time of contact) tc is the instant of time when the distance between two bodies is zero. If we can write a function for the distance between two bodies along time, what we want to find is the smallest root of this function. Thus, the time of impact computation is a root-finding problem.

For the time of impact computation, we consider the state (position and orientation) of the body in the previous step at time ti-1, and in the current step at time ti. To make things simpler, we assume linear motion between the steps.

Let’s simplify the problem by assuming the shapes are circles. For two circles C1 and C2, with radius r1 and r2, where their center of mass x1 and x2 coincide with their centroid (i.e., they naturally rotate about their center of mass), we can write the distance function as:

Considering linear motion between steps, we can write the following function for the position of C1 from ti-1 to ti

It is a linear interpolation from x1(ti-1) to x1(ti). The same can be done for x2. For this interval we can write another distance function:

Set this expression equal to zero and you get a quadratic equation on t. The roots can be found directly using the quadratic formula. If the circles don’t intersect, the quadratic formula will not have a solution. If they do, it might result in one or two roots. If it has only one root, that value is the time of impact. If it has two roots, the smallest one is the time of impact and the other is the time when the circles separate. Note that the time of impact here is a number from 0 to 1. It is not a time in seconds; it is just a number we can use to interpolate the state of the bodies to the precise location where the collision happened.

Continuous Collision Detection for Non-Circles

Writing a distance function for other kinds of shapes is difficult, primarily because their distance depends on their orientations. For this reason, we generally use iterative algorithms that move the objects closer and closer on each iteration until they are close enough to be considered colliding.

The conservative advancement algorithm moves the bodies forward (and rotates them) iteratively. In each iteration it computes an upper bound for displacement. The original algorithm is presented in Brian Mirtich’s PhD Thesis (section 2.3.2), which considers the ballistic motion of bodies. This paper by Erwin Coumans (the author of the Bullet Physics Engine) presents a simpler version that uses constant linear and angular velocities.

The algorithm computes the closest points between shapes A and B, draws a vector from one point to the other, and projects the velocity on this vector to compute an upper bound for motion. It guarantees that no points on the body will move beyond this projection. Then it advances the bodies forward by this amount and repeats until the distance falls under a small tolerance value.

It may take too many iterations to converge in some cases, for example, when the angular velocity of one of the bodies is too high.

## Resolving Collisions

Once a collision has been detected, it is necessary to change the motions of the colliding objects in a realistic way, such as causing them to bounce off each other. In the next and final installment in this theories, we’ll discuss some popular methods for resolving collisions in video games.

## References

If you are interested in obtaining a deeper understanding about collision physics such as collision detection algorithms and techniques, the book Real-Time Collision Detection, by Christer Ericson, is a must-have.

Since collision detection relies heavily on geometry, Toptal’s article Computational Geometry in Python: From Theory to Application is an excellent introduction to the topic.

Nilson Souto, Brazil

member since January 28, 2013

Parse SDKMac OSC++Adobe PhotoshopObjective-CUI KitiOSOpenGLXcodeGitOpenGL ES+   more

Nilson started programming in C/C++ after playing video games for the first time at a young age. Over the last few years, he has worked on iOS applications mostly, and he's now focusing on game development, computer graphics, physics simulation, and vehicle simulation. He's also a 2D/3D technical artist. [click to continue...]

Hiring? Meet the Top 10 Freelance Game Developers for Hire in March 2018

There are many reasons you might want to create a custom physics engine: first, learning and honing your skills in mathematics, physics and programming are great reasons to attempt such a project; second, a custom physics engine can tackle any sort of technical effect the creator has the skill to create. In this article I would like to provide a solid introduction on how to create a custom physics engine entirely from scratch.

Physics provides a wonderful means for allowing a player to immerse themselves within a game. It makes sense that a mastery of a physics engine would be a powerful asset for any programmer to have at their disposal. Optimizations and specializations can be made at any time due to a deep understanding of the inner workings of the physics engine.

By the end of this tutorial the following topics will have been covered, in two dimensions:

• Simple collision detection
• Simple manifold generation
• Impulse resolution

Here's a quick demo:

Note: Although this tutorial is written using C++, you should be able to use the same techniques and concepts in almost any game development environment.

## Prerequisites

This article involves a fair amount of mathematics and geometry, and to a much lesser extent actual coding. A couple prerequisites for this article are:

• A basic understanding of simple vector math
• The ability to perform algebraic math

## Collision Detection

There are quite a few articles and tutorials throughout the internet, including here on Tuts+, that cover collision detection. Knowing this, I would like to run through the topic very quickly as this section is not the focus of this article.

### Axis Aligned Bounding Boxes

An Axis Aligned Bounding Box (AABB) is a box that has its four axes aligned with the coordinate system in which it resides. This means it is a box that cannot rotate, and is always squared off at 90 degrees (usually aligned with the screen). In general it is referred to as a "bounding box" because AABBs are used to bound other more complex shapes.

The AABB of a complex shape can be used as a simple test to see if more complex shapes inside the AABBs can possibly be intersecting. However in the case of most games the AABB is used as a fundamental shape, and does not actually bound anything else. The structure of your AABB is important. There are a few different ways to represent an AABB, however this is my favorite:

struct AABB { Vec2 min; Vec2 max; };

This form allows an AABB to be represented by two points. The min point represents the lower bounds of the x and y axis, and max represents the higher bounds - in other words, they represent the top left and bottom right corners. In order to tell whether two AABB shapes are intersecting you will need to have a basic understanding of the Separating Axis Theorem (SAT).

Here's a quick test taken from Real-Time Collision Detection by Christer Ericson, which makes use of the SAT:

bool AABBvsAABB( AABB a, AABB b ) { // Exit with no intersection if found separated along an axis if(a.max.x < b.min.x or a.min.x > b.max.x) return false if(a.max.y < b.min.y or a.min.y > b.max.y) return false // No separating axis found, therefor there is at least one overlapping axis return true }

### Circles

A circle is represented by a radius and point. Here is what your circle structure ought to look like:

struct Circle { float radius Vec position };

Testing for whether or not two circles intersect is very simple: take the radii of the two circles and add them together, then check to see if this sum is greater than the distance between the two circles.

An important optimization to make here is get rid of any need to use the square root operator:

float Distance( Vec2 a, Vec2 b ) { return sqrt( (a.x - b.x)^2 + (a.y - b.y)^2 ) } bool CirclevsCircleUnoptimized( Circle a, Circle b ) { float r = a.radius + b.radius return r < Distance( a.position, b.position ) } bool CirclevsCircleOptimized( Circle a, Circle b ) { float r = a.radius + b.radius r *= r return r < (a.x + b.x)^2 + (a.y + b.y)^2 }

In general multiplication is a much cheaper operation than taking the square root of a value.

## Impulse Resolution

Impulse resolution is a particular type of collision resolution strategy. Collision resolution is the act of taking two objects who are found to be intersecting and modifying them in such a way as to not allow them to remain intersecting.

In general an object within a physics engine has three main degrees of freedom (in two dimensions): movement in the xy plane and rotation. In this article we implicitly restrict rotation and use just AABBs and Circles, so the only degree of freedom we really need to consider is movement along the xy plane.

By resolving detected collisions we place a restriction upon movement such that objects cannot remain intersecting one another. The idea behind impulse resolution is to use an impulse (instantaneous change in velocity) to separate objects found colliding. In order to do this the mass, position, and velocity of each object must be taken into account somehow: we want large objects colliding with smaller ones to move a little bit during collision, and to send the small objects flying away. We also want objects with infinite mass to not move at all.

In order to achieve such effects and follow along with natural intuition of how objects behave we'll use rigid bodies and a fair bit of math. A rigid body is just a shape defined by the user (that is, by you, the developer) that is implicitly defined to be non-deformable. Both AABBs and Circles in this article are non-deformable, and will always be either an AABB or Circle. No squashing or stretching allowed.

Working with rigid bodies allows for a lot of math and derivations to be heavily simplified. This is why rigid bodies are commonly used in game simulations, and why we'll be using them in this article.

### Our Objects Collided - Now What?

Assuming we have two shapes found to be intersecting, how does one actually separate the two? Lets assume our collision detection provided us with two important pieces of information:

• Collision normal
• Penetration depth

In order to apply an impulse to both objects and move them apart, we need to know what direction to push them and by how much. The collision normal is the direction in which the impulse will be applied. The penetration depth (along with some other things) determine how large of an impulse will be used. This means the only value that needs to be solved for is the magnitude of our impulse.

Now let's go on a long trek to discover how we can solve for this impulse magnitude. We'll start with our two objects that have been found to be intersecting:

$V^{AB} = V^B - V^A$ Note that in order to create a vector from position A to position B, you must do: . $$V^{AB}$$ is the relative velocity from A to B. This equation ought to be expressed in terms of the collision normal $$n$$ - that is, we would like to know the relative velocity from A to B along the collision normal's direction:

$V^{AB} \cdot n = (V^B - V^A) \cdot n$

We are now making use of the dot product. The dot product is simple; it is the sum of component-wise products:

$V_1 = \begin{bmatrix}x_1 \\y_1\end{bmatrix}, V_2 = \begin{bmatrix}x_2 \\y_2\end{bmatrix} \\ V_1 \cdot V_2 = x_1 * x_2 + y_2 * y_2$

The next step is to introduce what is called the coefficient of restitution. Restitution is a term that means elasticity, or bounciness. Each object in your physics engine will have a restitution represented as a decimal value. However only one decimal value will be used during the impulse calculation.

To decide what restitution to use (denoted by $$e$$ for epsilon), you should always use the lowest restitution involved in the collision for intuitive results:

// Given two objects A and B e = min( A.restitution, B.restitution )

Once $$e$$ is acquired we can place it into our equation solving for the impulse magnitude.

Newton's Law of Restitution states the following:

$V' = e * V$

All this is saying is that the velocity after a collision is equal to the velocity before it, multiplied by some constant. This constant represents a "bounce factor". Knowing this, it becomes fairly simple to integrate restitution into our current derivation:

$V^{AB} \cdot n = -e * (V^B - V^A) \cdot n$

Notice how we introduced a negative sign here. In Newton's Law of Restitution, $$V'$$, the resulting vector after the bounce, is actually going in the opposite direction of V. So how do we represent opposite directions in our derivation? Introduce a negative sign.

So far so good. Now we need to be able to express these velocities while under the influence of an impulse. Here's a simple equation for modify a vector by some impulse scalar $$j$$ along a specific direction $$n$$:

$V' = V + j * n$

Hopefully the above equation makes sense, as it is very important to understand. We have a unit vector $$n$$ which represents a direction. We have a scalar $$j$$ which represents how long our $$n$$ vector will be. We then add our scaled $$n$$ vector to $$V$$ to result in $$V'$$. This is just adding one vector onto another, and we can use this small equation to apply an impulse of one vector to another.

There's a little more work to be done here. Formally, an impulse is defined as a change in momentum. Momentum is . Knowing this, we can represent an impulse as it is formally defined like so:

$Impulse = mass * Velocity \\ Velocity = \frac{Impulse}{mass} \therefore V' = V + \frac{j * n}{mass}$

The three dots in a small triangle ($$\therefore$$) can be read as "therefore". It is used to show that the thing beforehand can be used to conclude that whatever comes next is true.

Good progress has been made so far! However we need to be able to express an impulse using $$j$$ in terms of two different objects. During a collision with object A and B, A will be pushed in the opposite direction of B:

$V'^A = V^A + \frac{j * n}{mass^A} \\ V'^B = V^B - \frac{j * n}{mass^B}$

These two equations will push A away from B along the direction unit vector $$n$$ by impulse scalar (magnitude of $$n$$) $$j$$.

All that is now required is to merge Equations 8 and 5. Our resulting equation will look something like this:

$(V^A - V^V + \frac{j * n}{mass^A} + \frac{j * n}{mass^B}) * n = -e * (V^B - V^A) \cdot n \\ \therefore \\ (V^A - V^V + \frac{j * n}{mass^A} + \frac{j * n}{mass^B}) * n + e * (V^B - V^A) \cdot n = 0$

If you recall, the original goal was to isolate our magnitude. This is because we know what direction to resolve the collision in (assumed given by the collision detection), and only have left to solve the magnitude of this direction. The magnitude which is an unknown in our case is $$j$$; we must isolate $$j$$ and solve for it.

$(V^B - V^A) \cdot n + j * (\frac{j * n}{mass^A} + \frac{j * n}{mass^B}) * n + e * (V^B - V^A) \cdot n = 0 \\ \therefore \\ (1 + e)((V^B - V^A) \cdot n) + j * (\frac{j * n}{mass^A} + \frac{j * n}{mass^B}) * n = 0 \\ \therefore \\ j = \frac{-(1 + e)((V^B - V^A) \cdot n)}{\frac{1}{mass^A} + \frac{1}{mass^B}}$

Whew! That was a fair bit of math! It is all over for now though. It's important to notice that in the final version of Equation 10 we have $$j$$ on the left (our magnitude) and everything on the right is all known. This means we can write a few lines of code to solve for our impulse scalar $$j$$. And boy is the code a lot more readable than mathematical notation!

void ResolveCollision( Object A, Object B ) { // Calculate relative velocity Vec2 rv = B.velocity - A.velocity // Calculate relative velocity in terms of the normal direction float velAlongNormal = DotProduct( rv, normal ) // Do not resolve if velocities are separating if(velAlongNormal > 0) return; // Calculate restitution float e = min( A.restitution, B.restitution) // Calculate impulse scalar float j = -(1 + e) * velAlongNormal j /= 1 / A.mass + 1 / B.mass // Apply impulse Vec2 impulse = j * normal A.velocity -= 1 / A.mass * impulse B.velocity += 1 / B.mass * impulse }

There are a few key things to note in the above code sample. The first thing is the check on Line 10, . This check is very important; it makes sure that you only resolve a collision if the objects are moving towards each other.

If objects are moving away from one another we want to do nothing. This will prevent objects that shouldn't actually be considered colliding from resolving away from one another. This is important for creating a simulation that follows human intuition on what should happen during object interaction.

The second thing to note is that inverse mass is computed multiple times for no reason. It is best to just store your inverse mass within each object and pre-compute it one time:

A.inv_mass = 1 / A.mass
Many physics engines do not actually store raw mass. Physics engines often times store inverse mass and inverse mass alone. It just so happens that most math involving mass is in the form of .

The last thing to note is that we intelligently distribute our impulse scalar $$j$$ over the two objects. We want small objects to bounce off of big objects with a large portion of $$j$$, and the big objects to have their velocities modified by a very small portion of $$j$$.

In order to do this you could do:

float mass_sum = A.mass + B.mass float ratio = A.mass / mass_sum A.velocity -= ratio * impulse ratio = B.mass / mass_sum B.velocity += ratio * impulse

It is important to realize that the above code is equivalent to the sample function from before. Like stated before, inverse masses are quite useful in a physics engine.

### Sinking Objects

If we go ahead and use the code we have so far, objects will run into each other and bounce off. This is great, although what happens if one of the objects has infinite mass? Well we need a good way to represent infinite mass within our simulation.

I suggest using zero as infinite mass - although if we try to compute inverse mass of an object with zero we will have a division by zero. The workaround to this is to do the following when computing inverse mass:

if(A.mass == 0) A.inv_mass = 0 else A.inv_mass = 1 / A.mass

A value of zero will result in proper calculations during impulse resolution. This is still okay. The problem of sinking objects arises when something starts sinking into another object due to gravity. Perhaps something with low restitution hits a wall with infinite mass and begins to sink.

This sinking is due to floating point errors. During each floating point calculation a small floating point error is introduced due to hardware. (For more information, Google [Floating point error IEEE754].) Over time this error accumulates in positional error, causing objects to sink into one another.

In order to correct this error it must be accounted for. To correct this positional error I will show you a method called linear projection. Linear projection reduces the penetration of two objects by a small percentage, and this is performed after the impulse is applied. Positional correction is very simple: move each object along the collision normal $$n$$ by a percentage of the penetration depth:

void PositionalCorrection( Object A, Object B ) { const float percent = 0.2 // usually 20% to 80% Vec2 correction = penetrationDepth / (A.inv_mass + B.inv_mass)) * percent * n A.position -= A.inv_mass * correction B.position += B.inv_mass * correction }

Note that we scale the by the total mass of the system. This will give a positional correction proportional to how much mass we are dealing with. Small objects push away faster than heavier objects.

There is a slight problem with this implementation: if we are always resolving our positional error then objects will jitter back and forth while they rest upon one another. In order to prevent this some slack must be given. We only perform positional correction if the penetration is above some arbitrary threshold, referred to as "slop":

void PositionalCorrection( Object A, Object B ) { const float percent = 0.2 // usually 20% to 80% const float slop = 0.01 // usually 0.01 to 0.1 Vec2 correction = max( penetration - k_slop, 0.0f ) / (A.inv_mass + B.inv_mass)) * percent * n A.position -= A.inv_mass * correction B.position += B.inv_mass * correction }

This allows objects to penetrate ever so slightly without the position correction kicking in.

## Simple Manifold Generation

The last topic to cover in this article is simple manifold generation. A manifold in mathematical terms is something along the lines of "a collection of points that represents an area in space". However, when I refer to the term manifold I am referring to a small object that contains information about a collision between two objects.

Here is a typical manifold setup:

struct Manifold { Object *A; Object *B; float penetration; Vec2 normal; };

During collision detection, both penetration and the collision normal should be computed. In order to find this info the original collision detection algorithms from the top of this article must be extended.

### Circle vs Circle

Lets start with the simplest collision algorithm: Circle vs Circle. This test is mostly trivial. Can you imagine what the direction to resolve the collision will be in? It is the vector from Circle A to Circle B. This can be obtained by subtracting B's position from A's.

Penetration depth is related to the Circles' radii and distance from one another. The overlap of the Circles can be computed by subtracting the summed radii by the distance from each object.

Here is a full sample algorithm for generating the manifold of a Circle vs Circle collision:

bool CirclevsCircle( Manifold *m ) { // Setup a couple pointers to each object Object *A = m->A; Object *B = m->B; // Vector from A to B Vec2 n = B->pos - A->pos float r = A->radius + B->radius r *= r if(n.LengthSquared( ) > r) return false // Circles have collided, now compute manifold float d = n.Length( ) // perform actual sqrt // If distance between circles is not zero if(d != 0) { // Distance is difference between radius and distance m->penetration = r - d // Utilize our d since we performed sqrt on it already within Length( ) // Points from A to B, and is a unit vector c->normal = t / d return true } // Circles are on same position else { // Choose random (but consistent) values c->penetration = A->radius c->normal = Vec( 1, 0 ) return true } }

The most notable things here are: we do not perform any square roots until necessary (objects are found to be colliding), and we check to make sure the circles are not on the same exact position. If they are on the same position our distance would be zero, and we must avoid division by zero when we compute .

### AABB vs AABB

The AABB to AABB test is a little more complex than Circle vs Circle. The collision normal will not be the vector from A to B, but will be a face normal. An AABB is a box with four faces. Each face has a normal. This normal represents a unit vector that is perpendicular to the face.

Examine the general equation of a line in 2D:

$ax + by + c = 0 \\ normal = \begin{bmatrix}a \\b\end{bmatrix}$

In the above equation, and are the normal vector for a line, and the vector is assumed to be normalized (length of vector is zero). Again, our collision normal (direction to resolve the collision) will be in the direction of one of the face normals.

Do you know what represents in the general equation of a line? is the distance from the origin. This is very useful for testing to see if a point is on one side of a line or another, as you will see in the next article.

Now all that is needed is to find out which face is colliding on one of the objects with the other object, and we have our normal. However sometimes multiple faces of two AABBs can intersect, such as when two corners intersect each other. This means we must find the axis of least penetration.

Here is a full algorithm for AABB to AABB manifold generation and collision detection:

bool AABBvsAABB( Manifold *m ) { // Setup a couple pointers to each object Object *A = m->A Object *B = m->B // Vector from A to B Vec2 n = B->pos - A->pos AABB abox = A->aabb AABB bbox = B->aabb // Calculate half extents along x axis for each object float a_extent = (abox.max.x - abox.min.x) / 2 float b_extent = (bbox.max.x - bbox.min.x) / 2 // Calculate overlap on x axis float x_overlap = a_extent + b_extent - abs( n.x ) // SAT test on x axis if(x_overlap > 0) { // Calculate half extents along x axis for each object float a_extent = (abox.max.y - abox.min.y) / 2 float b_extent = (bbox.max.y - bbox.min.y) / 2 // Calculate overlap on y axis float y_overlap = a_extent + b_extent - abs( n.y ) // SAT test on y axis if(y_overlap > 0) { // Find out which axis is axis of least penetration if(x_overlap > y_overlap) { // Point towards B knowing that n points from A to B if(n.x < 0) m->normal = Vec2( -1, 0 ) else m->normal = Vec2( 0, 0 ) m->penetration = x_overlap return true } else { // Point toward B knowing that n points from A to B if(n.y < 0) m->normal = Vec2( 0, -1 ) else m->normal = Vec2( 0, 1 ) m->penetration = y_overlap return true } } } }

### Circle vs AABB

The last test I will cover is the Circle vs AABB test. The idea here is to compute the closest point on the AABB to the Circle; from there the test devolves into something similar to the Circle vs Circle test. Once the closest point is computed and a collision is detected the normal is the direction of the closest point to the circle's center. The penetration depth is the difference between the distance of the closest point to the circle and the radius of the circle.

There is one tricky special case; if the circle's center is within the AABB then the center of the circle needs to be clipped to the closest edge of the AABB, and the normal needs to be flipped.

bool AABBvsCircle( Manifold *m ) { // Setup a couple pointers to each object Object *A = m->A Object *B = m->B // Vector from A to B Vec2 n = B->pos - A->pos // Closest point on A to center of B Vec2 closest = n // Calculate half extents along each axis float x_extent = (A->aabb.max.x - A->aabb.min.x) / 2 float y_extent = (A->aabb.max.y - A->aabb.min.y) / 2 // Clamp point to edges of the AABB closest.x = Clamp( -x_extent, x_extent, closest.x ) closest.y = Clamp( -y_extent, y_extent, closest.y ) bool inside = false // Circle is inside the AABB, so we need to clamp the circle's center // to the closest edge if(n == closest) { inside = true // Find closest axis if(abs( n.x ) > abs( n.y )) { // Clamp to closest extent if(closest.x > 0) closest.x = x_extent else closest.x = -x_extent } // y axis is shorter else { // Clamp to closest extent if(closest.y > 0) closest.y = y_extent else closest.y = -y_extent } } Vec2 normal = n - closest real d = normal.LengthSquared( ) real r = B->radius // Early out of the radius is shorter than distance to closest point and // Circle not inside the AABB if(d > r * r && !inside) return false // Avoided sqrt until we needed d = sqrt( d ) // Collision normal needs to be flipped to point outside if circle was // inside the AABB if(inside) { m->normal = -n m->penetration = r - d } else { m->normal = n m->penetration = r - d } return true }

## Conclusion

Hopefully by now you've learned a thing or two about physics simulation. This tutorial is enough to let you set up a simple custom physics engine made entirely from scratch. In the next part, we will be covering all the necessary extensions that all physics engines require, including:

• Contact pair sorting and culling