These tutorials focus mainly on OpenGL, Win32 programming and the ODE physics engine. OpenGL has moved on to great heights and I don't cover the newest features but cover all of the basic concepts you will need with working example programs.

 

Working with the Win32 API is a great way to get to the heart of Windows and is just as relevant today as ever before. Whereas ODE has been marginalized as hardware accelerated physics becomes more common.

 

Games and graphics utilities can be made quickly and easily using game engines like Unity so this and Linux development in general will be the focus of my next tutorials.    

  

 

And off we go...

 

Lets make things a little easier for our first example and just create a single object, a box, that bounces on a flat surface. We can leave the subject of joints (such as a hinge joint) and the finer details to future tutorials. To begin with we need to include the ode header file.


#inlcude <ode/ode.h>

And for the purpose of this example we will use a few global variables.


MATRIX GeomMatrix;

MyObject Object;

dWorldID World;

dSpaceID Space;

dJointGroupID contactgroup;

The first variable called GeomMatrix is a matrix used during the rendering stage to convert the 3x3 rotation matrix and position vector used in ODE to one we can use with OpenGL. I have used my own matrix class in this example but an array of twelve float values will do. I will discuss this further and give you the conversion function later in the tutorial.

 

The variable below called MyObject is simply a structure containing two object IDs. These object IDs are used throughout the ODE library to uniquely reference the many different objects created in the simulation.


#define GEOMSPERBODY 1  // maximum number of geometries per body



struct MyObject

{

    dBodyID Body;  // the dynamics body

    dGeomID Geom[GEOMSPERBODY];  // geometries representing this body

};

The first object ID in the structure is an identifier for this objects dynamics body. The dynamics body of an object includes information about its position, velocity, etc. The second object ID in the structure points to an array of the basic geometries, such as rectangles or spheres, that make up this object. In this example there is only one 'geom' in the group (the box) but most of the objects that you will create will be a collection of many different primitive shapes. Together the dynamics body and its geoms make up a single rigid body within the simulation.

 

The first three variables above are also object IDs for the following types of objects:

 

  • the dWorld object is a container for all the rigid bodies and their joints.
  • the dSpace object is similar to the world container except that it applies to collision instead of dynamics.
  • the dJointGroup object is a collection of joints.

 

You may be wondering why we need the dJointGroup object at all if we aren't going to be joining any objects together. In physics simulations all contacts between rigid bodies (or the ground) are calculated using temporary contact joints. When a box falls on a plane surface up to four contact joints, one for each corner of the box, may be created.

 

Initializing ODE

 

Now that we have covered the global variables and their uses we can move on to initializing the ODE simulation. And to keep things in order let's create a function called InitODE.


void InitODE()

