Guns, Cars, and Sports

I watched most of Sony's E3 presentation the other day. They were showing previews of their games lineup for the coming year. If you like guns, cars, or sports, Sony's got you covered. Especially if you like guns.

Presentations like this always give me a sinking feeling. There is a huge, huge gap between the kinds of games that could get made and the ones that do get made.

Rigid Body Dynamics: Friction

Friction is tricky. It's generally formulated as a force, integrated over time. In my program as it stands now, only gravity is integrated; the rest of the simulation is handled via impulses between frames.

Numerical integration methods make simplifying assumptions about how forces will behave over the frame duration. The simplest is to assume forces won't change. For a uniform gravity field this is not a bad assumption. For friction, though, it's not a good assumption. Friction forces should never cause objects to reverse direction, which is entirely possible if you assume the friction doesn't change during a frame.

I am going to try and work friction into my impulse model. Brian Mirtich's thesis covers this in quite a lot of detail. Basically I will assume that friction doesn't change during a frame, but clamp it so that it never reverses object directions or adds energy.

This week I did a tiny step toward that. When a collision happens, I compute an impulse that cancels all relative motion at the collision point. This is equivalent to a completely inelastic collision (because motion in the normal direction is canceled) combined with infinite static friction (because all motion in the tangential direction is canceled). I decided to try this out because it's sort of the maximum friction possible; sliding friction and bouncing will involve reducing it somehow.

To come up with the proper impulse, I wrote an expression for how the relative velocity changes as a function of the (unknown) impulse, then set it equal to the negative of the relative velocity to be canceled:



In the equation p is the unknown impulse. The r constants are the positions of the contact point with respect to the two bodies' centers of mass, but turned 90 degrees counterclockwise. I put the little T symbol on them to indicate that. Basically, these are the direction vectors for affecting the objects' rotational velocity. They scale up as the contact point gets further from the center of mass, since you get better leverage the further out you get. The M and I constants are the mass and rotational inertia tensor, respectively; in 2D they're both scalars.

This is a vector equation, so it works out to two scalar equations in two unknowns (the components of p). I used Cramer's Rule to solve the system of equations for the impulse.

With this friction, wheels roll, but they can never skid. I think I'll start to put together a vehicle this week. This will involve making some constraints to hold the wheels and the body together. I'll solve these the same way I solve nonpenetration constraints; by fixing up after integration.

Rigid Body Dynamics, part 3

(Part 1 and part 2 here.)

Flat-shaded discs aren't very good for seeing rotation, so I gave them a radial pattern:



Immediately, I saw something odd: occasionally a colliding disc would start spinning, sometimes at high speed. Since collisions are currently frictionless, this should not happen.

After poking at it for an hour I pinpointed the problem. I was collecting all contacts in a list before handling any of them. However, the collision handler moves objects apart when they are interpenetrating. The contact record stores a world-space contact position, which is incorrect if the objects in question have been moved by an earlier contact resolution. In frictionless disc collisions, the collision impulse should always be aimed through the center of the disc, but due to prior movement of the disc it would sometimes be off-center, which imparted a spin to the disc.

For the moment, then, I am back to finding and handling collisions one at a time.

I'm still daunted a bit by the quadratic programming espoused by some of the papers I've been reading. I've turned up a few newer papers, and it seems like people are going in the direction I have, which is to say: integrate without worrying about constraints, then iteratively fix up problems, rather than trying to come up with a giant system of equations (or inequalities) and solve them wholesale. I don't know if that's because it's better or just easier, though.

I'm currently trying to see if I can add friction to my current model, without going to the giant-system-of-equations approach. I need friction in order to implement things like wheels. Friction is usually modeled as two different kinds: resistance to motion while in motion (dynamic friction), and resistance to motion while stationary (static friction). A rolling wheel makes use of static friction since the part touching the ground is stationary with respect to it. A skidding wheel experiences dynamic friction instead.

Rigid Body Dynamics, part 2

