A Better Approach to Plane Intersection

One interesting problem, arises in many geometric algorithms is the intersection of two planes. For sure, this problem has a well studied textbook answer, which I don't find particularly elegant. This is because, the solution is depends on picking an arbitrary point, which lacks geometric intuition. Ideally, we'd like this point to have some meaning, such as being close to the planes, or the line or etc. For that, in this post, I'd like to remind you of a solution by John Krumm, which remains to be unnoticed by many. While the plane-plane intersection problem is well studied, for our purposes, we would like to summarize a different, more effective approach. The textbook solution to the problem depends on picking an arbitrary point, which might lack geometric intuition. Ideally, we'd like this point to have some meaning, such as being close to the points, planes, or the line or etc. Let \(\mathbf{p}=\{p_x,p_y,p_z\}\) and \(\mathbf{n}=\{n_x,n_y,n_z\}\) compose the plane \(\mathbf{P}=\{\mathbf{n},\mathbf{p}\}\). Let there be two planes \(P_1\) and \(P_2\), for which we'd like to compute the intersecting line \(\mathbf{l}\). Let's first review the standard solution. It's trivial to compute the direction as the cross-product: \begin{equation} \mathbf{l}_d=\mathbf{n}_1 \times \mathbf{n}_2 \end{equation} We then pick an arbitrary point (origin, or the midpoint) \(\mathbf{p}_0\), plug it into both the equations and solve for the unknowns. At this stage there is no meaning of \(\mathbf{p}_0\). In the following, we'll pick this point to be the midpoint of the input points. Basically, we desire that the resulting point \(\mathbf{p}\) is as close to the chosen point \(\mathbf{p}_0\) as possible. Thus, we could write a distance : \begin{equation} \lVert \mathbf{p}-\mathbf{p}_0 \rVert = (p_x-p_{0x})^2 + (p_y-p_{0y})^2 + (p_z-p_{0z})^2 \end{equation} Incorporating the other points in the similar fashion, and writing this constraint using Lagrange multipliers into an objective function results in: \begin{equation} J=\lVert \mathbf{p}-\mathbf{p}_0 \rVert+\lambda(\mathbf{p}-\mathbf{p}_1)^2 + \mu(\mathbf{p}-\mathbf{p}_2)^2 \end{equation} Using the standard Lagrange framework (omitting the details), one establishes a nice matrix, in the form: \begin{equation} M= \left[ {\begin{array}{ccccc} 2 & 0 & 0 & n_{1x} & n_{2x}\\ 0 & 2 & 0 & n_{1y} & n_{2y}\\ 0 & 0 & 2 & n_{1z} & n_{2z} \\ n_{1x} & n_{1y} & n_{1z} & 0 & 0\\n_{2x} & n_{2y} & n_{2z} & 0 & 0 \end{array} } \right] \end{equation} This matrix can now be used in a system of linear equations: \begin{equation} \mathbf{M}\left[ {\begin{array}{c} p_x \\ p_y \\ p_z \\ \lambda \\ \mu \end{array} } \right] = \left[ {\begin{array}{c} 2p_{0x} \\ 2p_{0y} \\ 2p_{0z} \\ \mathbf{p}_1 \cdot \mathbf{n}_1 \\ \mathbf{p}_2 \cdot \mathbf{n}_2 \end{array} } \right] \end{equation} to solve for the unknown point, \(\mathbf{p}\), as well as the Lagrange multipliers, \(\{\lambda, \mu\}\). While the multipliers are not of particular interest they would be interesting for understanding configuration of points, or for different parametrizations. Here is the MATLAB script:

% Robust plane to plane intersection
% Plane is parameterized by a point and a normal.
function [p, n] = intersect_planes(p1, n1, p2, n2, p0)

M = [2 0 0 n1(1) n2(1)
     0 2 0 n1(2) n2(2)
     0 0 2 n1(3) n2(3)
     n1(1) n1(2) n1(3) 0 0
     n2(1) n2(2) n2(3) 0 0
    ];

b4 = p1(1).*n1(1) + p1(2).*n1(2) + p1(3).*n1(3);
b5 = p2(1).*n2(1) + p2(2).*n2(2) + p2(3).*n2(3);
b = [2*p0(1) ; 2*p0(2) ; 2*p0(3); b4 ; b5];

x = M\b;
p = x(1:3);
n = cross(n1, n2);
end