Ageia paper on position based physics
Re: Ageia paper on position based physics
Yeah, I mostly want to try my formula alongside the PBD one, and then see what R/Q/e do. The difference in convergence speed especially is something I want to finally prove for myself, so far it's possible that it's just the result of an improper implementation of the PBD formulas.
Re: Ageia paper on position based physics
Here's a quick test; for some reason (e) makes things unstable (adds energy), quite possibly an implementation error. The PBD formula is definitely better than the one I've been using, and using Q the Kalman method converges even better.
I still need to add some mouseinteraction and figure out what's going wrong with (e).
Anyway, I'll upload my code, it's actionscript3. Thanks very much Antonio for inventing this! I'm glad I'll be able to get rid of the "cardinality" formula from my code.
I still need to add some mouseinteraction and figure out what's going wrong with (e).
Anyway, I'll upload my code, it's actionscript3. Thanks very much Antonio for inventing this! I'm glad I'll be able to get rid of the "cardinality" formula from my code.
 Attachments

 kalmanTest.zip
 (47.7 KiB) Downloaded 680 times

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
i think you should thank Kalman!:)raigan2 wrote:Here's a quick test; for some reason (e) makes things unstable (adds energy), quite possibly an implementation error. The PBD formula is definitely better than the one I've been using, and using Q the Kalman method converges even better.
I still need to add some mouseinteraction and figure out what's going wrong with (e).
Anyway, I'll upload my code, it's actionscript3. Thanks very much Antonio for inventing this! I'm glad I'll be able to get rid of the "cardinality" formula from my code.
anyway in my tests, (e) reduced the error even further, so there is probably an error somewhere. there is a numerically more stable version of (e) you may want to try:
e) W' = (I  K*J)*W*(I  K*J)_t
however it's much more likely that there is a bug somwhere, my test wasn't general enough, or convergence depends on Q and R.
as i have almost nothing installed on my laptop, what do i need to install in order to view your code?
cheers,
Antonio
PS: i had a quick glance at the code, i am not comfortable with the language used however after equation (e) W' is not a diagonal matrix anymore so it should be stored as a full matrix. is that the case at the moment? as (1expression) is evaluated inside the loop, either you are computing only the diagonal of W' or the I(identity matrix) has been replaced with a matrix(bug) having 1 in every entry.Either ways seem wrong;)
Re: Ageia paper on position based physics
My mistake  I thought W would stay as a diagonal matrix.
I may very well have made other mistakes, when I worked through it on paper it seemed that all the terms were either scalars (R, K*J, J*W*J_t, C), column/row matrices (X,J,K, W*J_t, ) or diagonal (I,W,Q,S). Some of the column matrices were columns of scalars, others columns of vectors.. i suppose this might be where I screwed up  I just assumed that all of the matrix operations would reduce to dotproduct or scaling. To be honest it was a bit confusing what should be a row and what a column..
Also, it's not clear how W can become nondiagonal: isn't A*B diagonal for any A, so long as B is diagonal? And since W is diagonal to start with..
You'll need FlashDevelop (http://www.flashdevelop.org/) to compile the project.
I may very well have made other mistakes, when I worked through it on paper it seemed that all the terms were either scalars (R, K*J, J*W*J_t, C), column/row matrices (X,J,K, W*J_t, ) or diagonal (I,W,Q,S). Some of the column matrices were columns of scalars, others columns of vectors.. i suppose this might be where I screwed up  I just assumed that all of the matrix operations would reduce to dotproduct or scaling. To be honest it was a bit confusing what should be a row and what a column..
Also, it's not clear how W can become nondiagonal: isn't A*B diagonal for any A, so long as B is diagonal? And since W is diagonal to start with..
You'll need FlashDevelop (http://www.flashdevelop.org/) to compile the project.

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
no, set B to identity A*I = Araigan2 wrote: isn't A*B diagonal for any A, so long as B is diagonal?
a*B is diagonal when a is a scalar and B is diagonal.
i think you tried to simplify it a bit before getting it to work, i just used a generic mxn matrix library for my tests as i wanted to exclude any other source of error.
if by setting Q and R = 0 and without using (e) you get a result equivalent to the PBD paper formulas(either on paper or numerically) then everything is fine, you can just concentrate on fixing (e).
cheers,
Antonio
Re: Ageia paper on position based physics
Argh, you're right! I'm a bit rustyAntonioMartini wrote: i think you tried to simplify it a bit before getting it to work, i just used a generic mxn matrix library for my tests as i wanted to exclude any other source of error.
if by setting Q and R = 0 and without using (e) you get a result equivalent to the PBD paper formulas(either on paper or numerically) then everything is fine, you can just concentrate on fixing (e).
cheers,
Antonio
I did in fact get identical results for PBD and Q=R=0 Kalman, that's promising  my main problem is that i have no matrix library.. it seemed like the matrices were all sparse, so I just implemented them using arrays of scalars or vec2s.
I think I'm a bit confused about row and column vectors; K is definitely a column vector, as is X, but is J a row vector? If so, this would mean that K*J is an outer product, which would explain my problems as I assumed it was an inner product. Initially I just assumed everything what whatever was necessary to turn all of the multiplications into inner products

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
yes you are correct J is a row vector if you have only one constraint.raigan2 wrote:AntonioMartini wrote:I think I'm a bit confused about row and column vectors; K is definitely a column vector, as is X, but is J a row vector? If so, this would mean that K*J is an outer product, which would explain my problems as I assumed it was an inner product. Initially I just assumed everything what whatever was necessary to turn all of the multiplications into inner products
in general the number of rows of J is the same as the number of contraints, the number of columns is the same as the number of degrees of freedom, 6 for a edge constraint or 2 row vectors of size 3 if written in block form. I think i have written the Jacobian for the edge constraint in the other thread.
cheers,
Antonio

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
i have just managed to run the code you posted, Q in fact seems to make a noticeable difference.raigan2 wrote:Here's a quick test; for some reason (e) makes things unstable (adds energy), quite possibly an implementation error. The PBD formula is definitely better than the one I've been using, and using Q the Kalman method converges even better.
cheers,
Antonio

 Posts: 3
 Joined: Sun Nov 11, 2007 8:58 pm