My rigid-body dynamics program now handles discs in addition to boxes. Since there is still no friction, there is no way yet to affect a disc's spin.



Boxes and discs share much in common: they have mass, rotational inertia, momentum, position, orientation, color, and elasticity. They differ only in their geometry, which affects how they are drawn and how they collide.

I created a simple class hierarchy to represent boxes and discs. Both of them derive from a class called Body which contains all of the common state. The Box class adds width and height, while the Disc class has a radius instead. Here is how the classes stand now:

class Box;
class Disc;

class Body
{
public:

    virtual void draw() const = 0;

    virtual void integrate(float dt);

    virtual void compute_contacts(Body &, Contacts &) = 0;
    virtual void compute_contacts(Box &,  Contacts &) = 0;
    virtual void compute_contacts(Disc &, Contacts &) = 0;

    vec2 pt_velocity(vec2 p) const;

    vec2 pos; // position
    vec2 vel; // linear velocity
    float angle; // heading
    float angle_vel; // angular velocity
    float inv_mass; // inverse mass
    float inv_rot; // inverse of rotational inertia tensor
    float life; // remaining lifetime, in seconds
    float bounciness; // 0 to 1
    vec3 color;
};

class Box : public Body
{
public:
    virtual void draw() const;

    virtual void integrate(float dt);

    virtual void compute_contacts(Body &, Contacts &);
    virtual void compute_contacts(Box &,  Contacts &);
    virtual void compute_contacts(Disc &, Contacts &);

    vec2 radius; // half width, height in local space
    vec2 x_axis; // computed, cached axes
    vec2 y_axis;
};

class Disc : public Body
{
public:
    virtual void draw() const;

    virtual void compute_contacts(Body &, Contacts &);
    virtual void compute_contacts(Box &,  Contacts &);
    virtual void compute_contacts(Disc &, Contacts &);

    float radius;
};

