Email Updates RSS Subscribe
Line

This blog is created and maintained by the technical team at Hook in an effort to preserve and share the insights and experience gained during the research and testing phases of our development process. Often, much of this information is lost or hidden once a project is completed. These articles aim to revisit, expand and/or review the concepts that seem worth exploring further. The site also serves as a platform for releasing tools developed internally to help streamline ad development.

Launch
Line

Hook is a digital production company that develops interactive content for industry leading agencies and their brands. For more information visit www.byhook.com.

Line

Particle-based Rigid Bodies using Shape Matching

Line
Posted on June 29th, 2010 by Brad
Line

While hopping between projects the idea of a particle-based approach to rigid-body kinematics soon became my white whale, obsession ensued and I quickly began filling enough notebooks to fill a John Doe apartment. While the mechanics of rigid-bodies can be fairly simply described and implemented with linear and angular properties, my goal with ‘fuzzy physics’ has been to at least test out the possibilities of a different approach to physics-based movement, one a little more accessible to someone that is programmer by trade but animator by craft (in other words, can’t do a Jacobi rotation to save his life)

The main challenge of rigid-body kinematics with particles is properly representing rotational dynamics, as particles by definition can only move linearly. That linear movement needs to then be translated to the other particles rigidly bound to the acting particle (of course, it is more likely they will ALL be moving and influencing each other). But to consider just the single particle in 2D (who cares about 3D, really?), the old lever & fulcrum model can’t be beat.

What this illustrates is that the movement of a particle can be broken down into components parallel and perpendicular to its relationship to another point in space.  The perpendicular component is what imparts the angular movement around the reference point while the parallel is translational.  Most importantly, rotational movement requires a point of reference, and a description of the particle relative to that point.  That point of reference for rigid-bodies is usually the center of mass, which in the classic representation becomes the single frame of reference for the body’s movement.  Not so simple with particles flying independently all over the place.

While much of the solution seems inherent to the problem described above, of course it gets more complicated.  What I am after is something that works with the position based dynamics I’ve been playing with, something that binds multiple particles together while conveying plausible rigid-body dynamics.  And to slather on the hurt, one that can have particle effectors added on the fly for collision handling and other fun features.  This is in keeping with the collision handling as described by Jakobsen (half-way down), which was the basis of my recent post dealing with more simplified ‘stick’ bodies.

The descent into rigid-body hell

I want to give a brief overview of the techniques I explored on the way to something that actually worked, as there were a fair number of them with varying success, but feel free to spare yourself more rambling and jump down to the good stuff.

Feeling optimistic about the stick results from my last post, I went straight for a 2-Dimensionalized version of the tetrahedron approach as shown in Advanced Character Physics, this is what was used for my quick haXe to HTML5 demo.  While in those simplified terms it got the job done, that method (or at least, my implementation of it) relied too heavily on distance constraint relaxation to maintain rigidity.  Much like the ‘stick’ physics post, the triangle was deformed by every correction.  Under minimal conditions this was acceptable, the 3 distance constraints relaxed over several iterations, but under more drastic movements the triangle particles would flip or invert.  While distance constraints are great for various things, they started to feel like Band-Aids in this case.

Second, I tried out a modification of the stick model, in which particle effectors could be attached anywhere, even outside of the line between the two particles.  Using some barycentric trickery, a center of mass was found as well as a correlative point along the stick, by which both points would be moved while maintaining distance.  This worked great for a single effector and could likely be used to simulate a more robust IK system, but when multiple effectors were introduced so were strange ghost-forces resulting in some really cool perpetual chaotic pendulums… just not quite the predictable system desired.

Last up was a 3-particle model.  This used a similar technique to above but had the added benefit of a well defined orientation vector and system for defining internal points.  Each particle effector’s movement created a virtual configuration of the model, which when blended together created an ideal position and orientation depending on embedded particle movement.  The idea here was that you could peg a single particle effector at the center of mass and just let it go, then add/remove particles as need dictated.  The first part worked great, the second suffered from more ghost forces.  Most likely this was due to the lack of a real center of mass, the planned emergence of a virtual one never materialized.

