1

I wrote a code testing for the intersection of two line segments in the plane. I won't bother you with all the details.

The code takes two line segments, each described by two end-points, then fits each segment to a line by fitting the a and b in y = a*x + b. Then it finds the intersection point of the two lines by x = (b2 - b1) / (a2 - a1). Last, it tests if the intersection point x is contained within the two line segments.

The relevant part looks like this:

# line parameterization by a = Delta y / Delta x, b = y - a*x
a1 = (line1.edge2.y - line1.edge1.y) / (line1.edge2.x - line1.edge1.x)
b1 = line1.edge1.y - a1 * line1.edge1.x
a2 = (line2.edge2.y - line2.edge1.y) / (line2.edge2.x - line2.edge1.x)
b2 = line2.edge1.y - a2 * line2.edge1.x
# The intersection's x
x = - (b2 - b1) / (a2 - a1)
# If the intersection x is within the interval of each segment
# then there is an intersection
if (isininterval(x, line1.edge1.x, line1.edge2.x) and
    isininterval(x, line2.edge1.x, line2.edge2.x)):
    return True
else:
    return False

For brevity I dropped a lot of tests handling specific cases like when the edges are parallel to each other (a1==a2), when they are on the same line, when an edge is of length 0, when the edge is along the vertical axis (then a becomes infinite) etc.

The function isininterval is simply

def isininterval(x0, x1, x2):
    """Tests if x0 is in the interval x1 to x2"""
    if x1 <= x0 <= x2 or x2 <= x0 <= x1:
        return True
    else:
        return False

Now the question: I find that due to roundoff errors, the test will give false results when the intersection coincides with the segment edge.

For example, with line1 between (0,0) and (3,5) and line2 between (3,5) and (7,1) the resulting intersection x is 2.9999999999999996, which gives the false answer. Should have been 3.

Can you please suggest a solution?

Jonathan Leffler
  • 730,956
  • 141
  • 904
  • 1,278
Aguy
  • 7,851
  • 5
  • 31
  • 58
  • 1
    Use the decimal class: `from decimal import Decimal`. http://stackoverflow.com/questions/2986150/python-floating-number . Also your `isininterval` function could simply return: `return x1 <= x0 <= x2 or x2 <= x0 <= x1` – Bahrom Jul 04 '16 at 14:21
  • The standard way to write Boolean conditions involving floating point numbers is to include some error tolerance in the condition. For example, replace `x1 <= x0 <= x1` by `x1 - eps <= x0 <= x2+ eps` where `eps` is something like 0.000001 – John Coleman Jul 04 '16 at 14:27
  • In the first place, do you have a good reason to worry about 2.9999999999999996 not being 3 ? –  Jul 04 '16 at 14:30
  • @YvesDaoust, well, my only worry is that it fails the condition... – Aguy Jul 04 '16 at 14:47
  • Are there good reasons for vertices of one segment to fall exactly on the other segment ? Are the segments isolated or do they form chains or graphs ? –  Jul 04 '16 at 14:51
  • @YvesDaost - They can be e.g. edges of a closed polygon. The test can ask which other edges intersect with a specific edge. I guess I could do a special case for adjacent edges, but I wanted something more elegant... – Aguy Jul 04 '16 at 15:09

1 Answers1

5

This is a problem/feature of floating point arithmetic. There are ways to minimise the error, by ordering instructions certain ways, but in the end, you will get approximate answers, because you're representing possibly infinite numbers with a finite number of bits.

You need to define whatever functions you build in such a way that they are tolerant of these errors. Looking at your example, the difference between the "proper" value and what you got is of the order 1e-16 - extremely extremely low.

With inequality and especially equality, it's worth it to relax the constraints for exact/bitewise matching. For example, if you wanted to test that x == 3, you would write this as abs(x - 3) < EPSILON, where EPSILON = 1e-6 or EPSILON = 1e-9. Basically, the difference between what you want to have and what you have is smaller than a very small value. Ditto, for inequality, you might test that 3 - EPSILON <= x or x <= 3 + EPSILON.

Horia Coman
  • 8,681
  • 2
  • 23
  • 25
  • Is there a machine defined epsilon? Is it universal for all machines? – Aguy Jul 04 '16 at 15:12
  • 1
    There is sys.float_info.epsilon from the sys module, but that might be a little bit on the lower side, because it is defined, technically as the difference between 1 and the smallest number larger than 1 which is represented as a float. So there aren't any floats between [1, 1+eps] more or less. However, this is quite application specific, and you should definitely feel free to experiment. – Horia Coman Jul 04 '16 at 15:15
  • Actually, you can provably mitigate this "problem/feature" by writing some *really* careful code. I believe the technique is due to Jonathan Shewchuk. See http://www.cs.cmu.edu/~quake/robust.html. – tmyklebu Jul 04 '16 at 21:47
  • 1
    You probably want to use a 'relative diff' such as `abs(x - target) / max(abs(x), abs(target)) < epsilon` so that if the numbers are very small (1E-32), you don't accept values that are actually massively different, and similarly don't run into problems with values that are very large (1E+32). – Jonathan Leffler Jul 04 '16 at 22:35