I use the function leastsq from scipy.optimize to fit sphere coordinates and radius from 3D coordinates.
So my code looks like this :
def distance(pc,point):
xc,yc,zc,rd = pc
x ,y ,z = point
return np.sqrt((xc-x)**2+(yc-y)**2+(zc-z)**2)
def sphere_params(coords):
from scipy import optimize
err = lambda pc,point : distance(pc,point) - pc[3]
pc = [0, 0, 0, 1]
pc, success = optimize.leastsq(err, pc[:], args=(coords,))
return pc
(Built thanks to : how do I fit 3D data.)
I started working with the variable coords as a list of tuples (each tuple being an x,y,z coordinate):
>> coords
>> [(0,0,0),(0,0,1),(-1,0,0),(0.57,0.57,0.57),...,(1,0,0),(0,1,0)]
Which lead me to an error :
>> pc = sphere_params(coords)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/michel/anaconda/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 374, in leastsq
raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
TypeError: Improper input: N=4 must not exceed M=3
Where N is the number of parameters stored in pc, and M the number of data points. Which makes it look like I haven't given enough data points while my list coords actually regroups 351 tuples versus 4 parameters in pc !
From what I read in minipack the actual culprit seems to be this line (from _check_func()) :
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
Unless i'm mistaken, in my case it translates into
res = atleast_1d(distance(*(pc[:len(pc)],) + args)
But I'm having a terrible time trying to understand what this mean alongs with the rest of the _check_func() function.
I ended up changing coords into an array before giving it as an argument to sphere_param() : coords = np.asarray(coords).T and it started working just fine. I would really like to understand why the data format was giving me trouble though !
In advance, many thanks for your answers!
EDIT : I notice my use of coords for the "distance" and "err" functions was really unwise and misleading, it wasn't so in my original code so it was not the core of the problem. Now make more sense.
Your err function must take the full list of coords and return a full list of distances. leastsq will then take the list of errors, square and sum them, and minimize that squared sum.
There are also distance functions in scipy.spatial.distance, so I would recommend that:
from scipy.spatial.distance import cdist
from scipy.optimize import leastsq
def distance_cdist(pc, coords):
return cdist([pc], coords).squeeze()
def distance_norm(pc, points):
""" pc must be shape (D+1,) array
points can be (N, D) or (D,) array """
c = np.asarray(pc[:3])
points = np.atleast_2d(points)
return np.linalg.norm(points-c, axis=1)
def sphere_params(coords):
err = lambda pc, coords: distance(pc[:3], coords) - pc[3]
pc = [0, 0, 0, 1]
pc, success = leastsq(err, pc, args=(coords,))
return pc
coords = [(0,0,0),(0,0,1),(-1,0,0),(0.57,0.57,0.57),(1,0,0),(0,1,0)]
sphere_params(coords)