This page is aimed at developers of games or physics based real-time simulations. Some background in physics and mathematics is recommended. The core code examples on this page are written in very simple python, which requires little background in software development to understand.
Some of the graphical demonstrations require additional extension modules. Where required these packages are mentioned again.
- PyOpenGL Bindings to OpenGL and associated rendering packages
- Matplotlib Matlab like plotting utilities
The following definitions are used within this article
- Rigid Body is a fixed size body that cannot be deformed. A rigid body is typically represented as a mesh with a position and an orientation.
- State In the study of dynamical systems, a physical system is typically modeled as a differential equation. The variable in such an equation is often called the state of the system. Throughout this article 'state' is assumed to be a vector. For two dimensional games this may include x and y position and velocities.
- State Transition Equation describes a differential equation which calculates the rate of change of the state vector, from the state and input vector, at the current time. Integrating the result of the state transition equation over a short timestep allows prediction of the state in the future.
Physics Simulation Loop
The following steps outline a basic approach to physics simulation.
while not finished // stage 1 for body in bodies body.applyForcesAndTorques() // stage 2 for body in bodies body.update( timeStep ) // stage 3 detectAndResolveCollisions( bodies )
The first part of the loop involves the application of all relevant forces and torques to the bodies in the simulation. Typically this includes things like:
- user applied forces (for example a space ship would have engine forces)
The second stage involves the updating the internal state of each rigid body. Typically this is described as a differential equation. At this stage the forces from the first step, together with the mass and moment of inertia of the body are used to calculate the new velocities, which are in turn used to calculate the linear and angular positions.
One of the most common approaches for preforming very simple physics updates is the Euler Integrator. See below for a discussion.
Collision detection and response are used to model situations in which multiple rigid bodies begin to interact. Collisions can be seen as mathematically as discontinuities in the differential equations describing the motion of the bodies. This means that the state transition equations can no longer fully predict the responses of the bodies. To ensure that our simulation is a reasonable approximation we need to identify when and where collisions happen, and to modify the response of the bodies in such situations.
Example 1 - Simulation
The following example demonstrates the basic principles of collision free simulation. The following code demonstrates how velocity and position of a ball change under the influence of gravity.
The equations for this system can be derived from Newtons laws of motion, assuming that acceleration is constant over a given timestep of length and the ball has mass m.
Discrete time displacement updates:
Discrete time velocity update
The state of a ball can be represented by a matrix containing the x and y components of displacement and velocity.
The state transition equation when written in full:
# ball.py from numpy import * class Ball: def __init__(self): self.fx = 0 self.fy = 0 # current state # vx,vy,x,y self.x = array([0,0,0,0]) self.m = 1.0 def clearForceAccumulator(self): self.fx = 0 self.fy = 0 def applyForce(self, f): fx,fy = f self.fx += fx self.fy += fy def stateTransitionEquation(self,x): # rate of change of state (xdot) # vx_dot,vy_dot,x_dot,y_dot xdot = array([0.0,0.0,0.0,0.0]) xdot = self.fx / self.m # derivative of v is a, a = f/m xdot = self.fy / self.m # derivative of v is a, a = f/m xdot = x # derivative of x is v (from last) xdot = x # derivative of x is v (from last) return xdot def eulerIntegrate(self,x,xdot,dt): # state at the next timestep xnext = x + xdot*dt return xnext def update(self,dt): xdot = self.stateTransitionEquation(self.x) xnext = self.eulerIntegrate(self.x,xdot,dt) self.x = xnext
A basic test application for the ball physics code might look like this:
from ball import * b = Ball() b.x = [0,0,0,10] # stationary, 10 meters up g = (0,-9.81) # gravitational constant simTime = 0 endTime = 1.0 timeStep = 0.01 while simTime < endTime: print simTime, b.x b.clearForceAccumulator() b.applyForce( g ) b.update( timeStep ) simTime += timeStep
Example 2 - Validation of Simulation
Although the results of example 1 can be validated by hand, it can be useful to plot the results, and ensure that they match up with the expected results from the original equation. This can be done with any plotting package, however here we use matplotlib.
from matplotlib.numerix import * import pylab from ball import * b = Ball() b.x = [0,0,0,10] # stationary, 10 meters up g = (0,-9.81) # gravitational constant dataLog =  simTime = 0 endTime = 1.0 timeStep = 0.01 while simTime < endTime: dataLog.append( (simTime, b.x) ) b.clearForceAccumulator() b.applyForce( g ) b.update( timeStep ) simTime += timeStep # convert the state data to a better form for plotting time = array([ log for log in dataLog ]) vx = array([ log for log in dataLog ]) vy = array([ log for log in dataLog ]) x = array([ log for log in dataLog ]) y = array([ log for log in dataLog ]) # clear any existing graphs (required for use in the interpreter) pylab.cla() # plot the position and velocity graphs pylab.plot(time,y,'r.',label='y position (m)') pylab.plot(time,vy,'g.',label='y velocity (m.s-1)') # setup the axes and their labels gca = pylab.gca() gca.set_title('Y position and velocity versus time') gca.set_xlabel('time (seconds)') gca.set_xlim([0,1.0]) gca.set_ylim([-10,10.0]) # plot the legend pylab.legend() # save to file pylab.savefig('results_vs_time.png',dpi=75)
Example 3 - Visualisation
Extending the ball dynamics of the previous examples, we can show a type of 2d fireworks using particle systems. The sample code here has been written using PyOpenGL.
Particles systems make an ideal starting point for examining physics simulation because particles do not rotate, they only translate. Particle systems are typically used as special effects elements in many games and movies. Amongst the best known effects are fireworks, smoke, rain and sparks.
See the NeHe OpenGL Particle System tutorial for more information on getting started.
from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * import sys import random import traceback from ball import * def randomState(): vx = random.random()*2.0-1.0 vy = random.random()*4.0+8.0 return array([vx,vy,0,0]) def randomBalls(n): balls = [ Ball() for idx in range(n)] for ball in balls: ball.x = randomState() ball.color = [ random.random(),random.random(),random.random() ] return balls class ParticleDemo: def __init__(self): self.balls = randomBalls(100) self.g = (0,-9.81) # gravitational constant self.simTime = 0 self.endTime = 1.0 self.timeStep = 0.01 def step(self): """ Do a physics update for one timestep """ for b in self.balls: b.clearForceAccumulator() b.applyForce( self.g ) b.update( self.timeStep ) self.simTime += self.timeStep def render(self): glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glPointSize(4.0) self.step() glBegin(GL_POINTS) for b in self.balls: glColor(*(b.color)) glVertex(b.x,b.x,0.0) glEnd() glutSwapBuffers() def initialiseView(self,w,h): # Select the background clear color glClearColor( 0.0, 0.0, 0.0, 0.0 ) #initialise viewing values glMatrixMode(GL_PROJECTION) glLoadIdentity() glOrtho(-w,w,-h,h,-1.0,1.0) def onKey(self,k,x,y): if ord(k) == 27: # Escape sys.exit(0) elif k == ' ': self.balls = randomBalls(100) else: return glutPostRedisplay() def onIdle(self): glutPostRedisplay() # Global variable for maintaining the state of # the particle system p = ParticleDemo() # glut callback function to redirect to the real # function in the ParticleDemo class. These # functions also contain additional error handling # code to print exception traces because the current # PyOpenGL doesn't print stack traces def idle(): try: p.onIdle() except: traceback.print_exc() raise def key(k,x,y): try: p.onKey(k,x,y) except: traceback.print_exc() raise def display(): try: p.render() except: traceback.print_exc() raise glutInit(sys.argv) glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH) glutInitWindowSize(640, 480) glutCreateWindow("Swin Particle Sim") p.initialiseView(10.0,10.0) glutIdleFunc(idle) glutDisplayFunc(display) glutKeyboardFunc(key) glutMainLoop()
Collision detection is the process of detecting interpenetration of two rigid bodies. There are a wide range of approaches to solving the collision detection problem, many of which are suitable in different circumstances. The accuracy of collision detection is often dependant on:
- The speed of the bodies involved in the collision
- The time between collision tests
- Whether the geometries of the bodies involved in the collision are convex or concave
- The number of bodies that can collide
- How the geometry of the colliding bodies is being managed. Typically three dimensional shapes are approximated using Triangle Meshes, although representations like Constructive Solid Geometry (CSG) can be used where higher precision is required. Two dimensional shapes can be represented either as 2d triangle meshes, or using two dimensional geometric primitives.
Collision detection can be a relatively slow process when a large number of bodies need to be processed - each body needs to be tested against every other body O(n squared). This problem quickly becomes intractable when many triangle meshes need to be tested for collision (each triangle in a given mesh needs to be tested against each triangle in every other mesh). In an effort to reduce the overhead of collision detection a one or more preliminary tests can be used to eliminate some bodies (or parts of bodies) from the collision detection process. These tests include
- Bounding spheres
- Axis Aligned Bounding Boxes
- Object Aligned Bounding Boxes.
It is worth noting that generalised collision detection is still an open area of research - particularly with regard to culling, robustness and speed. Some of the best open source packages for collision detection are provided by Team Gamma at the Department of Computer Science University of North Carolina (UNC). Significant packages include
Collision detection approaches
- At each time step each object is tested for intersection with every other object. If a collision occurs, time is rewound and a binary search is conducted to find the time at which the first collision occurs. This approach can be slow depending on the number of bodies and the frequency of collision.
- A popular alternative is to use the Separating Axis Theorem. The developers of the flash game 'N' have several easy to understand tutorials and diagrams.
Collision detection also has applications in
- Frustrum culling
- Path planning
Example 4 - Static Collisions
The following example shows a simple approach to detecting collisions between static (unmoving) circles and halfplanes. This example demonstrates how the collision detection process is tightly coupled to the geometric representation of the colliding bodies, and that basic mathematical concepts like the dot product and vector projection can be used to write simple and easy to read code. This example also demonstrates the importance of rigourous testing of collision detection routines.
# phyutils.py import numpy import math def magnitude(v): """ Euclidean length of a vector sqrt(x*x+y*y+z*z) """ return math.sqrt(numpy.dot(v,v)) def normalise(v): """ returns a normalised copy of the vector v. undefined if the length of v is zero """ return v / magnitude(v)
from numpy import * import math from physutils import * # defines an infinite plane passing through a point with # an arbitrary normal vector class HalfPlane: def __init__(self, origin, normal): self.origin = origin self.normal = normalise(normal) class Circle: def __init__(self, center, radius): self.center = center self.radius = radius def circleIntersectsHalfPlane(halfPlane,circle): vectorOC = circle.center-halfPlane.origin vectorOCPlaneNormalComponent = dot(vectorOC,halfPlane.normal) return (vectorOCPlaneNormalComponent - circle.radius) < 0 def testCircleHalfPlaneIntersection(): def createCircle(y): return Circle( array([0.0,y]), 1.0 ) # define the ground halfplane to pass through the # origin, 'upwards' is positive y ground = HalfPlane( array([0,0]), array([0,1]) ) # clearly in assert( circleIntersectsHalfPlane(ground, createCircle(2.0)) == False) # boundary conditions assert( circleIntersectsHalfPlane(ground, createCircle(1.001)) == False) assert( circleIntersectsHalfPlane(ground, createCircle(1.0)) == False) assert( circleIntersectsHalfPlane(ground, createCircle(0.009)) == True) # clearly out assert( circleIntersectsHalfPlane(ground, createCircle(0.0)) == True) # potential corner case assert( circleIntersectsHalfPlane(ground, createCircle(-1.0)) == True) if __name__ == "__main__": testCircleHalfPlaneIntersection()
To test the collision code in a more realistic environment, the following code demonstrates how static collisions can be added into an OpenGL test environment
from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * import sys import random import traceback from physutils import * RED = (1.0,0.0,0.0) GREEN = (0.0,1.0,0.0) def circleIntersectsHalfPlane(halfPlane,circle): vectorOC = circle.center-halfPlane.origin vectorOCPlaneNormalComponent = dot(vectorOC,halfPlane.normal) return (vectorOCPlaneNormalComponent - circle.radius) < 0 class Shape: def __init__(self): self.color = WHITE def getColor(self): return self.color def setColor(self,color): print color self.color = color def draw(self): pass # defines an infinite plane passing through a point with # an arbitrary normal vector. We use an additional piece # of domain specific knowledge (the size of the environment) # when drawing the line outward from the origin perpendicular # to the normal class HalfPlane(Shape): def __init__(self, origin, normal, widthHint= 20.0): Shape.__init__(self) self.setColor(RED) self.origin = origin self.normal = normalise(normal) self.perp = array( [self.normal,self.normal ]) self.widthHint = widthHint def draw(self): p1 = self.origin + self.perp*self.widthHint p2 = self.origin - self.perp*self.widthHint glColor(*self.color) glBegin(GL_LINES) glVertex(p1,p1,0.0) glVertex(p2,p2,0.0) glEnd() class Circle(Shape): def __init__(self, center, radius): Shape.__init__(self) self.setColor(GREEN) self.center = center self.radius = radius self.steps = 20 def draw(self): glColor(*self.color) glBegin(GL_LINE_LOOP) for i in range(self.steps): theta = (math.pi*2.0*i)/self.steps x = math.cos( theta ) * self.radius + self.center y = math.sin( theta ) * self.radius + self.center glVertex(x,y,0) glEnd() class CollisionDemo: def __init__(self,circle, plane): self.circle = circle self.plane = plane def render(self): glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glPointSize(4.0) if circleIntersectsHalfPlane(plane,circle): self.circle.setColor(RED) self.plane.setColor(RED) else: self.circle.setColor(GREEN) self.plane.setColor(GREEN) self.circle.draw() self.plane.draw() glutSwapBuffers() def initialiseView(self,w,h): # Select the background clear color glClearColor( 0.0, 0.0, 0.0, 0.0 ) #initialise viewing values glMatrixMode(GL_PROJECTION) glLoadIdentity() glOrtho(-w,w,-h,h,-1.0,1.0) def onKey(self,k,x,y): if ord(k) == 27: # Escape sys.exit(0) else: return def onSpecialKey(self,k,x,y): if k == GLUT_KEY_UP: self.circle.center += array([0.0,0.5]) elif k == GLUT_KEY_DOWN: self.circle.center += array([0.0,-0.5]) elif k == GLUT_KEY_LEFT: self.circle.center += array([-0.5,0.0]) elif k == GLUT_KEY_RIGHT: self.circle.center += array([0.5,0.0]) else: return def onIdle(self): glutPostRedisplay() # Global variable for maintaining the state of # the particle system origin = array([0.0,0.0]) up = array([0.0,1.0]) circle = Circle(origin+up,3.0) plane = HalfPlane(origin,up) demo = CollisionDemo(circle,plane) # glut callback function to redirect to the real # function in the ParticleDemo class. These # functions also contain additional error handling # code to print exception traces because the current # PyOpenGL doesn't print stack traces def idle(): try: demo.onIdle() except: traceback.print_exc() raise def key(k,x,y): try: demo.onKey(k,x,y) except: traceback.print_exc() raise def specialKey(k,x,y): try: demo.onSpecialKey(k,x,y) except: traceback.print_exc() raise def display(): try: demo.render() except: traceback.print_exc() raise glutInit(sys.argv) glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH) glutInitWindowSize(400, 400) glutCreateWindow("SwinBrain Static Collisions") demo.initialiseView(10.0,10.0) glutIdleFunc(idle) glutDisplayFunc(display) glutKeyboardFunc(key) glutSpecialFunc(specialKey) glutMainLoop()
The image below shows the results of the OpenGL collision demonstration for two different cases. Red shapes denote collision states, whilst the green shapes denote collision free states.
Example 5 - Dynamic Collisions
Accurately detecting dynamic collisions is significantly more challenging than detecting static collisions. To perform ideal collision detection time would need to be treated as a continuous variable, however limited resources mean that a discrete approximation to time needs to be used. The choice of a discretization method can have a significant impact on the time required for simulation - which has an impact on the number of physics bodies which can be managed in real time.
Additionally because a discretization of time is used, it is possible that very fast moving bodies (for example, bullets) may not be accurately tested for collision. One approach to minimising this type of effect is to be able to predict the maximum displacements of any physics body for a timestep. If collisions are likely to be missed then it may be necessary to reduce the timestep.
This example demonstrates how to detect the time (and consequently the location) of all the bodies involved in the collision. This is done using a binary search approach to identify the first time that is collision free (to some reasonable limit, epsilon).
More information is coming on this topic...
Collision response is the process of modifying the states of rigid bodies involved in collisions. Generally this is carried out in two stages - separation of penetrating bodies, and updating the state vectors of the bodies to reflect the collision.
Collision response is heavily influenced by the material properties of the colliding objects. This can lead to a range of situations, the two extremes are elastic collisions (conservation of kinetic energy) and inelastic collisions (conservation of momentum). The responses generated by different types of materials can be approximated by using the velocity equations based on the Coefficient of restitution
There are several major approaches to resolving collisions, which each have advantages and disadvantages, these include
- Impulse methods assume collision occurs over a very short period of time, and that the result of the collision is a change in velocity for each of the body.
- Penalty Forces construct a set of contact points, interpenetrating bodies are forced apart using spring like equations through the contact points.
Although individual up until now only rigid bodies have been considered it is possible to produce effects where multiple rigid bodies are connected together. If this level of physics is required it may be worth considering third party physics engines like ODE.
The following are a list of the most common integrators, and some guidelines for their use:
- Euler Integrator Although Euler integrators are simple to implement, they are only appropriate for simulating equations which are not stiff, when short timesteps are used and when cumulative errors are not significant. Gaffer On Games has a good explanation of why not to use Euler integrators.
- Runge-Kutta RK4 (Fourth order Runge-Kutta integrator) should be the integrator of choice for most types of simulations systems.
- Verlet and Leapfrog integrators are most commonly used for use in systems with large numbers of interconnected particles (good examples include cloth simulation and ragdoll physics), where the computational overhead of a higher order integrator or smaller timesteps is not feasible.
A wide range of open and closed source physics libraries are used in game development. A partial list of the most relevant engines can be found at gamedev.net physic forum
Good starting points for experimentation into game physics in two dimensions include
- Artilery simulations like Scorched Earth.
- Pinball simulations
- Any game played on a billiards table. This includes 8-ball (pool), 9-ball, billiards and snooker.
Nearly all of the physics tutorials currently available refer back to the Baraff and Witkins papers which are a valuable reference, but also quite readable.
Chris Heckers Rigid Body Dynamics Information has a checklist which describes the necessary steps for working towards a full physics engine.
Ragdoll Physics is the simulation of a constrained system of rigid bodies, typically bonemeshes.
Jeff Lander's Game Developer Magazine articles are good demonstrations of how to create basic physics simulations. Of particular interest are the spring systems and cloth simulation.