I saw a video online recently of several cubes flying into the scene and stacking up to form a larger cube. I thought it was a pretty simple yet powerful effect. The effect was achieved using Particle Flow, so the movement of the particles are entirely animated and not in real time. After seeing this, I immediately wanted to try it, but with a more general solution. So, off I went to write a physics system for my particle engine.
Particle System
I began this project with an old particle system I started working on to dig into DX9. It’s gone through a few iterations, some architectural and some performance oriented. This version is entirely object oriented, where each emitter, particle, camera, input, etc. are classes. Particle is an abstract class, inherited by a Cube particle in these examples. Most of the details regarding the particular graphics API being used have been abstracted into the respective classes, and the scene can be described without using D3D function calls.
As a performance enhancement, this particle system uses hardware instancing. All particle specific information is transferred to the graphics card via a vertex buffer, with all position, color, and other particle logic performed in the vertex shader. The next major revision will include physics calculations to be performed on the GPU.
Spring Physics
Now, there are a lot of challenges in getting a particle rocketing through 3D space to land in exactly the right position with physics alone. There are a few options. With physics, we can exponentially increase a force on the particle in the direction of the destination with relation to the distance from the destination. Another possibility is merging a pure physics system with a pure spline animation to allow a random initial velocity, but a guaranteed stopping point.
My first step however, is developing some way to attract the particles to their destination. I decided to start with a spring system. Although the system will exert forces on the particles to aim the particles at the destination, it is very unlikely that they will pass directly through a predetermined point, and instead orbit around the destination.
For those interested, here’s the nice little gem used to compute the spring force. I modified this algorithm from a cloth simulator used in the book “Physics for Game Developers”.
D3DXVECTOR3 Collector::GetForce(Particle particle)
{
D3DXVECTOR3 forceVector = GetVector(particle); // From particle to dest
float length = D3DXVec3Length(&forceVector);
float k = 1; // tensile constant
float d = 1; // damping constant
float dot = D3DXVec3Dot(&particle.m_velocity, &forceVector);
D3DXVECTOR3 force = (k * length + d * (dot / length)) * (forceVector / length);
return force;
}
Orbits and Decay
Attempt number 2 at getting particles to stick to a collector: apply gravitational pull to the collector. The ideal is simple; the closer a particle gets to the collector, the greater the force of gravity acting on the particle towards that collector. In this case, I used the inverse of the distance squared. This produces an “Ideal” orbital pattern as shown below:
D3DXVECTOR3 Collector::GetIdealOrbitalForce(const D3DXVECTOR3& magVector,
const float& magLength)
{
float grav = 1000.0f;
D3DXVECTOR3 force = magVector * grav;
if(magLength > 1)
{
force /= (magLength * magLength);
}
return force;
}
This ideal orbital pattern has the obvious problem of not actually attracting the particles into the center, and instead just pulling the particles into an orbit around the collector. To give more realistic physical behavior, we would like to model drag the particle encounters by entering the “atmosphere” of the collector. The thickness of the atmosphere was modeled as the inverse of the distance squared. Tweaking the constants a bit yields the following results:
D3DXVECTOR3 Collector::GetOrbitalDecayForce(const D3DXVECTOR3& magVector,
const float& magLength,
const D3DXVECTOR3& curVelocity,
const D3DXVECTOR3& curPos,
const D3DXVECTOR3& collectorPos,
const float& dt)
{
D3DXVECTOR3 force = GetIdealOrbitalForce(magVector, magLength);
D3DXVECTOR3 newVelocity = curVelocity + force * dt;
if (dt > 0)
{
D3DXVECTOR3 newPos = curPos + newVelocity * dt;
D3DXVECTOR3 magVector2 = collectorPos - newPos;
float magLength2 = D3DXVec3Length(&magVector2);
float damping = 4.0f;
D3DXVECTOR3 decayVelocity = ((magVector2 / magLength2) - (magVector / magLength)) / dt * damping;
if(magLength2 > 1)
{
decayVelocity /= (magLength2 * magLength2);
}
force += decayVelocity/dt;
}
return force;
}
Adding this behavior with spring physics we get:
This all adds to make a pretty interesting effect. I hadn’t optimized the code which models the orbital decay, so I was curious to see what the limitations of the particle system were after performing the physics computations. The physics for each particle is computed on the host, with the final position stored in the instance buffer for each cube to be drawn on the GPU. I’m interested to see what kind of speedup is achievable by performing the physics on the GPU, but I’ll save that for another day. For now, I was able to achieve 70 FPS on my GTX260 with 20,000 cubes, each with spring physics and orbital decay. Results speak for themselves:
After writing these physics calculations, I’m beginning to believe that a purely physically based solution to get the kind of flight paths I am interested in does not exist. So, the third iteration of this particle system will use a spline hybrid to guide the particles along a path shaped by the physics calculations, but also guaranteeing the particles will reach their destination within a fixed amount of time.
Collector Corner Cases
Well, it turns out I was able to guarantee the destination of particles by a deadline. This requires a little bit of cheating; as the lifetime of the particle advances, the velocity vector dictated by the physics is blended with a velocity vector of the same magnitude pointing directly at the destination. This, in combination with a maximum velocity can virtually guarantee that the particle will hit its destination by a specific time. If the velocity vector is blended to point directly at the destination by half the lifetime, and the forces acting on the particle only contribute to increase the velocity towards the destination regardless of its position, then the particle is guaranteed to hit the destination by the end of its life.
The following code accumulates the forces acting on the particle, blends the velocity vector with a vector pointing directly at the destination based on lifetime percentage, and detects if the particle passed through the destination in order to “stick” it to the destination. The code isn’t optimized, and is a bit messy, so take it with a grain of salt. It’s also not quite as elegant as I would like. Notice that within a distance of 1, the velocity vector points directly into the collector. This can be considered the event horizon; to save me the headache of dealing with corner cases when calculating forces given a magVector length < 1.
void Collector::UpdateParticle(float dt, Particle& particle, UINT index)
{
D3DXVECTOR3 magVector = GetVector(particle, index);
float magLength = D3DXVec3Length(&magVector);
if (magLength >= 1)
{
D3DXVECTOR3 force = D3DXVECTOR3(0.0f, 0.0f ,0.0f);
force += GetSpringForce(magVector, magLength, particle.m_velocity);
force += GetOrbitalDecayForce(magVector, magLength, particle.m_velocity,
particle.m_pos, GetDestination(index), dt);
particle.m_velocity += force*dt;
// Calculate blended vector
float velLength = D3DXVec3Length(&particle.m_velocity);
float perc = ((float)particle.m_timeElapsed / ((float)particle.m_lifetime)/2);
perc = min(perc, 1.0f);
perc = pow(perc, 1);
D3DXVECTOR3 vec1 = particle.m_velocity / velLength * (1.0f - perc);
D3DXVECTOR3 vec2 = magVector / magLength * perc;
D3DXVECTOR3 vec3 = vec1 + vec2;
D3DXVec3Normalize(&vec3, &vec3);
velLength = min(velLength, MAX_VELOCITY);
particle.m_velocity = vec3 * velLength;
particle.m_pos += particle.m_velocity * dt;
}
else
{
D3DXVec3Normalize(&magVector, &magVector);
float velLength = D3DXVec3Length(&particle.m_velocity);
velLength = min(velLength, MAX_VELOCITY);
particle.m_velocity = magVector * velLength;
particle.m_pos += particle.m_velocity * dt;
}
// Fix position if passing through particle
D3DXVECTOR3 newMagVector = GetVector(particle, index);
D3DXVec3Normalize(&newMagVector, &newMagVector);
float dot = D3DXVec3Dot(&magVector, &newMagVector);
if ( dot <= (-1.0f + EPSILON) && dot >= -1.0f)
{
particle.m_pos = GetDestination(index);
}
}
So, with this new code, I decided to make a nice little demo. Consider this a sneak peak at the various shapes and animations possible by associating individual particles at different collectors:
Shapes
This update adds the ability to create shapes using a particle collector. The physics for each particle is calculated independently for individual destinations; they do not need to have the same destination. This allows me to specify any number of destinations, which can be moved independently. Because the physics system takes care of tweening between positions, the collector destinations can jump from place to place, and the physics will result in a fluid animation.
There was a big feature implemented in this update as well. The last videos I captured of the particle collector had a particularly chaotic nature, because the framerate was tied to the simulation timestep (update and render once per frame), and because the framerate drops during video capture, the numerical integrator in my update phase became far less accurate. There is an excellent article by Glenn Fiedler at his site, gafferongames.com. The code below demonstrates how Glenn’s example fixed my timestep issues, but all credit should really go to him; it’s simply included here to provide another example.
// global declarations
ULONG timeNow; // the timestamp of the current frame
ULONG timePrev; // the timestamp of the previous frame
ULONG timeFrame; // the time delta for the current frame
ULONG timeSum; // accumulation of time
const ULONG timeDelta = 8; // the constant update time delta
const ULONG timeMax = 16; // the constant maximum time delta
// function prototypes
void InitScene(); // initialize the scene
void Update(float dt); // updates the simulation
void RenderFrame(); // renders a single frame
void DestroyScene(); // free memory allocated by the scene
// the entry point for any Windows program
int WINAPI WinMain(HINSTANCE hInstance,
HINSTANCE hPrevInstance,
LPSTR lpCmdLine,
int nCmdShow)
{
timePrev = timeGetTime();
timeNow = 0;
timeSum = 0;
InitScene();
MSG msg;
while(TRUE)
{
while(PeekMessage(&msg, NULL, 0, 0, PM_REMOVE))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
if(msg.message == WM_QUIT)
{
break;
}
timeNow = timeGetTime();
timeFrame = timeNow - timePrev;
timePrev = timeNow;
// practical limitation, don't let updates overcome framerate
if (timeFrame > timeMax)
{
timeFrame = timeMax;
}
timeSum += timeFrame;
// update at a rate of timeDiff, regardless of framerate
while (timeSum >= timeDelta)
{
Update(timeDelta * .001f);
timeSum -= timeDelta;
}
// render once per frame
RenderFrame();
}
DestroyScene();
return msg.wParam;
}
There were also a number of crucial bugfixes and optimizations I made. My last post has been updated to reflect the bug fixes. A particular performance enhancement was the calculation of the orbital decay force. I wasn’t happy with the previous code, which looked something like this:
D3DXVECTOR3 Collector::GetOrbitalDecayForce(const D3DXVECTOR3& magVector,
const float& magLength,
const D3DXVECTOR3& curVelocity,
const D3DXVECTOR3& curPos,
const D3DXVECTOR3& collectorPos,
const float& dt)
{
D3DXVECTOR3 force = GetIdealOrbitalForce(magVector, magLength);
D3DXVECTOR3 newVelocity = curVelocity + force * dt;
if (dt > 0)
{
D3DXVECTOR3 newPos = curPos + newVelocity * dt;
D3DXVECTOR3 magVector2 = collectorPos - newPos;
float magLength2 = D3DXVec3Length(&magVector2);
float damping = 4.0f;
D3DXVECTOR3 decayVelocity = ((magVector2 / magLength2) - (magVector / magLength)) / dt * damping;
if(magLength2 > 1)
{
decayVelocity /= (magLength2 * magLength2);
}
force += decayVelocity/dt;
}
return force;
}
The problem with the above code is that it computes the position from the force and velocity of the particle, calculates a vector from the previous position and the new position, generates a new velocity vector, and divides by the timestep to get back to a force. Well, I coded it this way just to validate my initial thought process. Now that I know it can achieve the sort of effect I’m looking for, I whittled it down to something like this:
inline D3DXVECTOR3 Collector::GetOrbitalDecayForce(const D3DXVECTOR3& magVector, const float& magLength,
const D3DXVECTOR3& velVector, const float& velLength)
{
D3DXVECTOR3 force = GetIdealOrbitalForce(magVector, magLength);
D3DXVECTOR3 magVector2 = (magVector) - (velVector);
float damping = 100.0f;
force += magVector2 * velLength * damping / (magLength * magLength);
return force;
}
Note that the vectors are now passed in normalized, and their length is included. I made this change to each of my force computations, since each force function would eventually calculate the vector length and normalize. The magVector always points at the collector, and the current velocity vector indicates the direction the particle is moving. We wish to apply some force acting against the particle as it tries to move laterally throughout its orbit of the collector. The simple vector subtraction yields a vector that will act against the lateral movement of the particle, as well as contribute to the force aimed directly at the collector. While this may not be entirely physically accurate, it does yield the same effect, with far less work. It’s good to remind yourself, sometimes you don’t need perfect physical accuracy, an approximation may be good enough.