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?

-

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.

-
may I know how do you determine $\rho$? –  Graviton Jul 21 '11 at 4:37
@Graviton: One has to minimize with respect to the four variables $a_k$ and $\rho$. The condition $a_1^2+a_2^2+a_3^2=1$ is a constraint that has to be dealt with by means of a Lagrange multiplier. –  Christian Blatter Jul 21 '11 at 8:37

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.

-

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.

-
Or before doing that computation determine using row-reduction which pairs of differences are linearly independent. Can you elaborate on this? And how does this avoid the numerical issue, i.e., when 2 lines are almost collinear, how to avoid using them to construct non-collinear vector? –  Graviton Jul 20 '11 at 13:23
You probably don't need to consider all pairs. Taking $e_i = P_i - P_{i+1}$ (and cyclic) consider only cross products $e_1 \times e_j$ for $j>1$. They should either be zero or in the same direction if the points are coplanar. –  Alice Jul 20 '11 at 14:00