{

    // Create a new, empty world and assign its ID number to World. Most applications will only need one world.

    World = dWorldCreate();



    // Create a new collision space and assign its ID number to Space, passing 0 instead of an existing dSpaceID.

       // There are three different types of collision spaces we could create here depending on the number of objects 

       // in the world but dSimpleSpaceCreate is fine for a small number of objects. If there were more objects we

       // would be using dHashSpaceCreate or dQuadTreeSpaceCreate (look these up in the ODE docs)	

    Space = dSimpleSpaceCreate(0);



    // Create a joint group object and assign its ID number to contactgroup. dJointGroupCreate used to have a 

       // max_size parameter but it is no longer used so we just pass 0 as its argument.

    contactgroup = dJointGroupCreate(0);



    // Create a ground plane in our collision space by passing Space as the first argument to dCreatePlane.

       // The next four parameters are the planes normal (a, b, c) and distance (d) according to the plane

       // equation a*x+b*y+c*z=d and must have length 1

    dCreatePlane(Space, 0, 1, 0, 0);



    // Now we set the gravity vector for our world by passing World as the first argument to dWorldSetGravity.

       // Earth's gravity vector would be (0, -9.81, 0) assuming that +Y is up. I found that a lighter gravity looked

       // more realistic in this case.	

    dWorldSetGravity(World, 0, -1.0, 0);



    // These next two functions control how much error correcting and constraint force mixing occurs in the world.

       // Don't worry about these for now as they are set to the default values and we could happily delete them from

       // this example. Different values, however, can drastically change the behaviour of the objects colliding, so

       // I suggest you look up the full info on them in the ODE docs.

    dWorldSetERP(World, 0.2);

    dWorldSetCFM(World, 1e-5);

   	

    // This function sets the velocity that interpenetrating objects will separate at. The default value is infinity.

    dWorldSetContactMaxCorrectingVel(World, 0.9);

    // This function sets the depth of the surface layer around the world objects. Contacts are allowed to sink into

       // each other up to this depth. Setting it to a small value reduces the amount of jittering between contacting

       // objects, the default value is 0. 	

    dWorldSetContactSurfaceLayer(World, 0.001);



    // To save some CPU time we set the auto disable flag to 1. This means that objects that have come to rest (based

       // on their current linear and angular velocity) will no longer participate in the simulation, unless acted upon

       // by a moving object. If you do not want to use this feature then set the flag to 0. You can also manually enable

       // or disable objects using dBodyEnable and dBodyDisable, see the docs for more info on this.

    dWorldSetAutoDisableFlag(World, 1);



    // This brings us to the end of the world settings, now we have to initialize the objects themselves.



    // Create a new body for our object in the world and get its ID.

    Object.Body = dBodyCreate(World);



    // Next we set the position of the new body

    dBodySetPosition(Object.Body, 0, 10, -5);



    // Here I have set the initial linear velocity to stationary and let gravity do the work, but you can experiment

       // with the velocity vector to change the starting behaviour. You can also set the rotational velocity for the new

       // body using dBodySetAngularVel which takes the same parameters.

    VECTOR tempVect(0.0, 0.0, 0.0);

    dBodySetLinearVel(Object.Body, tempVect.x, tempVect.y, tempVect.z);



    // To start the object with a different rotation each time the program runs we create a new matrix called R and use

       // the function dRFromAxisAndAngle to create a random initial rotation before passing this matrix to dBodySetRotation.

    dMatrix3 R;

    dRFromAxisAndAngle(R, dRandReal() * 2.0 - 1.0,

                          dRandReal() * 2.0 - 1.0,

                          dRandReal() * 2.0 - 1.0,

                          dRandReal() * 10.0 - 5.0);

    dBodySetRotation(Object.Body, R);



    // At this point we could add our own user data using dBodySetData but in this example it isn't used.

    size_t i = 0;

    dBodySetData(Object.Body, (void*)i);

    

    // Now we need to create a box mass to go with our geom. First we create a new dMass structure (the internals

       // of which aren't important at the moment) then create an array of 3 float (dReal) values and set them

       // to the side lengths of our box along the x, y and z axes. We then pass the both of these to dMassSetBox with a

       // pre-defined DENSITY value of 0.5 in this case.

    dMass m;

    dReal sides[3];

    sides[0] = 2.0;

    sides[1] = 2.0;

    sides[2] = 2.0;

    dMassSetBox(&m, DENSITY, sides[0], sides[1], sides[2]);



    // We can then apply this mass to our objects body.

    dBodySetMass(Object.Body, &m);



    // Here we create the actual geom object using dCreateBox. Note that this also adds the geom to our 

       // collision space and sets the size of the geom to that of our box mass.

    Object.Geom[0] = dCreateBox(Space, sides[0], sides[1], sides[2]);

    // And lastly we want to associate the body with the geom using dGeomSetBody. Setting a body on a geom automatically

       // combines the position vector and rotation matrix of the body and geom so that setting the position or orientation

       // of one will set the value for both objects. The ODE docs have a lot more to say about the geom functions.

    dGeomSetBody(Object.Geom[0], Object.Body);

}

That's it for the initialization stage for ODE, but I might just jump ahead and discuss the clean up stage as it is very short. Then we will delve into the actual simulation loop and look how ODE handles collisions.

 

Cleaning up


void CloseODE()

{

    // Destroy all joints in our joint group

    dJointGroupDestroy(contactgroup);



    // Destroy the collision space. When a space is destroyed, and its cleanup mode is 1 (the default) 

       // then all the geoms in that space are automatically destroyed as well.

    dSpaceDestroy(Space);

    // Destroy the world and everything in it. This includes all bodies and all joints that are not part of a joint group.

    dWorldDestroy(World);

}

The simulation loop

 

To update the simulation we call SimLoop each frame. It short, it calculates which geoms are colliding before advancing the simulation one step and then rendering the geoms.


void SimLoop()

{

    // dSpaceCollide determines which pairs of geoms in the space we pass to it may potentially intersect. We must also

       // pass the address of a callback function that we will provide. The callback function is responsible for

       // determining which of the potential intersections are actual collisions before adding the collision joints to our 

       // joint group called contactgroup, this gives us the chance to set the behaviour of these joints before adding them

       // to the group. The second parameter is a pointer to any data that we may want to pass to our callback routine.

       // We will cover the details of the nearCallback routine in the next section.

    dSpaceCollide(Space, 0, &nearCallback);



    // Now we advance the simulation by calling dWorldQuickStep. This is a faster version of dWorldStep but it is also

       // slightly less accurate. As well as the World object ID we also pass a step size value. In each step the simulation

       // is updated by a certain number of smaller steps or iterations. The default number of iterations is 20 but you can

       // change this by calling dWorldSetQuickStepNumIterations.

    dWorldQuickStep(World, 0.05);



    // Remove all temporary collision joints now that the world has been stepped

    dJointGroupEmpty(contactgroup);



    // And we finish by calling DrawGeom which renders the objects according to their type or class

    DrawGeom(Object.Geom[0], 0, 0, 0);

}

