generalrelativity

In the previous article we implemented a particle constrained to a curve in two-dimensional space. We’ll expand on the concepts presented there and go through the process of creating modular constraints which can be easily interconnected. You’ll want to have some experience with 2D rigid body dynamics as we’ll be leaving behind the simplicity of particles for now. If “rigid body dynamics” sounds foreign to you, I recommend starting with the first 2 articles by Chris Hecker here.

As a fun little introduction to what we’ll be solving, look at the multi-body simulation below. The wheel has a torque applied and each body experiences a gravitational force. 1 ground constraint, 3 revolute joints and a prismatic joint constrain the system. Click the checkbox to toggle torque.

Demo: Roll over to simulate slider-crank.

The reason I wanted to write this article is because I found an imbalance between theory and implementation in the decent amount of freely available literature on the topic. The major rift exists, I think, because actual implementation involves radical optimizations that transform the tidy mathematical nature of the theory. In the concrete process and examples below, we’ll avoid this deviation with the hope that the material presented is a bit more accessible.

Global Form of the Constraint Problem

When we solved the particle-to-curve problem, it was sufficient to write the equations of motion in terms of the parabolic constraint. In the case that multiple constraints are acting upon a physical element, we need to write the equations of motion in terms of all constraints, because solving one may violate another. We could decide on a system configuration and solve each constraint force to derive the equations of motion, but that would be a lot of work; and not just once, but for each unique configuration.

Instead, we can construct a linear system of the form Ax = b, which can be solved for x generically. This allows us to piece together constraints anyway we like, which is a much more manageable way of setting up a simulation than sketching out the equations of motion for a giant, complex configuration. Now all we have to do is create this system, which basically just involves a lot of organization.

First, we need to compile the state of the system into a vector, q:

where n is the number of bodies. Each body is defined by an index, its position along the world x and y-axes and its orientation (θ). Next, we compile all constraints imposed on the system into the global constraint vector, C:

where s is the number of constraint equations. Everything we need to construct the aforementioned linear system is either directly derived from q and C, or ordered to match these vectors. Now we’ll see what we’re really working toward (see Witkin for derivation):

Equation 1:

This looks a little scary at first, but it’s actually not that bad. M^-1 is a diagonal matrix of dimension 3n, where each element on the diagonal is the inverse intertial term associated with the element at the corresponding index in q:

Where m_k is the mass of body k and I_k is the moment of inertia about the center of mass of body k. Looking at q, we see that mass is associated with translation and moment of inertia is associated with rotation, which is what we know from the Newton-Euler equations.

F_ext is the global vector of all external forces and torques exerted on the bodies:

which is of dimension 3n, also. Here, τ_k is the torque acting on body k. We could now write the unconstrained equations of motion as:

which should make sense of all the purposeful ordering we’ve seen so far.

Next, we’ll define the only radically new term introduced in Equation 1, J. This is known as the Jacobian of C and is an s-by-3n matrix, being composed of 1 row per constraint function (element of C) and 3 columns per body (x, y and θ). Row i is the gradient of the ith constraint function with respect to q:

Since a constraint generally involves at most 2 bodies, the gradient will be composed of mostly zeroes. The global Jacobian is:

This all might seem a bit abstract right now, but we’ll see that these terms are easy to construct in practice when we derive a few constraints. The time derivative of the Jacobian is built the same way, but with the velocity constraint vector:

The coefficients ks and kd of equation 1 are spring and drag terms, respectively, and used to correct error. Because each constraint equation should equal zero when satisfied, the amount of error is defined in the evaluation of the constraint functions.

Finally, we show how Equation 1 relates to the general form of a linear system mentioned above:

A is an s-by-s matrix, and x and b are vectors of dimension s. Solving for λ leads to the global constraint force vector:

Finally, the equations of motion for the constrained system can be written abstractly as:

Now we have a method for formulating constraints which can participate in the global system:

1. Write the positional constraint, C, as a function of q.
2. Take the derivative of C to solve for the velocity constraint.
3. Solve for the Jacobian of C.
4. Solve for the derivative of the Jacobian with respect to time.

Computer Implementation

public class Constraint
{
 
	//index in global constraint vector
	public var index:uint;
 
	//indicates the number of constraint functions
	//this is equal to the number of system DOF it removes
	internal function get numJacobianRows() : uint
	{
		//abstract
		return 0;
	}
 
	//populate the system matrices and vectors
	internal function solve( J:MatrixND, JDot:MatrixND, C:Vector. ) : void
	{
		//abstract
	}
 
}

Note that while we use C’ to solve for J’, we don’t actually need to explicitly contribute to the velocity constraint because C’ = Jq’ and can be solved at the global level without additional input from the individual constraints. Although this adds computational overhead, I like the generality achieved. Further, because we can alternatively solve for J’ by directly taking the derivative of J with respect to time, solving for C’ becomes seemingly unnecessary altogether. That said, I’ve found that the intuition-confirming statement a velocity constraint makes is worth the little bit of extra differentiation.

Now we’ll see how the global system is constructed and solved. Below is the relevant portion of the ConstraintSolver class, minus some declarations for conciseness’ sake. It allows each constraint to evaluate its very own special part of the structure via the above interface:

