[SciPy-user] Fitting sphere to 3d data points

[David Goldsmith]

James Vincent wrote:
> Thanks for all the input. I _think_ I've got it. (emphasis added)
Have you tested it w/ points known to lie on a sphere, e.g., (+/-1, 0, 
0), (0,+/-1,0), (0,0,+/-1), which should of course yield a=b=c=0, r=1 
with zero residual error?  One reason why I ask is that, when fitting 
data to a polynomial model (a sphere is the locus of the solution set of 
the polynomial A(x^2 + y^2 + z^2) + Bx + Cy + Dz -1 = 0) I usually like 
to take advantage of the fact that a polynomial is a linear function of 
the "monomial" pieces, in this case x^2 + y^2 + z^2, x, y, & z.  Since 
in general one can't solve a polynomial for what one wants 
formulaically, this is the only general way to linear least squares fit 
data to a polynomial; of course in your quadratic case, one can solve 
for any one of the data variables in terms of the others (as you have 
done), but to do so you have to introduce the square-root function, 
which numerically speaking, isn't quite as exact as simply squaring.  In 
other words, I would respectfully suggest that you fit the four 
quantities (x^2 + y^2 + z^2), x, y, & z to the model A(x^2 + y^2 + z^2) 
+ Bx + Cy + Dz -1 = 0 and calculate your a, b, c, and r from the A, B, 
C, & D.  (And even if you choose not to do that, at least test your 
algorithm with some pseudo-data chosen to yield known exact results.)

DG

PS: I've had to fit data to a circle before (using IDL); if anyone else 
thinks they might have to do something like this in the future, I could 
contribute "under-test" polyFit, circleFit, and sphereFit functions, not 
immediately 'cause I'm under deadline, but in the not too distant future.
> This works:
>
>
> def resSphere(p,x,y,z):
>    """ residuals from sphere fit """
>   
>    a,b,c,r = p                             # a,b,c are center x,y,c 
> coords to be fit, r is the radius to be fit
>    distance = sqrt( (x-a)**2 + (y-b)**2 + (z-c)**2 )
>    err = distance - r                 # err is distance from input 
> point to current fitted surface 
>
>    return err
>
> params = [0.,0.,0.,0.]
> myResult = leastsq(resSphere, params, args=(myX,myY,myZ) )
> print myResult[0]
>
>
>
>
> On Jan 25, 2007, at 10:47 AM, David Douard wrote:
>
>> On Thu, Jan 25, 2007 at 09:31:10AM -0500, David Huard wrote:
>>> Hi James,
>>>
>>> As a first guess, I'd say the center of the sphere is simply the mean of
>>> your data points, if they're all weighted equally. 
>>
>> Hello,
>>
>> I would have rather said that you need to find the point that minimize
>> the distance to the normals of the triangles you have from your data
>> points (not sure this is really meaningful...).
>> If the points really are on a sphere, all the normal will cut on one
>> point. If not, there really is a minimization problema to solve.
>>
>> David
>>
>>
>>
>>> With only one parameter
>>> left to fit, it should be easy enough. However, you may want to look 
>>> at the
>>> paper:
>>>
>>> Werman, Michael and Keren, Daniel
>>> A Bayesian method for fitting parametric and nonparametric models to 
>>> noisy
>>> data
>>> Ieee Transactions on Pattern Analysis and Machine Intelligence, 23, 
>>> 2001.
>>>
>>> They write that the Mean Square Error approach overestimates the 
>>> radius in
>>> the case of circles. They don't talk about the 3D case, but I'd guess
>>> similar problems arise. They provide a method to fit parametric 
>>> shapes with
>>> some robustness to data errors.
>>>
>>> Cheers,
>>>
>>> David
>>>
>>>
>>>
>>> 2007/1/25, James Vincent <jjv5 at nih.gov <mailto:jjv5 at nih.gov>>:
>>>>
>>>> Hello,
>>>> Is it possible to fit a sphere to 3D data points using
>>>> scipy.optimize.leastsq? I'd like to minimize the residual for the 
>>>> distance
>>>> from the actual x,y,z point and the fitted sphere surface. I can 
>>>> see how to
>>>> minimize for z, but that's not really what I'm looking for. Is there a
>>>> better way to do this? Thanks for any help.
>>>>
>>>> params = a,b,c and r
>>>> a,b,c are the fitted center point of the sphere, r is the radius
>>>>
>>>> err = distance-to-center - radius
>>>> err = sqrt( x-a)**2 + (y-b)**2 + (z-c)**2) - r
>>>>
>>>>
>>>>
>>>> ----
>>>> James J. Vincent, Ph.D.
>>>> National Cancer Institute
>>>> National Institutes of Health
>>>> Laboratory of Molecular Biology
>>>> Building 37, Room 5120
>>>> 37 Convent Drive, MSC 4264
>>>> Bethesda, MD 20892 USA
>>>>
>>>> 301-451-8755
>>>> jjv5 at nih.gov <mailto:jjv5 at nih.gov>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> SciPy-user mailing list
>>>> SciPy-user at scipy.org <mailto:SciPy-user at scipy.org>
>>>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>>>
>>>>
>>>>
>>
>>> _______________________________________________
>>> SciPy-user mailing list
>>> SciPy-user at scipy.org <mailto:SciPy-user at scipy.org>
>>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>
>>
>> -- 
>> David Douard                             LOGILAB, Paris (France)
>> Formations Python, Zope, Plone, Debian : http://www.logilab.fr/formations
>> Développement logiciel sur mesure :      http://www.logilab.fr/services
>> Informatique scientifique :              http://www.logilab.fr/science
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user at scipy.org <mailto:SciPy-user at scipy.org>
>> http://projects.scipy.org/mailman/listinfo/scipy-user
>
> ----
> James J. Vincent, Ph.D.
> National Cancer Institute
> National Institutes of Health
> Laboratory of Molecular Biology
> Building 37, Room 5120
> 37 Convent Drive, MSC 4264
> Bethesda, MD 20892 USA
>
> 301-451-8755
> jjv5 at nih.gov <mailto:jjv5 at nih.gov>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>