And here's our first update.

This version contains a couple more Vector3 macros, and a small extension to the Integrate method.
After integrating the Position, we also integrate the Velocity.
This new step requires that we have computed the Linear Force of each Body (see ComputeForces), which in turn requires that we know the MASS.
This is our earliest introduction to Mass, later we'll learn how to calculate it properly, but for now we'll just assume that the Mass is always 1.0f
I've also introduced 'inverse mass', which is simply 1/Mass.
If a physics formula requires that we divide something by Mass, we can instead multiply it by 1/Mass, thus avoiding the slower division operation, and speeding up the math... integrating velocity is one such formula where this small optimization can be made, and makes a good example.

So now, assuming that Gravity is turned on, not only will Position change, but Velocity will too... gravity will accelerate the body toward the ground.

I've also taken the step of rewriting my old Timer code as a new OA32 object.
The next update will include this file, I just want to check it for bugs first.

EDIT: I've tested the Timer code, I've embedded a timer in the Simulator, I've moved the CreateSphere method from CollisionBody into Simulator, I've added Start and Stop methods to the Simulator that control a thread which is dedicated to driving the simulator.
Now the user doesn't need to worry about the Timer code, and the physics stuff is completely decoupled from the GUI and Rendering code (if any).
We should see big performance gains on multi-cpu systems, and it will still work on single-core machines at the cost of an extra thread.
Best of all, its more user-friendly.
Posted on 2008-12-15 07:15:58 by Homer
Here's our first working proof of the code so far.
Included is the new Timer object (yes, I found bugs, and rewrote it)... if you look, you'll notice that I used the "inversion" trick to avoid divide operations in the timer (the Period field is actually 1/TimerFrequency).

As mentioned, the simulator class now contains its own driving thread, which means that the simulation will execute concurrently with any other code that you might have...  and if you have more than one logical CPU, you'll see a big performance improvement when the simulator is placed under load, and/or we're doing expensive stuff in other threads at the same time, like rendering 3D graphics. This version also embeds its own Timer, and is very easy to use since the user doesn't need to worry about time.

The demo provided will add a Sphere object to the simulator, beginning at Position <0,0,0> and Velocity <0,0,0> with Gravity enabled.
The demo will emit messages to DebugCenter (do you have this installed?) which clearly shows the object falling in the Y-MINUS axis (since Y-PLUS is UP)... its velocity increases over time as it accelerates under the force of gravity, and its positional changes become greater in response.

Just nice to actually be able to see 'something' :)
Posted on 2008-12-16 08:09:09 by Homer
The current code uses a default Mass of 1.0 ...
It's time to learn another simple physics formula, this one is for calculating Mass from Volume and Density:

Mass = Volume * Density

Pretty easy, eh? Well, it's easy if we know how to calculate the Volume of a given geometry.
For a Box, the volume would be: Volume = Length * Width * Height
But we're working (at least for now) only with Spheres.

Assuming that the Sphere is SOLID, we can calculate the Volume of a Sphere like this:
Volume = 4/3 * PI * R^3

And if the Sphere is HOLLOW, we just calculate the major and minor volumes of the inner and outer walls of the sphere, and subtract the smaller value from the larger one.
Given R1 = outer radius of sphere and R2 = inner radius (ie, thickness = R1-R2) :
Volume = (4/3 * PI * R1^3) - (4/3 * PI * R2^3)

So now we're able to find the volume of solid and hollow spheres, what about the Density?
We can't actually calculate it, because it depends totally on the material that the Sphere is made of.
Basically, early scientists decided that water had a density of 1.0, and compared other materials to it.
Things that are lighter than water (and thus float, regardless of shape) have a density less than one.
And things that are heavier than water (and sink unless they displace a large volume, like boats do) have a density of greater than one.
Each kind of material has a specific Density, we can create a table of densities for common materials.
Here's a few:

Balsa Wood = 0.17
Willow        = 0.42
Oak            = 0.75

Water       = 1.0

Aluminium = 2.7
Iron          = 7.8
Copper      = 8.6
Gold         = 19.3