public class ConstraintSolver
{
 
	private var numJacobianRows:uint; //global # of constraint functions
	private var ode:DynamicsODE; 
 
	public function addConstraint( constraint:Constraint, initialize:Boolean = false ) : void
	{
		constraint.index = numJacobianRows;
		numJacobianRows += constraint.numJacobianRows;
		constraints.push( constraint );
		if( initialize ) initSystem();
	}
 
	public function initSystem() : void
	{
		W = new MatrixND( ode.length, ode.length, false );
		W.diagonal = ode.invMass;
		J = new MatrixND( numJacobianRows, W.m, false );
		JDot = new MatrixND( numJacobianRows, W.m, false );
		C = new Vector.( numJacobianRows, true );
		b = new Vector.( numJacobianRows, true );
		lambda = new Vector.( numJacobianRows, true );
	}
 
	public function solve() : void
	{
 
		for each( var constraint:Constraint in _constraints )
		{
			constraint.solve( J, JDot, C );
		}
 
		var JW:MatrixND = J.times( W );
		var JT:MatrixND = J.transpose;
		var JWJT:MatrixND = JW.times( JT );
 
		//J'q'
		var JDotqDot:Vector. = JDot.timesVector( ode.v );
		var JWFext:Vector. = JW.timesVector( ode.f );
 
		//C' = Jq'
		var CDot:Vector. = J.timesVector( ode.v );
 
		//b = -J'q' - JWQ - Cks - C'kd
		for( var i:int = 0; i < numJacobianRows; i++ )
		{
			b[ i ] = -JDotqDot[ i ] - JWFext[ i ] - C[ i ] * ks - CDot[ i ] * kd;
		}
 
		//solve the system
		JWJT.solveAnalytical( lambda, b );
 
		//constraint force = JTx
		var Fc:Vector. = JT.timesVector( lambda );
 
		//add the constraint force to the external forces
		for( i = 0; i < ode.length; i++ )
		{
			ode.f[ i ] += Fc[ i ];
		}
 
	}
 
}

With this structure in mind, we’ll now go through the process of deriving and implementing a couple constraints.

Ground Constraint

The first constraint we’ll derive with this method is the ground constraint, which keeps a body fixed in space. The body is not permitted to translate or rotate. First, we write the positional contraint:

where c_x is where the body is fixed along the global x-axis, c_y along the y, and c_θ the fixed orientation of the body. Note that 0 is the zero vector, [0,0,0]. Next, we solve the velocity constraint by taking the derivative of C_ground with respect to time:

which simply and intuitively states that linear and angular velocity are zero. This is what we expect considering the constraint prevents the body from experiencing any motion. Next, we derive J:

And, finally, the derivative of J w.r.t. time:

That was not so bad at all! The Jacobian can be a bit off-putting at times and make the problem seem more complicated than it is–just remember that we’re simply creating a matrix where each row is the gradient of a constraint function. In practice, you wouldn’t want to simulate the ground constraint, but it’s a good exercise because of how simple the math is. Other joints will have more complex positional constraints, but the underlying framework for formulating these quantities is the same.

The GroundConstraint class:

public class GroundConstraint extends Constraint
{
 
	public var body:Body2D;
 
	private var U:Point;
	private var theta:Number;
 
	public function GroundConstraint( body:Body2D )
	{
		super();
 
		this.body = body;
 
		U = body.position.clone();
		theta = body.theta;
 
	}
 
	override internal function get numJacobianRows() : uint
	{
		//removes all DOF
		return 3;
	}
 
	override internal function solve( J:MatrixND, JDot:MatrixND, C:Vector. ) : void
	{
 
		J.A[ index ][ body.index 	] = 1;
		J.A[ index ][ body.index + 1 	] = 0;
		J.A[ index ][ body.index + 2 	] = 0;
 
		J.A[ index + 1 ][ body.index 		] = 0;
		J.A[ index + 1 ][ body.index + 1 	] = 1;
		J.A[ index + 1 ][ body.index + 2 	] = 0;
 
		J.A[ index + 2 ][ body.index 		] = 0;
		J.A[ index + 2 ][ body.index + 1 	] = 0;
		J.A[ index + 2 ][ body.index + 2 	] = 1;
 
		JDot.A[ index ][ body.index 		] = 0;
		JDot.A[ index ][ body.index + 1 	] = 0;
		JDot.A[ index ][ body.index + 2 	] = 0;
 
		JDot.A[ index + 1 ][ body.index 	] = 0;
		JDot.A[ index + 1 ][ body.index + 1 	] = 0;
		JDot.A[ index + 1 ][ body.index + 2 	] = 0;
 
		JDot.A[ index + 2 ][ body.index 	] = 0;
		JDot.A[ index + 2 ][ body.index + 1 	] = 0;
		JDot.A[ index + 2 ][ body.index + 2 	] = 0;
 
		//positional error
		C[ index 	] = body.x - U.x;
		C[ index + 1 	] = body.y - U.y;
		C[ index + 2 	] = body.theta - theta;
 
	}
 
}

Here we finally see how the indexing of the constraint and the body play into the global vectors and matrices.

Revolute Joint

