How to find line intersections with single line equation array in a pythonic way

Question

What would the most pythonic way be, without using for loops, of finding line intersection points in an array comprised of m,c values?

lines=np.array([m0,c0],
               [m1,c1],
               [m2,c2],
               ....)

achieving the desired result with for loops would consist of something like:

for i in lines:
 for n in lines:
   np.linalg.solve(i, n)

Doesnt seem like the first argument to solve will be a square matrix in this setup; so how is that supposed to work? — Eelco Hoogendoorn
– Eelco Hoogendoorn, Commented Jul 18, 2017 at 9:32
instead of the nested for-loop, you can use itertools.combinations to generate the combinations — Maarten Fabré
– Maarten Fabré, Commented Jul 18, 2017 at 10:35
@WillemVanOnsem, yes the m is the slope, and select any two lines in the matrix and the matrix will become square, and thanks will look into [link]docs.python.org/3/library/itertools.html#itertools.combinations — Alex
– Alex, Commented Jul 18, 2017 at 12:56

MB-F · Accepted Answer · 2017-07-18 11:37:38Z

5

The equation for the intersection of two lines y1 = a1*x + b1 and y2 = a2*x + b2 is x = (b2 - b1) / (a1 - a2).

By making use of broadcasting it is easy to compute all intersections between any number of lines:

import numpy as np    

# lines of the form y = a * x + b
# with lines = [[a0, b0], ..., [aN, bN]]
lines = np.array([[1, 0], [0.5, 0], [-1, 3], [1, 2]])

slopes = lines[:, 0]  # array with slopes (shape [N])
slopes = slopes[:, np.newaxis]  # column vector (shape [N, 1])

offsets = lines[:, 1]  # array with offsets (shape [N])
offsets = offsets[:, np.newaxis]  # column vector (shape [N, 1])

# x-coordinates of intersections
xi = (offsets - offsets.T) / (slopes.T - slopes) 

# y-coordinates of intersections
yi = xi * slopes + offsets

This works by appling the element-wise - operator to a column vector of shape [N, 1] and it's transpose of shape [1, N]. The vectors are broadcast to a matrix of shape [N, N].

The final result are two symmetric matrices xi and yi. Each entry xi[m, n] is the intersection of lines m and n. nan means the lines are identical (they intersect in every point). inf means the lines do not intersect.

Let's show off the result:

#visualize the result
import matplotlib.pyplot as plt
for l in lines:
    x = np.array([-5, 5])
    plt.plot(x, x * l[0] + l[1], label='{}x + {}'.format(l[0], l[1]))
for x, y in zip(xi, yi)    :
    plt.plot(x, y, 'ko')
plt.legend()
plt.show()

edited Jul 18, 2017 at 11:37

answered Jul 18, 2017 at 10:39

MB-F

23.8k5 gold badges71 silver badges127 bronze badges

Sign up to request clarification or add additional context in comments.

6 Comments

Daniel F Over a year ago

Might be important to highlight that when you do slopes = lines[:, [0]] (as opposed to lines[:, 0]) you're creating a 2D slice instead of the more common 1D slice. This lets you broadcast using .T instead of the more verbose (but also more common) [None, :]. Not everyone can parse fancy indexing in their head like that :)

MB-F Over a year ago

@DanielF good point :) Are you sure something like lines[:, 0][:, np.newaxis] or lines[:, np.newaxis, 0] is more readable?

Daniel F Over a year ago

It certainly isn't more readable in 2D, but it is more explicit, and therefore more readable when you get beyond 2 dimensions (in my opinion). For 2D your method is much more readable, but the fancy indexing used to set it up may be confusing for beginners.

Daniel F Over a year ago

Not to confuse - normally I wouldn't construct the lines or offsets arrays to be 2D, but would do xi = (offsets - offsets[none, :]) / (slopes[none, :] - slopes) - making xi 2d explicitly during the broadcast step. This is less readable than your method, but becomes more readable (imo) as dimensionality gets larger.

MB-F Over a year ago

I agree. I would even write (offsets[:, None] - offsets[None, :]) to make it super explicit.

|

Collectives™ on Stack Overflow

How to find line intersections with single line equation array in a pythonic way

1 Answer 1

6 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

6 Comments

Your Answer

Sign up or log in

Post as a guest

Related