Using this knowledge, we can write a general-purpose SetMass method for all bodies:

;This method creates a new Sphere collision body and adds it to the simulator
Method Simulator.AddCollisionSphere, uses esi, fOuterRadius:real4, fMaterialDensity:real4, IsHollow:BOOL, fInnerRadius:real4
LOCAL fVolume:real8
    SetObject esi
    New CollisionBody,Init, esi, SHAPEID_SPHERE
    m2m .CollisionBody.fRadius,  fRadius,edx
    m2m .CollisionBody.Density,  fMaterialDensity,edx
    m2m .CollisionBody.IsHollow, IsHollow, edx
    ;Add the new sphere to the simulator
    OCall .Bodies::Collection.Insert, eax
    ;Calculate VOLUME
    mov edx,4
    fildReg edx
    mov edx,3
    fildReg edx
    fdiv            ;4/3
    fmul            ;* pi
    fld  fOuterRadius
    fmul fOuterRadius
    fmul fOuterRadius
    fmul            ;* r^3
    .if IsHollow==TRUE
        mov edx,4
        fildReg edx
        mov edx,3
        fildReg edx
        fdiv            ;4/3
        fmul            ;* pi
        fld  fOuterRadius
        fmul fOuterRadius
        fmul fOuterRadius
        fmul            ;* r^3
        fsub            ;find the difference
    fstp fVolume
    OCall eax::CollisionBody.SetMass, fMaterialDensity, fVolume

We've had to add an extra parameter to the AddCollisionSphere for the 'inner radius', but its only valid if the sphere is hollow (IsHollow = TRUE).

Given that we 'know' the density of the material that a sphere is made of, we now have code that can accurately calculate the Mass, whether it is hollow or solid.

Later we'll look at 'buoyancy' - we'll be able to determine whether a hollow sphere will float or sink in water (and other liquids) based on its 'displacement' and its mass.

Posted on 2008-12-17 00:19:25 by Homer
Now I'm not quite sure which way to take this thread next.
There's basically two ways to go: we can introduce rotation and then collision, in which case rotation must immediately be part of the collision response discussion.
Alternatively, we can introduce basic linear collision and response, then add rotation into the fray.

Which path should we take?

EDIT : That was a question, you can put your own opinion here.
Posted on 2008-12-17 05:52:22 by Homer
rotation :)
Posted on 2008-12-17 09:58:02 by Ultrano
OK, we'll do it the hard way then :P

The first thing we need to know : Linear and Angular calculatations are completely separate.
We'll need to write Angular counterparts to all our Linear code.

We've described Linear Velocity as being 'the rate of change of Position' - literally, the rate at which the Position is changing, or 'how fast its moving'.
Now we introduce Angular Velocity, which I'll defined as 'the rate of change of Orientation' - literally, the rate at which the body's Orientation is changing, or 'how fast its spinning'.

In the ComputeForces method, we calculated Linear Force from Linear Velocity.
For rotation, we calculate Angular Force (also known as Torque) from the Angular Velocity.

The formula for calculating Torque is essentially the same as its linear counterpart.
linear force += velocity * linear damping
angular force += angular velocity * angular damping

So the changes to the ComputeForces method are really easy.
However, the changes to the Integrate method are not ALL so simple... but we'll deal with each problem as we come to it.

In the Integrate method, after the 'Linear' code we already have, the first thing we need to add is some code to calculate the ANGULAR MOMENTUM. The reason we need Angular Momentum is that at the very end of the Integrate method, we'll use it to extract a new value for the Angular Velocity.

New AngularMomentum = Old AngularMomentum + DeltaTime * Old Torque

OK, so far, so good.
The remainder of the Integrate method contains some fairly complex maths, so before we dive head-first into that, let's pause briefly and talk about a couple of important things.
I've mentioned previously that I'm using a 3x3 Matrix to describe Orientation.
One of the reasons for this is that we'll be using a 4x4 Matrix to render our object on the screen, and the Mat33 Orientation matrix fits neatly in the top left corner of the 4x4 transformation matrix (the remaining cells are used to describe position and scale).
Now I'd like you to imagine that you are in outer space, playing a game of pool : you are poking a ball with a stick in 'mid-air' , with no gravity, no other forces.
When we apply an external force to a body in a line that passes through the center of mass, the linear velocity of the object will change, but its angular velocity won't ... the way it is spinning is not affected.
If we poke the ball 'off-center', we'll cause its spin to change as well as its linear velocity.
And if we poke it with a 'glancing blow' on its very edge, its spin will change considerable, while linear velocity is unaffected.