A revolute joint, in two dimensions, is like paper pinned up on a wall. The paper is prevented from falling to the ground, but can spin freely. To formulate the joint constraint, we need a mathematical representation of a point on a rigid body.

Point on Rigid Body

First, we’ll define the terms used to describe rigid body state. Just as with a particle, x is the body’s position vector:

A 2D body is also characterized by its orientation about the axis poking out of the screen. This is represented as the 2D rotation matrix, A:

Which is clearly parametrized by the angle, θ. Given these representations, we can solve for the world-space position of any point, P, that is defined in the body’s local coordinate frame. We designate this term r:

Equation 2:

Where u is the vector pointing from the body’s origin to P. Equation 2 can be decomposed from vector form into 2 scalar equations:

Figure 1: Rigid Body Transfromation

Figure 1 shows a rigid body’s local coordinate system and that same system transformed into world-space. Note that in the body’s local coordinate frame, all points and vectors are constants, whereas in the global frame, they are functions of time. Equation 2 can be written more explicitly to reflect this:

Although we’ll prefer the abbreviated form of equation 2, it’s important to remember which terms are functions of time, as we differentiate this function below.

Figure 2: Revolute Joint

The revolute joint pins 1 body to another at a shared point in world space. Because we can represent this point in terms of each body as shown in equation 2, the positional constraint is easy to write:

which states that the distance between the 2 bodies at these respective points is zero. Differentiating with respect to time leads to the velocity constraint:

indicating that the relative velocity of the bodies at these points is zero. Next, we evaluate the Jacobian of C_revolute:

And, finally, the time derivative of the Jacobian:

The RevoluteJoint class:

public class RevoluteJoint extends Constraint implements IRenderable
{
 
	public var b1:Body2D;
	public var b2:Body2D;
 
	private var u1:Point;
	private var u2:Point;
 
	public function RevoluteJoint( b1:Body2D, b2:Body2D, U:Point )
	{
 
		super();
 
		this.b1 = b1;
		this.b2 = b2;
 
                //convert U to b1's local coordinate frame
		var T:Matrix = b1.A.clone();
		T.invert();
		u1 = T.transformPoint( U.subtract( b1.position ) );
 
                //convert U to b2's local coordinate frame
		T = b2.A.clone();
		T.invert();
		u2 = T.transformPoint( U.subtract( b2.position ) );
 
	}
 
	override internal function get numJacobianRows() : uint
	{
		//removes 2 DOF
		return 2;
	}
 
	override internal function solve( J:MatrixND, JDot:MatrixND, C:Vector. ) : void
	{
 
		var uxc1:Number = u1.x * b1.A.a;	var uxc2:Number = u2.x * b2.A.a;
		var uyc1:Number = u1.y * b1.A.a;	var uyc2:Number = u2.y * b2.A.a;
		var uxs1:Number = u1.x * b1.A.b;	var uxs2:Number = u2.x * b2.A.b;
		var uys1:Number = u1.y * b1.A.b;	var uys2:Number = u2.y * b2.A.b;
 
		//J
		J.A[ index ][ b1.index 		] = 1;
		J.A[ index ][ b1.index + 1 	] = 0;
		J.A[ index ][ b1.index + 2 	] = -uyc1 - uxs1;
 
		J.A[ index + 1 ][ b1.index 		] = 0;
		J.A[ index + 1 ][ b1.index + 1 	] = 1;
		J.A[ index + 1 ][ b1.index + 2 	] = uxc1 - uys1;
 
		J.A[ index ][ b2.index 		] = -1;
		J.A[ index ][ b2.index + 1 	] = 0;
		J.A[ index ][ b2.index + 2 	] = uyc2 + uxs2;
 
		J.A[ index + 1 ][ b2.index 		] = 0;
		J.A[ index + 1 ][ b2.index + 1 	] = -1;
		J.A[ index + 1 ][ b2.index + 2 	] = -uxc2 + uys2;
 
		//J'
		JDot.A[ index ][ b1.index 		] = 0;
		JDot.A[ index ][ b1.index + 1 	] = 0;
		JDot.A[ index ][ b1.index + 2 	] = b1.av * ( -uxc1 + uys1 );
 
		JDot.A[ index + 1 ][ b1.index 	] = 0;
		JDot.A[ index + 1 ][ b1.index + 1 	] = 0;
		JDot.A[ index + 1 ][ b1.index + 2 	] = b1.av * ( -uyc1 - uxs1 );
 
		JDot.A[ index ][ b2.index 	] = 0;
		JDot.A[ index ][ b2.index + 1 	] = 0;
		JDot.A[ index ][ b2.index + 2 	] = b2.av * ( uxc2 - uys2 );
 
		JDot.A[ index + 1 ][ b2.index 	] = 0;
		JDot.A[ index + 1 ][ b2.index + 1 	] = 0;
		JDot.A[ index + 1 ][ b2.index + 2 	] = b2.av * ( uyc2 + uxs2 );
 
		//positional error
		C[ index 	] = b1.x + uxc1 - uys1 - b2.x - uxc2 + uys2;
		C[ index + 1 	] = b1.y + uxs1 + uyc1 - b2.y - uxs2 - uyc2;
 
	}
 
}

Considerations

