好的,我通过将Scipy的凸包算法与此网站上找到的一些Python函数相结合,找到了我的解决方案。
假设您得到了名为x、y和z的numpy向量的X坐标、Y坐标和Z坐标。这对我有效:
from scipy.spatial
import ConvexHull, convex_hull_plot_2d
import numpy as np
from numpy.linalg import eig, inv
def ls_ellipsoid(xx,yy,zz):
x = xx[:,np.newaxis]
y = yy[:,np.newaxis]
z = zz[:,np.newaxis]
J = np.hstack((x*x,y*y,z*z,x*y,x*z,y*z, x, y, z))
K = np.ones_like(x)
JT=J.transpose()
JTJ = np.dot(JT,J)
InvJTJ=np.linalg.inv(JTJ);
ABC= np.dot(InvJTJ, np.dot(JT,K))
eansa=np.append(ABC,-1)
return (eansa)
def polyToParams3D(vec,printMe):
if printMe: print('\npolynomial\n',vec)
Amat=np.array(
[
[ vec[0], vec[3]/2.0, vec[4]/2.0, vec[6]/2.0 ],
[ vec[3]/2.0, vec[1], vec[5]/2.0, vec[7]/2.0 ],
[ vec[4]/2.0, vec[5]/2.0, vec[2], vec[8]/2.0 ],
[ vec[6]/2.0, vec[7]/2.0, vec[8]/2.0, vec[9] ]
])
if printMe: print('\nAlgebraic form of polynomial\n',Amat)
A3=Amat[0:3,0:3]
A3inv=inv(A3)
ofs=vec[6:9]/2.0
center=-np.dot(A3inv,ofs)
if printMe: print('\nCenter at:',center)
Tofs=np.eye(4)
Tofs[3,0:3]=center
R = np.dot(Tofs,np.dot(Amat,Tofs.T))
if printMe: print('\nAlgebraic form translated to center\n',R,'\n')
R3=R[0:3,0:3]
R3test=R3/R3[0,0]
s1=-R[3, 3]
R3S=R3/s1
(el,ec)=eig(R3S)
recip=1.0/np.abs(el)
axes=np.sqrt(recip)
if printMe: print('\nAxes are\n',axes ,'\n')
inve=inv(ec)
if printMe: print('\nRotation matrix\n',inve)
return (center,axes,inve)
surface = np.stack((conf.x,conf.y,conf.z), axis=-1)
hullV = ConvexHull(surface)
lH = len(hullV.vertices)
hull = np.zeros((lH,3))
for i in range(len(hullV.vertices)):
hull[i] = surface[hullV.vertices[i]]
hull = np.transpose(hull)
eansa = ls_ellipsoid(hull[0],hull[1],hull[2])
print("coefficients:" , eansa)
center,axes,inve = polyToParams3D(eansa,False)
print("center:" , center)
print("axes:" , axes)
print("rotationMatrix:", inve)