Estimate center and radius of a sphere from points on the surface

Question

If we assume that our data points were sampled from the surface of a sphere (with some perturbation), how can we recover the center of that sphere?

In my searching, I found papers on something labeled "spherical regression", but it didn't quite seem like that was doing the same thing. Maybe I just didn't understand it.

Is there a straightforward formula, similar to linear regression, that finds a sphere center point and radius that minimize the sum-squared distance of a set of data points from the surface of the sphere?

---

Edit 1:

We can assume that the noise will be 2 or 3 orders of magnitude smaller than the radius of the sphere and uniformly spherically Gaussian. However, the samples themselves will definitely not be drawn uniformly from the sphere's surface, but will likely be clustered in a few patches on the surface, likely all within one hemisphere. A solution that works for data in $\mathbb R^3$ is fine, but a general solution for arbitrary dimensionality is great too.

---

Edit 2:

What are the chances that I might get a sensible answer if I were to use linear regression, $y = X\beta + \epsilon$, in the 7 dimensional space pretending that the squared components are independent from the other parameters:

$\begin{align} X &= \begin{array}{ccccccc}[-2x& -2y&-2z&1&1&1&-1]\end{array}\\ \beta &= \begin{array}{ccccccc}[x_0 & y_0 & z_0 & x_0^2 & y_0^2 & z_0^2 & r^2]'\end{array}\\ y &= x^2+y^2+z^2\end{align}$

At best, I suppose that my error metric will be a bit wacky. At worst the solution won't be even close to consistent.
...or that's silly because with four identical columns, we get a singular matrix when we try to do regression.

---

Edit 3:

So, it looks like these are my options:

1. Non-linear numerical optimization using some cost function: $f(x_0,y_0,z_0,r|X) = \frac{1}{2}\sum_{i=1}^n \left(r - \sqrt{(x_i-x_0)^2+(y_i-y_0)^2+(z_i-z_0)^2}\right)^2$
2. Hough-transform: discretize the plausible space or possible centers and radii around the data points. Each point casts a vote for the potential centers that it could be part of at each specific radius discretization. Most votes wins. This might be okay if there were potentially an unknown number of spheres, but with just one it's a messy solution.
3. Randomly (or systematically) select groups of 4 points and analytically compute the center. Reject the sampling if ill-conditioned (points are nearly co-planar). Reject outliers and find the mean center. From that we can find the mean radius.

Does anyone have a better method?

Answer 1

Here is some R code that shows one approach using least squares:

# set parameters

mu.x <- 8
mu.y <- 13
mu.z <- 20
mu.r <- 5
sigma <- 0.5

# create data
tmp <- matrix(rnorm(300), ncol=3)
tmp <- tmp/apply(tmp,1,function(x) sqrt(sum(x^2)))

r <- rnorm(100, mu.r, sigma)

tmp2 <- tmp*r

x <- tmp2[,1] + mu.x
y <- tmp2[,2] + mu.y
z <- tmp2[,3] + mu.z


# function to minimize
tmpfun <- function(pars) {
    x.center <- pars[1]
    y.center <- pars[2]
    z.center <- pars[3]
    rhat <- pars[4]

    r <- sqrt( (x-x.center)^2 + (y-y.center)^2 + (z-z.center)^2 )
    sum( (r-rhat)^2 )
}

# run optim
out <- optim( c(mean(x),mean(y),mean(z),diff(range(x))/2), tmpfun )
out


# now try a hemisphere (harder problem)

tmp <- matrix(rnorm(300), ncol=3)
tmp[,1] <- abs(tmp[,1])
tmp <- tmp/apply(tmp,1,function(x) sqrt(sum(x^2)))

r <- rnorm(100, mu.r, sigma)

tmp2 <- tmp*r

x <- tmp2[,1] + mu.x
y <- tmp2[,2] + mu.y
z <- tmp2[,3] + mu.z

out <- optim( c(mean(x),mean(y),mean(z),diff(range(y))/2), tmpfun )
out

If you don't use R then you should still be able to follow the logic and traslate it into another language.

Technically the radius parameter should be bounded by 0, but if the variability is small relative to the true radius then the unbounded method should work fine, or optim has options for doing the bounded optimization, (or you could just do the absolute value of the radius in the function to minimize).

Answer 2

You may be interested by the best fit d-dimensional sphere, i.e. minimizing the variance of the population of the squared distances to the center; it has a simple analytical solution (matrix calculus): see the appendix of the open access paper of Cerisier et al. in J. Comput. Biol. 24(11), 1134-1137 (2017), https://doi.org/10.1089/cmb.2017.0061 It works when the data points are weighted (it works even for continuous distributions; as a by-product, when d=1, a well-known inequality is retrieved: the kurtosis is always greater than the squared skewness plus 1).