Estimer le centre et le rayon d'une sphère à partir de points sur la surface

Si nous supposons que nos points de données ont été échantillonnés à partir de la surface d'une sphère (avec une certaine perturbation), comment pouvons-nous récupérer le centre de cette sphère?

Dans ma recherche, j'ai trouvé des articles sur quelque chose qui s'appelait "régression sphérique", mais il ne semblait pas que cela faisait la même chose. Peut-être que je ne l'ai tout simplement pas compris.

Existe-t-il une formule simple, similaire à la régression linéaire, qui trouve un point central et un rayon de sphère qui minimise la distance au carré d'un ensemble de points de données de la surface de la sphère?

Modifier 1:

On peut supposer que le bruit sera de 2 ou 3 ordres de grandeur plus petit que le rayon de la sphère et uniformément sphérique gaussien. Cependant, les échantillons eux-mêmes ne seront certainement pas tirés uniformément de la surface de la sphère, mais seront probablement regroupés en quelques patchs sur la surface, probablement tous dans un hémisphère. Une solution qui fonctionne pour les données dans est très bien, mais une solution générale pour la dimensionnalité arbitraire est également excellente.$\mathbb R^3$

Modifier 2:

Quelles sont les chances que j'obtienne une réponse raisonnable si je devais utiliser une régression linéaire, , dans l'espace à 7 dimensions en faisant comme si les composantes au carré étaient indépendantes des autres paramètres:$y = X\beta + \epsilon$

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

Au mieux, je suppose que ma métrique d'erreur sera un peu farfelue. Au pire, la solution ne sera même pas proche de la cohérence.
... ou c'est idiot car avec quatre colonnes identiques, nous obtenons une matrice singulière lorsque nous essayons de faire une régression.

Modifier 3:

Il semble donc que ce sont mes options:

1. Optimisation numérique non linéaire utilisant une fonction de coût:$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: discrétise l'espace plausible ou les centres et rayons possibles autour des points de données. Chaque point jette un vote pour les centres potentiels dont il pourrait faire partie à chaque discrétisation de rayon spécifique. La plupart des votes l'emportent. Cela pourrait être correct s'il y avait potentiellement un nombre inconnu de sphères, mais avec une seule, c'est une solution compliquée.
3. Sélectionnez au hasard (ou systématiquement) des groupes de 4 points et calculez analytiquement le centre . Rejetez l'échantillonnage s'il est mal conditionné (les points sont presque coplanaires). Rejeter les valeurs aberrantes et trouver le centre moyen. De là, nous pouvons trouver le rayon moyen.

Quelqu'un at-il une meilleure méthode?

Voici un Rcode qui montre une approche utilisant les moindres carrés:

# 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 vous ne l'utilisez pas, Rvous devriez toujours pouvoir suivre la logique et la traduire dans une autre langue.

Techniquement, le paramètre de rayon doit être limité par 0, mais si la variabilité est petite par rapport au vrai rayon, alors la méthode illimitée devrait fonctionner correctement, ou optim a des options pour faire l'optimisation limitée, (ou vous pouvez simplement faire la valeur absolue de la rayon dans la fonction pour minimiser).