We’ve followed the mathematical formulas closely, which has left us with a computationally inefficient application. Hopefully this trade has been for a solid understanding that lays the foundation for the reader to implement a system of constraints. As a next step, I recommend reading Erin Catto’s paper, Iterative Dynamics with Temporal Coherence, which highlights several optimizations in both time and memory.

Just as was the case with the particle-to-curve simulation, an RK4 integrator was necessary to not accumulate an insurmountable amount of error. The choice of the spring and drag coefficients becomes less trivial than expected and aid in the degradation of accuracy.

In later articles, we’ll switch from acceleration to a velocity-based constraint solver. This is simpler, more efficient and capable of handling a wider variety of constraint problems.

drew ActionScript 3.0, Mathematics, Physics

This is the first in a series on simulating constraints; in this case a particle is constrained to move along a curve in planar space. To follow along, you’ll need to know a little physics, calculus and some algebra, and hopefully you’ve implemented a simple particle system before!

I’ve chosen a parabolic curve. We’re actually going to solve the problem in 2 different ways, so we’ll be seeing 2 separate representations of the curve, the first being the implicit form:

Which looks like this for x {-200…200}:

The first method we’ll use for solving the constraint is described for a circle in this paper by Andrew Witkin. There are some convenient aspects of his example, and so it may not be entirely apparent how to constrain the particle to our parabola. I’ll follow Witkin’s format to show how his formula works in this context.

First, you write a positional constraint which is just our curve equation:

The reason the curve is the constraint is because it implies that when the function’s value at (x,y) is zero, then [x,y] is located on the parabola. Next, we define the velocity constraint, which is the derivative of the positional constraint with respect to time:

If the velocity constraint is satisfied, then we’ll always end up with a valid position. Finally, we find the derivative again with respect to time. This is the acceleration constraint function:

As was the case with velocity → position, if the acceleration constraint is satisfied, so will the velocity. This means that, provided we’re starting with a valid position and velocity, we only ever need to enforce the acceleration constraint on the system.

Let’s do a quick review of Newton’s 2nd law:

Which states that force is equal to mass times acceleration. We’ll be making use of the fact that F represents the sum of all forces acting on the particle. So:

Solving for acceleration gives:

The subscripts ext and c stand for external and constraint respectively. It’s worth pointing out for clarity’s sake that when a term is bold it represents a vector, whereas italicized variables are scalars.

Now we can replace the appearances of x and y’s 2nd derivatives, breaking F into its component parts to allow the constraint force to appear explicitly:

Equation 1

The 2 unknowns now appear, but only in the 1 equation, so we can’t solve yet. If we consider a particle attached to a frictionless wire, it should be apparent that the wire, which is the source of the constraint force, can do no work and so it can not add nor remove energy from the system. The kinetic energy of the system is defined as:

To gain anything useful out of this, we have to consider what direction a force can push without changing T. From examining the derivative of this function with respect to time:

it should be apparent that the constraint force has an important relationship with the velocity vector, specifically that their dot product should equal zero (or else T would change as a direct result of the constraint force). We can use this knowledge to write the constraint force components into a 2nd equation:

Equation 2

This makes sense not just in terms of energy conservation, but also based on the constraint function chain described above. We already know our particle has a legal velocity, and so we don’t want to add any force that would ruin the circle of trust–the constraint force should only pull the external forces onto the acceleration curve.

I first solved for the x constraint component in terms of y from equation 2:

Equation 3

Solving for y’ in the velocity constraint and substituting for x’ in equation 3 reduces the variable count:

Plugging the right-hand side in for F_cx in equation 1 and solving for F_cy yields:

We could explicitly solve for F_cx too, but that’s a non-trivial number of multiplications and a divide. Luckily, we don’t have to because this is for a computational simulation; after F_cy is known, the solution to F_cx is given in equation 3.

Plugged into F=ma, we have the equations of motion for a particle constrained to the parabola x² + 100y = 0.

Parametrized Form

Now we’ll solve the problem a different way, using a reformulation of Newton’s equations called Lagrangian mechanics. Lagrangian mechanics involves finding the number of degrees of freedom (DOF) in a system, then expressing that system in terms of generalized coordinates, one for each DOF. In our example we have a single DOF, even though our particle exists in a 2-dimensional space. First, we phrase our parabola parametrically:

Equation 4

Here, q is our lone generalized coordinate. Noting that x is a vector, the transformation from the generalized coordinate to our particle’s position is given in the parametric equation (4):

The Lagrangian equations of motion for a particle experiencing only conservative forces, as is the case with a particle constrained to a wire, is:

Equation 5

David Eberly gives a full derivation of equation 5 from Newton’s 2nd law in his excellent, if often daunting, Game Physics. F, again, represents all externally applied conservative forces and T is the kinetic energy of the system, which is a function of q and q’. Before we can define T for this system, we have to derive dx/dq and x’ because kinetic energy is dependent upon velocity.

and

Now we can define the kinetic energy of the system as a function of our generalized coordinate because x’ can be expressed entirely in terms of q and its first derivative:

Having determined the kinetic energy function, we can define the rest of the terms in equation 5. We’ll see that its first term includes q’s second derivative, which is the last piece needed to advance the simulation.

