CODEWITHSUNDEEP
numpyscipy
Adobe
Adobe

Reputation: 13467

how do I fit 3D data

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.

Upvotes: 5

Views: 3564

Answers (1)

Jaime
Jaime

Reputation: 67417

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])

Upvotes: 7

Related Questions