margin of error trying to check if a point lays on line

Question

Working with floating point values it is as easy as breathing to run on approximation errors by comparing quantities which should be the same. I want to know if there is a way built in some MSDN (or even external) library for c# to ignore the problem.

An example could be: more than comparing 2 float values like this

if(myVector3.X == anotherVector3.X)

I would appreciate something like this

if(myVector3.X.isInTheNeighbourhood(anotherVector3.X))

This is not well-written, I know. That is just to simplify the explaination. What I am exactly doing is checking if a point (Vector3) lays on line segment. So basically the calculations I make are nothing more than

(x - x1)/(x2 - x1) = (y - y1)/(y2 - y1) = (z - z1)/(z2 - z1)

But these values won't be always the same, so I need to write down some code which includes a sort of tolerance a sort of Neighbourhood mathematical concept to accept values close enought to the line.

Hope that I made myself clear. Has anyone a solution for this?

Answer 1

I'm proposing the use of an exact predicate. This is an alternative approach to what you actually asked, but it might be worth considering.

Suppose your three points really are at locations indicated by their respective double precision coordinates. A simple double precision computation like the one suggested in your question will likely still return a wrong result. However, there are techniques to obtain exact results to detect these cases. One technique is by turning all numbers into arbitrary precision floating point (or integer) numbers and do most of the computation using integers under the hood. Another, which should work faster on modern hardware, expresses intermediate results as sums of doubles. E.g. when computing a+b, you obtain two resulting numbers, the first is the sum as you'd usually compute it, the other is a correction term to note down the error. In many cases, the bigger terms are sufficient to make a choice, which leads to the concept of adaptive precision.

All of this, including the application to geometric predicates, has been outlined nicely by Jonathan Richard Shewchuk in his page Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry. He has a paper on it, and some C code which should be possible to adapt to C#. Or perhaps to compile in C and link to C#, thus forming a mixed language project. Note however that it makes some assumptions on how the compiler treats floating point computations. In particular, you have to be careful that certain intermediate results don't have excess precision, like the 80-bit numbers stored in a 80387-style floating point unit. I've looked into this myself recently, and the easiest solution might be asking the compiler to use SSE instructions instead of x87 ones.

In your case, you are asking whether a point p lies on the segment connecting p1 and p2. One predicate which would be very much in the spirit of that paper would be the position of a fourth point in relation to a plane spanned by three others: it is either above, below or within that plane, and the orient3d predicate can tell you which. So one possible approach would be taking four arbitrary points q1 through q4 in general position, i.e. not all in a single plane, and no three on a single line. Then you could check whether p,p1,p2,qi are coplanar (sign is zero) for all i ∈ {1,2,3,4}. If they are, then p does at least lie on the line spanned by p1 and p2. Next you could check the orientations of p,p1,qi,qj for different i,j until you find a non-zero sign, then you could see whether p,p2,qi,qj has a different sign. If it does, then p is indeed between p1 and p2, hence on the same line. If you find no i,j such that p,p1,qi,qj is non-zero, then p=p1. Likewise if you have found one non-zero sign, but the corresponding sign of p,p2,qi,qj is zero, then p=p2. It is up to you whether you want to include the endpoints of your line segment. This paragraph might not be the most elegant way to do this, but it leverages the existing orient3d implementation, so it might be easier to use than writing a new predicate from scratch.

Note that in most cases, a point will not exactly lie on a line segment, due to rounding of the point coordinates themselves. The above will only help you to reliably detect those rare cases when it does, so you can deal with them. It may also allow you to make consistent choices for other predicates as well. In the world of CGAL, this approach would be termed “exact predicates, inexact constructions”, since the predicates can be computed reliably, but the geometric objects on which they operate are still subject to approximation. If you really need points that reliably lie on the line segments, and not just close by, then an exact construction approach using some exact number type would be preferable. Or you go with the approaches the other answers suggest.

Answer 2

You need to calculate the distance of the point to the line. That's simple mathematics. Then you decide how far of a distance "close enough" is in your case.

Answer 3

This problem is described as floating-point tolerance in this article and it points out the importance of measuring relative tolerance rather than absolute when comparing large values: http://realtimecollisiondetection.net/blog/?p=89

I've never had a situation where large floating point values are possible so I've always hard-coded a magic value into my comparisons. Eg:

if (Math.Abs(value1, value2) < 0.05f)
{
   ...
}

Which is "bad" but it works so long as value1 and value2 can't get too big.

In your case you really want to calculate the distance of the point from the line and then check that distance is small enough, rather than testing that the point is exactly on the line. I am rubbish at math myself so don't understand this but let me copy someone else's answer for calculating the shortest distance between a 3D point and a 3D line:

Vector3 u = new Vector3(x2 - x1, y2 - y1, z2 - z1);
Vector3 pq = new Vector3(x3 - x1, y3 - y1, z3 - z1);
float distance = Vector3.Cross(pq, u).Length / u.Length;

3D Perpendicular Point on Line From 3D point