Ageia paper on position based physics
Regularized methods
Someone might be interested in looking at this:
http://www.vrlab.umu.se/research/flexible.shtml
And specifically at:
Interactive simulation of elastic deformable materials at
http://www.vrlab.umu.se/research/flexible/ISEDM.pdf
This implementation has not been optimized for speed, but mainly for evaluating the method, constitutive equations etc. The solver can be optimized way beyond the current implementation.
/Kenneth
http://www.vrlab.umu.se/research/flexible.shtml
And specifically at:
Interactive simulation of elastic deformable materials at
http://www.vrlab.umu.se/research/flexible/ISEDM.pdf
This implementation has not been optimized for speed, but mainly for evaluating the method, constitutive equations etc. The solver can be optimized way beyond the current implementation.
/Kenneth
Re: Regularized methods
I looked at it. As they note, it's really two parts: one is a more physically accurate way to model elastic tetrahedrons, the other is a new integration scheme.KenB wrote: And specifically at:
Interactive simulation of elastic deformable materials at
http://www.vrlab.umu.se/research/flexible/ISEDM.pdf
I'm still trying to fully understand exactly what they're doing in the new integration scheme and wondering if it can be applied successfully in a more general case, for example an iterative solver.
I'm working on one in 2D, however I've recently run into a problem: http://www.continuousphysics.com/Bullet ... php?t=1152jiangwei wrote:Is there any open source physics engine which used position based method?
thanks
The solver is excellent  superstable just like Jakobsen's, but much more applicable as you can extend it to rigid bodies, motors, etc..

 Posts: 2
 Joined: Tue Oct 02, 2007 12:54 am
Re: Ageia paper on position based physics
Hello
I'm interested in implementing this paper for a cloth simulation I'm working on. The problem is that the "solver" section is fairly vague. Well, it's not vague, but it implies knowledge about GaussSeidel which  even after a bit of research  I cannot find a clear example of its usage. I studied numerical methods a bit in college, but I'm not familiar with all the intimate details of course.
I'm confused as to whether or not I'm supposed to be creating a giant matrix and solving that, or if GaussSeidel is just a nice name for a relatively simple practice (with a strong mathematical concept "behind the scenes").
I'm mostly looking to implement the cloth solution from this paper. Any direction, code samples, etc. should be useful. If someone were to walk me through an iteration of a basic 2x2 piece of cloth, I'd be very thankful
So say I have edge constraints:
C1(p2, p3), C2(p3, p1), C3(p1, p4) C4(p4, p2)
Should I also have diagonal constraints? (e.g. C5(p1, p2) & C6(p4, p3)) Reading the ppaer, I believe I should have C5(p1, p2) but not the other diagonal since it's not a triangle edge. Am I reading that correctly? Why not have the other diagonal?
And I know I'll have one bend constraints Cbend(p1, p2, p3, p4) between two two triangles /\(p1, p3, p2), /\(p1, p2, p4).
So, I have those 57 constraints.. How do I solve this thing?
Am I reading too heavily into this? I have a feeling I am.
Am I supposed to just go from constraint C1 through CN and "apply" (project) the constraint to p1p4.
So C1 updates p2 & p3, C2 updates p3 & p1 (using the results from C1), and so on? In the case of cloth, this would be extremely similar to the Jacobson method except adding bend and collision constraints?
Is the real novelty here for cloth is that we are using predicted positions to modify the velocity ala the Nonconvex Rigid Bodies w/ Stacking paper (albeit in a nonimpulse way)?
I'm new to a lot of this, so any help or guidance would be appreciated. While I value the mathematical theory behind things, all I've been able to find is theory so some concreteness would help!
Thanks,
john
I'm interested in implementing this paper for a cloth simulation I'm working on. The problem is that the "solver" section is fairly vague. Well, it's not vague, but it implies knowledge about GaussSeidel which  even after a bit of research  I cannot find a clear example of its usage. I studied numerical methods a bit in college, but I'm not familiar with all the intimate details of course.
I'm confused as to whether or not I'm supposed to be creating a giant matrix and solving that, or if GaussSeidel is just a nice name for a relatively simple practice (with a strong mathematical concept "behind the scenes").
I'm mostly looking to implement the cloth solution from this paper. Any direction, code samples, etc. should be useful. If someone were to walk me through an iteration of a basic 2x2 piece of cloth, I'd be very thankful
Code: Select all
p2  p4
 \ 
 \ 
