Sympy: Intersection of segments

Viewed 715

I need to compute intersection of segments, I do so using sympy and its intersection function. It seems to work well when coordinates are given directly, but in my case I construct endpoints iteratively. As an example, let us compute the intersection of segments s0 (blue) and s1 (red) in the following figure:

drawing

This can be done by directly like so:

import sympy
from sympy.geometry import *
from numpy import sqrt, cos, sin

PI = float(sympy.pi)

p0 = Point(0.0, 0.0, evaluate=False)
p1 = Point(1.0, 0, evaluate=False)
p2 = Point(0.5, sqrt(3)/2.0, evaluate=False)
p3 = Point(0.5, -sqrt(3)/2.0, evaluate=False)

s0 = Segment(p0, p1)
s1 = Segment(p2, p3)

print(intersection(s0, s1))
# [Point2D(0.500000000000000, 0)]

Or by starting from p0 and constructing the other points (named q1, q2, q3 for later comparison) by adding the corresponding vectors:

angle = 0.0
q0 = Point(0.0, 0.0, evaluate=False)
q1 = q0 + Point(cos(angle), sin(angle), evaluate=False)
angle = 2.0*PI/3.0
q2 = q1 + Point(cos(angle), sin(angle), evaluate=False)
angle = -0.5*PI
q3 = q2 + Point(sqrt(3)*cos(angle), sqrt(3)*sin(angle), evaluate=False)

r0 = Segment(q0, q1)
r1 = Segment(q2, q3)

print(intersection(r0, r1))
# []

Which fails (the intersection returned is empty).

We can check that p0 = q0, etc., at least numerically:

print("|p0-q0| = ", p0.distance(q0)) # |p0-q0| =  0
print("|p1-q1| = ", p1.distance(q1)) # |p1-q1| =  0
print("|p2-q2| = ", p2.distance(q2)) # |p2-q2| =  2.48253415324727e-16
print("|p3-q3| = ", p3.distance(q3)) # |p3-q3| =  3.51083346857670e-16

I want this to be reasonably fast, so I'm using evaluate=False to only deal with floats, so I cannot make sense of the difference in results. Also note that the cos, sin and sqrt functions are those of the numpy package, so that we (should?) avoid lazy evaluation shenanigans. Am I missing something else?

Edit: this is using python 3.7.4 and sympy 1.4

1 Answers

This looks like a bug in SymPy. You can use srepr() to show the internal representation of s0, s1, r0, r1 and notice they are very similar in structure.

The problem is that your coordinates are almost but not exactly equal. So, this is a more numerical than a symbolic problem.

Note that SymPy shines in symbolic calculations and gives beautiful exact results. For problems with this complexity you really don't win much time letting SymPy do the calculations numerically. My advice would be to either leave out the numerical calculations, or to also do the intersection calculation inside NumPy

Related