## Orthogonal least squares fitting of a sphere/2

As said in my previous post, to obtain an orthogonal least squares fitting of a sphere to a cloud of points $P_i = \left[ x_i \hspace{1 mm} y_i \hspace{1 mm} z_i \right], i=0 \cdots m-1$ one should minimize the function

$E(x_c, y_c, z_c, r) = \sum_{i=0}^{m-1} e_i^2 = \\ \\ \hspace{5 mm} \sum_{i=0}^{m-1} (\sqrt{(x_i - x_c)^2+(y_i - y_c)^2 + (z_i - z_c)^2} - r)^2$

Dave Eberly calls this the energy function, probably as a metaphorical reference to the minimum total potential energy principle.

Setting $L_i (x_c, y_c, z_c) = \sqrt{(x_i - x_c)^2 + (y_i - y_c)^2 + (z_i - z_c)^2}$

and considering that

$\dfrac{\partial L_i}{\partial x_c} = \dfrac{x_c - x_i}{L_i}$

$\dfrac{\partial L_i}{\partial y_c} = \dfrac{y_c - y_i}{L_i}$

$\dfrac{\partial L_i}{\partial z_c} = \dfrac{z_c - z_i}{L_i}$

$\dfrac{\partial L_i}{\partial r} = 0$

one can write the components of the gradient as

$\dfrac{\partial E}{\partial x_c} = 2 m x_c - 2 \sum_{i=0}^{m-1} x_i - 2 r \sum_{i=0}^{m-1} \dfrac{x_c - x_i}{L_i}$

$\dfrac{\partial E}{\partial y_c} = 2 m y_c - 2 \sum_{i=0}^{m-1} y_i - 2 r \sum_{i=0}^{m-1} \dfrac{y_c - y_i}{L_i}$

$\dfrac{\partial E}{\partial z_c} = 2 m z_c - 2 \sum_{i=0}^{m-1} z_i - 2 r \sum_{i=0}^{m-1} \dfrac{z_c - z_i}{L_i}$

$\dfrac{\partial E}{\partial r} = 2 m r - 2 \sum_{i=0}^{m-1} L_i$

Each $L_i$ in the denominators approaches zero when the centre gets close to  $P_i$, what should not happen if a sufficiently good initial approximation is chosen.

In the minimum point one expects a null gradient. Eberly suggests to solve the equations $\dfrac{\partial E}{\partial x_c} = 0$, $\dfrac{\partial E}{\partial y_c} = 0$, $\dfrac{\partial E}{\partial z_c} = 0$ and $\dfrac{\partial E}{\partial r} = 0$ with a fixed point iteration, and says:

Warning. I have not analyzed the convergence properties of this algorithm. In a few experiments it seems to converge just fine

In my experience, fixed point iteration resulted very slow. I had better results with Levenberg-Marquardt.