instantaneous collision resolution

This page tackles the problem of how to resolve the collision of two rigid bodies. My first encounter with this question was while designing a computer game in 2004. For over three years, I did not come across a satisfactory treatment neither in books, nor in the Internet. During my semester in Kiruna, I finally had the opportunity to thoroughly investigate the problem myself.
The input to the problem are:

Initial states (0. and 1. order) of the bodies at the moment of contact,
point, and normal of contact.

Without loss of generality, we assume:

The center of total mass is at position 0, and
the total linear momentum (impulse) is the null vector.

We are interested in

resolving the collision within a single instance of time
the velocity, and the angular velocity after the moment of contact.

We derive a mathematical solution to the problem step by step. To perform simulations and to visualize the objects in motion, we utilize a computer program, which is available for download. The implementation uses double precision, and simple Euler integration to solve the differential equations of rotation.

Rigid Body Collision Resolution (C++) *	ribocore.zip	230 kB
Download this article	rigid_body_collis...	320 kB

* The program checks the math, and produces the animations on this website.

It might be the warriors who get the glory,
but it's the engineers who build societies.
B'Elanna Torres

Towards a mathematical solution

This section suggests forces and torques that - when applied to the bodies during collision - will preserve total linear momentum as well as total angular momentum. In this model, we allow the bodies to intrude each other. However, the forces and torques applied during collision make the bodies separate again.
In the illustrations, the ellipsoids visualize the inertia tensors of the bodies. The ellipsoids are not identical to the exterior of the body.

billiard balls
baseball
stuntmen

Most likely, the following notions are familiar to the reader.

The rigid bodies are enumerated by

. Each body has

Recall the assumptions stated above. The center of total mass is at position 0, i.e.

and the total linear momentum (impulse) is the null vector:

We treat the collection of bodies as a closed system. Within this system, the total linear momentum, and the total angular momentum are perserved, i.e. constant vectors, at all times.

The total linear momentum is the sum of the linear momentums of each body

The total angular momentum is the sum of the angular momentums of each body with respect to the center of total mass

Since in our world, the total momentums are constant, the derivative with respect to time is the null vector. I.e. we demand

In the above derivation, we assume that all quantities depend on time, except mass and the inertia tensors. The introduced variables represent the following:

while

is the angular velocity in matrix form

We are interested in forces and torques applied to the bodies during the period of contact.

We specialize the formulas on two rigid bodies, i.e. we assume

. We introduce a few more symbols that characterize the contact between the two bodies:

If bodies are intruding each other, then let

be proportional to the distance of intrusion. If the bodies are not in contact, we set

. We suggest to apply the following forces and torques

We show that under these conditions both, the total linear momentum and the total angular momentum are preserved. Conservation of total linear momentum follows because

The total angular momentum does not change because

These formulas produce the animations shown earlier.

Man muss wissen, bis wohin man zu weit gehen kann.
Jean Cocteau

There is no time for procrastination.

In the introduction, we have stated the ultimate goal of this project: Given the linear and angular velocities just prior to the contact, we would like to have a formula that computes the linear and angular velocities just after the contact. This section suggests a solution to this problem. The solution is consistent with the definitions and derivations carried out earlier.
In the illustrations, the ellipsoids visualize the inertia tensors of the bodies. The ellipsoids are not identical to the exterior of the body.

billiard balls
baseball
stuntmen

I hope you enjoy the upcoming formulas as much as I do.

First, we introduce a few more symbols

We suggest: The linear and angular velocities just prior to the contact relate to the linear and angular velocities just after the contact in the following way

We show that under these conditions both, the total linear momentum and the total angular momentum are preserved. Conservation of total linear momentum follows because

The total angular momentum after collision is

From that, we subtract the total angular momentum before collision

and yield a difference of

At this point,

is the only remaining unknown. As we will show next,

is uniquely determined under the assumption that the two bodies collide elastically.

In an elastic collision, the total energy before the collision equals the total energy after the collision, i.e.

Previously, we have suggested the relations

where

Now, determining

so that the total energy is preserved reduces to solve a quadratic polynomial. Evidently,

is one solution to conserve total energy, which corresponds to "no collision". However, we are interested in

Note, the numerator of the fraction is the scalar product between the velocity of contact and the surface normal

, where

Basically, we are done at this point. The latter value of

generates the three animations you see above. We conclude with some remarks:

white spheres indicate location of previous contact

The forces and torques stated in the first section produce the same outcome if the collision process is modelled stiffly enough.
Taking a fraction (for instance 95%) of will damp the collision, which might result in more realistic animations. Different fractions may correspond to different (non rigid) materials. The total linear momentum as well as the total angular momentum are preserved nevertheless.
More general, both momentums are conserved for any vector , which is a convenient gateway for adapting the formula to less-elastic collisions.