In order to convert forces into changes in angular velocity, we employ a mathematical tool called INERTIA TENSOR.
This is another 3x3 Matrix, however it does NOT represent an orientation.
It describes the way that Mass is distributed across the volume of the body, which is a function of its SHAPE.
Just as for Density, we can create a table of inertia tensors for common shapes, however we can also calculate the inertia tensor for ARBITRARY shapes - I won't be getting into that now, and I won't talk much about the inertia tensor, only that its a Mat33, and that its content is "known" for common shapes like Sphere.

In my next post, we'll write the remaining code for the Integrate method, which will mention the inertia tensor.
Posted on 2008-12-18 00:09:22 by Homer

So, we've just calculated the relatively easy 'instantaneous angular momentum'.
The reason we can't simply track this value is that (for anything other than a sphere) it is constantly changing based on the shape of the object and its current axis of rotation... but calculating it wasn't hard, was it?

Now we're ready to find the new Orientation.
To do that, we'll need to apply the angular velocity for some Time ... deltaTime.
We'll be using another mathematical tool based on the Mat33 called a 'skew'symmetric matrix', which is built from our angular velocity vector.
That matrix will be multiplied by the Body's OLD Orientation, and scaled by DeltaTime, resulting in an 'offset rotation matrix' which describes the CHANGE in orientation over the given time period.
Finally, we add to it the original Orientation matrix, resulting in a NEW Orientation matrix...

OffsetOrientation = (DeltaTime * SkewSymmetric(Source.AngularVelocity) * OLD Orientation)
NEW Orientation = OLD Orientation + OffsetOrientation

We can create a SkewSymmetric Mat33 from a Vec3 using the STAR Matrix Operator.
I'll provide a Mat33_Star macro for this purpose.

So now we have our NEW ORIENTATION! Woohoo! But we're not done yet.
Due to the imprecise nature of floating point math, we will find that the above operations will 'deform' our Orientation Matrix... after a few iterations, it will no longer represent a rotation matrix, and that's bad.
So the very next thing we'll do is 'squash' our Orientation matrix back into shape, forcing it to be a nice rotation matrix. This is known as 'orthonormalizing', and I'll provide a Mat33_OrthoNormalize macro for you to use.

So we have a new orientation matrix, and its a nice rotation matrix, are we there yet? NO! Quiet in the back seat or I'll pull over and make you walk!

The Inertia Tensor we stored in our Body is given in 'Body Space', so that it stays legal no matter what orientation the body has... and we take the liberty of INVERTING this matrix to avoid Divisions...
Our next step is to transform the Inverse BodySpace Inertia Tensor into WorldSpace using the new orientation matrix we just calculated. To do this, we use a 'Similarity Transform' as follows:

InverseWorldInertiaTensor = New Orientation * InverseBodyInertiaTensor * Transpose(New Orientation)

Finally, we can extract the new Angular Velocity from the tensor we just calculated, and the angular momentum we calculated earlier:

New AngularVelocity = InverseWorldInertiaTensor * New AngularMomentum

Whew !! We now have all the stuff we need for both the CalculateForces and Integrate methods.
Furthermore, we have everything in place that we need to apply external forces to our body and see its rotation changing correctly, so when we introduce collisions and friction, our body will behave as it should.

I told you rotation wasn't as simple, but we've completely covered the stuff we needed to, and we can put it behind us and move on (just tell yourself it was a bad dream).
Actually, we've just covered the most complicated stuff in the physics simulator - from now on, everything is pretty easy!

