# How to find a 2D basis within a 3D plane - direct matrix method?

I have a plane equation in 3D, in the form $Ax+By+Cz+D=0$ (or equivalently, $\textbf{x}\cdot\textbf{n} = \textbf{a}\cdot\textbf{n}$), where $\textbf{n}=\left[A\:B\:C\right]^T$ is the plane normal, and $D=-\textbf{a}\cdot\textbf{n}$ with $\textbf{a}$ as a point in the plane, hence $D$ is negative perpendicular distance from the plane to the origin. I have some points in 3D space that I can project onto the plane, and I wish to express these points as 2D vectors within a local 2D coordinate system in the plane. The 2D coordinate system should have orthogonal axes, so I guess this is a case of finding a 2D orthonormal basis within the plane?

There are obviously an infinity of choices for the origin of the coordinate system (within the plane), and the in-plane $x$ axis $\textbf{i}$ may be chosen to be any vector perpendicular to the plane normal. The 3D unit vector of the in-plane $y$ axis $\textbf{j}$ could then be computed as the cross product of the $x$ axis 3D unit vector and the plane normal. One algorithm for choosing these could be:

• Set origin as projection of $\left[0\:0\:0\right]^T$ on plane
• Compute 2D $x$ axis unit vector $\textbf{i}$ direction as $\left[1\:0\:0\right]^T\times\textbf{n}$ (then normalize if nonzero)
• If this is zero (i.e. $\textbf{n}$ is also $\left[1\:0\:0\right]^T$) then use $\textbf{i}$ direction as $\left[0\:1\:0\right]^T\times\textbf{n}$ instead (then normalize)
• Compute $\textbf{j}=\textbf{i}\times\textbf{n}$

However, this all seems a bit hacky (especially the testing for normal along the 3D $x$ axis, which would need to deal with the near-parallel case on a fixed-precision computer, which is where I'll be doing these sums).

I'm sure there should be a nice, numerically stable, matrix-based method to find a suitable $\textbf{i}$ and $\textbf{j}$ basis within the plane. My question is: what might this matrix method be (in terms of my plane equation), and could you explain why it works?

Thanks,

Eric

-
If I understand you correctly the offset $D$ has no impact on the basis you want to find. So we can just choose $D=0$ for the time being. Then you can choose basis of the subspace given by the plane and apply Gram-Schmidt. –  Simon Markett May 31 '12 at 15:46
If this helps I can write some more details,well also if not ;) –  Simon Markett May 31 '12 at 15:50
The only effect the offset $D$ has is on the origin point for the basis - choosing the direction vectors $\textbf{i}$ and $\textbf{j}$ only depends on the normal, so yes, we can assume $D=0$. Please could you give more details on Gram-Schmidt, specifically what the nicest way of computing the two 3D unit vectors for the new coordinate axes using $A$, $B$ and $C$ is? –  Eric Greveson May 31 '12 at 16:13

As we discussed in the comments we can set $D=0$ for the time being, then the plane is not longer affine but a honest to god subspace of $\mathbb R^3$. There we can find two vectors which span this subspace. You want this basis to look nice (orthonormal) but the trick is that Gram-Schmidt allows us to deal with this later.

First we just have to find any two vectors which span then plane $Ax+By+Cz=0$. In other words we just need two linearly independent solutions of this equation. Assume wlog $A\neq 0$ (rearrange if necessary), then two solutions would be for example $(B/A,-1, 0)$ and $(C/A,0,-1)$. These two vectors span the plane. Then you can just apply Gram-Schmidt to these vectors to get a orthonormal basis.

It seems that you are looking for a way to implement this quickly. I am no expert on this but I know for sure that many languages have ready to use packages for Gram-Schmidt.

Edit: Maybe Gram-Schmidt is a really big word for such a low dimensional case. In fact what you do is the following. Call our vectors $b_1$ and $b_2$, then compute $b_1'=\frac {b_1}{|b_1|}$, i.e. make it into a length one vector. Then compute $b_2'=b_2-\left< b_1',b_2\right>b_1'$, this vector is now orthogonal to $b_1'$ and finally take $b_2''=\frac{b_2'}{|b_2'|}$. The vectors $b_1'$ and $b_2''$ are now an orthonormal basis.

-
Simon, thanks for the explanation. Gram-Schmidt therefore sounds very similar to the algorithm I described - the test for A not equal to zero, coming up with an initial vector in the plane based on this, etc. Thinking about it, I was hoping for an alternative method to "automagically" deal with such cases (e.g. if $A$ is very small), but I suppose as the axes are unconstrained due to the degree of freedom about the normal, this won't happen. Essentially I was hoping that there might be a method that uses $A$, $B$, and $C$, and no if-tests. –  Eric Greveson May 31 '12 at 17:06
Well you could sort $A,B,C$ wrt to their absolute value first and take the biggest one. –  Simon Markett May 31 '12 at 17:10
Or you could take the three vectors $(B,-A,0)$, $(C,0,-A)$ and $(0,C,-B)$ and find a different method to decide which one of those is superfluous (usually any one of the three, if any of the parameters is 0 you have limited choice) –  Simon Markett May 31 '12 at 17:16

Your original approach is fine, except for one small improvement: instead of choosing $\textbf{i}$ by default to obtain the first basis vector, you should choose the one of $\{\textbf{i}, \textbf{j}, \textbf{k}\}$ whose angle with $\textbf{n}$ is closest to 90°. This will be the one whose dot product with $\textbf{n}=(n_1, n_2, n_3)$ is smallest (by absolute value). Since $\textbf{i}\cdot \textbf{n} = n_1$ etc., simply choose the one that corresponds to the smallest $|n_i|$ (e.g. choose $\textbf{j}$ if $|n_2|$ is smallest). This guarantees a minimum angle (ca. 54°) between the two vectors and is therefore numerically robust.

-