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