Estimar el centro y el radio de una esfera a partir de puntos de la superficie

Question

Si suponemos que nuestros puntos de datos fueron muestreados de la superficie de una esfera (con alguna perturbación), ¿cómo podemos recuperar el centro de esa esfera?

En mi búsqueda, encontré documentos sobre algo llamado "regresión esférica", pero no me pareció que hiciera lo mismo. Tal vez no lo entendí.

¿Existe una fórmula directa, similar a la de la regresión lineal, que encuentre un punto central y un radio de la esfera que minimicen la distancia de la suma cuadrada de un conjunto de puntos de datos a la superficie de la esfera?

---

Editar 1:

Podemos suponer que el ruido será 2 o 3 órdenes de magnitud menor que el radio de la esfera y uniformemente esférico gaussiano. Sin embargo, las muestras no se extraerán de manera uniforme de la superficie de la esfera, sino que se agruparán en algunas zonas de la superficie, probablemente en un hemisferio. Una solución que funciona para los datos en $\mathbb R^3$ está bien, pero una solución general para una dimensionalidad arbitraria también es genial.

---

Editar 2:

¿Qué posibilidades tengo de obtener una respuesta sensata si utilizo la regresión lineal? $y = X\beta + \epsilon$ en el espacio de 7 dimensiones pretendiendo que los componentes al cuadrado sean independientes de los otros parámetros:

$\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}$

En el mejor de los casos, supongo que mi métrica de errores será un poco disparatada. En el peor de los casos, la solución no será ni de lejos consistente.
...o eso es una tontería porque con cuatro columnas idénticas, obtenemos una matriz singular cuando intentamos hacer la regresión.

---

Edita 3:

Así que parece que estas son mis opciones:

1. Optimización numérica no lineal utilizando alguna función de coste: $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. Transformación de Hough: discretizar el espacio plausible o los posibles centros y radios alrededor de los puntos de datos. Cada punto emite un voto para los posibles centros de los que podría formar parte en cada discretización de radio específica. El mayor número de votos gana. Esto podría estar bien si hubiera un número potencialmente desconocido de esferas, pero con sólo una es una solución desordenada.
3. Seleccione al azar (o sistemáticamente) grupos de 4 puntos y calcular analíticamente el centro . Rechazar el muestreo si está mal condicionado (los puntos son casi coplanares). Rechazar los valores atípicos y encontrar el centro medio. A partir de ahí podemos encontrar el radio medio.

¿Alguien tiene un método mejor?

Answer 1

Aquí hay algunos R código que muestra un enfoque que utiliza los mínimos cuadrados:

# 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

Si no utiliza R entonces deberías ser capaz de seguir la lógica y trasladarla a otro idioma.

Técnicamente, el parámetro del radio debería estar acotado por 0, pero si la variabilidad es pequeña en relación con el radio verdadero, entonces el método no acotado debería funcionar bien, o la optimización tiene opciones para hacer la optimización acotada, (o podría simplemente hacer el valor absoluto del radio en la función a minimizar).

Answer 2

Puede interesarle la esfera d-dimensional mejor ajustada, es decir, la minimización de la varianza de la población de las distancias al cuadrado al centro; tiene una solución analítica sencilla (cálculo matricial): véase el apéndice del artículo de acceso abierto de Cerisier et al. en J. Comput. Biol. 24(11), 1134-1137 (2017), https://doi.org/10.1089/cmb.2017.0061 Funciona cuando los puntos de datos están ponderados (funciona incluso para distribuciones continuas; como subproducto, cuando d=1, se recupera una desigualdad bien conocida: la curtosis es siempre mayor que la asimetría al cuadrado más 1).