Finding the line along the intersection of two planes

Question

I am trying to draw the line formed by the intersections of two planes in 3D, but I am having trouble understanding the math, which has been explained here and here.

I tried to figure it out myself, but the closest that I got to a solution was a vector pointing along the same direction as the intersection line, by using the cross product of the normals of the planes. I have no idea how to find a point on the intersection line, any point would do. I think that this method is a dead end. Here is a screenshot of this attempt:

I tried to use the solution mentioned in this question, but it has a dead link to the original explanation, and the equation didn't work for me (it has unbalanced parentheses, which I tried to correct below).

var planeA = new THREE.Plane((new THREE.Vector3(0, 0, 1)).normalize(), 100);
var planeB = new THREE.Plane((new THREE.Vector3(1, 1, 1)).normalize(), -100);

var x1 = planeA.normal.x,
    y1 = planeA.normal.y,
    z1 = planeA.normal.z,
    d1 = planeA.constant;

var x2 = planeB.normal.x,
    y2 = planeB.normal.y,
    z2 = planeB.normal.z,
    d2 = planeB.constant;

var point1 = new THREE.Vector3();
point1.x = 0;
point1.z = (y2 / y1) * (d1 - d2) / (z2 - z1 * y2 / y1);
point1.y = (-z1 * point1.z - d1) / y1;

var point2 = new THREE.Vector3();
point2.x = 1;
point2.z = (y2 / y1) * (x1 * point2.x + d1) - (x2 * point2.x - d2) / (z2 - z1 * y2 / y1);
point2.y = (-z1 * point2.z - x1 * point2.x - d1) / y1;

console.log(point1, point2);

output:

THREE.Vector3 {x: -1, y: NaN, z: NaN, …}
THREE.Vector3 {x: 1, y: Infinity, z: -Infinity, …}

expected output:

• A point along the intersection where x = 0, and
• Another point on the same line where x = 1

If someone could point me to a good explanation of how this is supposed to work, or an example of a plane-plane intersection algorithm, I would be grateful.

Answer 1

Here is an implementation of a solution for plane-plane intersections described at http://geomalgorithms.com/a05-_intersect-1.html . Essentially, you first use the cross product of the normals of the planes to find the direction of a line in both planes. Secondly, you use some algebra on the implicit equation of the planes (P . n + d = 0 where P is some point on the plane, n is the normal and d is the plane constant) to solve for a point which is on the intersection of the planes and also on one of the x=0, y=0 or z=0 planes. The solution is then the line described by a point and a vector. I was using three.js version 79

/*

Algorithm taken from http://geomalgorithms.com/a05-_intersect-1.html. See the
section 'Intersection of 2 Planes' and specifically the subsection
(A) Direct Linear Equation

*/
function intersectPlanes(p1, p2) {

  // the cross product gives us the direction of the line at the intersection
  // of the two planes, and gives us an easy way to check if the two planes
  // are parallel - the cross product will have zero magnitude
  var direction = new THREE.Vector3().crossVectors(p1.normal, p2.normal)
  var magnitude = direction.distanceTo(new THREE.Vector3(0, 0, 0))
  if (magnitude === 0) {
    return null
  }

  // now find a point on the intersection. We use the 'Direct Linear Equation'
  // method described in the linked page, and we choose which coordinate
  // to set as zero by seeing which has the largest absolute value in the
  // directional vector

  var X = Math.abs(direction.x)
  var Y = Math.abs(direction.y)
  var Z = Math.abs(direction.z)

  var point

  if (Z >= X && Z >= Y) {
    point = solveIntersectingPoint('z', 'x', 'y', p1, p2)
  } else if (Y >= Z && Y >= X){
    point = solveIntersectingPoint('y', 'z', 'x', p1, p2)
  } else {
    point = solveIntersectingPoint('x', 'y', 'z', p1, p2)
  }

  return [point, direction]
}


/*

This method helps finding a point on the intersection between two planes.
Depending on the orientation of the planes, the problem could solve for the
zero point on either the x, y or z axis

*/
function solveIntersectingPoint(zeroCoord, A, B, p1, p2){
    var a1 = p1.normal[A]
    var b1 = p1.normal[B]
    var d1 = p1.constant

    var a2 = p2.normal[A]
    var b2 = p2.normal[B]
    var d2 = p2.constant

    var A0 = ((b2 * d1) - (b1 * d2)) / ((a1 * b2 - a2 * b1))
    var B0 = ((a1 * d2) - (a2 * d1)) / ((a1 * b2 - a2 * b1))

    var point = new THREE.Vector3()
    point[zeroCoord] = 0
    point[A] = A0
    point[B] = B0

    return point
}


