Labs.byHook 

Scripts, Tools & Methods Developed at Hook 
While hopping between projects the idea of a particlebased approach to rigidbody 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 rigidbodies 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 physicsbased 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 rigidbody 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 rigidbodies 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 rigidbody 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 (halfway down), which was the basis of my recent post dealing with more simplified ‘stick’ bodies.
The descent into rigidbody 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 2Dimensionalized 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 BandAids 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 ghostforces resulting in some really cool perpetual chaotic pendulums… just not quite the predictable system desired.
Last up was a 3particle 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 rigidbody 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 timestepped displacement (p, the mirror of q) from the relocated center of mass. In my implementation I am using Verlet style integration and running the code as a positionconstraint vs. modifying velocity via Euler integration. This choice was made to keep things in line with previous work, though the more forcebased approach is something worth exploring, especially for deformable bodies.
The following example shows the rigidbody 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 mousewheel) 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 timestep mode, the colored arrows are shaded green or red to indicate the direction of rotational influence of that individual particle.
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 psuedocode:
// 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 timestep, 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 rigidbody system particles[i].pos += (g  particles[i].pos) * k; }
Just plug that into the relaxation loop of the particle system and magic ensues:
The rotational values above are based on this rule of the dot product of two vectors a and b:
or:
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:
or:
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 framerate 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.
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
@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.
Hi Brad,
I’ve solved the problem in the mean time. Explanation and proofofconcept programs here:
http://www.gamedev.net/community/forums/topic.asp?topic_id=579582
More is coming up.
Cheers,
Mike
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 2degreeoffreedom 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 rigidbody 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
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
Hi Brad, how’s it going with your alternate approach?
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
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!
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
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
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 shapematching 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!
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 lightbulb 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 Jacobimagic 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 upness of the current system may be leading to a less stable system. Also some more standard box siming, 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
Wow, I would *love* it if you could explain how you arrived at your method for calculating the rotation of the shapematching 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 calculateframeorientation 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
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
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…
Ah yes, system.mass, not right, durned pasteover 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?
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.
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?
@niall: Wow, thanks for reading so closely. I’m firing my proofreaders, are you looking for work?
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 rigidbody 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
@Sven: Sorry for any misunderstanding there, I was just throwing down psuedocode 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 ultranerdy vector implementation debates any time. Cheerz.
@Jeff: TY, compliments from the grave are always welcome. I’m sure the coming zombiepocalypse will suite you well. Best.
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
[...] This post was mentioned on Twitter by karannnnnnnnnnn3. karannnnnnnnnnn3 said: Particlebased Rigid Bodies using Shape Matching: While hopping between projects the idea of a particlebased appr… http://bit.ly/9Pvg2j [...]
Great post man! too bad my noggin is so small.