Plane-plane intersection in python

Question

I need to calculate intersection of two planes in form of AX+BY+CZ+D=0 and get a line in form of two (x,y,z) points. I know how to do the math, but I want to avoid inventing a bicycle and use something effective and tested. Is there any library which already implements this? Tried to search opencv and google, but with no success.

"Questions asking us to recommend or find a book, tool, software library, tutorial or other off-site resource are off-topic for Stack Overflow as they tend to attract opinionated answers and spam. Instead, describe the problem and what has been done so far to solve it." All that being said, perhaps numpy, scipy, or something alike has something suitable. — Ilja Everilä, Jan 06 '18 at 11:27
Also check out http://www.sagemath.org/ and the packages it relies on. But, this is off topic for SO. — trincot, Jan 06 '18 at 11:36
Have you looked into [`sympy`](http://docs.sympy.org/latest/modules/geometry/plane.html) library? — Mr. T, Jan 06 '18 at 11:49
Add a random 3rd equation, which will together with the other 2 yield a point when using numpy.linag.inv. Take that equation away and add another random 3rd equation, which will together with the other 2 yield another point. The line you're looking for runs through these 3 points. Graphically you intersect 2 random planes with your intersection line. — Jacques de Hooge, Jan 06 '18 at 13:01
@JacquesdeHooge: taking random equations is a kind of Russian roulette because you can get close to degeneracies. Much better to choose the planes smartly, as a function of the given coefficients. — , Jan 06 '18 at 16:06
I believe this should be reopened, as if it were phrased "how to do that with numpy" or similar, it would be a great question that is quite hard to find an answer to online. — Gulzar, Jul 24 '21 at 13:01

f5r5e5d · Answer 1 · 2018-01-07T00:19:07.017

my numpy solution:

def plane_intersect(a, b):
    """
    a, b   4-tuples/lists
           Ax + By +Cz + D = 0
           A,B,C,D in order  

    output: 2 points on line of intersection, np.arrays, shape (3,)
    """
    a_vec, b_vec = np.array(a[:3]), np.array(b[:3])

    aXb_vec = np.cross(a_vec, b_vec)

    A = np.array([a_vec, b_vec, aXb_vec])
    d = np.array([-a[3], -b[3], 0.]).reshape(3,1)

# could add np.linalg.det(A) == 0 test to prevent linalg.solve throwing error

    p_inter = np.linalg.solve(A, d).T

    return p_inter[0], (p_inter + aXb_vec)[0]


a, b = (1, -1, 0, 2), (-1, -1, 1, 3)
plane_intersect(a, b)

Out[583]: (array([ 0.,  2., -1.]), array([-1.,  1., -3.]))

a test, subs points back in:

p1, p2 = plane_intersect(a, b)
a_vec, b_vec = np.array(a[:3]), np.array(b[:3])

(np.dot(p1, a_vec), np.dot(p2, a_vec), np.dot(p1, b_vec), np.dot(p2, b_vec))
Out[585]: (-2.0, -2.0, -3.0, -3.0)

score 0 · Accepted Answer · answered Jan 06 '18 at 18:47

Finally I've reused sympy library, by converting Ax+By+Cz+D=0 equation to (n,pt) formulation:

def planeCoeffToPoint(n,d):
    nabs = [abs(i) for i in n]
    i=nabs.index(max(nabs))
    r=[0,0,0]
    r[i]=-d/n[i]
    return r

import sympy as sp
n1=(A1,B1,C1)
D1=...
n2=(A2,B2,C2)
D2=...
pt1=planeCoeffToPoint(n1,D1)
pt2=planeCoeffToPoint(n2,D2)
pl1 = sp.Plane(pt1, normal_vector=n1)
pl2 = sp.Plane(pt2, normal_vector=n2)
r=pl1.intersection(pl2)
rp1=r[0].points[0]
rp2=r[0].points[1]
[float(rp1[0]), float(rp1[1]), float(rp1[2])] # first point on line
[float(rp2[0]), float(rp2[1]), float(rp2[2])] # second point on line

You have to call float() explicitly, becuse sympy might return Zero or Rational object, representing float.

Plane-plane intersection in python

2 Answers2