col. prediction at 2nd order

Started by
17 comments, last by Eelco 19 years, 6 months ago
I have studied the boring underlying maths required to know in which conditions quadratic approximation of motions would be enough to deal with the prediction of contact points in time and space. The result tends to prove that physics based on predictive algos can be efficient and robust. It's obvious that any given metric precision of quadratic approximations can be achieved with enough rate of update (physics updates frequency). It's also intuitive that it depends on the maximum accelerations, angular speeds and radii of the objects. Here is the formula I have obtained. I spare you the details, unless one requests them. This is the metric error in world coordinates for a solid translating (no acceleration) and rotating at angular speed w and with a max radius R around its gravity center : d = R*((w*t)^3)/6! For R=1m, w=2PI/s, d=1mm, t=1/40 s. This means that a physics update of 40 FPS is enough for 1mm precision, 400FPS for 1 micro meter precision, that is below the floating point precision (relative to meters). When considering two objects rotating errors cumulate. The same if one considers the movement of one in the motion frame of the other since metrics are conserved. Thus if both objects have roughly the same R and w params we need 80FPS for 1mm precision and 800FPS for floating point precision or more. It's also easy to cumulate the metric error due to the accelerations relative to the previous problem. Since err = 1/2 a*t*t. For t=1/40 s, err = a/3200, if a=10m/s*s (1G), err=3mm. For typical real life accelerations, this means this error is slightly preponderant. For 400FPS, we have err=30 micro meters. Now here is a how quadratics could be used for quasi immediate contact point resolution. Imagine we're in the context of 2 convexes, an algo uses the separating plane theorem. At t0=0 (start of the current frame), the currently cached plane is one support plane Ha of object A, the plane of face Fa. Thus I deal with Face/Vert case for simplicity. We work in the frame of A. We use the polynomial matrix (M0,M1,M2) quadratics. Say, the nearest vertex Vb of object B has 4 adjacent edges thus vertex neighbours. We seek the first vertex hitting the plane, the time being t. That's very easy and root finding inside a [0,1] (by scaling time) interval (thus optimizable) three outcomes are possible : a) t> tmax, the end of the time interval of research (maybe the end of the current frame. Then no collision (this frame). b) The first vertex that crosses the plane is not the nearest vertex at t0. c) The nearest vertex is not in the Voronoi cell of Fa. d) else the collision is found. In cases c and d a new maximum separation axis must be found. t will advance very quickly and d will be < epsilon very quickly if a contact exists. I have dug the details, and the algo will not stall. In my opinion a contact is found in 1 step mainly, and 3 steps at max, very rarely. No contact will be missed.
"Coding math tricks in asm is more fun than Java"
Advertisement
interesting stuff. makes sense, accuracy with t^(order+1)

an inituative and at first glance good maximum interval size for swepth spheres is 1/(wb-wa), ie the period of their relative rotation speed, possibly multiplied by a constant to tweak accuracy. however, your intergrationstep is going to be smaller than this value for about all cases anyway.

accuracy indeed doesnt seem a problem then.

i understand how to apply all this to spheres: evaluate position and speed at tmin, and position at tmax, and you can use these three constraints to obtain a quadratic (on the interval 0-1 for simplicity indeed).

then:
if at both times relative position < 0, collision along the whole interval.
if one relative position < 0, there is one root: use the quadratic to find it
if both are positive (most likely), get the determinant to see if there are roots, and find them if so.


however, how would you use a quadratic to sweep edge-edge or vertex-plane? youd need a way to evaluate their positions at a time t (ok easy) but also one derivative. the formula for the distance between a rotating plane and rotating point is probably a bitch, and its derivative even worse. the same goes for edges. id be interested to know what you think is a good way to handle this.
why was i not logged in?
@Eelco

interesting stuff. makes sense, accuracy with t^(order+1)

Yes, t^3 indeed for 2nd order, and I think your formula also works at any degree. This comes from the third order term of sinus, cosinus does not have one. The theory of alternate limited series is the heart of this mathematic demonstration. At the 3rd order it would be O(t^4/4!), a coefficient that comes from cosinus.

This all means that you really have to compute the appropriate sampling frequency, adapt to each collision pair or to the whole context if your speeds and radii can be bounded.


... an inituative and at first glance good maximum interval size for swepth spheres is 1/(wb-wa) ...


- wb-wa only works if the axes of rotation are the same.
- In case each rotations are around the center of mass=center of sphere, then wa and wb do not even count except for the collision response. Else, right the sampling frequency it linear to the absolute (I forgot to mention) rotation frequency.



However, your intergrationstep is going to be smaller than this value for about all cases anyway.

