Position correction tutorial
Position correction tutorial
I read the position correction introduction : http://www.gphysics.com/archives/35
They really help a lot. But I don't know what exactly the algorithms are.
Could you please help me provide some more detailed materials introducing these algorithms?
Thanks a lot!
They really help a lot. But I don't know what exactly the algorithms are.
Could you please help me provide some more detailed materials introducing these algorithms?
Thanks a lot!
Re: Position correction tutorial
I have sort of wondered the same thing, as I tried to implement the Full NGS position correction but apparently failed.
If you write out the whole velocity correction equation first, you've got:
Vcorrection = M^1 * J^T * (J * M^1 * J^T)^1 * J * V,
where:
M^1 is the inverted mass matrix
J is the Jacobian (and J^T is the same transposed), and
V is the original velocities
I might be missing a negation since I usually have to subtract the Vcorrection from the original velocities.
Anyway, I *think* you're supposed to plug in some sort of position errors somewhere, and you'll get position corrections out the other side. I could not get this to work, however. I tried using the original velocities plus position errors, and sometimes the result got me moderately closer to the correct answer, but often it just sent me into the weeds.
So I'm pretty sure I was not doing it the correct way. Perhaps I was confused because there's no proper way to come up with an orientation error in 3 components to match up with the 3 components of the rotation. Or perhaps you have to measure the position error in the specific directions of the impulse.
I ended up correcting position errors another way that was a bit more intuitive to me. I would measure the effective mass in the direction of the error, then I applied some sort of instantaneous position correction (which I call teleportation) in the same direction. Even when I had to do two separate corrections for the position and orientation (for joints like hinges), it typically corrected 95+% of the error in one pass. And it was probably faster than recomputing a Jacobian for each position correction pass.
However, if someone can explain how to use position errors in an equation using the Jacobians, I'd love to hear it and try it out.
If you write out the whole velocity correction equation first, you've got:
Vcorrection = M^1 * J^T * (J * M^1 * J^T)^1 * J * V,
where:
M^1 is the inverted mass matrix
J is the Jacobian (and J^T is the same transposed), and
V is the original velocities
I might be missing a negation since I usually have to subtract the Vcorrection from the original velocities.
Anyway, I *think* you're supposed to plug in some sort of position errors somewhere, and you'll get position corrections out the other side. I could not get this to work, however. I tried using the original velocities plus position errors, and sometimes the result got me moderately closer to the correct answer, but often it just sent me into the weeds.
So I'm pretty sure I was not doing it the correct way. Perhaps I was confused because there's no proper way to come up with an orientation error in 3 components to match up with the 3 components of the rotation. Or perhaps you have to measure the position error in the specific directions of the impulse.
I ended up correcting position errors another way that was a bit more intuitive to me. I would measure the effective mass in the direction of the error, then I applied some sort of instantaneous position correction (which I call teleportation) in the same direction. Even when I had to do two separate corrections for the position and orientation (for joints like hinges), it typically corrected 95+% of the error in one pass. And it was probably faster than recomputing a Jacobian for each position correction pass.
However, if someone can explain how to use position errors in an equation using the Jacobians, I'd love to hear it and try it out.

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
Re: Position correction tutorial
Yeah, J has to be the Jacobian of your position constraint.
We want to solve:
M * (v2  v1) = J^T * lambda
J * v2 + C = 0
Where C is the current position error.
We set v1 = 0 because there is no momentum in position correction.
Solve:
v2 = inv(M) * J^T * lambda
Substitute:
J * inv(M) * J^T * lambda + C = 0
Define:
K = J * inv(M) * J^T
Solve:
lambda = inv(K) * C
Pseudo Integrate:
x2 = x1 + v2
q2 = q1 + 0.5 * quat(omega2) * q1 (pseudo quaternion integration)
q2 = normalize(q2)
We have to recompute J, M, and C after each position update.
We want to solve:
M * (v2  v1) = J^T * lambda
J * v2 + C = 0
Where C is the current position error.
We set v1 = 0 because there is no momentum in position correction.
Solve:
v2 = inv(M) * J^T * lambda
Substitute:
J * inv(M) * J^T * lambda + C = 0
Define:
K = J * inv(M) * J^T
Solve:
lambda = inv(K) * C
Pseudo Integrate:
x2 = x1 + v2
q2 = q1 + 0.5 * quat(omega2) * q1 (pseudo quaternion integration)
q2 = normalize(q2)
We have to recompute J, M, and C after each position update.

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Position correction tutorial
You can also look at the Muller paper about position based dynamics. It has a pretty good explanation as well. You only need to transpose the gradient dC/dx and treat it as a matrix. Note that the Jacobian basically has the gradients in its rows and that you can write the dot product also as matrix multiplication. Basically: a * b = transpose( a ) * b. This is a useful trick. Matrix multiplication is associative. E.g. you need this when building the Jacobian of a springdamper.
Example: (a * b) * c != a * (b * c) where a, b and c are vectors and * denotes the dot product
This other gotcha is the pseudo integration of the angular velocity. If you feel unsure with this just use a function that integrates the angular velocity and pass a timestep of dt = 1.0.
The PhD of Cline explains the concept of postprojection as well. I think Erwin has it somewhere on his FTP server. Maybe he can add/link it to his reference material.
HTH,
Dirk
Example: (a * b) * c != a * (b * c) where a, b and c are vectors and * denotes the dot product
This other gotcha is the pseudo integration of the angular velocity. If you feel unsure with this just use a function that integrates the angular velocity and pass a timestep of dt = 1.0.
The PhD of Cline explains the concept of postprojection as well. I think Erwin has it somewhere on his FTP server. Maybe he can add/link it to his reference material.
HTH,
Dirk
Re: Position correction tutorial
Is the paper "PostStabilization for Rigid Body Simulation with Contact and Constraints"?
Re: Position correction tutorial
Hi Erin:
those your method mentioned above a NGS method or a Pseudo Velocities?
those your method mentioned above a NGS method or a Pseudo Velocities?
We want to solve:
M * (v2  v1) = J^T * lambda
J * v2 + C = 0
Where C is the current position error.
We set v1 = 0 because there is no momentum in position correction.
Solve:
v2 = inv(M) * J^T * lambda
Substitute:
J * inv(M) * J^T * lambda + C = 0
Define:
K = J * inv(M) * J^T
Solve:
lambda = inv(K) * C

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Position correction tutorial
It depends on how you see it. Basically it is NGS. The gotcha to understand is that when you solve for the position error directly (and not in terms of a velocity bias) the Jacobians and inertia change after each solved constraint. The pseudo velocities keep the Jacobians from the beginning of the time frame, but don't add the bias to the velocity but bypass them. This way it doesn't add to the momentum. Position correction with NGS is similar, but now you need to solve a nonlinear system using an outer Newton loop and an inner GaussSeidel loop. These two loops are usually merged into one.
Re: Position correction tutorial
Thanks Dirk Gregorius:
So in pseudo velocities method,the Solver is very similar like the PGS Solver when we solver constraints J * inv(M) * J^T * lambda + C = 0, the J * inv(M) * J^T and C remains const in all the iterators,right?
But as you mentioned in NGS Method,
while(iteratornum < totalnum)
{
for(every constraint row)
{
Recalculate Jacobian ;
Recalculate PositionError;
solve the i'th lambda;
update the corresponding body's position
}
}
Could you give more information about the
Thanks a lot!
So in pseudo velocities method,the Solver is very similar like the PGS Solver when we solver constraints J * inv(M) * J^T * lambda + C = 0, the J * inv(M) * J^T and C remains const in all the iterators,right?
But as you mentioned in NGS Method,
,So is that mean when we are solving the i'th row of the constraint ,the corresponding Jacobian Row in J * inv(M) * J^T Matrix need to be Recaculate due to current mass and position? And whether the position error C need to be recalculate too?the Jacobians and inertia change after each solved constraint
I don't quite understand this, what is the code like when using outer Newton loop and an inner GaussSeidel loop and merge them:will the solver like this:but now you need to solve a nonlinear system using an outer Newton loop and an inner GaussSeidel loop. These two loops are usually merged into one
while(iteratornum < totalnum)
{
for(every constraint row)
{
Recalculate Jacobian ;
Recalculate PositionError;
solve the i'th lambda;
update the corresponding body's position
}
}
Could you give more information about the
or some paper or reference that I can read and understand it more clearlyouter Newton loop and an inner GaussSeidel loop
Thanks a lot!

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Position correction tutorial
What you have posted is correct. You need to recalculate the Jacobian, inertia and position error in each constraint!
Regarding the inner and outer loop stuff. Imagine solving a nonlinear system of equations using Newton method (see Generalizations here: http://en.wikipedia.org/wiki/Newton%27s_method). You basically solve:
for ( n Newton iterations)
{
J*W*JT * lambda = C
}
J*W*JT * lambda = C is a linearized system of equations using Taylor expansion. So you need to rebuild it for each iteration (Jacobian, inertia and position error change). Now you can solve this inner linearized system using an iterative method as well > Gauss Seidel. So you get:
for ( n Newton iterations)
{
for ( m GaussSeidel iterations )
{
J*W*JT * lambda = C
}
}
Now you can kind of "merge" these two iterative loops which is basically done in the way we already discussed. You rebuild the Jacobian, inertia and position error every constraint. So we can say that the outer loop takes care of the nonlinear problem and the inner loop solves the linearized system. In the method suggested by Erin (which is similar to Jacobsen Cloth paper for example) you iteratively solve the constraints as usual, but you take care of the nonlinearity by constantly updating the Jacobian, inertia and position error.
HTH,
Dirk
Regarding the inner and outer loop stuff. Imagine solving a nonlinear system of equations using Newton method (see Generalizations here: http://en.wikipedia.org/wiki/Newton%27s_method). You basically solve:
for ( n Newton iterations)
{
J*W*JT * lambda = C
}
J*W*JT * lambda = C is a linearized system of equations using Taylor expansion. So you need to rebuild it for each iteration (Jacobian, inertia and position error change). Now you can solve this inner linearized system using an iterative method as well > Gauss Seidel. So you get:
for ( n Newton iterations)
{
for ( m GaussSeidel iterations )
{
J*W*JT * lambda = C
}
}
Now you can kind of "merge" these two iterative loops which is basically done in the way we already discussed. You rebuild the Jacobian, inertia and position error every constraint. So we can say that the outer loop takes care of the nonlinear problem and the inner loop solves the linearized system. In the method suggested by Erin (which is similar to Jacobsen Cloth paper for example) you iteratively solve the constraints as usual, but you take care of the nonlinearity by constantly updating the Jacobian, inertia and position error.
HTH,
Dirk

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Position correction tutorial
Regarding a reference the "Position Based Dynamics" paper by Matthias Mueller is pretty decent. Note that the gradients dC/dx are the rows of the Jacobians. So you only need to transpose the gradients and you end in the matrix form we use here.
Re: Position correction tutorial
Thanks Dirk:
I will take a close look for the Newton method first:)
I will take a close look for the Newton method first:)
Re: Position correction tutorial
Hi Dirk:
Sorry for disturb you again,but I occasionnaly read the topic
http://www.bulletphysics.com/Bullet/php ... ?f=4&t=878
which is also talk about the position correction and I am now a little confused because in that topic Jan blender also mentioned his linear position correction method ,which you replyed that "I don't see a big difference between your approximation and the Newton solver"
I am confused
1:) if Jan blender 's method has the same effect as the Newton Solver,why we need Newton Solver?Because blender 's method only need to rebuild the Jacobian once,but the Newton Method Need to Rebuild Jacobian every time when one constraint is solved.
2:)I am confused about the difference between the Blender's method and the split impluse method ,it seems they are the same expect the calculation of the position error(Blender use integret to recalculate position error and split impluse use linear method?)
3:)I want to know the paper that introduce the Blender's method ,Could you tell me which paper of Blender contain the method
Thanks
Sorry for disturb you again,but I occasionnaly read the topic
http://www.bulletphysics.com/Bullet/php ... ?f=4&t=878
which is also talk about the position correction and I am now a little confused because in that topic Jan blender also mentioned his linear position correction method ,which you replyed that "I don't see a big difference between your approximation and the Newton solver"
I am confused
1:) if Jan blender 's method has the same effect as the Newton Solver,why we need Newton Solver?Because blender 's method only need to rebuild the Jacobian once,but the Newton Method Need to Rebuild Jacobian every time when one constraint is solved.
2:)I am confused about the difference between the Blender's method and the split impluse method ,it seems they are the same expect the calculation of the position error(Blender use integret to recalculate position error and split impluse use linear method?)
3:)I want to know the paper that introduce the Blender's method ,Could you tell me which paper of Blender contain the method
Thanks

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Position correction tutorial
You are absolutely correct about the Jacobian. Bender's method only uses the Jacobian's from the beginning of the time step while a correct Newton solver rebuilds them for every constraint. Think of it practically. If the position error is small then you can get away with this. In games you use incredible big time steps (3060 Hz), you have large angular velocities (which break linearizations) and you use an iterative solver (which means you have only an approximation of the correct solution). All these factors introduce large position errors into system and now you have a root finding solver like Newton which needs to converge against the global solution. To my knowledge Bender used much smaller timesteps (100200 Hz) and much more iterations (untill the error falls below some epsilon threshold). This let's him get away with his approximation. My personal opinion is that this method is not suited for realtime applications like games.
What I meant in my other post was that my intuition tells me that Benders method and the NGS are "kind of " equivalent, but Bender's method is more restricted since he doesn't update the positions directly and therefore cannot update the Jacobians as easy as the NGS. Of course he could use candidate positions, but you see hos hacky this becomes.
Maybe a mixture of Jan Benders and the original NGS method would make sense:
 Update positions
 Correct position
 Update velocities
 Correct velocities
These is a little better than the Jacobsen/Mueller method since you now project the velocities onto the velocity constraint manifold instead of using the difference between the old and the new position. On the downside you don't have symplectic Euler anymore (see Erin's latest GDC presentation on integrators about this topic if you are interested).
I recommend reading all papers by Jan since one can learn a lot from them.
Dirk
What I meant in my other post was that my intuition tells me that Benders method and the NGS are "kind of " equivalent, but Bender's method is more restricted since he doesn't update the positions directly and therefore cannot update the Jacobians as easy as the NGS. Of course he could use candidate positions, but you see hos hacky this becomes.
Maybe a mixture of Jan Benders and the original NGS method would make sense:
 Update positions
 Correct position
 Update velocities
 Correct velocities
These is a little better than the Jacobsen/Mueller method since you now project the velocities onto the velocity constraint manifold instead of using the difference between the old and the new position. On the downside you don't have symplectic Euler anymore (see Erin's latest GDC presentation on integrators about this topic if you are interested).
I recommend reading all papers by Jan since one can learn a lot from them.
Dirk
Re: Position correction tutorial
Thanks a lot,I will read these paper you recommanded