Coming up next, we'll spend some time looking at collision detection, we'll implement that code in a demo, and then we'll move on to collision response without Friction, and finally, we'll implement a friction model that will allow our bodies to realistically interact with the surfaces they encounter, causing things like 'rolling' to 'magically happen'.
Posted on 2008-12-18 00:42:25 by Homer
In regard to NETWORKED Physics, it's worth noting at this point that we have everything we need for CLIENTSIDE physics...

The simulation is executed on both the Client and Server.
But the Client only integrates bodies, it does not perform collision testing, or collision resolving.... it simply performs 'Dead Reckoning' based on the last known physics state of each simulated object.

The Server implements a full-featured simulator which does 'everything', and periodically sends state information to the Clients, who make 'corrections' (as directed by the Server) to their simulation.

If it wasn't for the Server, the Client would be able to walk through walls etc.

Doing things this way makes a lot of sense - usually, the Server doesn't have a 3D rendering interface, it doesn't render anything at all, it's just a networking program crossed with a physics simulator.

The client is busy doing a lot of 3D rendering, by NOT performing collision testing, the client can concentrate on doing its main job of drawing eyecandy, and just integrate the bodies it simulates until the Server tells it 'hey, this thing just changed directions due to a collision'.

Posted on 2008-12-18 03:21:56 by Homer
Attached is a complete update, containing code for pretty much everything I've talked about so far.
Of particular note is the new Mat33 macro set for 3x3 matrix math operations.

The demo doesn't appear to do anything more than it was before, because I haven't asked it to - but all the new code is being called, and the application is not crashing and burning :P

Over 400 views in 10 days, that's just under two hits per hour, 24 hours a day - someone out there is watching..
Posted on 2008-12-18 07:44:12 by Homer

It's time to begin talking about collision detection - I'll treat collision response as a separate topic.

There's two basic kinds of collision testing, we can call them 'swept and instantaneous', or we can call them 'broadphase and narrowphase'. Let's begin by looking at the narrow phase (instantaneous) testing, which is quite simply testing for INTERSECTION of two entities at a specific moment in time.

For two spheres, the instantaneous test involves calculating the distance between the centers of the two spheres, and comparing that distance with the sum of the radii of both spheres...
If the distance is less than the sum of the radii, then the spheres are penetrating one another, ie, intersecting.
If the distance is equal to the summed radii, then the spheres are colliding (touching one another).
If the distance is greater than the summed radii, then the spheres are clear of one another.

What if we wanted to check for intersection of a sphere and a flat surface?
All flat surfaces lay apon a 3D PLANE.
Planes always have a 'surface normal', which tells us which direction in 3D space the plane is facing.
There is a well-known test to find the distance between a POINT and a PLANE... we calculate the Dot Product of the Normal and the Point - the result is the Signed distance, and negative values indicate which SIDE of the Plane the point lays. This Point/Plane test can be extended to deal with Spheres, we must take the radius into account - we can either modify the distance result by the radius, or we can shift the plane normal by the radius before we calculate the distance. Either way, we can determine whether the plane intersects the sphere (and if so, by how much) or whether the sphere is totally behind or in front of the plane (and find the distance between them).

Posted on 2008-12-18 21:35:04 by Homer
And what if we want to check for intersection (perform instantaneous collision detection) of a sphere and a triangle?
First, we find the Plane of the triangle, and perform the sphere/plane test, and if the sphere is colliding, penetrating, or BEHIND the plane (anything but CLEAR of it), we find the point of intersection on the plane, and then test if that point is within the bounds of the Triangle ...Here's two possible methods for doing that .. "Same Side" and "Barycentric" techniques for determining if a point is within a triangle:

Now we've talked about a few instantaneous intersection tests, we can move on to swept tests.
I've already provided the sourcecode for an implementation of the sphere/sphere sweep test.
You can find algorithms and pseudocode for sphere/sphere and sphere/plane sweep tests here:

Instantaneous tests are only good for testing if theres a collision at a specific instant in time.
Swept tests are different because they test a whole slice of time - we can define the start time and the end time, and check all the time inbetween with a single test by using 'quadratic equations'. Of course, this linear solution is not accurate over large slices of time, we're talking about a maximum of (PhysicsTimeStep) seconds of time.
So we're likely to be passing 0.100 as the time-range for these tests. Swept tests are particularly good at detecting collisions of small and fast-moving entities such as bullets.

