Epsilons in floating point math
cbloom
7-19-04
Whenever you do floating point math, you generate inaccurate results
(actually that's not strictly true, but I'm not going to get into all
the details of numerical accuracy here; there are other articles on that
if you care). For example, whenever you add floating point numbers, some
error might be introduced. If you later compare that sum to the known
exact sum, there will be inaccuracies. The big problem arises when you
want to generate a boolean query on some innaccurate math. A standard
and good solution is to use a fuzzy "epsilon" range which changes your
boolean result to a trinary result - true, false, and maybe. There are
two things that people often do wrong - the definition of epsilon, and
the handling of the "maybe" case.
First of all, many people think that an "epsilon" fuzzy query is just a
fudge factor or hack to cover problems. That is not true if done
correctly. In fact, floating point math without epsilons is generally
just not correct and not robust. With epsilons, your algorithm can be
completely robust and well defined and consistent. It may not be exactly
"correct" in the sense that it may produce slightly different results
than if you had done exact arithmetic, but it is perfectly well-defined
and self-consistent.
Now, for concreteness I'm going to talk about a very specific problem
for the rest of the article. This stuff applies to all types of FPU
math, but there's a very common case in game development, which is
collision detection and response. I'm going to talk about a moving point
particle, which may collide with some triangle soup.
The basic way the movement is done is like this :
At time T the particle has some position P, velocity V, and acceleration A
We wish to evolve time by a delta "d"
Compute the desired end-point, Q = P + Vd+Ad*d/2
This is done with floating point math, so Q is innacurate, but that doesn't matter
Now we try to do the move P->Q using our collision test
Our collision test returns a boolean - whether you hit something or not.
If you did hit something, it also give you the point you can move up to,
and the surface normal of what you hit. We'll have an addendum about how
to respond to the collision. The collision test here is a line segment
vs. triangle test.
The collision test must be "epsilon robust", in the sense that it must
never allow the particle to penetrate solid space by more than a
distance of Epsilon. Note that we are now talking about a very specific
concrete value Epsilon here - and it has *units* of distance. You can
never combine values that have different units, so our value Epsilon
should only be used in formulas involving distance (doing that sort of
thing is what made that mars explorer miss its target). There are two
primary queries that will stress the epsilon tolerance - a test where
the line segment query goes through an edge of the triangle, and a test
where the start point is nearly on the face of the triangle.
We only gaurantee correct results on a closed mesh, so to handle the
first test, we just need to make sure that a line segment that goes
through an edge is either counted as being hitting one triangle face or
on the other (two triangle faces are incident on an edge). Whether or
not the edge actually has two faces is irrelevant. Also, the test should
not return true for line segments that hit the plane of the triangle
farther than Epsilon from the edge.
The second case (the point starts nearly on the face) comes up often.
The primary way it arises is after a collision, the particle is moved up
to being in contact with the triangle face. Note that "in contact" here
is a fuzzy term defined by Epsilons. Once it is moved to being "in
contact", you will try to move again from that point. These cases must
be handled correctly such that the particle cannot slip inside the body.
Now, there are many ways to handle these things, but I'll describe one
specifically.
We are going to return more collisions than exact math would. That is,
any ray which actually collides (in exact math) should produce a
collision. Any ray which is farther than epsilon away from colliding
should report no collision. Any ray that's within Epsilon may or may not
report a collision.
First, we consider every triangle to have an Epsilon (which is a
distance) region around it. Triangles have a front and back, and a point
is classified as either on the "front", "back", or "on" the triangle
face, which means inside the epsilon region. Our collision is
directional, that is - segments only colliding when moving into faces.
Any point which starts in the "on" region of the face is not allowed to
get any more embedded - it must move strictly *away* from the face.
Points that start in the "back" region are considered not colliding -
they can move in or out of the face. Points that start in the "front"
region may move anywhere away from the face. Let me make this precise
with pseudo-code.
N = world space normal of the triangle
T0,T1,T2 = coordinates of the triangle (world space)
segment end-points are P and Q
D(P) is the signed distance of point P to the plane of the triangle
D(P) = P*N - T0*N
if D(P) > Epsilon (E) then P is on front
else if D(P) < -E then P is behind
else P is "on" the triangle
if Q is on front -> no collision
if Q is "on" the face ->
if P is "on" the face -> allow it if D(Q) >= D(P) , otherwise return collision at P
if P is in front or behind the face -> allow the move, no collision
if Q is behind ->
if P is behind -> allow the move
if P is "on" the face -> return collision at P
if P is in "front" -> do the triangle test (more details later)
One thing we notice is that if you have polygonal features on the order
of the size of Epsilon, you will get very wrong results. This is the
limit that the Epsilon tolerance puts on your world. The exact value of
Epsilon depends on the number of bits of precision in your math, what
math you do exactly, and how big your world is.
The actual triangle test is not that interesting. We use the
Moller-Haines ray test with epsilons; this is well described various
places, including on my site. For concreteness, the code is -
static bool intersect_triangle(P,Q,T0,T1,T2)
{
/* begin calculating determinant - also used to calculate U parameter */
Vec3 dir = Q - P;
Vec3 pvec = dir ^ (T2-T0); // ^ means cross product
/* if determinant is near zero, ray lies in plane of triangle */
float det = fabsf( (T1-T0) * pvec ); // * = dot product
/* calculate distance from vert0 to ray origin */
Vec3 tvec = P - T0;
/* calculate U parameter and test bounds */
float u = tvec * pvec; // = Vec3U::TripleProduct(dir,edge2,tvec);
if ( u < -INTERSECT_EPSILON || u > det + INTERSECT_EPSILON)
return false;
/* prepare to test V parameter */
/* calculate V parameter and test bounds */
float v = dir * (tvec ^ (T1-T0));
if ( v < -INTERSECT_EPSILON || (u + v) > det + INTERSECT_EPSILON)
return false;
return true;
}
intersect_triangle returns true or false for whether or not a collision
occurred. Note that here we're using INTERSECT_EPSILON as a tolerance
rather than our Epsilon. This is quite a different beast -
INTERSECT_EPSILON has units of distance cubed. It is a strict expansion
- our only goal with it is to make sure that we return true for any ray
that would have collided (or grazed an edge) in exact arithmetic. It
must also be small enough that it doesn't pick up rays farther than
Epsilon away. We could be more precise with this, but we don't want to
give up speed, and this is good enough (see appendix for figuring out an
appropriate value for INTERSECT_EPSILON).
Now, if intersect_triangle returns true, we have a collision, and we
must return a hit point. We do a simple plane-distance computation :
Hit = D(P) * Q - D(Q) * P / ( D(P) - D(Q) )
note that D(P) > Epsilon and D(Q) < -Epsilon. For our system to be
consistent and robust, this hit point must be strictly in the Epsilon
region for the face. That is :
D(Hit) = Hit*N - T0*N
require abs(D(Hit)) <= Epsilon. This kind of self-consistency check is
very valuable, because you can assert() with it to verify your system is
working right.
That's all there is to it! You now have a robust consistent collision
detection system. We'll now go over some of the details that we skipped
over.
Appendix 1 : Collision Response
Once you get a collision, you must respond to it. This doesn't really
need to be done very carefully, because it won't lead to major bugs if
you get it wrong, but we'll describe one reasonably way to do it. First,
we go ahead and move our particle up to the hit position, Hit. Then, we
take the previous velocity and clamp it out of the contact normals. I'm
not going to go into that here because it's rather complex to get it
right with multiple contact surfaces, etc. Remember that "contact" here
is defined as being in the Epsilon region. One thing we do is create a
velocity which is definitely *away* from the face. We clamp the velocity
out of the normal, apply "friction" to reduce the magnitude of the
velocity, and then give an extra little boost of delta away from the
surface to make sure we are moving away. Note that this boost is NOT
needed for the system to be robust, and we aren't ever moving the
particle without collision checking! This is NOT a lift-off the surface
normal, that would cause bugs. We are just changing the velocity to make
sure we make a velocity that will be allowed to move, so that we can hit
the surface and slide smoothly without sticking.
Appendix 2 : Coordinate Transformations
One common problem is caused by coordinate space issues. Your collision
geometry might be in some object space, and your particle might be in
world space. To do the collision, you need to transform the query into
local space, do the collision, then transform it back. The main thing
this does is just introduce a bit more inaccuracy into the tests,
because the coordinates get distorted a bit by the transforms. Also, if
the transform has any scale or shear you get additional problems. The
Epsilon is defined in world space, so if you have scale in the
transform, and then do the work in local space, you may be using Epsilon
in the wrong coordinate system. One way to handle this is to define a
limit on the range of scales allowed. For example, you might require the
scale to be in the range 1/5 - 5. Then, you have to choose an Epsilon
and constrain your world such that any tolerance in the range Epsilon/5
- 5*Epsilon works.
Appendix 3 : How big my Epsilon ?
Epsilon must be big enough so that the consistency condition abs(D(Hit))
<= Epsilon is satisfied for all queries. You want the smallest epsilon
you can have that meets that constraint. Normal floats have a 24 bit
mantissa. When you do any float addition, you may have round-off error
that makes the result off by roughly 2^-24 * result. Our result is
scaled by the position values. If our world is strictly required to be
in a box of world size W (each coordinate in -W to W), then the maximum
error is 2^-24 * W. Thus Epsilon must be at least >= 2^-24 * W. If
you're doing coordinate transforms, that may scale your error up by some
amount, so you'll need a bigger epsilon. In general something like
2^-22*W is reasonable. If you allow scaled transforms, it needs to be
something like 2^-22*W*MAX_SCALE.
Appendix 4 : Constraints on your geometry
Your geometry must fit within the world-size box (-W to W). It must also
not have any features smaller than Epsilon. For example, if you have two
faces with opposite facing closer than Epsilon to each other, you can
get points stuck there, since a point between those two faces is "on"
both, and won't be allowed to move in either direction. Sliver triangles
are problematic mainly because they make the computation of the triangle
normal difficult, and they create a large dependence on the choice of
the "base" vertex T0. For example, consider a triangle {A,B,C} where the
points B and C are very close together, and both are very far from A.
This makes (B-C) much more accurate than (A-B) and (A-C). Now if you try
to compute a triangle normal using (A-B) X (A-C) , you will get a very
innacurate normal; if you use (A-B) x (B-C), it will be very different.
This mainly affects the barycentric computation and INTERSECT_EPSILON.
Exact study of the affect of slivers is beyond the scope of this
article.