9

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!

amoeba
  • 93,463
  • 28
  • 275
  • 317
Suzy
  • 93
  • 1
  • 3
  • what's a midplane? – MichaelChirico Jul 27 '15 at 17:05
  • 1
    Apologies, midplane is the middle plane of the cranium in this case, also known as the sagittal plane which cuts the skull into two symmetrical halves (left and right) – Suzy Jul 29 '15 at 08:48
  • sounds like machine learning is particularly well suited to this problem – MichaelChirico Jul 29 '15 at 11:51
  • @Suzy, please take a look at this thread http://stats.stackexchange.com/questions/13152. It should answer your question, but if it doesn't or you are still confused, please be more precise about where you get confused. – amoeba Jul 30 '15 at 15:07
  • I have had a look at that thread in the past but it's not that helpful as the examples seem to deal with 2D points/lines, so I wouldn't know how to fit it to my 3D data. Could you explain further? –  Aug 03 '15 at 09:58
  • Yes, but you need to specify more exactly what you want to get. Let's say you have a $27\times 3$ as an input -- right? In what form do you want to obtain the output? A plane can be specified in many different ways. Do you need any particular way? E.g. if you run PCA on your data matrix, then the first two eigenvectors will span your plane. Is it enough for you, or do you need some other representation? – amoeba Aug 03 '15 at 12:55
  • 1
    The format of the output doesn't matter hugely, as long as I can work out the equation for the 3D plane in the following format: ax + by + cz = d. I have read about the first two eigenvectors being related and can run a PCA, however I don't understand what to do with the eigenvectors once I have obtained them, and how they can be used to work out the equation of the plane. Thanks in advance. – Suzy Aug 06 '15 at 10:23
  • Suzy, a tip for the future: put `@amoeba` somewhere in your comment and then I will receive it as a notification in my inbox. I noticed your last comment only by chance. I've posted an answer. – amoeba Aug 06 '15 at 15:15
  • Hi @Suzy! Have you had a chance to look at my answer below? Does it answer your question? Don't hesitate to ask if something remains unclear. – amoeba Aug 18 '15 at 15:41
  • 1
    Sorry for the delayed response @amoeba, this did answer my question and I have now tried out the method and gotten it to work! Many thanks! – Suzy Sep 02 '15 at 11:57
  • @Suzy, I am glad that it worked for you. If you think that your question is now answered, please "accept" my answer by clicking a tick sign to the left of it. You should also upvote all posts that you find helpful, but you might need to wait until you have reputation of 15 to be able to do it. Thanks. – amoeba Sep 02 '15 at 13:05

1 Answers1

9

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:

PCA plane

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}

amoeba
  • 93,463
  • 28
  • 275
  • 317