All hope was pretty much lost, but there was still the idea of ‘shape matching’ from another paper Muller et al. This one I saved for a last resort because, well, the math in it is just plain SCARY.  But, sometimes you just gotta punch that fear in the throat and go for it, after all, this is 2D, and angles in 2D are EASY.

Shape matching

After taking the time to work this one out on paper, I was treated with one of those wonderful ‘AH HA’ moments that really makes you feel like a total idiot for overcomplicating every experiment you did beforehand.  While the paper is primarily about deformable objects, the initial description of the technique is fairly general and as stated by the author’s themselves, is well suited for rigid-body behavior.  The gist of the method is storing a configuration of the particles you want to bind together.  For each particle you store a displacement vector (q) from the center of mass, which is analogous to the lever/fulcrum diagram linked above, and provides exactly enough information to describe the particles default relationship to the center of mass.  An optimal position and orientation of the ‘body’ can be determined based on any future configuration (movement) of the bound particles by comparing that stored displacement vector to the time-stepped displacement (p, the mirror of q) from the re-located center of mass. In my implementation I am using Verlet style integration and running the code as a position-constraint vs. modifying velocity via Euler integration.  This choice was made to keep things in line with previous work, though the more force-based approach is something worth exploring, especially for deformable bodies.

The following example shows the rigid-body configuration as well a ‘step in time’ illustration of what goes into the correction on each step, toggled using the space bar.  In the configuration view you can create and delete particles and arrange the placement and mass distribution (using the mouse-wheel) of each, then in the ‘time step’ view you can adjust the placement of each particle as they would move during any time step to see the resulting translation and rotation.  While in time-step mode, the colored arrows are shaded green or red to indicate the direction of rotational influence of that individual particle.

Get Adobe Flash player

You should see that the movement is largely dependent on the placement of the center of mass, while ‘heavier’ points are more effective in imparting movement than ‘lighter’ ones.

The code consists of two main functions, a configuration routine and a relaxation routine.  The configuration routine is simply resetting the center of mass and the displacement vectors (q), in psuedo-code:

// initialize the center of mass and system mass to 0
centerOfMass.reset (0, 0);
systemMass = 0;
 
// go through the particles and add their positions, weighted by mass
// this will give 'heavier' particles greater influence over the movement of the system
 
for (i = numParticles; --i >= 0; ) {
   centerOfMass += particle[i].pos * particle[i].mass;
   systemMass += particle[i].mass;
}
 
// divide by the total system mass to normalize out the particle weighting from above
 
centerOfMass /= systemMass;
 
// with the new system of mass we need to now recalculate the relative position of each particle
// to the center of mass
 
var qs:Array = [];
 
for (i = 0; i < numParticles; i++) {
   qs[i] = particles[i].pos - centerOfMass;
}

Every time we make a change to the system, be it a mass or placement adjustment or by adding or removing a particle, we need to reconfigure it. Nothing really fancy there, thats just how you find a center of mass. For each iteration of the solver, in each time-step, this ‘relaxation’ routine is run to bind the particles together and distribute the movement of the individual particles to each other:

// as above, we need to find our current center of mass
 
var currentCOM:Vector2 = new Vector2 (0, 0);
 
for (i = numParticles; --i >= 0; ) {
 
    currentCOM += particles[i].pos * particles[i].mass;
 
}
 
currentCOM /= system.mass;
 
// with the center of mass found, we can go through each particle
// and calculate its current displacement from the center of mass to
// calculate a new orientation in the form of a total angular displacement (1 dimensional)
 
var angleDelta:Number = 0;
 
for (i = numParticles; --i >= 0; ) {
 
    var p:Vector2 = particle[i].pos - currentCOM;
    var q:Vector2 = qs[i];                                    // this is the stored default displacement from the COM, our 'home' position if you will
 
    // this used the dot and cross product values to find the current
    // angle, a more detailed explanation on this to come
    var cos:Number = p.x * q.x + p.y * q.y;
    var sin:Number = p.y * q.x - p.x * q.y;
 
    // using those two values we can calculate an angle, weighted by the particle's mass
    // this coupled with the shifted center of mass should give the feel of diminishing rotational effect based on particle mass and distance
    angleDelta += Math.atan2 (sin, cos) * particle[i].mass;
 
}
 