p3  p1
C1(p2, p3), C2(p3, p1), C3(p1, p4) C4(p4, p2)
Should I also have diagonal constraints? (e.g. C5(p1, p2) & C6(p4, p3)) Reading the ppaer, I believe I should have C5(p1, p2) but not the other diagonal since it's not a triangle edge. Am I reading that correctly? Why not have the other diagonal?
And I know I'll have one bend constraints Cbend(p1, p2, p3, p4) between two two triangles /\(p1, p3, p2), /\(p1, p2, p4).
So, I have those 57 constraints.. How do I solve this thing?
Am I reading too heavily into this? I have a feeling I am.
Am I supposed to just go from constraint C1 through CN and "apply" (project) the constraint to p1p4.
So C1 updates p2 & p3, C2 updates p3 & p1 (using the results from C1), and so on? In the case of cloth, this would be extremely similar to the Jacobson method except adding bend and collision constraints?
Is the real novelty here for cloth is that we are using predicted positions to modify the velocity ala the Nonconvex Rigid Bodies w/ Stacking paper (albeit in a nonimpulse way)?
I'm new to a lot of this, so any help or guidance would be appreciated. While I value the mathematical theory behind things, all I've been able to find is theory so some concreteness would help!
Thanks,
john
Re: Ageia paper on position based physics
You're correct  you just iterate over all constraints, projecting onebyone. The way that positions are updated (using "current" and "predicted") are in fact identical to Jakobsen ("oldpos" and "newpos") if you work it through on paper.JohnJensen wrote: Am I supposed to just go from constraint C1 through CN and "apply" (project) the constraint to p1p4.
So C1 updates p2 & p3, C2 updates p3 & p1 (using the results from C1), and so on? In the case of cloth, this would be extremely similar to the Jacobson method except adding bend and collision constraints?
Is the real novelty here for cloth is that we are using predicted positions to modify the velocity ala the Nonconvex Rigid Bodies w/ Stacking paper (albeit in a nonimpulse way)?
The main benefit of this paper is that it shows you how to properly extend Jakobsen's method to any sort of constraint (provided it's described as a function of particle positions). Personally I struggled for a long time with getting angle constraints to be stable, and the solver presented by Muller works excellently for this (in 2d at any rate).

 Posts: 2
 Joined: Tue Oct 02, 2007 12:54 am