To produce a final demonstration, we use the collision detection library RAPID.

Wenn du Zeit gewinnen willst,
musst du Zeit verlieren.
Jean-Jacques Rousseau

Software and further reading

An open source physics engine that I have checked out myself, in the Newton Game Dynamics by Julio Jerez and Alain Suero. The software comes with a demo that explores the broad spectrum of applications. I find the results to be efficient, stable, and realistic. Moreover, the authors C++ style is just perfect.
The Open Dynamics Engine (ODE) originally by Russell Smith seems to be a popular C++ library for developers in gaming, but also in robotics. The library is suitable for simulating articulated rigid body structures. However, ODE emphasizes speed and stability over physical accuracy.
Computer graphics researchers from Stanford University came up with a successful strategy to animate and collide innumerable rigid bodies. The paper Nonconvex Rigid Bodies with Stacking describes how to achieve breathtaking results. An open source implementation of their method is written by Danny Chapman. To resolve a multi-contact multi-body collision, the authors propose to simulate the scenario a number of times: The sequence of pairwise collisions are resolved with different ordering and the outcomes are averaged.

Simulation of Ideal Gas Particles

Steven Hochstadt (sdh78)
Eric Chen (cc459)

Introduction:

Our project was simulating an ideal gas on a VGA monitor. The idea is to have the monitor display moving white pixels that represent particles in an ideal gas and perform the proper calculations when they collide with each other or with the edge of the screen. To achieve this goal, we built a state machine on the DE2 development board that keeps track of the particle coordinates and velocities in m4k blocks, and updated each particles coordinates and redrew them on the screen once per frame. Since an ideal gas is generally considered a large mass of free moving particles, a fluid simulation based on an ideal gas should have a large density of simulated particles. This proves to be a difficult but attractive challenge that we determined would be a rewarding experience to pursue. Also, the resulting visual phenomenon from such a simulation is interesting to study and can provide an educational insight into the general fluid dynamics of an ideal gas.

High Level Design:

Before providing details about our design, we will present a brief outline on the key properties of an ideal gas. The first property relies on the law of conservation of momentum, which states that the total net momentum in a closed system must remain constant unless imposed upon by a source outside the system. In an ideal gas, we can assume that all particles have the same mass, which means that sum of all particle velocities in a closed gaseous system will remain a constant. This also means that in any given collision between two particles, the net velocity of the two particles must be the same before and after the collision. There are two types of collisions: elastic, and inelastic. In an elastic collision, all energy in the system is conserved as momentum. In an inelastic collision, some of this energy may be dissipated as heat, which results in less energy in the momentum and an overall decrease in the average magnitude of velocity. Collisions in a gas can be considered to be approximately elastic, which simplifies the math involved. The change in velocity between two particles, (assumed to be spherical balls), in an elastic collision between can be calculated from the projection of the relative velocities onto the normal of the collision site:

The vector vi is the change in the velocity vector of ball i, rij is the relative position vector between balls i and j, and vij is the relative velocity vector between balls i and j. Intuitively, it makes sense for the collision to be a projection along the normal of the collision site because the net force on each ball would be along that normal.
For our design, we weighed the advantages and disadvantages of developing a hardware solution versus a software solution. Due to having to keep track of all the particles (or balls) in memory, the software solution sounded attractive because we would not have to worry about the memory overhead. However, in order to find the collisions, the coordinates of all balls have to be compared with the coordinates all other balls. To keep track of just 500 balls, we have to make (500*501/2 comparisons) * (60 frames per second) = 7,515,000 comparisons per second. Since every instruction in a software solution uses at least one clock cycle, a single comparison may take a few clock cycles, which means the number of instructions required per second will increase to a lot more than 7 million. In hardware, we can use a state machine that performs multiple tasks at the same time, which saves a lot of clock cycles. Another advantage of hardware is the ability to pipeline tasks that take multiple clock cycles. Overall, the hardware solution allows for many more tasks to be completed within one clock cycle than the software solution. Therefore, our entire design is written in verilog.
Our main design is a statemachine that is controlled by the vga control clk (VGA_CTRL_CLK):

The four main states are checkboundaries, checkcollisions, collideotherball, removeballs, and drawballs. The removeballs and drawballs states are the most straightforward. The removeballs state clears the SRAM and transitions to the drawballs state, which loops through all the m4k block addresses and draws each ball to its new location in SRAM. There is a oneshot state that waits until the VGA is not being used before transitioning to the removeballs state. The checkboundaries state checks whether each ball has reached the edge of the screen and changes the ball�s position and velocity accordingly. If the ball is not at a boundary, the statemachine transitions to the checkcollisions state, which is the most interesting state. For each ball, the checkcollisions state loops through all the other ball coordinates and if it finds a collision, it updates the current ball�s position and velocity accordingly and transitions to the collideotherball state, which updates the update�s the other ball�s position accordingly. The checkcollisions and collideotherball states utilize a helper module called collisioncalculator, which performs a pipelined version of the projection equation given earlier.

