2013-04-03 37 views
5

Tôi có một danh sách các điểm 3D, và tôi muốn để phù hợp với một lĩnh vực:làm cách nào để phù hợp với dữ liệu 3D

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

Vì vậy, tôi nghĩ, tôi muốn bày tỏ z, và phù hợp với dữ liệu 2D với 4 tham số (x0, y0, z0 và R):

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

Dưới đây là một mã (nó là một phần của dự án lớn hơn):

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

tôi nhận được lỗi:

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

Đây là liên kết đến data.npy.

Trả lời

7

Bạn cần phải xác định đúng của bạn fitfunc:

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

Tôi không nghĩ rằng cách tiếp cận của bạn là rất mạnh mẽ, bởi vì khi bạn dành sqrt có hai giải pháp, một dương, một tiêu cực, và bạn chỉ xem xét tích cực. Vì vậy, trừ khi tất cả các điểm của bạn nằm ở nửa trên của hình cầu, cách tiếp cận của bạn sẽ không hoạt động. Có thể tốt hơn là tạo r fitfunc của bạn:

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

Tại sao bạn 'coords.T'? – Adobe

+1

@Adobe Bởi vì 'coords' có hình dạng' (100, 3) ', nhưng' coords.T' có hình dạng '(3, 100)', và do đó nó unpacks một hàng vào mỗi biến khi bạn làm 'x, y, z = coords.T'. – Jaime

+0

Ồ, thật tuyệt. Tôi đã làm 'zCoords = Tọa độ [:, 2]' cho điều mô phỏng. – Adobe

Các vấn đề liên quan