Fitting a plane to points using SVD

Question

I am trying to find a plane in 3D space that best fits a number of points. I want to do this using SVD. To calculate the SVD:

1. Subtract the centroid of the points from each point.
2. Put the points in an mx3 matrix.
3. Calculate the SVD (e.g. [U, S, V] = SVD(A)).

The last column of V, (e.g. V(:,3)), is supposed to be a normal vector to the plane. While the other two columns of V (e.g. V(:,1) and V(:,2)) are vectors parallel to the plane (and orthogonal to each other). I want to find the equation of the plane in ax+by+cz+d=0 form. The last column of V (e.g. V(:,3)) gives "a", "b", and "c", however, in order to find "d", I need a point on the plane to plug in and solve for d. The problem is that I don't know what are valid points to use to plug in.

My question is: does the centroid of the points necessarily lie on the fitted plane? If so, then it's easy to just plug in the centroid values in the equation (along with the from the norm) and solve for "d". Otherwise, how can I calculate "d" in the above equation? The matrix U apparently gives the point values but I don't understand which values to take.

Answer 1

Suggested by @Ailurus (thanks!!!)

Excerpted from: PRINCIPAL AXES AND BEST-FIT PLANES, WITH APPLICATIONS, by Christopher M. Brown, University of Rochester.

Consider the problem of finding a plane which “best fits” a swarm of weighted points. If the points are in n-space, the plane is a hyperplane; we will refer to it as a plane. Represent k-weighted points in n-space by row n-vectors x$_{\mathrm{i}}$, i=1, 2, ..., k; let the weight of the ith point be w$_{\mathrm{i}}$. Represent an n-1-dimensional subspace Π of n-space (a hyperplane) by a unit n-vector $\vec {\boldsymbol{z}}$ normal to Π and a point $\vec {\boldsymbol{v}}$ in Π. In the "sequel" (the following equations), all summations run from 1 to k.

The signed perpendicular distance from x$_{\mathrm{i}}$ to the plane (Π) is: \begin{equation*} \mathrm{d}_{\mathrm{i}Π }=\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right)\cdot \vec {\boldsymbol{z}}. \end{equation*}

The error measure we wish to minimize is the sum over all points of the square of this distance times the weight (mass) of the point, i.e. \begin{equation*} \begin{array}{c} e=\sum _{i}\left(\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right)\cdot \vec {\boldsymbol{z}}\right)^{2}w_{i}=\sum _{i}\vec {\boldsymbol{z}}\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right)^{T}\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right)\vec {\boldsymbol{z}}^{T}w_{i}=\vec {\boldsymbol{z}}\boldsymbol{M}\vec {\boldsymbol{z}}^{T}. \left(1\right) \end{array} \end{equation*} The equation (1), ${\boldsymbol{M}}$ is a real, symmetric n x n matrix, sometimes called the "scatter matrix" of the points. \begin{equation*} \begin{array}{c} \boldsymbol{M}=\sum _{i}w_{i}\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right)^{T}\left(\vec {\boldsymbol{x}}_{i}-\vec {\boldsymbol{v}}\right). \left(2\right) \end{array} \end{equation*} First, we show that the best-fit plane passes through the center of mass (C. of M.) of the points. (The answer to the original question.)

$\underline{Proposition 1}.$ For ${e}$ of equation (1) to be minimized, the plane must pass through the C. of M. of the point swarm. Thus, in equation (1), ${e}$ may attain a minimum when $\vec {\boldsymbol{v}}$ is the C. of M.

$\underline{Proof}$: \begin{equation*} e=\sum _{i}w_{i}\vec {\boldsymbol{z}}\left(\vec {\boldsymbol{x}}_{i}\cdot \vec {\boldsymbol{x}}_{i}\right)\vec {\boldsymbol{z}}^{T}-2\sum _{i}w_{i}\vec {\boldsymbol{z}}\left(\vec {\boldsymbol{x}}_{i}\cdot \vec {\boldsymbol{v}}\right)\vec {\boldsymbol{z}}^{T}+\sum _{i}w_{i}\vec {\boldsymbol{z}}\left(\vec {\boldsymbol{v}}\cdot \vec {\boldsymbol{v}}\right)\vec {\boldsymbol{z}}^{T}. \end{equation*} Since $\vec {\boldsymbol{z}}\vec {\boldsymbol{z}}^{T}=1$ by definition, \begin{equation*} \frac{\partial e}{\partial \vec {\boldsymbol{v}}}=-2\sum _{i}w_{i}\vec {\boldsymbol{x}}_{i}+2\vec {\boldsymbol{v}}\sum _{i}w_{i}; \end{equation*} Setting $\frac{\partial e}{\partial \vec {\boldsymbol{v}}}=0$ implies \begin{equation*} \vec {\boldsymbol{v}}=\frac{\sum _{i}w_{i}\vec {\boldsymbol{x}}_{i}}{\sum _{i}w_{i}}, \end{equation*} which is the center of mass.

So, it is possible to find a point in the best-fit plane; the plane would be determined completely if a normal vector for it could be obtained. (The original paper explains how to do this in the subsequent propositions/proofs.)

Answer 2

Let $c$ denote the centroid of the points $x_1, x_2,\ldots x_n$.

In the first step, you translated the centroid $c$ to the origin, and obtained new points $x_i'=x_i-c$. If you do not believe me that the origin is the centroid of $x'_i$, then you should just compute it yourself to see.

In the last step the plane normal to your vector contains the origin. The normal plane, call it $P'$, is the best fit through $x'_1,\ldots, x'_n$, and it contains the origin, since it is defined as the plane perpendicular to a vector, and every such plane contains the origin.

Now, translating everything by a vector is an isometry. This is true both for $x\mapsto x-c$ and $x\mapsto x+c$. They preserve all distances, and therefore preserve all the nice properties of a plane fitted to points.

Since $P'$ you find in step 3 is optimally fitted to the translated points, we may use the inverse translation by translating everything back with $x\mapsto x+c$ to find the best fit for the original points. Every point $p$ on the plane from step 3 goes to $p+c$. This produces a new plane $P$ which fits the original $x_1,x_2,\ldots x_n$.

Since the centroid (the origin, $(0,0,0)$) of $x'_1, \ldots, x'_n$ was on the unshifted plane, the shifted origin $(0,0,0)+c=c$ is on the fitted plane of the original points.