Lastly, F, in this case, is just the gravitational force, which we define as the vector mg[0,1]:

Next we plug all of this back into equation 5, which then defines the Lagrangian equations of motion for a particle constrained to the parametrized parabola (q,-q²/100):

Equation 6

Solving for acceleration yields:

Now that we’ve solved for acceleration, we can advance our 2-dimensional simulation by numerically integrating our 1-dimensional differential equations. Below is a demo, showing both methods running side-by-side. If you watch long enough, you’ll see the implicit form (left) eventually leaves the parabola as error accumulates. In contrast, the parametric form will always be on the parabola, but energy will dissipate. These errors are associated with numerical integration inaccuracies and floating point imprecision.

Roll over to run simulation. Left is the implicit form, right is the parametric

Considerations

Each method has advantages and disadvantages. Both offer a relatively easy-to-follow format. Explicitly solving for a constraint force as we did with the implicit curve is more intuitive to me, but the parametric form is more computationally efficient, and the algebra isn’t as messy. A downside to the parametric form is that it isn’t a generic solution–adding new constraints to the system means solving a totally new set of equations of motion.

Below are the separate implementations that each run in a generic dynamics harness. You can see that there is more involved computationally in the implicit form:

Implicit implementation

var particle:Particle = Particle( p[ 0 ] );
particle.exertForce( G ); //apply gravity
 
var m:Number = particle.mass;
 
var X:Vector3D = particle.position;
var V:Vector3D = particle.velocity;
 
var Fext:Vector3D = particle.force; //generically grab all applied forces
 
var Fc:Vector3D = new Vector3D();
 
//solve for the constraint force
Fc.y = ( -50 * m * V.x * V.x - 100 * X.x * Fext.x - 2500 * Fext.y ) / ( X.x * X.x + 2500 );
Fc.x = ( Fc.y * X.x ) / 50;
 
particle.exertForce( Fc ); //exert the constraint force on the particle

Parametric implementation

var q:Number = x[ 0 ];
var qPrime:Number = v[ 0 ];
 
a[ 0 ] = ( -q * qPrime * qPrime - 50 * G * q ) / ( 2500 + q * q );

As mentioned above, both of these are prone to computational error. In the demo, I’m using an RK4 integrator. Using a Euler integration scheme results in visible error within the first second of simulation. Although RK4 does a good job of reducing error, it isn’t perfect. Ideally, an error term is introduced to keep drift in check.

Move onto part 2, where we cover joints and global constraints.

drew ActionScript 3.0, Mathematics, Physics

I was tidying up my laptop today and found the slides from a presentation I gave at Schematic a couple years ago on Bezier curves. I can’t find the demos unfortunately, but the slides are pretty funny on their own.

If I find the demos, I’ll post them here.

drew ActionScript 3.0, General, Mathematics

I’m working on a game at work where a ball interacts physically with a series of capsules. The system was originally written to support circle vs line collision tests. It was then decided that it would be capsules instead of lines.

This post will show how finding the closest point on a line segment to a point can be easily extended to determine if and where a circle and line segment are intersecting; and how that test can be extended again to identify the point of intersection between a circle and capsule.

First, we need a test for the closest point, P, on a line segment, L, to point X.

We can define L as its endpoints A and B:

L = AB

We can derive the vector, v, as that which points from A to B:

v = AB = B - A

Next, we define a vector, w, pointing from A to our point, X.

w = AX = X - A

Now, we’ve expressed the relationship between X and L in terms of vectors. By defining w and v each pointing from A, we’re working in a 2D Cartesian coordinate frame where A is the origin. Because of this, we can take advantage of some geometrical properties of vectors that will allow us to find the closest point on L to X. Specifically, we’re going to project w onto v:

proj(w,v) = dot(w,v) / dot(v,v)

where dot is the dot product:

dot(w,v) = w.x * v.x + w.y * v.y

When envisioning projection, I find it helpful to think of vectors as having a spatial domain. If you compute the dot product of 2 vectors, it represents how much 1 vector exists in the domain of the other. The dot product of a vector and itself is a value which is a maximum of that vector’s domain. Anything greater than this value exceeds the domain. The dot product of a vector and one perpendicular to it is zero, and defines the minimum of that domain (exclusive). So the domain of a vector, u, in practice, ends up being spanned by any scalar value, q, which satisfies the constraint:

0 < q <= dot(u,u)

Now let’s apply this geometric intuition to our problem of finding the closest point on a line segment. What happens if the vector w points in an opposite direction than v? dot(w,v) will be negative. Therefore w cannot exist in the domain of v. If this is the case, it should be evident that A will be the closest point on L to X. Next, let’s consider the easy-to-imagine scenario of w exceeding the maximum of v’s domain, by defining w as pointing in the same direction as v, but having twice its length. Clearly, B in this case will be the closest point to X.

Example 1. Closest point on segment to user’s mouse:

0 <= t <= 1

Where t is the projection. By clamping t to this range, and using it to scale the vector v, we find the point on v closest to X. Because we’re working in the space of A, we need to transform back to world coordinates, by adding A:

P = A + v * clamp( proj(w,v), 0, 1 )

