Linear least squares fitting of a circumference in 3D

In my last post I described how it is possible to find the centre and radius of a sphere given four or more points on its surface and not belonging to the same plane.

What happens if the points do belong to the same plane?

Let me assume that the points belong to a circumference on the sphere (the intersection of the sphere with the plane). Be the points \left[ x_i y_i z_i \right], i=1 \cdots m, with m \geq 3. The linear system

(1) \alpha x_i + \beta y_i + \gamma z_i + \delta = \epsilon_i

in the unknowns \alpha, \beta, \gamma, \delta (see my previous post) is underdetermined (rank 3) and does not individuate one single sphere, but rather a pencil of spheres sharing the common circumference.

So in this conditions the method is not useful for finding a sphere, but it can be used to find the circumference in 3D containing the points: here is how.

The m x 4 matrix of the system, A=\left[ \begin{array}{cccc} x_1 & y_1 & z_1 & 1 \\ \cdots & \cdots & \cdots & \cdots \\ x_m & y_m & z_m & 1\end{array} \right], can be decomposed in singular values as

A = U \cdot W \cdot V^T,

where U (m x 4) is the matrix having the left singular vectors as columns, W (diagonal 4 x 4) is the matrix of the singular values w_{1 1}, \cdots w_{4 4} and V (4 x 4) is the matrix having the right singular vectors as columns.

Being 3 the rank of the matrix, one of the singular values is zero (within machine precision) and the corresponding right singular vector spans the nullspace of A.

As all the points lie on the circumference, the second hand vector \left[ \epsilon_1 \cdots \epsilon_m\right] ^T lies in the range of A and a solution can be found as

(2) \left[ \alpha^{(0)} \hspace{1 mm} \beta^{(0)} \hspace{1 mm} \gamma^{(0)} \hspace{1 mm} \delta^{(0)}\right] ^T = V W' U^T \left[ \epsilon_1 \cdots \epsilon_m\right] ^T

remembering that V^{-1} = V^T and U^{-1} = U^T; W' is the diagonal matrix having w'_{i i} = \dfrac{1}{w_{i i}} if w_{i i} > 0 and w'_{i i} = 0 if w_{i i} = 0.

As always with underdetermined system, there is an infinity of solutions, obtained by adding to (2) any vector in the nullspace. In this case, the nullspace is spanned by the right singular vector corresponding to the null singular value, say v = \left[ v_{\alpha} \hspace{1 mm} v_{\beta} \hspace{1 mm} v_{\gamma} \hspace{1 mm} v_{\delta} \right]^T.

So \left[ \alpha \hspace{1 mm} \beta \hspace{1 mm} \gamma \hspace{1 mm} \delta \right]^T = \left[ \alpha^{(0)} \hspace{1 mm} \beta^{(0)} \hspace{1 mm} \gamma^{(0)} \hspace{1 mm} \delta^{(0)}\right]^T + \tau \left[ v_{\alpha} \hspace{1 mm} v_{\beta} \hspace{1 mm} v_{\gamma} \hspace{1 mm} v_{\delta}  \right]^T is the parametric equation of the pencil of spheres sharing the circumference one is looking for.

Now, the circumference will be an equator of the smallest sphere in the pencil. So, remembering from my previous post that \alpha = - 2 x_c\beta = - 2 y_c, \gamma = - 2 z_c and \delta = x_c^2 + y_c^2 + z_c^2 - r^2,

\min_{\tau} r^2 = \min_{\tau} x_c^2 + y_c^2 + z_c^2 - \delta = \min_{\tau} \dfrac{1}{4} (\alpha^2 + \beta^2 + \gamma^2) - \delta

= \min_{\tau} \dfrac{1}{4} ((\alpha^{(0)} + \tau v_{\alpha})^2 + (\beta^{(0)} + \tau v_{\beta})^2 + (\gamma^{(0)} + \tau v_{\gamma})^2) - \delta^{(0)} - \tau v_{\delta}

By equating to zero the derivative of the expression to be minimized one gets:

\tau = \dfrac{2 v_{\delta} - \alpha^{(0)} v_{\alpha} - \beta^{(0)} v_{\beta} - \gamma^{(0)} v_{\gamma}}{v_{\alpha}^2 + v_{\beta}^2 + v_{\gamma}^2}

whence the coordinates x_c, y_c, z_c of the centre and the radius r of the circumference can be obtained, while the the plane of the circumference is orthogonal to \left[ v_{\alpha} \hspace{1 mm} v_{\beta} \hspace{1 mm} v_{\gamma} \hspace{1 mm} \right]^T.

Advertisements
This entry was posted in Uncategorized and tagged , , , , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s