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:
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
