Come posso adattare i dati 3D

Question

Ho una lista di punti 3D, e voglio montare la a una sfera:

R^2 = (x-x0)^2 + (y-y0)^2 + (z-z0)^2 

Così ho pensato, mi esprimo z, e in forma di dati 2D con 4 parametri (x0, y0, Z0 e R):

z = sqrt(R^2 - (x-x0)^2 - (y-y0)^2) + z0 

Ecco un codice (si tratta di una parte del più ampio progetto):

#!/usr/bin/python 

from scipy import * 
from scipy.optimize import leastsq 

Coordinates = load("data.npy") 

xyCoords = Coordinates[:, [0, 1]] 
zCoords = Coordinates[:, 2] 

p0 = [149.33499, 148.95999, -218.84893225608857, 285.72893713890107] 

fitfunc = lambda p, x: sqrt(p[3]**2 - (x[0] - p[0])**2 - (x[1] - p[1])**2) + x[2] 
errfunc = lambda p, x, z: fitfunc(p, x) - z 
p1, flag = leastsq(errfunc, p0, args=(xyCoords, zCoords)) 

print p1 

ottengo l'errore:

012.

ValueError: operands could not be broadcast together with shapes (2) (1404) 

Ecco un collegamento a data.npy.

Answer 1

È necessario definire correttamente il vostro fitfunc:

fitfunc = lambda p, x: sqrt(p[3]**2 - (x[:, 0] - p[0])**2 - (x[:, 1] - p[1])**2) + p[2] 

Non credo che il vostro approccio è molto robusto, perché quando si prende il sqrt ci sono due soluzioni, una positiva, una negativa, e si stiamo considerando solo il positivo. Quindi, a meno che tutti i tuoi punti siano nella metà superiore di una sfera, il tuo approccio non funzionerà. Probabilmente è meglio rendere r il tuo fitfunc:

import numpy as np 
from scipy.optimize import leastsq 

# test data: sphere centered at 'center' of radius 'R' 
center = (np.random.rand(3) - 0.5) * 200 
R = np.random.rand(1) * 100 
coords = np.random.rand(100, 3) - 0.5 
coords /= np.sqrt(np.sum(coords**2, axis=1))[:, None] 
coords *= R 
coords += center 

p0 = [0, 0, 0, 1] 

def fitfunc(p, coords): 
    x0, y0, z0, R = p 
    x, y, z = coords.T 
    return np.sqrt((x-x0)**2 + (y-y0)**2 + (z-z0)**2) 

errfunc = lambda p, x: fitfunc(p, x) - p[3] 

p1, flag = leastsq(errfunc, p0, args=(coords,)) 

>>> center 
array([-39.77447344, -69.89096249, 44.46437355]) 
>>> R 
array([ 69.87797469]) 
>>> p1 
array([-39.77447344, -69.89096249, 44.46437355, 69.87797469])