5

我有一个 3D 点列表,我想将其拟合到一个球体:

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

所以我想,我会表达 z,并用 4 个参数(x0、y0、z0 和 R)拟合 2D 数据:

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

这是一个代码(它是更大项目的一部分):

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

我得到错误:

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

这是data.npy的链接。

4

1 回答 1

7

您需要正确定义您的fitfunc

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

我不认为您的方法非常稳健,因为当您采用时,sqrt有两种解决方案,一种是正面的,一种是负面的,您只考虑正面的。因此,除非您的所有点都在球体的上半部分,否则您的方法将行不通。让r你的 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])
于 2013-04-03T11:38:01.110 回答