Re: Ageia paper on position based physics
Hello, I'm new to these boards...
After spending a couple of years implementing a physics engine based on the Barraff papers (it took me a really long time to get rid of numerical instabilities), I'm trying to implement the position based solution presented in this paper and Jakobsen's original paper.
I'm mostly interested in applying it to rigid bodies.
About that, Jakobsen mentions
"Here, we have collided a single rigid body with an immovable world. The above method generalizes to handle collisions of several rigid bodies. The collisions are processed for one pair of bodies at a time. Instead of moving only p, in this case both p and q are moved towards each other."
but he doesn't explain by how much should each rigid bodies should be moved to correct the constraint (somewhere between p and q).
For example if a node (at p) of body_1 penetrates a triangle of body_2 (p1, p2, p3), at point q, we should correct p by moving it along the normal of p1_2_3, toward q, and we should move p2, p3, p4 as well (from q toward p, correcting each node p1, p2, p3 proportionally). But where between p and q should be final correction?
Is this something that I can automatically derive with the constraint gradient technique?
After spending a couple of years implementing a physics engine based on the Barraff papers (it took me a really long time to get rid of numerical instabilities), I'm trying to implement the position based solution presented in this paper and Jakobsen's original paper.
I'm mostly interested in applying it to rigid bodies.
About that, Jakobsen mentions
"Here, we have collided a single rigid body with an immovable world. The above method generalizes to handle collisions of several rigid bodies. The collisions are processed for one pair of bodies at a time. Instead of moving only p, in this case both p and q are moved towards each other."
but he doesn't explain by how much should each rigid bodies should be moved to correct the constraint (somewhere between p and q).
For example if a node (at p) of body_1 penetrates a triangle of body_2 (p1, p2, p3), at point q, we should correct p by moving it along the normal of p1_2_3, toward q, and we should move p2, p3, p4 as well (from q toward p, correcting each node p1, p2, p3 proportionally). But where between p and q should be final correction?
Is this something that I can automatically derive with the constraint gradient technique?
Re: Ageia paper on position based physics
If you read the "Position Based Dynamics", they actually present the solution for collision between a particle and an edge/triangle, that should do exactly what you want.

 Posts: 3
 Joined: Sun Nov 11, 2007 8:58 pm
Re: Ageia paper on position based physics
ok, I guess you mean finding how to project
c(q,p1,p2,p3)=(qp1).[(p2p1)X(p3p1)]
thanks
c(q,p1,p2,p3)=(qp1).[(p2p1)X(p3p1)]
thanks
Re: Ageia paper on position based physics
If you need the actual position where the two contact points will be, the only way to get this (as far as I can tell) is to use the projection formula to project the point. I may be missing something however..

 Posts: 3
 Joined: Sun Nov 11, 2007 8:58 pm
