Ageia paper on position based physics

raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

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.
raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

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 mouse-interaction 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
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Ageia paper on position based physics

Post by Antonio Martini »

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 mouse-interaction 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 think you should thank Kalman!:)

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 (1-expression) 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;)
raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

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 non-diagonal: 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.
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Ageia paper on position based physics

Post by Antonio Martini »

raigan2 wrote: isn't A*B diagonal for any A, so long as B is diagonal?
no, set B to identity A*I = A
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
raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

AntonioMartini 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
Argh, you're right! I'm a bit rusty :)

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 ;)
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Ageia paper on position based physics

Post by Antonio Martini »

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 ;)
yes you are correct J is a row vector if you have only one constraint.

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
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Ageia paper on position based physics

Post by Antonio Martini »

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 have just managed to run the code you posted, Q in fact seems to make a noticeable difference.

cheers,
Antonio
flyingwaffle
Posts: 3
Joined: Sun Nov 11, 2007 8:58 pm

Re: Ageia paper on position based physics

Post by flyingwaffle »

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?
raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

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.
flyingwaffle
Posts: 3
Joined: Sun Nov 11, 2007 8:58 pm

Re: Ageia paper on position based physics

Post by flyingwaffle »

ok, I guess you mean finding how to project
c(q,p1,p2,p3)=(q-p1).[(p2-p1)X(p3-p1)]
thanks
raigan2
Posts: 212
Joined: Sat Aug 19, 2006 11:52 pm

Re: Ageia paper on position based physics

Post by raigan2 »

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..
flyingwaffle
Posts: 3
Joined: Sun Nov 11, 2007 8:58 pm

Re: Ageia paper on position based physics

Post by flyingwaffle »

Yea, I hadn't realized at first that projecting using the Constraint gradient formula was an approximation (we're doing a newton-raphson 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 ) . (( p2-p1)X(p3-p1))
or rather its 2D version
C = (q-p1) . (1z X (p2-p1)) (1z coming out of the page)
Image

I get
grad_q (C) = 1z X (p2-p1) (q moves down in the direction of the normal of p1p2)
grad_p1 (C) = 1z X (q-p2) (p1 moves in the direction of the normal of qp2)
grad_p2 (C) = 1z X (p1-q) (p2 moves in the direction the normal of qp2)
the projections are (unless I made some error... I wish :P)
dq = - (q-p1).(1z X p12) (1z X p12)
-------------------------------
|p12|^2 + |qp1|^2 + |qp2|^2

dp2 = - (q-p1).(1z X p12) (1z X (qp1))
-------------------------------
|p12|^2 + |qp1|^2 + |qp2|^2

dp1 = - (q-p1).(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 ).

Image

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)
Image

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 p1-p and p2-p):

Image

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.
strangelet
Posts: 3
Joined: Sun Mar 25, 2012 12:54 pm

Conservation of linear and angular momentum

Post by strangelet »

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.
Post Reply