1
votes

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?

1
Use the decimal class: from decimal import Decimal. stackoverflow.com/questions/2986150/python-floating-number . Also your isininterval function could simply return: return x1 <= x0 <= x2 or x2 <= x0 <= x1 - Bahrom
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
In the first place, do you have a good reason to worry about 2.9999999999999996 not being 3 ? - Yves Daoust
@YvesDaoust, well, my only worry is that it fails the condition... - Aguy
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 ? - Yves Daoust

1 Answers

5
votes

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.