Here’s the implementation in my geometry library (Vector2D is not a native data type), with the variables changed to match the description of the problem:

public static function closestPointLineSegment( X:Vector2D, A:Vector2D, B:Vector2D ) : Vector2D
{
    var v:Vector2D = B.minus( A );
    var w:Vector2D = X.minus( A );
    var wDotv:Number = w.dot( v );
    var t:Number = w.dot( v ) / v.dot( v );
    t = MathUtil.clamp( 0, 1, t );
    return A.plus( v.times( t ) );
}

Circle/Line Collision Test

Now that we’ve found the closest point on a line segment, L, to a point, X, determining intersection of a circle and line segment is as simple as checking whether the distance between X and P is less than the radius of the circle. Of course, X defines the center of the circle.

Example 2. Intersecting a circle and line segment:

Note that in the code below, I’m checking for distance by comparing the squared distance against the squared length of the circle’s radius. This is to avoid an expensive square root call.

public static function testCircleSegment( X:Vector2D, r:Number, A:Vector2D, B:Vector2D ) : Vector2D
{
    var P:Vector2D = closestPointLineSegment( X, A, B );
    P.minusEquals( X );
    if( P.dot( P ) <= r * r ) return P;
    return null;
}

Circle/Capsule Collision Test

Lastly, intersecting a circle and a capsule is as simple as inflating the circle’s radius by the radius of the capsule. You can see in the below example that we’re really just doing a circle/line test on an inflated circle.

Example 3. Intersecting a circle and capsule:

One extremely nice aspect of working with circles is that, not allowing the shapes to penetrate, circle’s can be touching another convex shape at a maximum of 1 point. Due in part to this convenience, the collision normal will always point from the closest point on the convex shape through the circle’s center.

public static function testCircleCapsule( X:Vector2D, r1:Number, A:Vector2D, B:Vector2D, r2:Number ) : Vector2D
{
    var P:Vector2D = closestPointLineSegment( X, A, B );
    P.minusEquals( X );
    if( P.dot( P ) <= ( r1 + r2 ) * ( r1 + r2 ) ) return P;
    return null;
}

drew ActionScript 3.0, Geometry, Mathematics

Here’s a demo of the auto-gen physics engine (SLRBS- Super Lightweight Rigid Body Simulator) I’ve been working on which allows you to select images to use as rigid bodies. First, the image is analyzed, ignoring transparent pixels, to determine its convex hull. Then mass/center of mass/inertia tensor are solved for. Finally, the rigid body is created and dropped into the world.

I’m sure that in playing with it, you’ll notice all sorts of wonky behavior… but I think it’s fun anyway. Let me know if there’s a class or type of shape/image that tends to be particularly misbehaved!

demo
image pack

drew ActionScript 3.0, Flash Player 10, Geometry, Mathematics, Physics

This was something I did a while ago and never got around to posting about. It’s a simulation of an incompressible fluid- when you render the density and apply temperature by adding an upward force, it looks like smoke!

It’s fairly computationally intensive. The state of the system has to be solved simultaneously. But because it uses an iterative solver to approximate a solution, it’s easy to tune for performance to realism.

Click and drag to knock the smoke around, release to add density (more smoke).

My implementation is based on this paper by Joe Stam. The leap from the Navier Stokes equations to this compact form is pretty radical. Mr. Stam’s a pretty smart guy.

drew ActionScript 3.0, Mathematics, Physics

For the past few months I’ve been working on a general purpose geometry library for representing your interactive elements in a more precise geometrical way. There are some pretty cool features like collision detection management, auto-generation of geometry, and an easy way to bind the model to your DisplayObject hierarchy (or not!). 

It will be released under the MIT license alongside SLRBS- a physics engine that uses the library. The absolute last thing the Flash world needs is another physics engine, but this has some really nice features for beginner developers, prototypers etc. And it does a great job of showing how to use the geometry library.

One of the most frequent questions or requests I get/got with FOAM was something along the lines of “how do I make my MovieClip a rigid body?” It can be hard from my position to know how to tell someone who probably uses the Flash authoring tool a proper way to handle the difference between the model and view. So with SLRBS, I set out to address this issue.

With SLRBS, you have a typical API to create physics simulation, but it also has some streamlined auto-creation methods that make setting it up a breeze. To ensure this, I even downloaded a copy of Flash CS4 to play with it!!! AAAAHH!!! wtf happened to Flash!? 

So, all you do is lay out your MovieClips on the stage, and call a single static method and voila! physics simulation! There’s some pretty cool stuff going on that determine how to handle overlapping DisplayObjects (each demo does it differently), and it’s just a single parameter.

Here’s a screenshot of the stage in a demo where a car is being created:

Car Demo

Because 2 of the tires overlap the car, they get pinned to the car at their center of mass.

This is the only ActionScript in the .fla:

1
2
3
4
5
6
7
8
9

import org.generalrelativity.slrbs.SLRBSimulator;
 
var simulator:SLRBSimulator = SLRBSimulator.createSLRBSimulatorFromContainer( this );
addEventListener( Event.ENTER_FRAME, onEnterFrame );
 
function onEnterFrame( event:Event ) : void
{
   simulator.step( 1 / 60 );
}