// as we did for the center of mass, divide the angleDelta by the system mass to normalize it
angleDelta /= systemMass;
 
// we now have the two values we need to calculate the goal positions of each particle:
// location of the center of mass and the orientation
// the goal position is then simply the default displacement vector (q) rotated by our orientation
// and added to the center of mass
 
// we need to first convert our angular displacement to the values of a transformation matrix
 
cos = Math.cos (angleDelta);
sin = Math.cos (angleDelta);
 
// to avoid using any function calls to a matrix object for what would be an easy 2x2 matrix operation we just use the following
// to find the target displacement vector
// x = cos * q.x - sin * q.y
// y = sin * q.x + cos * q.y
 
for (i = numParticles; --i >= 0; ) {
 
    // start with the goal displacement vector, using
 
    var g:Vector2 = new Vector2 (cos * qs[i].x - sin * qs[i].y, sin * q.x + cos * q.y);
    g += currentCOM;
 
    // now to just move the particle towards its goal
    // this can be adjusted by a strength variable k to create the appearance of a softer rigid-body system
 
    particles[i].pos += (g - particles[i].pos) * k;
 
}

Just plug that into the relaxation loop of the particle system and magic ensues:

Get Adobe Flash player

The rotational values above are based on this rule of the dot product of two vectors a and b:

a{circ}b={vert}a{vert}{vert}b{vert} cos theta or: cos theta={a{circ}b}/{{vert}a{vert}{vert}b{vert}}

To find the sine value, we use the same formula but with the perpendicular vector of b, or the cross product of a and b:

a*b={vert}a{vert}{vert}b{vert} sin theta or: sin theta={a*b}/{{vert}a{vert}{vert}b{vert}}

Now, since the Math.atan2 is used to find the angular displacement per particle based on a ratio of these two values, the length measurements can be cancelled out.

While I haven’t run this approach through any serious frame-rate annihilation tests, the results thus far have been very promising and surely warrant some more experimentation down the line.  One issue worth noting is that because the default ‘shape’ is stored without any sort of orientation information, angular movement limits become a little tricky.  But still, the absolute simplicity of the position based dynamics/particle approach is maintained while adding new features without just throwing in more distance constraints to keep things together.