Yes I think 400FPS for the physics can be achieved on any modern PC, specially if collision prediction accelerates things, which is the even more case for small time steps. 400 FPS means don't worry anywhere, we're beyond floating point precision. That's nearly magic, we're 'exact'.

... and you can use these three constraints to obtain a quadratic (on the interval 0-1 for simplicity indeed).


Yes rescaling the polynomial is a few SSE operations. Next it will be ez to cull much root computations very early as you said. I have written the code in the past. If you detail the computation more than you did, make quick study of the params and quadratic properties, you'll find how to exploit the [0,1] constraint first and very efficiently. You'll discover analytically a condition that means reject early if the spheres go in opposite directions. Then a second condition (for t=1) with a bit more computations. Computing the determinant will be in rare cases. And getting the root only once you are sure one contact exists in the interval.



However, how would you use a quadratic to sweep edge-edge or vertex-plane ?


I mentionned it in () and you probably failed to notice or did not see what it was. The collision prediction needs what I call a 2nd degree polynomial matrix for each mesh, that is three matrices. That will be the only additionnal information required between the physics engine and the collision prediction module.

These matrices are easilly derived from the body state informations at time t=0. During the canonical physics frame (time rescaled), you'll have for any point attached to the solid, P(t)/world = M(t)*P/local. The inaccuracy is described by the previous post. P(t) is thus a quadratic vector <=> position at t=0, velocity and double acceleration.

Face vs Vert
------------
the formula for the distance between a rotating plane and rotating point is probably a bitch

Now here is how you can do. Say we have cached a separating plane that is the support of a face (of A), that we have to check vs a vert (of B).

Compute M(t)b/a = tMa(t)*Mb(t). Keep only up to the second order coeffs for the reasons I mentionned earlier : "cumulating errors". Thus 1+2+3 = 6 matrix muls (use SIMD). Maybe quaternions and translations would help, or any other quicker method. But anyway it has to be done only once per candidate mesh-pair during a frame, thus not a prioritary optimization.

For edge/edge :
---------------
In fact we usually cache a separating plane statically attached to A (or B, whatever) in the frame of A. Then it's the same as previous. When a plane is invalidated, a new one is searched on adjacent d-facets, until we get close to a contact in space and time. This should cull most cases in one or two steps.

If the case is not culled, a true edge/edge must be done. There a contact is quasi certain (99%). A special computation is required. One edge is fixed since we're in the frame of A. It's a bit trickier, since the nearest points are moving on both meshes. Thus the root finding is of higher degree in theory.

There are special conditions to be sure a contact will happen in the Voronoi cell of the edge (thus remove roots on the infinite line that would not belong to the finite edge). This is tested quickly. If it fails a new nearest d-facet pair (thus separating plane) must be searched, at a given time (static algo, usually one step required). Then, when a contact is certain, a root extraction must be done using a Taylor algo.

Roughly it's a coplanar test between 3 vectors so it results in a 3x3 determinant and a degree 4 polynomial : (AB^AC(t))*AD(t). But this can be optimized with the useful time interval constraint. I think it can be optimized and made degree 2 or 3 root finding.
"Coding math tricks in asm is more fun than Java"
Yes, quadratic spline approximates cicrle pretty well.. it's kinda related to that spherical mirrors can be used in telescope instead of parabolic ones, even with required accuracy of light wavelength /8

Truing to derive edge-edge formula:
let edge 1 have ends A(t),B(t) and edge2 have ends A'(t),B'(t)