Here’s a second demo’s stage in the Flash authoring tool- this time a ragdoll. The overlap this time gets pinned at the center of overlap instead of the center of mass as in the car demo- this works well for creating ragdoll joints.

Ragdoll Demo

drew ActionScript 3.0, Flash Player 10, Mathematics, Physics

The next release of the Flash Player has support for 2.5D- which boils down to 3D properties, appropriate perspective, but no depth sorting. This method takes a DisplayObjectContainer, transforms its children into world space and sorts them along that z-axis. It’s sloppy because it’s DisplayObject-based and so there will be noticeable “pop” when a DisplayObject is sorted onto the top of the display list in close proximity to another. That said it should be perfect for carousels and other views that aren’t super tight.

Here’s a demo

And here’s the algorithm:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

public static function simpleZSort3DChildren( doc:DisplayObjectContainer, recurse:Boolean = true ) : void
{
 
    //transforms from local to world oordinate frame
    var transform:Matrix3D = doc.transform.getRelativeMatrix3D( doc.stage );
 
    var numChildren:int = doc.numChildren;
 
    //v = ( n * 3 )- (x,y,z) set for each child
    var vLength:int = numChildren * 3;
 
    var vLocal:Vector.<Number> = new Vector.<Number>( vLength, true );
    var vWorld:Vector.<Number> = new Vector.<Number>( vLength, true );
 
    //insertion point for child’s coordinates into state vector
    var vIndex:int = 0;
 
    for( var i:int = 0; i <numChildren; i++ )
    {
 
        var child:DisplayObject = doc.getChildAt( i );
        if( recurse && child is DisplayObjectContainer ) simpleZSort3DChildren( DisplayObjectContainer( child ), true );
 
        vLocal[ vIndex ] = child.x;
        vLocal[ vIndex + 1 ] = child.y;
        vLocal[ vIndex + 2 ] = child.z;
 
        vIndex += 3;
 
    }
 
    //populates vWorld with children coordinates in world space
    transform.transformVectors( vLocal, vWorld );
 
 
 
    //bubble sorts children along world z-axis
    for( i = numChildren - 1; i > 0; i-- )
    {
 
        var hasSwapped:Boolean = false;
 
        vIndex = 2;
 
        for( var j:int = 0; j < i; j++ )
        {
 
            //z value at that index for each child
            var z1:Number = vWorld[ vIndex ];
 
            vIndex += 3;
 
            var z2:Number = vWorld[ vIndex ];
 
            if( z2> z1 )
            {
 
                //swap
                doc.swapChildrenAt( j, j + 1 );
 
                //swap z values (don’t need to change x and y because they’re not used anymore)
                vWorld[ vIndex - 3 ] = z2;
                vWorld[ vIndex ] = z1;
 
                //mark as swapped
                hasSwapped = true;
 
            }
 
        }
 
        //if there was no swap, we don’t need to iterate again
        if( !hasSwapped ) return;
 
    }
 
 
}

drew ActionScript 3.0, Flash Player 10, Mathematics

This is the first in a series of quick tutorials on helpful maths, physics and algorithms, and their ActionScript implementations. Today we’re going to learn how to shoot a projectile out of the sky. Sweet.

DEMO

First, let’s define the algorithm’s structure. We’re going to pass the time to collision. The algorithm will return the muzzle velocity that our projectile requires to hit a moving target out of the air after this specified amount of time.

function solveMuzzleVelocity( timeToImpact:Number ) : Vector

Next, we have to define the equations of motion:

x1 is a particle’s position after t seconds. x0 is the original position and the dots represent the first and second positional derivatives (velocity and acceleration, respectively). This equation allows us to solve for position after any specified time (given constant acceleration). Referring to our algorithm’s structure, we know t (timeToImpact) to be defined. This means that we can determine where the target will be when we hit it, which gives us a point of impact:

//point of impact
var poi:Vector = target.position.plus( target.velocity.times( timeToImpact ) );
poi.plusEquals( new Vector( 0, Projectile.G * timeToImpact * timeToImpact * 0.5 ) );

It’s assumed that the only accelerative force is gravity (Projectile.G here), that we know the starting position of the target and our projectile, and the current velocity of the target.
Plugging this point of impact in for our projectile’s equation of motion (as x1), we eliminate all unknowns except the muzzle velocity. Through basic algebraic manipulation, we can isolate and solve for this unknown:

Finally, divide each side of the equation by t:

The algorithm will return this value. Here’s my implementation, in full:

1
2
3
4
5
6
7
8
9
10
11
12
13

private function solveMuzzleVelocity( timeToImpact:Number ) : Vector
{
 
    //point of impact
    var poi:Vector = target.position.plus( target.velocity.times( timeToImpact ) );
    poi.plusEquals( new Vector( 0, Projectile.G * timeToImpact * timeToImpact * 0.5 ) );
 
    var diff:Vector = poi.minus( projectile.position );
    diff.y -= Projectile.G * timeToImpact * timeToImpact * 0.5;
 
    return diff.dividedBy( timeToImpact );
 
}

Super simple! Now you can integrate dynamically and knock some shit out of the sky.

Here’s a verbose implementation that doesn’t use my Vector class and isn’t dependent on my application:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