Re: Ageia paper on position based physics
Yea, I hadn't realized at first that projecting using the Constraint gradient formula was an approximation (we're doing a newtonraphson step)
The main assumption is
C(p+Dp) ~ C(p) + Grad(C(p)).Dp = 0
so
applying
Dp_i =  C(p) Grad_i(C(p))/[Sum_i Grad_i(C)^2]
*once* isn't enough to satisfy the constraint (by an epsilon) if the penetration is significant (although in the case of the length constraint it works, it's all linear).
Is that correct?
I worked with the vertex/triangle collision constraint
C = ( q  p1 ) . (( p2p1)X(p3p1))
or rather its 2D version
C = (qp1) . (1z X (p2p1)) (1z coming out of the page)
I get
grad_q (C) = 1z X (p2p1) (q moves down in the direction of the normal of p1p2)
grad_p1 (C) = 1z X (qp2) (p1 moves in the direction of the normal of qp2)
grad_p2 (C) = 1z X (p1q) (p2 moves in the direction the normal of qp2)
the projections are (unless I made some error... I wish )
dq =  (qp1).(1z X p12) (1z X p12)

p12^2 + qp1^2 + qp2^2
dp2 =  (qp1).(1z X p12) (1z X (qp1))

p12^2 + qp1^2 + qp2^2
dp1 =  (qp1).(1z X p12) (1z X (p2q))

p12^2 + qp1^2 + qp2^2
if q is close to the segment p1p2, q moves along the normal of p1p2 and p1 & p2 almost move along the normal as well, but a bit towards q (1z x p2q ~ 1z X p1p2 ).
But if q, p1 and p2 form a equilateral triangle with side 1 (admittedly a very extreme case of collision).
each vertex moves toward the center of the triangle, by the same distance (sqrt(3)/6 along the normal)
And projecting repeatedly such a collision constraint doesn't help (I guess you'd have to apply it once, then project the length constraints, project the collision constraint again, etc...).
My point with all this is that it's quite different from what Jakobsen is doing, where he moves the vertices all in the same direction qp (with some weighing based on geometry (ratio p1p and p2p):
in the case q and p1p2 all belong to rigid non static bodies, he just says to move them along some point between q and p, but not where that point lies.
The main assumption is
C(p+Dp) ~ C(p) + Grad(C(p)).Dp = 0
so
applying
Dp_i =  C(p) Grad_i(C(p))/[Sum_i Grad_i(C)^2]
*once* isn't enough to satisfy the constraint (by an epsilon) if the penetration is significant (although in the case of the length constraint it works, it's all linear).
Is that correct?
I worked with the vertex/triangle collision constraint
C = ( q  p1 ) . (( p2p1)X(p3p1))
or rather its 2D version
C = (qp1) . (1z X (p2p1)) (1z coming out of the page)
I get
grad_q (C) = 1z X (p2p1) (q moves down in the direction of the normal of p1p2)
grad_p1 (C) = 1z X (qp2) (p1 moves in the direction of the normal of qp2)
grad_p2 (C) = 1z X (p1q) (p2 moves in the direction the normal of qp2)
the projections are (unless I made some error... I wish )
dq =  (qp1).(1z X p12) (1z X p12)

p12^2 + qp1^2 + qp2^2
dp2 =  (qp1).(1z X p12) (1z X (qp1))

p12^2 + qp1^2 + qp2^2
dp1 =  (qp1).(1z X p12) (1z X (p2q))

p12^2 + qp1^2 + qp2^2
if q is close to the segment p1p2, q moves along the normal of p1p2 and p1 & p2 almost move along the normal as well, but a bit towards q (1z x p2q ~ 1z X p1p2 ).
But if q, p1 and p2 form a equilateral triangle with side 1 (admittedly a very extreme case of collision).
each vertex moves toward the center of the triangle, by the same distance (sqrt(3)/6 along the normal)
And projecting repeatedly such a collision constraint doesn't help (I guess you'd have to apply it once, then project the length constraints, project the collision constraint again, etc...).
My point with all this is that it's quite different from what Jakobsen is doing, where he moves the vertices all in the same direction qp (with some weighing based on geometry (ratio p1p and p2p):
in the case q and p1p2 all belong to rigid non static bodies, he just says to move them along some point between q and p, but not where that point lies.

 Posts: 3
 Joined: Sun Mar 25, 2012 12:54 pm
Conservation of linear and angular momentum
Article from Müller about position based dynamics says that correcting/moving particles in direction of the gradient of the constraint function conserves linear and angular momentum. Could pls someone help me explain that? I don't see it.