Program/hardware design:

The entire project consists of the state machine, the collision calculator, and the m4kblocks. As stated before, our main goal is to simulate as many particles as we can. In order to do this, we had to tweak the different components so that as much as possible can be completed in a given number of clock cycles. One major strategy for freeing up clock cycles is to perform pipelining. Due to the enormous number of comparisons, the checkcollsions state was our focus for the pipelining. There were two major areas that were pipelined. The first area is the m4kblocks. Since reading and writing take 2 clock cycles each, we can make a two stage pipeline for reading and writing. The other major area is the collisioncalculator module. Due to the complexity of the projection calculation within the module, the propagation time for performing the calculation all at once limited our ability to increase the clock frequency. Changing the calculation into a four stage pipeline solved this problem enabling us to double the clock frequency from 25.2MHZ to 50.4MHZ. The four stage pipeline and the two stage pipeline were combined together within the checkcollisions state as a 6 stage pipeline. Unfortunately, we did not gain as many clock cycles as expected. Raising the frequency to any multiple of the vga control clock higher than 50.4MHZ results in a timing analysis failure.

Pipeline

Another area for optimization was within the collision calculation itself, which involves division by the magnitude squared of the relative position vector. In general, division is a resource intensive calculation but in our case, we can estimate the relative position vector. If we assume all balls have the same radius, then when they collide, they will always be approximately the same distance apart, which means that if we choose the radius to be a power of 2, we can divide by just performing a binary shift. The reason for them being approximately the same distance apart is because a collision is detected by checking if the distance between two particles is below a threshold. Unfortunately, this approximation turned out to be too crude to provide a decent simulation. If the particles are moving extremely fast relative to their radius, then the approximation becomes extremely inaccurate due to the distance in a collision being well below the specified threshold. On the other hand, our approximation of the square root for the distance calculation worked relatively well. The approximation of the sqrt(x^2 + y^2) where x > y is (x + �*y) and comes from the Taylor series expansion.
Another interesting challenge was to ensure that all collisions were actually detected. In order to do this, we could not allow any particle to move at a speed that is grater than its own radius per frame. Two particles on a collision course each traveling fasdter than 1 radius per frame may end up never falling below the collision threshold in any frame. Since our particles have a diameter of one pixel, they must move at fractions of a pixel per frame, which means the position and velocity representations within the m4k blocks have to have a fixed point representation that is more precise than the screen pixel coordinates. Along with this, we also want to be able to easily perform mathematics on the position and velocity representations. Therefore, we decided to have the x-position, y-position, x-velocity, and y-velocity of each ball to be represented by an 18 bit signed 10.8 fixed point number. This means we have one sign bit, 9 whole number bits, and 8 fractional bits. The 8 fractional bits means that a speed of half a pixel per second would be 2^7 = 128, which means there are 128 different speeds for each ball.

Results:

This project has been assessed according to accuracy, speed, and how well we met our goals. Unfortunately, our accuracy is much lower than expected due to the inaccuracy in the distance approximation in collisions. This inaccuracy caused the particles to lose energy over time due to dividing by a distance that is too big. The end result is that the elastic collisions ended up not really being elastic. It also causes particles to sometimes get �stuck� and clump together due to never being able to escape to a distance beyond the collision threshold. However, there are many methods for improving the accuracy and clumping including implementing a hardware divider or shifting the dot product by different amounts according what range the distance was in. We avoided the divider because of the complexity and timing factors involves but we did try the latter solution and, though still inaccurate, it resulted in a significant improvement. There is also a mysterious bug with particles on the top half of the screen being darker than particles on the bottom half. If we lowered the number of particles to a very small number (like 30), the ones on the top half of the screen would not show up at all. We believe this may be because clearing the old ball locations by clearing the entire screen does not fit within the small percentage of time when the VGA controller is not drawing. Overall, the visual result still has an educational value. We can change the speeds and locations of the particles and see how they interact on the screen. This is a simple change inside two different mif files for the velocity and the position. A small matlab program for generating random ball locations for pasting into the mif files is attatched here.

Conclusion:

Overall, the project is functional. It can simulate over 500 balls on the screen and they bounce off each other properly in most cases. However, the resulting simulation was not as accurate as we had hoped and we expected to simulate over 1000 balls. Having said that, if we were to do this project over again we would focus more on parallelization instead of serializing the process at a high frequency. We initially avoided this approach because of the complexity of communication between the parallelized components. We also thought that we would be able to use a higher clock frequency than just 50.4MHZ. Visually however, the project turned out to be an interesting simulation to explore.