Least-squares fitting of sphere

Hypotheses

This page explains how to fit a 3D sphere to a cloud of point by minimizing least squares errors.

The point cloud is given by $n$ points with coordinates $ {x_i, y_i, z_i} $. The aim is to estimate $ x_c $ , $ y_c $, $ z_c $ and $ r $, the parameters of the sphere that fit the best the points :

$ x_c $ is the x-coordinate of the center of the sphere
$ y_c $ is the y-coordinate of the center of the sphere
$ z_c $ is the z-coordinate of the center of the sphere
$ r $ is the radius of the sphere

The equation of the ideal sphere is given by:

$$ (x_i-x_c)^2 + (y_i-y_c)^2 + (z_i-z_c)^2 = r^2 $$

The matrix system

The previous equation is rewritten as:

$$ x_i^2 + x_c^2 - 2x_ix_c + y_i^2 + y_c^2 - 2y_iy_c + z_i^2 + z_c^2 - 2z_iz_c = r^2 $$

, then:

$$ 2x_cx_i + 2y_cy_i + 2z_cz_i + r^2 - x_c^2 - y_c^2 - z_c^2 = x_i^2 + y_i^2 + z_i^2 $$

and:

$$ ax_i + by_i + cz_i + d = x_i^2 + y_i^2 + z_i^2 $$

with:

$ a = 2x_c$
$ b = 2y_c$
$ c = 2z_c$
$ d = r^2 - x_c^2 - y_c^2 - z_c^2 $

The whole system (for all the points) can be rewritten as:

$$ \begin{array}{rcl} ax_1 + by_1 + cz_1 + d & = & x_1^2 + y_1^2 + z_1^2 \\ ax_2 + by_2 + cz_2 + d & = & x_2^2 + y_2^2 + z_2^2 \\ & ... & \\ ax_n + by_n + cz_n + d & = & x_n^2 + y_n^2 + z_n^2 \\ \end{array} $$

The matrix form of the system is given by:

$$ \left[ \begin{matrix} x_1 & y_1 & z_1 & 1 \\ x_2 & y_2 & z_2 & 1 \\ ... & ... & ... & ... \\ x_n & y_n & z_3 & 1 \\ \end{matrix} \right]. \left[ \begin{matrix} a \\ b \\ c \\ d \end{matrix} \right] = \left[ \begin{matrix} x_1^2 + y_1^2 + z_1^2 \\ x_2^2 + y_2^2 + z_2^2 \\ ... \\ x_n^2 + y_n^2 + z_n^2 \\ \end{matrix} \right] $$

Let's define $A$, $B$ and $X$:

$$ \begin{matrix} A=\left[ \begin{matrix} x_1 & y_1 & z_1 & 1 \\ x_2 & y_2 & z_1 & 1 \\ ... & ... & ... \\ x_n & y_n & z_1 & 1 \\ \end{matrix} \right] & B=\left[ \begin{matrix} x_1^2 + y_1^2 + z_1^2 \\ x_2^2 + y_2^2 + z_2^2 \\... \\ x_n^2 + y_n^2 + z_n^2 \\ \end{matrix} \right] & X=\left[ \begin{matrix} a \\ b \\ c \\ d \end{matrix} \right] \end{matrix} $$

The system is now given by:

$$ A.X=B $$

Solving the system

The optimal solution is given by:

$$ \hat{x}=A^{+}.B = A^{T}(A.A^{T})^{-1}.B $$

Where $ A^{+} $ is the pseudoinverse of $ A $. $ A^{+} $ can be computed thanks to the following formula :

$$ A^{+}=A^{T}(A.A^{T})^{-1} $$

As there has been a change of variables, it only remains to calculate $ x_c $ , $ y_c $ and $ r $:

$$ x_c = \dfrac{a}{2} $$ $$ y_c = \dfrac{b}{2} $$ $$ z_c = \dfrac{c}{2} $$ $$ r = \dfrac{ \sqrt{4d + a^2 + b^2 + c^2}}{2} $$

And you get the parameters of the sphere:

Source code

The following Matlab source code was used for drawing the above figure:

close all;
clear all;
clc;

%% Parameters of the system
%% This is the parameters we need to estimate
% Center
C = [30 , 20, 15];
% Radius
R = 12;

% Number of points
N = 1000;
% Noise
k = 0.002;

%% Create noisy sphere
alpha = 2*pi*rand(1,N);
beta  = 2*pi*rand(1,N);
noise = k*2*randn(1,N)-1;
Points = C + [ R*noise.*cos(alpha) ; R*noise.*sin(alpha).*cos(beta) ; R*noise.*sin(alpha).*sin(beta)  ]';

%% Draw the noisy sphere
plot3 (Points(:,1) , Points(:,2), Points(:,3), 'b.');
axis square equal;
grid on;
xlabel ('X');
ylabel ('Y');
zlabel ('Z');

%% Prepare matrices
A = [Points(:,1) Points(:,2) Points(:,3) ones(N,1)];
B = [Points(:,1).*Points(:,1) + Points(:,2).*Points(:,2) + Points(:,3).*Points(:,3)];

%% Least square approximation
X=pinv(A)*B;

%% Calculate sphere parameters
xc = X(1)/2
yc = X(2)/2
zc = X(3)/2
r = sqrt(4*X(4) + X(1)*X(1) + X(2)*X(2) + X(3)*X(3) )/2

%% Draw the fitted sphere
[X, Y, Z] = sphere;
X = xc + X*r;
Y = yc + Y*r;
Z = zc + Z*r;
hold on;
surf (X, Y, Z, 'FaceAlpha', 0.2, 'EdgeAlpha', 0.1);
legend ('Points', 'Least-square sphere');

Download

Matlab source code (example on this page) can be download here: