2

A lot of similar questions have already been asked and answered here, but after browsing them none quite solved my issue. I am looking for a reliable algorithm to find the intersection of two infinite lines specified by two points each. In my case, there are two complications:

  1. It is not guaranteed that the intersection point lies between the two vertices specifying each line (which means solutions like this one won't work), and
  2. the lines can be oriented arbitrarily in the 2-D plane, which means that solutions based on slope and intercept won't always work (since one of the lines can be oriented vertically, yielding infinite slope and intercept)

My current approach illustrates the shortomings of a slope-and-interception based approach. This issue could be partially circumvented by implementing a step which rotates the entire system so no line is vertical, but this doesn't seem very elegant. Do you know a better approach?

import numpy as np

# The intersection point of the example below should be (0,0)

# Vertices for the first line
p1_start    = np.asarray([-5,   0])
p1_end      = np.asarray([-3,   0])

# Vertices for the second line
p2_start    = np.asarray([0,    4])
p2_end      = np.asarray([0,    2])

# Calculate slope and intercept for the first line
m_1         = (p1_end[1]-p1_start[1])/(p1_end[0]-p1_start[0])
t_1         = p1_start[1] - m_1*p1_start[0]

# The slope and intercept are zero
print('First line')
print('slope = '+str(m_1))
print('intercept = '+str(t_1))

# Calculate slope and intercept for the second line
m_2         = (p2_end[1]-p2_start[1])/(p2_end[0]-p2_start[0])
t_2         = p2_start[1] - m_2*p2_start[0]

# The slope and intercept are infinite
print('Second line')
print('slope = '+str(m_2))
print('intercept = '+str(t_2))

# Find out where these points interset
# Doesn't work if one of the slopes is infinite
intersection_point_x = (t_2-t_1)/(m_1-m_2)
intersection_point_y = intersection_point_x*m_1 + t_1
print('Intersection point')
print(intersection_point_x)
print(intersection_point_y)
Reblochon Masque
  • 35,405
  • 10
  • 55
  • 80
J.Galt
  • 529
  • 3
  • 15
  • 1
    Are you sure e.g. https://stackoverflow.com/a/38112653/51685 doesn't do the trick? – AKX Oct 21 '20 at 12:21
  • 1
    I am afraid so, yes. The solutions in this thread seem to pertain only to line segments, i.e. finite lines. In my case, the lines specified by the points are infinite, so the intersection point may lie outside the segments the points span up. – J.Galt Oct 21 '20 at 12:28
  • 3
    The https://stackoverflow.com/a/565282/51685 answer has "If r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1" – removing the 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1 conditions should make that solution work for infinite lines. – AKX Oct 21 '20 at 12:34
  • Ah, yes, that seems to work. Thank you! I have updated the question above. – J.Galt Oct 21 '20 at 12:52
  • 1
    Pleas do not alter questions with answers; you can, if you like, post an answer to your own question, in the space reserved for answers. I rolled back your post to your original question. – Reblochon Masque Oct 21 '20 at 13:23
  • Ah, sorry - I wasn't aware that I can answer my own question this early. I'll post the solution AKX provided me with as a separate answer. Thanks for the heads-up! – J.Galt Oct 21 '20 at 14:07

2 Answers2

5

Following AKX's solution, I have adapted the solution reported in his linked thread into a brief Python snippet:

import numpy as np

# The intersection point of the example below should be (0,0)

# Vertices for the first line
p1_start    = np.asarray([-5,   0])
p1_end      = np.asarray([-3,   0])

# Vertices for the second line
p2_start    = np.asarray([0,    4])
p2_end      = np.asarray([0,    2])

p       = p1_start
r       = (p1_end-p1_start)

q       = p2_start
s       = (p2_end-p2_start)

t       = np.cross(q - p,s)/(np.cross(r,s))

# This is the intersection point
i       = p + t*r
J.Galt
  • 529
  • 3
  • 15
2

Why not use, sympy geometry module. below is the code.

import sympy as sy
import sympy.geometry as gm
sy.init_printing()
line1=gm.Line(gm.Point(1,1),gm.Point(8,5)) #Line1
line2=gm.Line(gm.Point(30,4),gm.Point(35,-4)) #Line2
 #These are two infinite lines defined by two points on the line
intersection=line1.intersection(line2)
print(intersection[0].evalf())

Lines intersection