1
$\begingroup$

There are numerous techniques to fit a sphere (with unknown centre and radius) through points in $R^3$, such that the fitted sphere passes through the points as closely as possible (in the least-squared error sense, for example), but I have a slightly different problem: instead of points on the sphere, I have points, each of which is a known (measured) distance away from the (unknown) sphere.

To formulate this problem precisely, let's call the four unknowns $cx, cy, cz$ and $r$ (being the x, y and z coordinates of the sphere's centre, and the radius of the sphere), and let's call the $i$th input point $P[i]$, which has coordinates $(P[i]_x, P[i]_y, P[i]_z)$ and distance from the sphere $P[i]_d$.

Then, I believe a solution to the following would be a good solution to the problem:

$min\sum_{i=0}^{n-1} (\sqrt{(P[i]_x - cx)^2 + (P[i]_y - cy)^2 + (P[i]_z - cz)^2} - r - P[i]_d)^2$

What I am wondering is:

  • Is this true? Can the above "minimum sum of squared errors" be expected to yield a reasonable fit?

  • If so, what is a reasonable way to solve this computationally, for cases where there are $\approx 100$ points?

I can compute the derivatives of the above expression with respect to the four unknowns $cx, cy, cz$ and $r$, and could attempt e.g. gradient descent, but is there a more efficient way to solve this problem? For example, can this be recast as a linear or quadratic program?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

After some false starts involving very crude gradient descent code, I think I can now answer this question:

  • Yes, the minimum sum of squared errors formulation does work for this problem, with one caveat: it is possible for the derivative to be undefined, due to a zero denominator. For example, this happened during testing when both the initial guess and the "true" center was the origin.

  • Despite the small number of variables (4), the only reasonable way to solve this problem numerically that I have found so far is to use a dedicated nonlinear least-squares solver (I am using Ceres) which solves cases with $\approx 100$ points, with a residual below $\approx 10^{-15}$ in about one millisecond. I imagine a general purpose nonlinear solver could also do this, but I haven't tried.

Improve this answer
$\endgroup$
1
  • $\begingroup$ If you are talking about the $\frac{c-p_i}{\|c-p_i\|}$ term making the gradient "undefined", I think you should evaluate that to zero if the denominator is zero. That is, if you get a NaN turn it into a 0. As far as gradient descent goes I suppose that it should work, but you may have to do a line search for the step size. But if you are unhappy with it you can look at conjugate gradients, damped Newton or Gauss-Newton/Levenberg-Marquardt. I wanted to make a python code to illustrate that but I just didn't find any time, sorry. $\endgroup$ Commented Dec 9, 2024 at 20:00
0
$\begingroup$

If the points $p_i$ are exactly $d_i$ distance away from the sphere then $\|p_i-c\|^2 - (R + d_i)^2 = 0$. Then you can write the objective:

\begin{align} E(c) &= \sum_{i=1}^n \|p_i-c\|^2 - (d_i + R)^2 \\ &= n\|c\|^2 - 2c^T \sum_{i=1}^n p_i - nR^2 - 2R\sum_{i=1}^nd_i + \sum_{i=1}^n \|p_i\|^2 - d_i^2. \end{align}

To solve $\min_{c} E(c)$ differentiate wrt $c$ and set to zero:

\begin{align} \nabla E(c) &= 2n c - 2 \sum_{i=1}^n p_i = 0 \implies c = \frac{\sum_{i=1}^n p_i}{n}. \end{align}

You get the barycenter.

Since $\|p_i - c\| = R + d_i$ you can write the (original) objective $$G(R,c) = \sum_{i=1}^n (R- (\|p_i -c \|-d_i))^2 = n R^2 - 2R (\|p_i -c \|-d_i) + (\|p_i -c \|-d_i)^2.$$

Differentiating wrt $R$ and setting to zero gives you: $$R = \frac{\sum_{i=1}^n \|p_i -c \|-d_i}{n},$$ where $c$ is the barycenter. Note that this split the problem into two subproblems.

Let's also look at the original formulation in order to figure out what we simplified by considering $E$ instead of $G$:

\begin{align} G(c,R) &= \sum_{i=1}^n (\|p_i -c \| - (R + d_i))^2 \\ &= \|p_i-c\|^2 - 2 \sum_{i=1}^n(R+d_i)\sqrt{c^2-2cp_i + p_i^2} + (R + d_i)^2. \end{align} The gradient wrt $c$ is: $$\nabla_c G(c,R) = 2nc - 2\sum_{i=1}^{n}p_i - 2\sum_{i=1}^n (R+d_i) \frac{c - p_i}{\|c-p_i\|}.$$ The derivative wrt $R$ is as derived previously: $$\frac{d}{dR} G(c,R) = 2nR - 2\sum_{i=1}^n \left(\|p_i-c\|-d_i\right).$$

The above is valid even if $\|p_i-c\| \ne R+ d_i$, e.g., maybe you don't have exact equality due to noise or inexact arithmetic. Since there is no closed form solution for $\nabla_c G = 0$ you can use gradient descent or Newton. The initial guess can be the solution from the first approach.

Since you have an explicit solution for $R$ you could also plug this into the expression $\nabla_c G$, and do away with $R$ entirely, which yields:

\begin{align} \nabla_c G(c,R^*) &= 2nc - 2\sum_{i=1}^{n}p_i - \frac{2}{n}\sum_{j=1}^n\sum_{i=1}^n (\|c-p_j\|-d_j+n d_i) \frac{c - p_i}{\|c-p_i\|}. \end{align}

You can then apply Newton just to the above set to zero, the initial guess being the barycenter.

As noted by Norman in the comment below, the first approach is not good enough, especially if the $p_i$ are not uniformly distributed around $c$.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I'm confused by the first formulation. Suppose there are no errors, and the sphere is centered on the origin with radius 1, and all the points are close to $(1,0,0)$ - the barycenter of these points will be close to $(1,0,0)$, not $(0,0,0)$. Regarding the original formation, I will attempt to implement some Newton iterations and report back. Hopefully it will work! $\endgroup$ Commented Dec 6, 2024 at 11:40
  • $\begingroup$ @Norman You're correct, it seems like the first approach is not a good substitute, at least if the points are not uniformly distributed around the center. I'll think about it and edit the answer if I figure something out. The second approach should still be valid, since it uses your formulation directly, except maybe the barycenter is not the best initial guess. On the other hand you have to solve nonlinear equations for it (or perform gradient descent). $\endgroup$ Commented Dec 6, 2024 at 11:47
  • $\begingroup$ $E(c)$ is not a good formulation because it is unbounded from below: If you let $R\to\infty$, you get $E(c,R)\to -\infty$. $\endgroup$ Commented Dec 6, 2024 at 20:30
  • $\begingroup$ @WolfgangBangerth That's why I explicitly wrote $E(c)$ and not $E(c,R)$, and also why I did not derive the expression for $R$ from it. $\endgroup$ Commented Dec 6, 2024 at 20:43
  • $\begingroup$ However, it is effectively equivalent to just $E(c) = \sum_i \|p_i - c\|^2$ which is why I got the barycenter. So it isn't really useful unless the points are uniformly distributed wrt the sphere or normally distributed. I'll edit that tomorrow. $\endgroup$ Commented Dec 6, 2024 at 20:53

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.