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
 A: 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.
A: 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.
