Fitting a plane to a set of points in 3D using PCA

Question

I am trying to estimate a midplane of a 3D model using the midpoints of paired landmarks, in order to reconstruct missing data (midplane refers here to the middle/saggital plane of the cranium which cuts the skull into two symmetrical halves, left and right).

I therefore need to estimate a plane from 27 points in 3D. I need to get the equation of the plane as $$ax+by+cz=d.$$

I have looked into orthogonal regression and principal component analysis (PCA) as methods, however I didn't take maths past A-levels and am struggling. I know I can supposedly use the eigenvectors to get the equation of the plane of best fit, but need someone to explain exactly how. I am using R for the PCA but am not great at R either.

Alternatively, if there's a better way to estimate the plane I would be glad to hear it!

Answer 1

When you perform principal component analysis (PCA) on your 27 points in 3D, you first subtract the mean vector $\mathbf m$ and then obtain three eigenvectors $\mathbf e_1, \mathbf e_2, \mathbf e_3$ of the covariance matrix. The first two eigenvectors (with two largest eigenvalues) span the plane that you want to find, so the geometric situation looks like that:

The question is: how to get from here to the equation of this plane in the form $$ax+by+cz=d.$$ This equation can be reformulated as follows: $\mathbf a \cdot \mathbf x = d$, where $\mathbf x$ is any point lying in the plane and $\mathbf a$ is a vector $(a,b,c)$. In other words, we need to find a vector $\mathbf a$ such that its dot product with any point in the plane gives the same constant value $d$.

From the picture above, we see that any point belonging to the plane can be written as $\mathbf x = \mathbf m+ g\mathbf e_1 + h\mathbf e_2$, where $g$ and $h$ are some real numbers. It follows that the dot product between $\mathbf e_3$ and $\mathbf x$ is given by $$\mathbf e_3 \cdot \mathbf x = \mathbf e_3 \cdot (\mathbf m+ g\mathbf e_1 + h\mathbf e_2) = \mathbf e_3 \cdot \mathbf m = \mathrm{const}.$$ So here we have it: we can take $\mathbf a = \mathbf e_3$ and $d = \mathbf e_3 \cdot \mathbf m$.

Putting it all together, the solution is: \begin{align} a &= e_{31}\\ b&=e_{32} \\ c&=e_{33} \\ d&=e_{31}m_1+e_{32}m_2 + e_{33}m_3. \end{align}