Line
24 Responses to “Particle-based Rigid Bodies using Shape Matching”
  1. Niall says:

    Hey Brad, unsure if you’re still alive!

    I’ve been revisiting this post and the method behind it, and running into some troubles. Running into a strange angle ‘snapping’ issue, I’ve run a variety of tests and can’t figure out the source of it.

    I’ve tried both methods of calculating angle differences, your way and using precalculated angles and finding the difference between them. Ignoring particles and directly setting positions, here’s what it looks like – http://filz.us/nrH
    ..and here it is moving by increments of .1*distance – http://filz.us/nrJ

    This is the bulk of my code, would you be able to take a look and provide some guidance? I’ve been staring at it for hours and asking friends to help and none of us can figure it out. http://pastebin.com/9g9eJqTV

  2. brad says:

    @kyle and mike: sorry for delay in responses to your earlier posts, our spam filter got a little overzealous.

    @kyle: Some of this behavior looks familiar, and there is definitely some weirdness going on. If you want to email me a link to the files (brad@byhook.com), I can take a peak and try and figure out what exactly the problem is.

    @mike: The feel of moving through a viscous fluid (air) was something I was going for here, and that was achieved through a basic dampening system as suggested by Jakobsen. While I’m not suggesting that the simple velocity dampening is the best way to go about it, I feel that it does the trick for basic motion pretty well.

  3. Mike says:

    Hi Brad,

    I’ve solved the problem in the mean time. Explanation and proof-of-concept programs here:

    http://www.gamedev.net/community/forums/topic.asp?topic_id=579582

    More is coming up.

    Cheers,
    Mike

  4. brad says:

    Hey Mike, we’ve been pretty busy here and I haven’t had the availability to devote enough brain bandwidth to it at home to really get a good start on it yet. But I can tell you my hypothesis going into experimentation. The problems I am mostly trying to address with the above solution are:

    1. The lack of orientation information on the body as a whole for movement restrictions and better shape recognition
    2. The cost of embedding new points/effectors on the fly

    So my intuition on that is leading me to try something that has a more rigidly defined origin and orientation and the particle effectors contribute to the movement of those two values based on a weighted 2-degree-of-freedom routine, where the rotational and translation effects can be independently weighted need be. I’m hoping running those effectors in sequence much like the standard constraints will then cause a rigid-body behavior to emerge. The major issue that I’ve been researching is how exactly to deal with the rotational interpolation in a way that is inexpensive enough to run numerous times per update.

    It’s a lot easier to explain in drawings

  5. kstyle says:

    Great Tutorial!

    I used the code provided and was able to reproduce your results fairly well. (Green dot is center of mass)
    http://www.truploader.com/view/586827

    However, when I tried to use the same set up with just 3 dots, it freaked out.
    http://www.truploader.com/view/176019

    Any idea what might be causing this? (I can show you class files and a .fla if it would help

  6. Mike says:

    Hi Brad, how’s it going with your alternate approach?

  7. Mike says:

    Hi Brad,

    The reason I think your model does not conserve momentum, is because it feels like your objects are moving through a viscous fluid. If f. inst. you make the soft body on a spring pendulate from side to side, it stops moving again within a few swconds. An energy / momentum preserving model would keep the body pendulating back and forth indefinitely.

    As for my soft body damping algorithm, I’ve made these small interactive open source apps, which are heavily expanded versions of Maciej Matyka’s pressurized soft body model (which I assume you know?), including among others angular springs, so you can give the body a distinct shape, and the mentioned global damping system for making very viscous bodys (notice how these do preserve momentum even though they are damped):

    http://www.jernmager.dk/stuff/cpp/soft_body_03.zip
    http://www.jernmager.dk/stuff/cpp/softbody_glfw.zip

    I would enjoy debating this in greater depth with you and anyone likeminded. May I suggest that we continue this discussion over at the gamedev.net forum? They have a nice math & physics section where this would fit in very neatly :-)

    Cheers,
    Mike

  8. brad says:

    Thanks Mike,

    While I agree with your sentiments on good and less good, I’m a little curious as to what leads you to say it doesn’t seem to conserve linear and angular momentum, ie. is it because the code itself doesn’t explicitly conserve or that it doesn’t ‘feel’ like it?

    I’ve started an alternate approach at this in with an orientation is stored directly, rather than the ‘matched shape’ always being in the same orientation. I think this might possibly solve some issues with angular momentum and movement in general, but it unfortunately complicates the calculations significantly.

    And please let me know how that global dampening algo goes, sounds interesting!

  9. Mike says:

    Hi Brad,

    This looks very interesting, and I’m definitely implementing it asap.I read the Muller et al. paper and – like you and raigan – found the math a little too obfuscated, especially since I also “just” wanted a 2d implementation.

    It just bother me your implementation does not seem to conserve linear and angular momentum. Imo this is what separates good and less good physics simulation. I once made a global damping algo which conserved momentum, and I think a modified version of it can solve the problem – will go and check…

    cheers,
    Mike

  10. Michelle says:

    Brad,

    Just to clarify what Robbie said above about the double spacing, it’s a preference not a rule. If you ever need to run something by me for grammatical proofreading, you know where to find me.

    Michelle

  11. raigan says:

    Well that was a hell of a brilliant insight, seriously.. that’s AWESOME. Or maybe it’s obvious and I’m just hellishly stupid; either way, when I tried to implement the paper I ended up with lots of seriously complicated matrix math building and inverting/sqrting the A matrix.

    > Or are you saying the difference between calculated atan2’s may produce more accurate results

    I don’t know if it would make any difference in terms of how it behaves; I just meant to point out that you could skip a step during the calculation of a particle’s angle, because what you care about is the delta angle (current angle relative to initial angle), not the delta vector (current vector relative to initial vector).

    AFAICT, this is what’s happening in your code:
    (let p be the current COM->particle vector and q be the initial COM->particle vector)
    1) calculate the delta vector d = (cos,sin); d is “p transformed into a coordinate system centered on q”, i.e you project p onto the basis vectors q and Perp(q)
    2) convert d from delta vector to delta angle

    But you could directly calculate the delta angle (which is what you really want) by operating on angles rather than vectors; if you let p be the current COM->particle _angle_ and q be the initial COM->particle _angle_, then the delta angle is just p – q and you don’t need to mess around with vectors at all.

    ..I hope this makes sense, what I’m trying to say is very simple but as usual my explanation seems very convoluted :(

    I’d be really interested to know how “rigid’ this shape-matching is for you; the problem I had with it was that collision and constraints were unavoidably soft: imagine you have a block made of 10 particles of equal mass, and one of the particles is colliding with the ground. When you correct the collision by moving that particle, the block itself only feels 1/10th the collision response that you applied to the particle, which means the collision is only reduced, not corrected fully; the more iterations of [solve collision constraints],[solve shape matching constraints] you do the harder the collision is, but it seems like it fundamentally doesn’t want to be rigid.

    We’re successfully using the “modified stick method” you talk about; I suspect that the ghost forces you saw may have been a result of programming errors and/or miscalculated Jacobians — when we started experimenting it took us three complete rewrites before things worked, because it was so easy to screw up the math! :)

  12. Brad says:

    Wow, thanks raigan!

    I wish I could give a detailed explanation of how I arrived at my conclusions, but it really just came down to a light-bulb icon appearing above my head as I was trying to work out Muller’s formulations from that paper. I was going to attempt a ‘flattened’ 3D approach, but all the examples of Jacobi rotations implemented into code seemed far too expensive just to end up getting rid of a good portion of the matrix. So instead of implementing verbatim, I did some reflection on the purpose of the Jacobi-magic and realized (or assumed) it was just to separate out the rotational transformation from the translational for each particle’s movement, which in 2D is a bit easier just to scissor out. So that I did, and ended up with this.

    The precalculated atan2 values seem like they would work, and I definitely started out storing a lot more information per particle, but because the dot product between q and p works out for the angle pretty well I just went with that. Or are you saying the difference between calculated atan2′s may produce more accurate results (there is definitely still some weirdness around the extremes in my implementation).

    One direction I would like to take this would be to actually store an orientation vector for the ‘body’ and recalculate the q’s per step, as I’m feeling the constant and default up-ness of the current system may be leading to a less stable system. Also some more standard box sim-ing, to see how that holds up.

    You’re right about that fulcrum business, I’ll scratch that out of there. There may have been some experimentation with factoring in the default distance along with the mass weighting, but I believe it’s the repositioning of the virtual center of mass that imparts that feel to it.

    Putting these things to words is always challenging :)

    @kevin: thanks, I feel the same way… often. you should check out raigan’s from above’s robot demo’s on vimeo

    @niall: sure, my aim is: beastskeels, we are in EST how ever that may differ from where you are

  13. raigan says:

    Wow, I would *love* it if you could explain how you arrived at your method for calculating the rotation of the shape-matching frame — Muller’s paper presents the really complex solution of first building the matrix A and then decomposing it to extract the rotational component R. This is how I implemented it when I tried this out in 2D, and it sucked, since calculating R involves finding the inverse of the square root of a matrix.

    You seem to have found a direct solution that’s much simpler: average together the “rotations of the particles”, i.e you’re calculating the angle from the original CoM->particle vector to the CoM->particle vector, and using this as the particle’s “vote” on what the frame’s orientation should be (just like the particle used its position as its “vote” for where the CoM/frame’s origin should be).

    Genius!!!

    If I understand the math properly, then you could equivalently do this:

    angleDelta += (Math.atan2 (p.y, p.x) – Math.atan2 (q.y, q.x)) * particle[i].mass;

    where the Math.atan2 (q.y, q.x) term could be precalculated, making it even simpler than it currently is.

    (Then again this might not work, since taking delta between two directions described in radians tends to suck — you might need to make sure the delta is always in [0,2*pi] or [-pi,pi] or whatever.)

    Your version also seems very stable while mine was quite “energetic”.. possibly this is related.

    One thing I don’t understand is this comment in your calculate-frame-orientation code: “this is in keeping with the principle that a force has diminishing rotational influence as it moves away from the fulcrum”

    AFAICT the distance to the fulcrum isn’t considered at all.. right? Atan2() doesn’t care about the length of the vector, only its direction. But I might be totally missing the point here…

    thanks,
    raigan

  14. Kevin says:

    Wow, you just blew my mind. This stuff is really cool I’m just learning AS3 and feel like I have acomplished a lot so far, but then I stumble across something like this and realize I have a long way to go. Ha ha!

    Keep up the cool experiments!

    Kevin

  15. Niall says:

    I’m afraid I don’t know what HoN stands for / is, so I’ll have to go with no for now.

    Tell me Brad, do you have an aim/skype/msn? I’d like to discuss this stuff a bit more, plus I feel a little guilty bombarding your blog with comments…

  16. brad says:

    Ah yes, system.mass, not right, durned paste-over typos.

    As far as the semantics of effect vs. affect, I see your point, but I’m going to go with effector here as something with ‘The power to produce an outcome or achieve a result’

    Do you play HoN by any chance?

  17. Robbie says:

    Brad! I also noticed that you use double spaces after your periods. Generally this was common practice in the days of typewriters and currently still acceptable with the use of monospaced fonts, but unnecessary otherwise.

  18. Niall says:

    Sure, heh.

    Another 2 things (I swear I’ll stop) – shouldn’t effector be affector? Effect is like a special effect, affect is to change something in a certain way.

    In your second block of code, I don’t know if it’s intentional, but you’re referencing systemMass as system.mass. Would it be less confusing (as should be the purpose of pseudocode) to reference it using the same name?

  19. brad says:

    @niall: Wow, thanks for reading so closely. I’m firing my proofreaders, are you looking for work?

  20. Niall says:

    Fantastic article. The math’s still a little over my head (terrible I know, but trig isn’t so familiar to the 13 year old mind), I’m trying to make sense of it though, heh.

    Just a couple of things I noticed:

    “decent into rigid-body hell” should be “descent into rigid body hell”
    And I’m not too sure, but on line 10 of your first block of code, shouldn’t “totalMass” be “systemMass”? You’re never updating systemMass so you’re dividing centerOfMass by 0 :P

  21. brad says:

    @Sven: Sorry for any misunderstanding there, I was just throwing down psuedo-code there for clarity’s sake, I wish there was operator overriding in AS3, that would be amazing. In my implementation of a vector class I made specific methods to be like an operator overload (plus, plusEquals, minus, minusEquals, etc. ), but I wanted to keep the code here free of my specific implementation for various reasons, mostly just to make the operations as clear as possible, as there seem to be a million flavors of vector class out there. Mine specifically has been slowly transforming as I use it for these experiments, and is actually a suite of 3 classes (2 utility ones), which makes me hesitant to ‘release’ them just yet. But I’m open for ultra-nerdy vector implementation debates any time. Cheerz.

    @Jeff: TY, compliments from the grave are always welcome. I’m sure the coming zombie-pocalypse will suite you well. Best.

  22. Sven says:

    Hi,

    while i must admit, i have my problems with Math, i am interested in this:

    var g:Vector2 = new Vector2 (…,…);
    g += currentCOM;

    What kind of class is Vector2? I wonder, how += works here.

    Thanks and cheers
    Sven

  23. [...] This post was mentioned on Twitter by karannnnnnnnnnn3. karannnnnnnnnnn3 said: Particle-based Rigid Bodies using Shape Matching: While hopping between projects the idea of a particle-based appr… http://bit.ly/9Pvg2j [...]

  24. Great post man! too bad my noggin is so small. :(


Leave a Reply to Brad

Click here to cancel reply.

*

Line
Line
Pony