Collision happen when
( (A(t)-B(t) )Cross( A'(t)-B'(t) ))Dot( A(t)-A'(t) )=0
(geometrically-derived in mind)

i'm using notation
A(t).x = t2*A.x2+t*A.x1+A.x0
and same for y and z.
so A is 3 quadratic polynomials for xyz, with "overloaded operator()"

Let
P=A-B
Q=A'-B'
R=A-A'

So let's rewrite equation component-by-component:
(P(t).y*Q(t).z-P(t).z*Q(t).y)*R(t).x + [same blabla for other components,y and z] = 0
P(t).y * Q(t).z * R(t).x - P(t).z * Q(t).y * R(t).x + [same blabla for other components] = 0

There [same blabla for other components] can be written by substituting z->x, y->z , x->y

It's looks like we will have 6th degree polynomial, not 4th , and also it looks like it will be very long to write.... Probably i started with bad collision equation.:(

6th degree polynomials can be solved numerically pretty well and it's not that slow.

Notation: "." mean the same as in C++, like if all variables is a classes/structs, and polynomials also have overloaded operator (t).

It's looks like we will have 6th degree polynomial, not 4th , and also it looks like it will be very long to write.... Probably i started with bad collision equation.:(

No it's OK. But in the case I described, that is compute in the frame of A, P=A-B is constant. So your own formula gives 4th degree too. I have explained why I can transform B into the frame of A and keep its 2nd order motion equations valid (somehow clamping the coefficients of higher degree).

Metrics approximations just cumulate in each transfo. So as long as the preconditions concerning wa, wb are true, the results concerning precision hold true.
"Coding math tricks in asm is more fun than Java"
Yes,you're right, i haven't noticed we can do things in one system.

So, as i understand, you interpolate matrix that transforms from one coordinate system to other? That's good idea!
strange. I was logged in before typing reply(i have avatars turned off so i always notice if i'm not logged in)
Quote:Original post by Charles B

... an inituative and at first glance good maximum interval size for swepth spheres is 1/(wb-wa) ...


- wb-wa only works if the axes of rotation are the same.

my formula was meant for the 2D case. in 3D finding a good metric for the number of subdivisions of your interval is probably more difficult, but for 2D, the above will work fine, since its a worst case formula.

Quote:
- In case each rotations are around the center of mass=center of sphere, then wa and wb do not even count except for the collision response. Else, right the sampling frequency it linear to the absolute (I forgot to mention) rotation frequency.

im not sure what youre saying here... the rotation velocity certainly is important for contact determination. if its several radians/frame, only one quadratic approximation certainly isnt going to suffice.

oh wait i think were talking about different things here.. i was thinking along the lines of as little physics updates as possible, ie as much as the framerate or something, and get the desired accuracy by subdividing the timeinterval where needed, whereas you seem talking about getting the desired accuracy by making more calls to UpdatePhysics.

the former is faster i think, and can handle unlimited rotationspeeds without loss of accuracy, by taking only smaller timesteps where its actually needed. it does require finding a good criterium on which to base your choice of timestep though. but i think its better than brute force, ie just making a huge number of calls to the physics update, and still introducing a limit on rotationspeed.

Quote:

However, how would you use a quadratic to sweep edge-edge or vertex-plane ?


I mentionned it in () and you probably failed to notice or did not see what it was. The collision prediction needs what I call a 2nd degree polynomial matrix for each mesh, that is three matrices. That will be the only additionnal information required between the physics engine and the collision prediction module.

These matrices are easilly derived from the body state informations at time t=0. During the canonical physics frame (time rescaled), you'll have for any point attached to the solid, P(t)/world = M(t)*P/local. The inaccuracy is described by the previous post. P(t) is thus a quadratic vector <=> position at t=0, velocity and double acceleration.

Face vs Vert
------------
the formula for the distance between a rotating plane and rotating point is probably a bitch

Now here is how you can do. Say we have cached a separating plane that is the support of a face (of A), that we have to check vs a vert (of B).

Compute M(t)b/a = tMa(t)*Mb(t). Keep only up to the second order coeffs for the reasons I mentionned earlier : "cumulating errors". Thus 1+2+3 = 6 matrix muls (use SIMD). Maybe quaternions and translations would help, or any other quicker method. But anyway it has to be done only once per candidate mesh-pair during a frame, thus not a prioritary optimization.

For edge/edge :
---------------
In fact we usually cache a separating plane statically attached to A (or B, whatever) in the frame of A. Then it's the same as previous. When a plane is invalidated, a new one is searched on adjacent d-facets, until we get close to a contact in space and time. This should cull most cases in one or two steps.

If the case is not culled, a true edge/edge must be done. There a contact is quasi certain (99%). A special computation is required. One edge is fixed since we're in the frame of A. It's a bit trickier, since the nearest points are moving on both meshes. Thus the root finding is of higher degree in theory.

There are special conditions to be sure a contact will happen in the Voronoi cell of the edge (thus remove roots on the infinite line that would not belong to the finite edge). This is tested quickly. If it fails a new nearest d-facet pair (thus separating plane) must be searched, at a given time (static algo, usually one step required). Then, when a contact is certain, a root extraction must be done using a Taylor algo.

Roughly it's a coplanar test between 3 vectors so it results in a 3x3 determinant and a degree 4 polynomial : (AB^AC(t))*AD(t). But this can be optimized with the useful time interval constraint. I think it can be optimized and made degree 2 or 3 root finding.


hmm youve kinda lost me here... however, i only today realized you dont really need to bother with derivatives of the distance function: simply calcing the distance on THREE timeintervals (does require one more state evaluation, quat->matrix (ouch) though) instead of using two and a derivative is much easier and just as accurate. then find roots and make sure that for each root found were inside the voronoi cell. quite easy actually. your method does sound faster though, since there is more potential for early rejection.
EDIT : something went wrong here too with the login.

Quote:
- In case each rotations are around the center of mass=center...
I'm not sure what youre saying here...
oh wait i think were talking about different things here...

I mentionned the case of a single swept sphere, not part of a tree, centered at it's mass center. Then due to spherical symetry collision prediction here would not care of the rotations. But that did not make much sense to point this trivial case out I admit. I hope you do not feel insulted ;)


Quote:
... by making more calls to UpdatePhysics.


Yes, it is actually required if you want to have a constant accuracy whatever the angular speed..

But it does not have to be a full integration necessarilly. To be more precise one can still subdivide the physics in two functionalities :
- The dynamics integration (may be 50FPS). With updated forces, forces that might depend on the first order body state, example : damping.
- The kinetics integration that supposes both linear accelerations and angular velocities (or accelerations) are constant during 1/fq.

What counts for me is the kinetic part. I want to be able to retrieve State(Body)(t) enough times. Depending on the implementation, using quaternions and slerps maybe, it should handle cycloid vertex movements accurately. (Imagine a dart thrown. You want to know which part hits the wall, on the blade end or not). I suppose these computations do not approximate sinus and cosinus. Else it's the problem of the integration, not the prediction.

I mentionned a context where dynamic parameters (accelerations and angular velocities) are bounded as an example. Then the physics fq (both dynamics and kinetics all at once) can be set once for all for a given precision, at compile time for instance. The prediction also works mainly at the same fq. The interval of search is the same as the physics frame except the tmax optimization I have mentionned in the swept sphere thread.

I also said that in the general case the (kinetic) fq could be locally adaptable. Think of typical engines, that manage multiple constraints and collisions they also request more physics updates from time to time. (And usually they do it the whole scene which causes huge stalls). Also consider that when you deal with contacts, you also have to call the physics too. There is nothing new here. Still using prediction instead of detection warrants that 400FPS is almost always sufficient for 24bits of precision ! (40FPS for 10bits). Plus in each interval, operations are done in quasi-linear time.

So let's take a context where the upper limitation of w is not known. A vertex path can then be a portion of cycloid during your 10 milliseconds (say fq=100Hz).

If you keep you physics (integration) frequency constant and only subdivide the prediction part. Then the quadratic approximation of the kinetics between two integration steps might be too unprecise for two given solids. As I explained it, with rotations, the metric error grows in t cube ! So if you don't call back more integration steps and retrieve 'fresh' body states, whatever the number of subdivisions you'll make in your collision prediction (your stance if I understood you well), this error exists. And it may well not be acceptable.

It's basically the conclusion of the math I have studied in details in the first post.

Quote:
... simply calcing the distance on THREE timeintervals (does require one more state evaluation, quat->matrix (ouch) though) ...


It this does not change the required physics fq for a given context and precision tolerance. So I'll compare the two solutions in the same interval.

This is a possibility, it's just a matter of implementation and interface between the physics and the collision prediction modules. But it also heavilly impacts on the underlying algorithms :

Here is what you do :
 M0             Mhalf           M1-|--------------|---------------|--> t 0              0.5             1


- 3 matrices required.
- You do not need to retrieve full body states. Just 3 frames cached, giving you any point at t=0, t=0.5 and t=1. You need a small computation to get the quadratic (easy). The problem to me is it requires 3 instants of physics (or just kinetic) evaluation.
- When you request P(t) you need 3 matrix*vert plus the construction of the vertex quadratic.

 M0                               +M1*t+M2*t*t-|------------------------------|--> t 0                              1

- 3 matrices required
- For the same time step, no body state is required in advance, which is certainly more convenient. A quadratic matrix just requires one body state at t=0.
- When you want P(t) you need three matrix*vert
- I can also get V(t) with two matrix*vert.
- Crucial : I can get the path of Pb(t) relative to the kinetic frame of the other solid (A). This is done with 3 matrix multiplications for the whole interval and once. Cf my solution for moving faces against moving vertices.

Both solutions are just equivalent in theory but the last crucial point. Yours gives the polynomial by 'control points' and mine by its coefficients more directly. It's mainly a question of convenience and implementation. It depends on your engine architecture. But think of the face-vert case again.


Last word :
On the largest scale, to me the goal of a physics engine mostly based on prediction is to have a computational complexity linear to the number of objects (inputs) AND linear to the discrete events it has to treat (outputs) : the changes of contact states, impulsions, start/ends of constraints, triggered events (ex : enter teleport). The goal is thus algorithmic optimality.

[Edited by - Charles B on September 22, 2004 7:36:27 AM]
"Coding math tricks in asm is more fun than Java"

This topic is closed to new replies.

Advertisement