8
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Notice that the two forms of your question are not equivalent: it is not necessarily the case that minimizing the sum of squares of distances from the surface gives the best estimates unless a strong assumption is made about the nature of the perturbations. It would therefore help to know more about how the perturbations occur (and how large they may be compared to the sphere's size). Also: in how many dimensions is your sphere? $\endgroup$ – whuber Jan 5 '12 at 15:49
  • $\begingroup$ @whuber I intended to define the best fit as the one that minimizes sum-squared distance of the data from the nearest point on the sphere's surface. I didn't think much about the assumptions that entails. I expect proportionally small errors; so perhaps the exact metric doesn't matter too much, although I would like to know what the function is minimizing. I've added more information about the noise to the question. $\endgroup$ – JCooper Jan 5 '12 at 16:18
  • $\begingroup$ @Max I did see that. But it's a site for a black-box commercial product. It's the actual formula I was interested in. It's starting to look like there's no closed-form solution and I'll have to use a numerical approach instead (which is what I assume the nlReg software is also doing). $\endgroup$ – JCooper Jan 5 '12 at 17:15
  • $\begingroup$ it seems like this could be a straight minimization problem with a nonlinear objective function (the one you mentioned above). if the errors are assumed to be gaussian, you need only compute the distributional parameters of the errors once you've found the center of the sphere minimizing the objective function. edit: i left the page open too long and didn't see your comment. we have the same idea. $\endgroup$ – assumednormal Jan 5 '12 at 17:17
  • 2
    $\begingroup$ re Edit 3: Given $(x_0, y_0, z_0)$, $r$ is easy to find. To obtain $(x_0, y_0, z_0)$, Newton's Method should converge rapidly from some reasonable starting value obtained as in (3). $\endgroup$ – whuber Jan 5 '12 at 20:32
2
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.