Re: Ageia paper on position based physics
Thanks for your help.. It's getting clearer by the minute.
So is the "solver" the section where he described constraint projection that conserves momentum? E.g, is the conservation the "condition" he needs to derive a numerically stable solver?
So what's the "GaussSeideltype iteration"? Is it "merely" the process of repeatedly applying the constraints?
What about diagonal constraints, are there one or two in the sample I presented above?
Sorry I'm so questionyquestiony  I am experimenting with this at the moment, I just want to make sure I'm understanding things.
Thanks!
So is the "solver" the section where he described constraint projection that conserves momentum? E.g, is the conservation the "condition" he needs to derive a numerically stable solver?
So what's the "GaussSeideltype iteration"? Is it "merely" the process of repeatedly applying the constraints?
What about diagonal constraints, are there one or two in the sample I presented above?
Sorry I'm so questionyquestiony  I am experimenting with this at the moment, I just want to make sure I'm understanding things.
Thanks!

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
a couple of alternative methods for solving the problem in a more general way are mentioned here:raigan2 wrote:The main benefit of this paper is that it shows you how to properly extend Jakobsen's method to any sort of constraint (provided it's described as a function of particle positions). Personally I struggled for a long time with getting angle constraints to be stable, and the solver presented by Muller works excellently for this (in 2d at any rate).
http://www.bulletphysics.com/Bullet/php ... f=4&t=1593
cheers,
Antonio
Re: Ageia paper on position based physics
Thanks  i've been reading the thread (and papers) with great interest. Sadly I've only taken one year of linear algebra so I'm not entirely confident about my understanding..
I'm especially interested in the equivalent rowbyrow/SCAAT version (as opposed to the bigmatrix approach) since it seems like it would end up looking almost the same as the PBD solver except that it would handle masses/weighting differently; currently I'm using an invented scaling term since the revised paper's formula led to very slow convergence.
I've been meaning to post about my simulator, as increasingly the main problem has been collision detection (or, updating jacobians/attach points per iteration). My current solution is very similar to the Stewart/Trinkle or Anitescu/Potra solvers that have been posted here from time to time (as least, as far as I understand)  unfortunately this results in a large number of constraints being generated.
I'm especially interested in the equivalent rowbyrow/SCAAT version (as opposed to the bigmatrix approach) since it seems like it would end up looking almost the same as the PBD solver except that it would handle masses/weighting differently; currently I'm using an invented scaling term since the revised paper's formula led to very slow convergence.
I've been meaning to post about my simulator, as increasingly the main problem has been collision detection (or, updating jacobians/attach points per iteration). My current solution is very similar to the Stewart/Trinkle or Anitescu/Potra solvers that have been posted here from time to time (as least, as far as I understand)  unfortunately this results in a large number of constraints being generated.

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
you can try to parametrize the convergence rate by redefining the mass weights the following way.raigan2 wrote:Thanks  i've been reading the thread (and papers) with great interest. Sadly I've only taken one year of linear algebra so I'm not entirely confident about my understanding..
I'm especially interested in the equivalent rowbyrow/SCAAT version (as opposed to the bigmatrix approach) since it seems like it would end up looking almost the same as the PBD solver except that it would handle masses/weighting differently; currently I'm using an invented scaling term since the revised paper's formula led to very slow convergence.
w1 = (m2+q)/(m1+m2+q+q)
w2 = (m1+q)/(m1+m2+q+q)
for q = 0 you get the usual weighting
for q = large number you get approximately 0.5 for both weights
q can be a global parameter.
the idea is to start with a low q and increase it at each iteration, so you get the mass difference effect and a faster convergence. This is what the Q covariance matrix in the Kalman filter does. So referring to my other post(i have updated it slightly) this is the same as not having equation (e) and just updating the inverse mass at every frame by eq(b). By including equation (e) the convergence is even faster. The total error may still be the same, however the maximum local error should be smaller depending on the choice of Q. So we dont get localised overstretched edges. My tests until now have been purely at numerical level(comparing algorithms by running them on the same data) with no graphical output however the weighting above make sense also at intuitive level.
cheers,
Antonio
Last edited by Antonio Martini on Thu Oct 25, 2007 12:11 am, edited 1 time in total.
Re: Ageia paper on position based physics
That's a really interesting idea!
I was referring to the application of the weights  using them to scale the projections. My (possibly mistaken) understanding of this process is that the weights determine how far we descend along each particle's gradient, acting like the damping parameter in a dampedGaussNewton optimization/leastsquaresfitting process.
The individual particle weights are as in the paper: w_i = 1/m_i
We project each particle along its gradient by a certain amount: dp_i = step_i*dC_i
where dp_i is the projection to be applied to particle i, and dC_i is the gradient of the constraint function wrt the particle.
The original paper suggests this formula for the step length (apologies for uglyness):
step_i = ((n*w_i)/Sum_j(w_j)) * (C / Sum_j( Dot( dC_j, dC_j) ))
Where C is the current value of the constraint function, and n is the "cardinality" of the constraint (the number of particles). This formula doesn't work, it leads to overshooting/oscillation when the mass ratios aren't 1:1.
The revised PBD paper gives the formula as:
step_i = w_i * (C / Sum_j( w_j * Dot( dC_j, dC_j) ))
I've tried implementing this three times, and the results have always been the same: very slow convergence.
I tried using the following variation, however it didn't work (things exploded):
step_i = w_i * (C / Sum_j( Dot( w_j*dC_j, w_j*dC_j) ))
The formula that I'm currently using is the result of a mistake, when implementing the original paper I tried scaling the gradients _before_ calculating their squared length: let s_i = (n*w_i)/Sum_j(w_j)
step_i = s_i * (C / Sum_j( Dot( s_j*dC_j, s_j*dC_j) ))
This formula provides the best scaling so far; it never overshoots (even with mass ratios of 1:100), and provides very fast convergence.
I'm sure that there is a more optimal/correct scaling factor, since the use of the cardinality term is totally adhoc and shouldn't be necessary. For instance, you should be able to consider each component of the particle position seperately (i.e convert the constraint function C(p1,p1) to one which takes scalars C(p1x,p1y,p2x,p2y)) and the results should be identical, however since there are 4 scalars rather than 2 vectors, this changes the cardinality and the resulting projections.. that doesn't seem right to me.
I may be incorrect as far as construing the PBD solver to be doing a leastsquaresfit/optimization process, the formulas simply seem very similar to Newtontype optimization. I would _really_ like to understand why my formula works so much better, and what the optimal formula is  this is what interests me most about your Kalman filter method, it seems like a more correct method (no magical "cardinality" values floating around). I haven't wrapped my head around it yet, unfortunately "filter" seems very similar to "kernel" in that it seems to be only an abstract engineering tool and not something that can be easily visualized!
I was referring to the application of the weights  using them to scale the projections. My (possibly mistaken) understanding of this process is that the weights determine how far we descend along each particle's gradient, acting like the damping parameter in a dampedGaussNewton optimization/leastsquaresfitting process.
The individual particle weights are as in the paper: w_i = 1/m_i
We project each particle along its gradient by a certain amount: dp_i = step_i*dC_i
where dp_i is the projection to be applied to particle i, and dC_i is the gradient of the constraint function wrt the particle.
The original paper suggests this formula for the step length (apologies for uglyness):
step_i = ((n*w_i)/Sum_j(w_j)) * (C / Sum_j( Dot( dC_j, dC_j) ))
Where C is the current value of the constraint function, and n is the "cardinality" of the constraint (the number of particles). This formula doesn't work, it leads to overshooting/oscillation when the mass ratios aren't 1:1.
The revised PBD paper gives the formula as:
step_i = w_i * (C / Sum_j( w_j * Dot( dC_j, dC_j) ))
I've tried implementing this three times, and the results have always been the same: very slow convergence.
I tried using the following variation, however it didn't work (things exploded):
step_i = w_i * (C / Sum_j( Dot( w_j*dC_j, w_j*dC_j) ))
The formula that I'm currently using is the result of a mistake, when implementing the original paper I tried scaling the gradients _before_ calculating their squared length: let s_i = (n*w_i)/Sum_j(w_j)
step_i = s_i * (C / Sum_j( Dot( s_j*dC_j, s_j*dC_j) ))
This formula provides the best scaling so far; it never overshoots (even with mass ratios of 1:100), and provides very fast convergence.
I'm sure that there is a more optimal/correct scaling factor, since the use of the cardinality term is totally adhoc and shouldn't be necessary. For instance, you should be able to consider each component of the particle position seperately (i.e convert the constraint function C(p1,p1) to one which takes scalars C(p1x,p1y,p2x,p2y)) and the results should be identical, however since there are 4 scalars rather than 2 vectors, this changes the cardinality and the resulting projections.. that doesn't seem right to me.
I may be incorrect as far as construing the PBD solver to be doing a leastsquaresfit/optimization process, the formulas simply seem very similar to Newtontype optimization. I would _really_ like to understand why my formula works so much better, and what the optimal formula is  this is what interests me most about your Kalman filter method, it seems like a more correct method (no magical "cardinality" values floating around). I haven't wrapped my head around it yet, unfortunately "filter" seems very similar to "kernel" in that it seems to be only an abstract engineering tool and not something that can be easily visualized!

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
equations (ae) of my post dont require any filter notion, you can just plug in the values and and even keep them in matrix form. If you get an improvement then they can be written in a much simpler form as many terms typically cancel out.raigan2 wrote: I may be incorrect as far as construing the PBD solver to be doing a leastsquaresfit/optimization process, the formulas simply seem very similar to Newtontype optimization. I would _really_ like to understand why my formula works so much better, and what the optimal formula is  this is what interests me most about your Kalman filter method, it seems like a more correct method (no magical "cardinality" values floating around). I haven't wrapped my head around it yet, unfortunately "filter" seems very similar to "kernel" in that it seems to be only an abstract engineering tool and not something that can be easily visualized!
for my tests i usually use newmat:
http://www.robertnz.net/nm_intro.htm
it's a very simple to use and flexible library. However it's terribly slow, when im happy with the results then i simplify the computations and move on other far more efficient libraries.
if you want to try the EKF approach just send me a private message so we can get it to work.
cheers,
Antonio
Last edited by Antonio Martini on Wed Oct 24, 2007 11:48 pm, edited 3 times in total.
Re: Ageia paper on position based physics
Thanks very much, I'm planning on working through this tomorrow, hopefully I'll be able to manage. I feel quite stupid compared to all of the experts on this board, I'm definitely more of a "technical artist" than a proper programmer in that my math background is weak.AntonioMartini wrote: equations (ae) of my post dont require any filter notion, you can just plug in the values and and even keep them in matrix form. If you get an improvement then they can be written in a much simpler form as many terms typically cancel out.
I asked Mr. Muller about this issue, and he said the revised version is working for him, and that my formula might ismply be introducing a successiveoverrelaxationtype of bias. However, I've received a private message from someone who was also having slow convergence with the revised formulas.Dirk Gregorius wrote:Strange that you have these issues. The formulas seem correct. I derived the formulas myself and in the new version it correct now. The derivation is extremely simple (using block vector notation):
The revised formula may very well be correct; the issue may be with my expectation of how it should perform  my formula might simply be more aggressively descending the gradient, and it turns out that for the types of constraints I'm using this bias never results in problems (overshooting/etc), thus I get accelerated convergence. Currently, I get very solid behaviour with a low number of iterations (4 iterations for this biped: http://www.metanetsoftware.com/robotolo ... ktest.html); you can violently swing it around and the bodies stay together  even two iterations is enough to prevent obvious visual gaps.
However with the revised formulas, 4 iterations feels very weak/loose  if you move one of the bodies across the screen, after the solve the remaining bodies will be spread out across the screen.. huge gaps. At around 1020 iterations it becomes solid  perhaps this is normal behaviour, I know box2D uses 10s of iterations. So perhaps the only problem is my high expectations of convergence; nevertheless, the fact that my current method behaves much better is something I'd like to understand! Also, since with the PBD approach you have to rebuild the gradients/jacobians at each iteration, keeping the iteration count low is desirable.
There is also the possibility that my code is flawed; once I have time to return to this I'll clean up and post some of the code.
Re: Ageia paper on position based physics
Working through ad (skipping e, and using Q=0, R=0) I get eq9 in the PBD paper; this is the same formula as Dirk's most recent post.
If I understand things correctly, my adhoc formula is equivalent to using a nonzero value for R; I'm going to try a sanity check implementation of a simple chain of jointed rigid bodies.
If I understand things correctly, my adhoc formula is equivalent to using a nonzero value for R; I'm going to try a sanity check implementation of a simple chain of jointed rigid bodies.

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
Re: Ageia paper on position based physics
if they are equivalent you should be able get almost identical numerical results. If thats the case you may want to test the effect of Q and (e) on the final result.raigan2 wrote:Working through ad (skipping e, and using Q=0, R=0) I get eq9 in the PBD paper; this is the same formula as Dirk's most recent post.
If I understand things correctly, my adhoc formula is equivalent to using a nonzero value for R; I'm going to try a sanity check implementation of a simple chain of jointed rigid bodies.
cheers,
Antonio