Find the Plane a 3D  (possibly degenerate) Face is Lying On Is there any sure fire way to find the plane which a 3D face ( defined by points $P_1$, $P_2$, $P_3$, $P_n$) is lying on? One might be tempted to say that obtain the following orthonormal  vectors
$$e_1=P_2-P_1$$
$$e_2=P_3-P_2$$
would do the job, but the problem is that $P_1$, $P_2$ and $P_3$ may all lie on a single line and hence the above vectors are not useful. 
Of course one can always loop through all the points until one finds two different sets of orthonormal vectors, but it is messy and subjected to numerical error when implemented on a computer.
Is there any other more elegant, sure fire way of finding the 3D plane on which the 3D face is lying on?
 A: If you know the polygon is a simple quadrilateral, then the diagonals $(P_1-P_3)$ and $(P_2-P_4)$ should give you non-collinear vectors.  (They might be close to collinear if your quadrilateral is very long and narrow, but in that case you'll have numerical issues anyway.)
For polygons with more sides, the worst case is that all but one of the points might lie on the same line, so in general you can't avoid looping through the points.  One possibility might be to take the cross product $(P_{k-1}-P_k)\times(P_{k+1}-P_k)$ of the adjacent sides for each vertex $P_k$ and sum them together: this should give you a sort of "consensus normal" of the plane the polygon lies on, hopefully with any numerical inaccuracies averaged out.
A: Here is another way to look at the problem. 
As the given 3D face may have $n>3$ vertices the numerically given points $P_i=(p_i^1, p_i^2, p_i^3)$ $\ (1\leq i\leq n)$ will not lie in a plane $a_1x^1+a_2x^2+a_3x^3=\rho$. Therefore finding a covector $a=(a_1, a_2, a_3)$ and a $\rho\geq 0$ such that $|a|^2=1$ and
$$\sum_{k=1}^3 a_k\ p_i^k =\rho\qquad(1\leq i\leq n)\qquad(*)$$
is actually an overdetermined system, and it will not be possible to fulfill all $n$ equations $(*)$ simultaneously. 
Now the error $\bigl|\sum_{k=1}^3 a_k\ p_i^k -\rho\bigr|$ is equal to the distance of the point $P_i$ from the plane $a_1x^1+a_2x^2+a_3x^3=\rho$. It is therefore reasonable to choose the $a_k$ (with $|a|=1$)  and $\rho\geq 0$ in such a way that
$$\sum_{i=1}^n\bigl(\sum_{k=1}^3 a_k\ p_i^k -\rho\bigr)^2$$
becomes minimal. This least squares principle is a standard tool in numerical linear algebra; see, e.g., Golub/van Loan: Matrix computations.
Entering the words "least squares hyperplane points" into Google one obtains a host of papers addressing exactly the problem posed here.
A: Take the cross product of pairs of differences $(P_i - P_j) \times (P_k - P_\ell)$ one at a time until you find a non-zero vector. Or before doing that computation determine using row-reduction which pairs of differences are linearly independent.
