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.