For moving spheres we see that instantaneous testing isn't generally useful to us : collision detection should always begin with a Swept test.
If the sweep test finds a collision, it will return the time-delta which represents the first moment of contact.
We can then wind the simulation forward to that moment in time, handle the collision event, and then continue the simulation from the time of impact.
Posted on 2008-12-19 17:54:01 by Homer
So - for two spheres. the goal of broad-phase collision testing is to determine if a collision is possible, and my implementation also tells us about the moment of impact, and where it occurred.
Narrow-phase collision is not necessary for Spheres, but other shapes will need them.
I'd rather talk about collision response and complete the engine for Spheres before talking about other shapes.

My collision testing/resolving will begin with testing all pairs of bodies, and making a list of pairs of colliding bodies, sorted by the time at which the collision occurs.
Then I will process this list in a loop:
1. integrate the system to that time,
2. we take the pair with the earliest collision time from the list, we try to resolve this collision
3. If a collision was resolved, we swap the states of all our bodies.
   (there can be 'non collisions' where a collision has become invalid)
3. we subtract the 'earliest collision time' from all remaining collision pairs in the list
4. we LOOP on step 2 until the list is empty.

Posted on 2008-12-20 00:25:05 by Homer
Hey Homer,

I'm quite curious how do you check for intersection between polygons...  ;)
Posted on 2008-12-20 03:59:32 by roticv
I didn't want to talk about more complex shapes yet - but since you insist..

Polygons begin with the most simple 3D polygon - the tetrahedron.
This is a triangular pyramid with four faces, and four vertices.
We can say a cube has six faces, and eight vertices.

Any 3D polygon can be described as a set of faces and a set of planes relative to the origin of that polygon (ie, the geometric description of the polygon is made in 'Body Space').
The reason that we choose the BodySpace representation is that it doesn't change while our Body is rotating.
We perform the collision test from one Body's Space , so that only the geometric features of the other Body need to be transformed before testing. Collision points are returned to WorldSpace, where the Collision Normal is found.

One way to check for intersection of two arbitrary polygons: we must compare the closest vertices of each body against the closest faces of the other body.  An improved method uses the Separating Axis Theorem, with Voronoi regions used to identify the closest features.
There's a totally new algorithm too, the name completely slips my mind, which doesn't require voronoi.

Anyway, the more complex a geometry becomes, the more expensive it is to collide with.
So beyond some limit, we should surround our complex 3D geometry with an imposter - a convex hull made from one or more geometric primitives which closely represents our more complex geometry.
There are algorithms which can construct a convex approximation of an arbitrary geometry which accept a 'cloud of points' as input.
For complex geometries, we should always perform a preliminary test against an imposter, and if that collides, we can perform full intersection testing using the real geometries.

Posted on 2008-12-20 19:05:17 by Homer
I now have another choice of direction - composite rigid body physics, or collision resolution.
Please vote for something :)
Posted on 2008-12-21 02:59:35 by Homer
Hope I at least got you thinking?
Posted on 2008-12-21 06:43:32 by Homer
I vote for collision resolution.
Posted on 2008-12-21 17:41:16 by ti_mo_n
I vote for collision resolution too.