Callback Function

 

All pairs of geoms that are potentially intersecting will be passed by dSpaceCollide to this function. Here we can determine which objects are actually colliding by making a call to dCollide and change the behavior of the joints before adding them to the joint group. The first parameter would be the user data pointer passed to dSpaceCollide if we had provided it. The second and third parameters are the two potentially intersecting geoms.


static void nearCallback (void *data, dGeomID o1, dGeomID o2)

{

    // Temporary index for each contact

    int i;



   // Get the dynamics body for each geom

    dBodyID b1 = dGeomGetBody(o1);

    dBodyID b2 = dGeomGetBody(o2);



    // Create an array of dContact objects to hold the contact joints

    dContact contact[MAX_CONTACTS];



    // Now we set the joint properties of each contact. Going into the full details here would require a tutorial of its

       // own. I'll just say that the members of the dContact structure control the joint behaviour, such as friction,

       // velocity and bounciness. See section 7.3.7 of the ODE manual and have fun experimenting to learn more.  

    for (i = 0; i < MAX_CONTACTS; i++)

    {

        contact[i].surface.mode = dContactBounce | dContactSoftCFM;

        contact[i].surface.mu = dInfinity;

        contact[i].surface.mu2 = 0;

        contact[i].surface.bounce = 0.01;

        contact[i].surface.bounce_vel = 0.1;

        contact[i].surface.soft_cfm = 0.01;

    }



    // Here we do the actual collision test by calling dCollide. It returns the number of actual contact points or zero

       // if there were none. As well as the geom IDs, max number of contacts we also pass the address of a dContactGeom

       // as the fourth parameter. dContactGeom is a substructure of a dContact object so we simply pass the address of

       // the first dContactGeom from our array of dContact objects and then pass the offset to the next dContactGeom

       // as the fifth paramater, which is the size of a dContact structure. That made sense didn't it?  

    if (int numc = dCollide(o1, o2, MAX_CONTACTS, &contact[0].geom, sizeof(dContact)))

    {

        // To add each contact point found to our joint group we call dJointCreateContact which is just one of the many

           // different joint types available.  

        for (i = 0; i < numc; i++)

        {

            // dJointCreateContact needs to know which world and joint group to work with as well as the dContact

               // object itself. It returns a new dJointID which we then use with dJointAttach to finally create the

               // temporary contact joint between the two geom bodies.

            dJointID c = dJointCreateContact(World, contactgroup, contact + i);

            dJointAttach(c, b1, b2);

        }

    }

}

An important thing to remember about the callback function is that you are not allowed to modify a collision space while that collision space is being processed with dSpaceCollide. For example, you can not add or remove geoms from a space, and you can not reposition the geoms within a space. Doing so will trigger a runtime error.

 

Drawing The Geoms

 

At the end of the SimLoop function, after advancing the simulation one step and then deleting the contact joints, we called a function called DrawGeom. This function serves as a generic rendering routine for all types of geoms, using a different rendering function for each class of geom. In keeping with the examples that came with the ODE library, the DrawGeom function takes four parameters. The first is the geom's ID, then its position vector, rotation matrix and lastly a flag used to toggle the rendering of the geoms axis aligned bounding box.


void DrawGeom (dGeomID g, const dReal *pos, const dReal *R, int show_aabb)

{

    // If the geom ID is missing then return immediately.

    if (!g)

        return;



    // If there was no position vector supplied then get the existing position.

    if (!pos)

        pos = dGeomGetPosition (g);



    // If there was no rotation matrix given then get the existing rotation.

    if (!R)

        R = dGeomGetRotation (g);



    // Get the geom's class type.

    int type = dGeomGetClass (g);



    if (type == dBoxClass)

    {

        // Create a temporary array of floats to hold the box dimensions.

        dReal sides[3];

        dGeomBoxGetLengths(g, sides);



        // Now to actually render the box we make a call to DrawBox, passing the geoms dimensions, position vector and

           // rotation matrix. And this function is the subject of our next discussion. 

        DrawBox(sides, pos, R);

    }

}

