CODEWITHSUNDEEP
pythonfloating-pointprecision
Aguy
Aguy

Reputation: 8059

Arithmetic accuracy issues while testing for intersection of two line segments

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?

Upvotes: 1

Views: 395

Answers (1)

Horia Coman
Horia Coman

Reputation: 8781

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.

Upvotes: 5

Related Questions