The idea of using convex hull to approximate the collision is a good idea. It seems extremely complicated and time consuming to check every vertex. Imagine 100 objects which are complex polygons with 1000 vertices each and checking the collision for everything.
Posted on 2008-12-22 09:55:35 by roticv
Depending on the implementation, I would use convex imposters only to determine if a collision of the high-resolution geometry is LIKELY - ie, as a preliminary test. If we see the convex imposter collide, then we can perform the most expensive test using the true geometry (for games, we're probably happy with the result of the imposter collision).

OK, now let's start talking about collision resolution.
I'm going to make the assumption that we have not only detected a collision will occur in the current timestep at some small future time delta, but we have used that time delta to integrate the simulation forwards to the exact moment of collision - the body or bodies must be in the colliding state where they are 'just touching'.
In this example, I'm going to discuss the most simple collision of a single spherical body and a fixed Plane.
So we have this sphere, and it's positioned so that it has just made contact with the plane.
In order to resolve a collision, we need two important pieces of information.
The first is the Collision Normal: for a collision of a body and a fixed plane, its simply the normal of the plane... for two bodies, we use the convention that it is the normal from body B back to body A.
The second piece of information we need is the 'Contact Point' - the place where the two entities actually made 'first contact'.
My solution of the collision is performed in 'Contact Space' - although this is not necessary for 'frictionless' collisions, it is VERY necessary when we introduce friction, so let's learn how it's done.
We're going to calculate the linear forces at the collision point which are in the direction of the collision normal, use that to calculate a Collision Impulse Per Unit Velocity, and then scale the impulse to achieve a Desired Post-Collision Velocity.

In the next post I'll cover all the steps we must follow in order to calculate a Collision Impulse in Contact Space, and transform it back to WorldSpace.

Posted on 2008-12-23 21:05:50 by Homer
For now, I'll assume that the Body is 'Clear' in its Old state, and 'Colliding' in its New state.

The first thing we need is a 3x3 Matrix for transforming (rotating) between worldspace and contact space.
Mat33_OrthoNormalBasis is a macro I'll be providing which requires only the Collision Normal to produce a matric we'll call ContactToWorld.

Just a quick word about Vectors in ContactSpace coordinates - the X axis represent the direction of the collision normal - the Y and Z axes represent the vectors orthogonal and normal to X, ie the 'sliding' vectors.. it will make more sense as we go.

The next thing we need is the position of the collision point with respect to the center of mass of each body.
RelativePosition = CollisionPoint - Body.NewState.CMPosition

Now I want you to imagine that the collision point is actually 'stuck' to the body - theres a particle on the surface of the body that happened to coincide with the collision point at the collision moment in time.
Now we're going to learn how to find the linear velocity of that particle on the surface of a body.
Given that the body is probably moving AND rotating, this is NOT the same as the linear velocity of the Center of Mass of that body !!! Furthermore, for my solution to work, I want that velocity in Contact Space.
We're going to need the RelativePosition and the ContactToWorld matrix.

First off, to find the velocity of the collision point, in world space:
Linear Velocity of Point in WorldSpace = CrossProduct(Angular velocity of body, RelativePosition) + Linear velocity of body

Now, we want to transform that from WorldSpace to ContactSpace:
    LinVel of Point in ContactSpace = Transpose(ContactToWorld) * Linear Velocity of Point in WorldSpace

To make that transformation a little cheaper, I've included a 'TransMult' macro (it performs the multiplication against what WOULD be the transpose of the input matrix).

We need to also find 'the change in linear velocity that was due to acceleration in this timestep", convert it to Contact coordinates, and zero the X axis of the resulting vector - we're not interested in the collision normal direction this time, but we're interested in acceleration along the 'friction' vectors - with sufficient friction, these will be eliminated - with NO friction, they must be present. We add the two vectors we produced, and the result is the Linear Velocity in ContactSpace which we wanted.
If this collision involves TWO bodies, we'll need to go through all of that for the second body as well.

We are now ready to calculate our 'Desired Delta Velocity' which we would like our body to have AFTER the collision, and using that, we can calculate an IMPULSE which would produce that change in velocity.

Please note that frictionless collision resolution can be a whole lot more simple than I am making it by using ContactSpace - but collisions with friction really require that we work in contact space, and they're a whole lot more interesting and realistic than non-friction collisions.
Imagine a ball falling onto a flat surface, and that the ball is spinning (has angular velocity).
With no friction, it will bounce straight back up every time... the fact that it is rotating will be ignored!
But with friction, it will go flying off in the direction that it was rotating, and its angular velocity will be reduced.. some of the energy it was using to rotate will be transformed into linear motion. That's much more realistic, isnt it?

The next post will take us through the process of calculating DesiredVelocityDelta (WorldSpace) and Impulse (ContactSpace).
Posted on 2008-12-23 23:58:01 by Homer