private function solveMuzzleVelocity( t:Number, xp:Point, xt:Point, vt:Point, pixelToMeterRatio:Number = 40 ) : Point
{
 
    const g:Number = 9.81;
    var vDotTerm:Number = g * pixelToMeterRatio * t * t * 0.5;
 
    var poiX:Number = xt.x + vt.x * t;
    var poiY:Number = xt.y + vt.y * t + vDotTerm;
 
    var diffX:Number = poiX - xp.x;
    var diffY:Number = poiY - xp.y - vDotTerm;
 
    return new Point( diffX / t, diffY / t );
 
}

t is time to impact, xp holds the x and y coordinates of the projectile’s starting position, xt is the target’s starting position, vt is the target’s starting velocity and pixelToMeterRatio defines how many pixels represent one meter.

drew ActionScript 3.0, Mathematics, Physics

I’m sorry that this took so long- especially those of you who’ve written me and patiently waited for my help with rendering… I simply don’t have as much time as I would like. Anyway- this will offer a solution for rendering your FOAM physics simulation with:

• Drawing API (SimpleFoamRenderer)
• BitmapData (BitmapFoamRenderer)
• DisplayObjects (DisplayObjectFoamRenderer)

Here’s a very silly demo- Mouse down anywhere to switch from the DisplayObject renderer to the Bitmap renderer.

demo
source

Note that these new renderers are meant as example only- I do think in most cases the developer will need a more specific setup than I could hope to generically offer. So consider this my implementation and at best a nudge in the right direction.

Hopefully both of these will be fairly self-explanatory, but I’ll briefly go over the DisplayObject renderer and what my thinking was when addressing this issue.

I wanted a means to render a DisplayObject (base class of Sprite, MovieClip etc.) that’s the visual representation of a RigidBody being simulated. I knew that I’d like to be able to pass either the DisplayObject’s Class or the instance itself. I also knew that I’d need an easy way to adjust the position of the DisplayObject- so first, I wrote a datatype to hold these properties- DisplayObjectData.

You’ll see that it has a setter for defining its displayObject property. This is what facilitates passing either a DisplayObject Class or a DisplayObject instance. The other properties:

• offsetX - the amount to offset the DisplayObject from its origin in the x direction
• offsetY - the amount to offset the DisplayObject from its origin in the y direction
• autoCenter - Boolean indicating whether to center the contents about the origin- Because of the way a RigidBody rotates about its center of mass, I figured this would be typical.
• hasBeenDisplayed - a helper Boolean for determining whether to perform the indicated transformations

In Rocket.as, you can see the way I’m embedding my visual assets:

[Embed(source="../../../bin/images/rocket/lunarPod.png")]
private var rocketSpriteClass:Class;
 
[Embed(source="../../../bin/images/rocket/asteroid.png")]
private var asteroidSpriteClass:Class;

If this is unfamiliar with you, simply note that it’s analogous to the classes you’re instantiating as MovieClips in Flash CS3. As you’ll see in a moment, it’s not necessary to pass the Class, we can pass the DisplayObject instance as well.

Here’s how I’m passing the class to FOAM to tie it to to a RigidBody. Note that the instance of DisplayObjectData simply defines Foam.addElement’s renderData argument.

foam.addElement( rocket, true, true, new DisplayObjectData( rocketSpriteClass, 0, -5 ) );

Compare this with how I’m creating the asteroid (I wanted to scale the asteroid, or else I’d have created it in the same way as the rocket- it worked out well for the purpose of demonstration):

var asteroidSprite:DisplayObject = new asteroidSpriteClass();
asteroidSprite.scaleX = asteroidSprite.scaleY = size / 64;
var asteroid:Circle = new Circle( sx, sy, size, size * 4, svx, svy, 0.5, 0.25, 0, -0.1 + Math.random() * 0.2 );
foam.addElement( asteroid, true, true, new DisplayObjectData( asteroidSprite ) );

There’s one last important thing to note in the Rocket example class. Because I don’t want to use the default FOAM renderer, I have to tell it which renderer to use:

foam.setRenderer( new DisplayObjectFoamRenderer() );

Now let’s move on to DisplayObjectFoamRenderer. Note that I’ve chosen to extend SimpleFoamRenderer so that I don’t have to rewrite drawing routines for objects I don’t have graphics for (such as the mouse dragger bungee). You’ll see my use of the DisplayObjectData properties in the addRenderable method.

So- ultimately all of this boils down to 3 lines of code which line up our DisplayObject with our RigidBody:

displayObject.x = ISimulatable( renderable.element ).x;
displayObject.y = ISimulatable( renderable.element ).y;
if( renderable.element is IBody ) displayObject.rotation = IBody( renderable.element ).q * 180 / Math.PI;

It places the DisplayObject at the x and y coordinates of the RigidBody’s center of mass. Because rotation is given in radians, we have to convert to an angle which is what DisplayObject.rotation expects.

And that’s pretty much it; for the most part a FOAM renderer is bookkeeping. It might suit you far better to just maintain reference to your DisplayObject and update it every frame just as I am in the above 3 lines of code… simply set render to false in Foam.addElement.

drew ActionScript 3.0, FOAM, Mathematics, Physics