how do I fit 3D data

Question

I have a list of 3D points, and I want to fit the to a sphere:

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

So I thought, I'd express z, and fit 2D data with 4 parameters (x0, y0, z0 and R):

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

Here's a code (it is a part of larger project):

#!/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

I get the error:

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

Here's a link to data.npy.

Answer 1

You need to properly define your fitfunc:

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

I don't think that your approach is very robust, because when you take the sqrt there are two solutions, one positive, one negative, and you are only considering the positive. So unless all your points are on the top half of a sphere, your approach will not work. It is probably better to make r your 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])