var planeA = new THREE.Plane((new THREE.Vector3(0, 0, 1)).normalize(), 100)
var planeB = new THREE.Plane((new THREE.Vector3(1, 1, 1)).normalize(), -100)

var [point, direction] = intersectPlanes(planeA, planeB)

Answer 2

When I have problems like this, I usually let a symbolic algebra package (Mathematica in this case) deal with it. After typing

In[1]:= n1={x1,y1,z1};n2={x2,y2,z2};p={x,y,z};

In[2]:= Solve[n1.p==d1&&n2.p==d2,p]    

and simplifying and substituting x=0 and x=1, I get

                d2 z1 - d1 z2        d2 y1 - d1 y2
Out[5]= {{{y -> -------------, z -> ----------------}}, 
                y2 z1 - y1 z2       -(y2 z1) + y1 z2

            d2 z1 - x2 z1 - d1 z2 + x1 z2
>    {{y -> -----------------------------, 
                    y2 z1 - y1 z2

            d2 y1 - x2 y1 + (-d1 + x1) y2
>      z -> -----------------------------}}}
                  -(y2 z1) + y1 z2

Answer 3

It is easy to let three.js solve this for you.

If you were to express your problem in matrix notation

m * x = v

Then the solution for x is

x = inverse( m ) * v

We'll use a 4x4 matrix for m, because three.js has an inverse() method for the Matrix4 class.

var x1 = 0,
    y1 = 0,
    z1 = 1,
    d1 = 100;

var x2 = 1,
    y2 = 1,
    z2 = 1,
    d2 = -100;

var c = 0; // the desired value for the x-coordinate

var v = new THREE.Vector4( d1, d2, c, 1 );

var m = new THREE.Matrix4( x1, y1, z1, 0, 
                           x2, y2, z2, 0,
                           1,  0,  0,  0,
                           0,  0,  0,  1
                         );

var minv = new THREE.Matrix4().getInverse( m );

v.applyMatrix4( minv );

console.log( v );

The x-component of v will be equal to c, as desired, and the y- and z-components will contain the values you are looking for. The w-component is irrelevalent.

Now, repeat for the next value of c, c = 1.

three.js r.58

Answer 4

Prerequisites

---

Recall that to represent a line we need a vector describing its direction and a point through which this line goes. This is called parameterized form:

line_point(t) = t * (point_2 - point_1) + point_1

where point_1 and point_2 are arbitrary points through which the line goes, and t is a scalar which parameterizes our line. Now we can find any point line_point(t) on the line if we put arbitrary t into the equation above.

NOTE: The term (point_2 - point_1) is nothing, but a vector describing the direction of our line, and the term point_1 is nothing, but a point through which our line goes (of course point_2) would also be fine to use too.

The Algorithm

---

1. Find the direction direction of the intersection line by taking cross product of plane normals, i.e. direction = cross(normal_1, normal_2).

2. Take any plane, for example the first one, and find any 2 distinct points on this plane: point_1 and point_2. If we assume that the plane equation is of the form a1 * x + b1 * y + c1 * z + d1 = 0, then to find 2 distinct points we could do the following:

   y1 = 1
   z1 = 0
   x1 = -(b1 + d1) / a1
   
   y2 = 0
   z2 = 1
   x2 = -(c1 + d1) / a1
   

   where point_1 = (x1, y1, z1) and point_2 = (x2, y2, z2).

3. Now that we have 2 points, we can construct the parameterized representation of the line lying on this first plane: line_point(t) = t * (point_2 - point_1) + point_1, where line_point(t) describes any point on this line, and t is just an input scalar (frequently called parameter).

4. Find the intersection point intersection_point of the line line_point(t) and the second plane a2 * x + b2 * y + c2 * z + d2 = 0 by using the standard line-plane intersection algorithm (pay attention to the Algebraic form section as this is all you need to implement line-plane intersection, if you haven't done so already).

5. Our intersection line is now found and can be constructed in parameterized form as usual: intersection_line_point(s) = s * direction + intersection_point, where intersection_line_point(s) describes any point on this intersection line, and s is parameter.

NOTE: I didn't read this algorithm anywhere, I've just devised it from the top of my head based on my knowledge of linear algebra. That doesn't mean that it doesn't work, but it might be possible that this algorithm can be optimized further.

Conditioning

---

When 2 normal vectors normal_1 and normal_2 are almost collinear this problem gets extremely ill-conditioned. Geometrically it means that the 2 planes are almost parallel to each other and determining the intersection line with acceptable precision becomes impossible in finite-precision arithmetic which is floating-point arithmetic in this case.