Anyone who has worked with OpenGL before will see that there is nothing out of the ordinary here. After we push the modelview matrix onto the stack we want to transform it by the position and rotation of the geom. We do this by calling a member function of my matrix class called ODEtoOGL. It takes the position vector and rotation matrix used by the ODE library and converts it to an OpenGL compatible 4x4 matrix (I will cover the inner workings of this conversion routine at the end of this tutorial.) Once we have a suitable transformation matrix we simply call glMultMatrixf to transform the current modelview matrix and then continue rendering the triangles of a box as normal.


void DrawBox(const float sides[3], const float pos[3], const float R[12])

{

    float mat_ambient[] = { 0.8, 0.8, 0.8, 1.0 };

    float mat_diffuse[] = { 0.8, 0.8, 0.8, 1.0 };

    glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);

    glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);



    glBindTexture(GL_TEXTURE_2D, texture[0].TexID);



    glPushMatrix();

    GeomMatrix.ODEtoOGL(pos, R);

    glMultMatrixf(GeomMatrix.Element);

    glBegin(GL_TRIANGLES);

        // Front Face

            glNormal3fv(&polygon[0].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[0].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[0].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[0].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[1].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[1].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[1].Vertex[2].x);

        // Back Face

            glNormal3fv(&polygon[2].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[2].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[2].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[2].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[3].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[3].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[3].Vertex[2].x);

        // Top Face

            glNormal3fv(&polygon[4].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[4].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[4].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[4].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[5].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[5].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[5].Vertex[2].x);

        // Bottom Face

            glNormal3fv(&polygon[6].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[6].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[6].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[6].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[7].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[7].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[7].Vertex[2].x);

        // Right face

            glNormal3fv(&polygon[8].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[8].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[8].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[8].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[9].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[9].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[9].Vertex[2].x);

        // Left Face

            glNormal3fv(&polygon[10].Vertex[0].nx);

        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[10].Vertex[0].x);

        glTexCoord2f(1.0f, 0.0f); glVertex3fv(&polygon[10].Vertex[1].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[10].Vertex[2].x);



        glTexCoord2f(0.0f, 0.0f); glVertex3fv(&polygon[11].Vertex[0].x);

        glTexCoord2f(1.0f, 1.0f); glVertex3fv(&polygon[11].Vertex[1].x);

        glTexCoord2f(0.0f, 1.0f); glVertex3fv(&polygon[11].Vertex[2].x);

    glEnd();

    glPopMatrix();

}

Note that in this example I have ignored the geom dimensions as the box is always a cube with side lengths of 2.

 

If you don't want to use my matrix class then the following function will do the job. Just enter an array of 16 floats as the first argument when calling the function and use the same array when calling glMultMatrixf.


void ODEtoOGL(float* M, const float* p, const float* R)

{

    M[0]  = R[0]; M[1]  = R[4]; M[2]  = R[8];  M[3]  = 0;

    M[4]  = R[1]; M[5]  = R[5]; M[6]  = R[9];  M[7]  = 0;

    M[8]  = R[2]; M[9]  = R[6]; M[10] = R[10]; M[11] = 0;

    M[12] = p[0]; M[13] = p[1]; M[14] = p[2];  M[15] = 1;

}

The Matrix Details

 

The ODE library uses a 3x3 matrix in mathematical notation, i.e. row first, column second. However, internally ODE stores its matrices as a 4x4 ordered matrix, padded with 0's. The good news is that the dBodyGetRotation function we used to retrieve the rotation matrix in DrawGeom above returns a 4x3 rotation matrix. So the ODEtoOGL function transposes the elements over to the OpenGL order (column first, row second) and plugs in the position vector into the 12th, 13th and 14th elements.




        | 0  1  2  3  |            | 0  4  8  12 |

        |             |            |             |

    R = | 4  5  6  7  |            | 1  5  9  13 |

        |             |        M = |             |

        | 8  9  10 11 |            | 2  6  10 14 |

                                   |             |

    p = | 0  1  2 |                | 3  7  11 15 |



And you're done!

 

You deserve a cuppa and a pat on the back if you were able to get most of that. We've come a long way from the days of manually rendering a polygon pixel by pixel along a scanline, now we can also add realistic rigid body physics to our programs.

 

You may have noticed that the gravity and density values being used in this example are a little strange. Higher values caused the box to penetrate too far into the ground surface or to fall right through the ground. I believe this is due to the slow message loop or other timing issues as it isn't as noticeable when using a separate rendering thread.

 

I hope you enjoyed the tutorial (in some intellectually masochistic way) and got a lot from it.

 

Until next time.

Alan Baylis

 

Do not read the following unless you are a search bot - I'm serious:

inertia matrix, integration, Runge-Kutta, linear and angular dampening, rotational torque, restitution, ordinary differential equations, rigid body dynamics