The draw method is overridden by Box and Disc to draw the appropriate geometry. integrate is implemented at the Body level since it is the same for both kinds of bodies. The Box version call's Body's, then recalculates its axes based on the new box orientation. (The axes are used a bunch during contact determination, so it's good to have them computed only once per frame.)

The first step of collision handling is to collect a list of contacts between pairs of objects. (Last week, I handled each collision as I found it, but I'm trying to move toward the creation of forces to handle resting contact, so I need the whole set of contacts at once.)

To compute contacts, the program has to consider pairs of objects and decide:

• If they are touching
• Where they are touching
• What the surface normal is at the contact
• How deeply they are interpenetrating

As in the tutorials, I have a Contact structure that contains this information:

struct Contact
{
    Body * a;
    Body * b;
    vec2 pos; // world-space contact position
    vec2 normal; // world-space contact normal (points away from body a)
    float depth; // interpenetration distance (positive when interpenetrating)
};

The central loop that compares each pair of bodies does not know whether a particular body is a box or a disc. I use a two-stage virtual method dispatch to resolve the types of the two bodies. Here's the top-level loop (at the moment I am not doing any culling by bounding-box or otherwise, so all pairs are considered):

static void find_contacts(const Bodies & bodies, Contacts & contacts)
{
    contacts.clear();

    Bodies::const_iterator body_end = bodies.end();
    Bodies::const_iterator body1 = bodies.begin();
    for (; body1 != body_end; ++body1)
    {
        Bodies::const_iterator body2 = body1;
        for (++body2; body2 != body_end; ++body2)
        {
            if ((*body1)->inv_mass == 0 && (*body1)->inv_rot == 0 &&
                (*body2)->inv_mass == 0 && (*body2)->inv_rot == 0)
                continue;

            (*body1)->compute_contacts(**body2, contacts);
        }
    }
}

This method is overridden for boxes and discs, so the type of body1 is established. These methods then call the appropriate methods on body2, passing the first object as an argument (but with known type now). Here's the code:

void Box::compute_contacts(Body & body, Contacts & contacts)
{
    body.compute_contacts(*this, contacts);
}

void Box::compute_contacts(Box & box, Contacts & contacts)
{
    box_box_contacts(*this, box, contacts);
}

void Box::compute_contacts(Disc & disc, Contacts & contacts)
{
    disc_box_contacts(disc, *this, contacts);
}

void Disc::compute_contacts(Body & body, Contacts & contacts)
{
    body.compute_contacts(*this, contacts);
}

void Disc::compute_contacts(Box & box, Contacts & contacts)
{
    disc_box_contacts(*this, box, contacts);
}

void Disc::compute_contacts(Disc & disc, Contacts & contacts)
{
    disc_disc_contacts(*this, disc, contacts);
}

Once both body types are known, the correct function for computing contacts can be called. It's one of these three:

void box_box_contacts(Box & box1, Box & box2, Contacts &);
void disc_box_contacts(Disc & disc, Box & box, Contacts &);
void disc_disc_contacts(Disc & disc1, Disc & disc2, Contacts &);

As you can see, this is sort of a hack to get the functionality of multimethods, which most object-oriented languages don't support. With multimethods, we could just write four functions that took all the permutations of box and disc

The mechanics of detecting contacts between objects is dealt with to some extent in the tutorials, although they don't cover round objects. Round objects are pretty straightforward to add, although I had to work at it a bit. Some time I will write an entry about sphere vs. convex polyhedron collision; it's an interesting problem.

My contact handling hasn't changed from last week. I start by moving the two objects apart so they don't interpenetrate; then apply impulses to them (if they are moving into each other) so they aren't moving toward each other.

This coming week I hope to look into friction and/or resting contact.

Rigid Body Dynamics

I am continuing to learn how to simulate rigid bodies. This is something I've worked at in desultory fashion for many years. Unfortunately, since I'm in the middle of learning it, my thoughts are not especially well-ordered.

As I said last week, I'm following the Witkin/Baraff tutorials. I also found an interesting-looking paper by François Faure, which looks like it might show a better way to compute contact forces in a game.

Here's a shot of how far I'd gotten a long time ago. I was working with balls, simplified by ignoring rotation. A 2D box contains marbles of different sizes. Rotating the box makes them roll and bounce from one end to the other:



Exact Collision Times vs. Post-Collision Fixup

The Baraff tutorials assume that integration proceeds from collision to collision. At each exact collision time, the colliding bodies are adjusted so they won't pass through each other, and then the integration is resumed with the new velocities and accelerations.

The problem with this is that the amount of work per frame can vary quite a bit. A step of integration costs the same no matter how much time you're covering, so the per-frame cost is going to be proportional to the number of collisions during that frame.

In games, the ideal simulation costs the exact same amount every frame no matter what. We're even willing to tolerate a less-accurate simulation to achieve this, since it lets us budget CPU and memory better.

I have been deviating from the tutorial in that I'm attempting to essentially treat all collisions as if they occur at the start of a frame. Integration proceeds in whole-frame increments. Of course this means that objects can wind up inter-penetrating each other at the start of the next frame. To correct for this I do physical position fixup.

When an inter-penetration is detected, the objects are moved apart the minimum distance necessary. The movement is weighted by mass, so the lighter object does a bigger share of the movement. In the case of a collision with an immobile object, the mobile object does all of the moving.

Moving one pair of objects can cause another pair of objects to inter-penetrate. In simple situations, things converge toward a solution where nothing interpenetrates over the course of the simulation. For objects that aren't too interdependent this works surprisingly well. Problem cases involve resting contact (a stack of crates, for instance). I hope to deal with resting contact as in the tutorial, by introducing additional forces to keep objects apart during the integration. We'll see how well that works.

Progress

This week, I've gotten 2D boxes working. The environment is made out of (tan) boxes with infinite mass, which makes them immovable. I've got a little chute into which I can drop random boxes:



There's no friction yet. At the moment, gravity is the only force being integrated, so I just integrate assuming constant acceleration during the frame. I may be adding springs or things that require more complex integration, but I'll cross that bridge when I get to it.

In 3D, you have to consider point/face contacts and edge/edge contacts. In 2D, there are only point/edge contacts, so things are a bit simpler. In 3D, the rotational inertial tensor is a matrix, while in 2D it is reduced to a single scalar value.

I'm still struggling to work out a good code architecture. Next up, I plan to generalize my colliding bodies so I can support both boxes and circles; I want to make a simple 2D vehicle, and circles should work for wheels. After that, friction and resting contact.

Everett Kaser's games

Everett Kaser is the master of deductive logic games. From his lair near Corvallis, Oregon he unleashes a steady stream of clever, addictive puzzlers.

If you've ever tried out one of those pen-and-paper puzzles that has three people, three houses, three jobs, and some clues like “The person who lives in the red house is not a carpenter” you know roughly what the games entail. Kaser has designed his games so that the clues are presented in graphical form, and has blended in visual puzzle elements.

For instance, in Honeycomb Hotel you are trying to determine which of several options belongs in each hexagon in the board; but you are also simultaneously trying to determine a Hamiltonian path that enters the board, visits each hex, and exits. Here's an example of a very simple puzzle with two clues:



Each hex contains its possible options; you can eliminate them as you work through the clues. The central hex's contents have been given. The left clue says that the mouse and the letter H are adjacent and connected by a path. The right clue says that the grasshopper and the letter Y are adjacent but separated by a wall.

There cannot be a path connecting the entrance and exit hexes, because then the path wouldn't visit any of the other hexes on the board. This rules out the possibility of mouse and H being on the right side of the board. Once they are placed, Y is known, which places grasshopper and the rest. Here's a partial solution; can you see how to finish it?



Baker Street, another personal favorite, replaces the Hamiltonian path with a tree and is on a board of mixed squares and hexagons, for more varied connectivity:



My wife and I have spent more time on these games than any other. Free demos are available for everything. Give them a try and see if they are to your taste.

2D/3D Vector Classes

Most of my projects end up using vector mathematics, so I've written a couple of classes for handling the 2D and 3D cases. They have evolved a bit over time but have now reached a point where they do most everything I need.

The classes are templatized on the component type, so that I can support vectors of float and vectors of double with the same piece of code. I even use vectors of int sometimes.

There are two different notations for referring to the components of two- and three-dimensional vectors. Sometimes people use the letters x, y, and z since the coordinate system axes are traditionally named the same way. Sometimes people use numeric subscripts instead.

Numeric subscripts are superior to letter subscripts in at least one situation: when you need to choose components programmatically. For instance, let's say you've got a plane in 3D space, and you want to project that plane to 2D (to triangulate a polygon in the plane, say). An easy way to do this is to determine which of the three components of the plane's normal has the largest absolute value, and then use the axes corresponding to the other two components as the 2D projection plane. This type of projection is cheap since it just involves selecting two of the three components of a 3D point; a more general projection would involve dot products against basis vectors, which is computationally more expensive.

Because of the scenario above, I've gone with numeric subscripts in my vector classes. Some vector classes use a union of an array and members named x, y, and z to support both notations; I didn't bother.

Functions that operate on vectors can be standalone, or defined as methods of the vector class. There isn't too much difference either way. I implemented the dot, cross, len, etc. functions as methods of the class because the a.dot(b) notation puts the dot between the operands, like you would do in written math. It also keeps those names out of the global namespace; they are inside the vector class.

The one situation where you can't use methods is when you have an operator that doesn't take a vector as its first argument, for instance if you want to be able to multiply a scalar times a vector with the scalar on the left. This type of function has to be standalone.

The < operator gives vectors an ordering, which allows them to be used as keys in a lookup data structure such as std::map.

The empty constructor allows an uninitialized vector to be constructed, which is nice when you have an if-statement that initializes a vector in one of two different ways.

The +=, *= and similar operators which mutate a vector all return a reference to the vector itself, mimicking the behavior of these operators for standard types such as int. It's not typical or recommended, but you can write a = b += 5;. The vector class supports this behavior in order to be as much like a standard numeric type as possible.

This code is completely standalone, with one exception: square root finding. I use the standard C library functions sqrt and sqrtf for finding square roots of 64- and 32-bit floating point values, respectively. In order to choose the right square-root function for the scalar type at hand, I use template specialization.

These vector classes don't do anything to make use of the SIMD instruction sets that most processors have these days. Using special machine instructions from within C++ is platform-specific, and I have wanted to keep my code simple. Using SIMD badly can result in code that is slower than plain floating-point, too.

Here's the code:

#ifndef VEC_H
#define VEC_H

#include <cmath>

// Square root template, specialized for float or double.

template<typename T> T sqrt_template(T);
template<> inline float sqrt_template<float>(float value)
{
    return std::sqrtf(value);
}
template<> inline double sqrt_template<double>(double value)
{
    return std::sqrt(value);
}

// A basic 2D vector.
// vec2 = float version
// dvec2 = double version

template<typename T> class vec2_template
{
public:
    vec2_template() {}
    vec2_template(T x, T y)
    {
        m_v[0] = x;
        m_v[1] = y;
    }

    bool operator < (const vec2_template & rhs) const
    {
        if (m_v[0] < rhs.m_v[0]) return true;
        if (rhs.m_v[0] < m_v[0]) return false;
        return m_v[1] < rhs.m_v[1];
    }

    bool operator == (const vec2_template & rhs) const
    {
        return m_v[0] == rhs.m_v[0] && m_v[1] == rhs.m_v[1];
    }

    bool operator != (const vec2_template & rhs) const
    {
        return m_v[0] != rhs.m_v[0] || m_v[1] != rhs.m_v[1];
    }

    T & operator [] (int dim)
    {
        return m_v[dim];
    }

    const T & operator [] (int dim) const
    {
        return m_v[dim];
    }

    vec2_template operator + (const vec2_template & rhs) const
    {
        return vec2_template(
            m_v[0] + rhs.m_v[0],
            m_v[1] + rhs.m_v[1]);
    }

    vec2_template operator - (const vec2_template & rhs) const
    {
        return vec2_template(
            m_v[0] - rhs.m_v[0],
            m_v[1] - rhs.m_v[1]);
    }

    vec2_template operator / (T f) const
    {
        return vec2_template(m_v[0] / f, m_v[1] / f);
    }

    vec2_template & operator += (const vec2_template & rhs)
    {
        m_v[0] += rhs.m_v[0];
        m_v[1] += rhs.m_v[1];
        return *this;
    }

    vec2_template & operator -= (const vec2_template & rhs)
    {
        m_v[0] -= rhs.m_v[0];
        m_v[1] -= rhs.m_v[1];
        return *this;
    }

    vec2_template & operator *= (T f)
    {
        m_v[0] *= f;
        m_v[1] *= f;
        return *this;
    }

    vec2_template & operator /= (T f)
    {
        m_v[0] /= f;
        m_v[1] /= f;
        return *this;
    }

    vec2_template operator - () const
    {
        return vec2_template(-m_v[0], -m_v[1]);
    }

    vec2_template operator * (T f) const
    {
        return vec2_template(m_v[0]*f, m_v[1]*f);
    }

    T dot(const vec2_template & rhs) const
    {
        return m_v[0]*rhs.m_v[0] + m_v[1]*rhs.m_v[1];
    }

    // "Perp dot" product
    T perpdot(const vec2_template & rhs) const
    {
        return m_v[0]*rhs.m_v[1] - m_v[1]*rhs.m_v[0];
    }

    T sqlen() const
    {
        return dot(*this);
    }

    T len() const
    {
        return sqrt_template<T>(sqlen());
    }

    vec2_template normalized() const
    {
        return *this / len();
    }

    void normalize()
    {
        *this /= len();
    }

private:
    T m_v[2];
};

typedef vec2_template<float> vec2;
typedef vec2_template<double> dvec2;

// A basic 3D vector.
// vec3 = float version
// dvec3 = double version

template<typename T> class vec3_template
{
public:
    vec3_template() {}
    vec3_template(T x, T y, T z)
    {
        m_v[0] = x;
        m_v[1] = y;
        m_v[2] = z;
    }

    bool operator < (const vec3_template & rhs) const
    {
        if (m_v[0] < rhs.m_v[0]) return true;
        if (rhs.m_v[0] < m_v[0]) return false;
        if (m_v[1] < rhs.m_v[1]) return true;
        if (rhs.m_v[1] < m_v[1]) return false;
        return m_v[2] < rhs.m_v[2];
    }

    bool operator == (const vec3_template & rhs) const
    {
        return
            m_v[0] == rhs.m_v[0] &&
            m_v[1] == rhs.m_v[1] &&
            m_v[2] == rhs.m_v[2];
    }

    bool operator != (const vec3_template & rhs) const
    {
        return
            m_v[0] != rhs.m_v[0] ||
            m_v[1] != rhs.m_v[1] ||
            m_v[2] != rhs.m_v[2];
    }

    T & operator [] (int dim)
    {
        return m_v[dim];
    }

    const T & operator [] (int dim) const
    {
        return m_v[dim];
    }

    vec3_template operator + (const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[0] + rhs.m_v[0],
            m_v[1] + rhs.m_v[1],
            m_v[2] + rhs.m_v[2]);
    }

    vec3_template operator - (const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[0] - rhs.m_v[0],
            m_v[1] - rhs.m_v[1],
            m_v[2] - rhs.m_v[2]);
    }

    vec3_template operator / (T f) const
    {
        return vec3_template(m_v[0] / f, m_v[1] / f, m_v[2] / f);
    }

    vec3_template & operator += (const vec3_template & rhs)
    {
        m_v[0] += rhs.m_v[0];
        m_v[1] += rhs.m_v[1];
        m_v[2] += rhs.m_v[2];
        return *this;
    }

    vec3_template & operator -= (const vec3_template & rhs)
    {
        m_v[0] -= rhs.m_v[0];
        m_v[1] -= rhs.m_v[1];
        m_v[2] -= rhs.m_v[2];
        return *this;
    }

    vec3_template & operator *= (T f)
    {
        m_v[0] *= f;
        m_v[1] *= f;
        m_v[2] *= f;
        return *this;
    }

    vec3_template & operator /= (T f)
    {
        m_v[0] /= f;
        m_v[1] /= f;
        m_v[2] /= f;
        return *this;
    }

    vec3_template operator - () const
    {
        return vec3_template(-m_v[0], -m_v[1], -m_v[2]);
    }

    vec3_template operator * (T f) const
    {
        return vec3_template(m_v[0]*f, m_v[1]*f, m_v[2]*f);
    }

    T dot(const vec3_template & rhs) const
    {
        return
            m_v[0]*rhs.m_v[0] +
            m_v[1]*rhs.m_v[1] +
            m_v[2]*rhs.m_v[2];
    }

    vec3_template cross(const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[1]*rhs.m_v[2] - m_v[2]*rhs.m_v[1],
            m_v[2]*rhs.m_v[0] - m_v[0]*rhs.m_v[2],
            m_v[0]*rhs.m_v[1] - m_v[1]*rhs.m_v[0]);
    }

    T sqlen() const
    {
        return dot(*this);
    }

    T len() const
    {
        return sqrt_template<T>(sqlen());
    }

    vec3_template normalized() const
    {
        return *this / len();
    }

    void normalize()
    {
        *this /= len();
    }

private:
    T m_v[3];
};

typedef vec3_template<float> vec3;
typedef vec3_template<double> dvec3;

#endif

At the moment I am taking another crack at learning how to implement rigid body dynamics. This is a fairly well-trod path. I'm working from the highly-influential tutorials by Andrew Witkin and David Baraff. In my case, I'm developing it in 2D, not necessarily because it's any simpler, but because that's what I need it for. I hope to begin writing about it next week.