Python：将 3D 椭球体 (oblate/prolate) 拟合到 3D 点
Python: fit 3D ellipsoid (oblate/prolate) to 3D points
亲爱的 Whosebug 用户，
我遇到如下问题：我想在我的 python 脚本中将 3D 椭圆体拟合到 3D 数据点。
起始数据是一组 x、y 和 z 坐标（笛卡尔坐标）。我想得到的是3D数据点凸包最佳拟合椭球定义方程中的a和c。
方程在正确旋转和平移的坐标系中是：
所以我最想做的任务是：
求3D数据点的凸包
将best-fit椭球拟合到凸包上得到a和c
您知道 Python 中是否有一些库允许使用最少的代码行来做到这一点吗？还是我必须用我有限的数学知识（在寻找最合适的椭球体时基本上等于零）对这些步骤中的每一个进行明确编码？
在此先感谢您的帮助，祝您有愉快的一天！
好的，我通过将 scipy 的凸包算法与 this website 上的一些 python 函数相结合找到了我的解决方案。
让我们假设您得到一个 x 坐标的 numpy 向量、一个 y 坐标的 numpy 向量和一个 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):
#finds best fit ellipsoid. Found at http://www.juddzone.com/ALGORITHMS/least_squares_3D_ellipsoid.html
#least squares fit to a 3D-ellipsoid
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz = 1
#
# Note that sometimes it is expressed as a solution to
# Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1
# where the last six terms have a factor of 2 in them
# This is in anticipation of forming a matrix with the polynomial coefficients.
# Those terms with factors of 2 are all off diagonal elements. These contribute
# two terms when multiplied out (symmetric) so would need to be divided by two
# change xx from vector of length N to Nx1 matrix so we can use hstack
x = xx[:,np.newaxis]
y = yy[:,np.newaxis]
z = zz[:,np.newaxis]
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz = 1
J = np.hstack((x*x,y*y,z*z,x*y,x*z,y*z, x, y, z))
K = np.ones_like(x) #column of ones
#np.hstack performs a loop over all samples and creates
#a row in J for each x,y,z sample:
# J[ix,0] = x[ix]*x[ix]
# J[ix,1] = y[ix]*y[ix]
# etc.
JT=J.transpose()
JTJ = np.dot(JT,J)
InvJTJ=np.linalg.inv(JTJ);
ABC= np.dot(InvJTJ, np.dot(JT,K))
# Rearrange, move the 1 to the other side
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz - 1 = 0
# or
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz + J = 0
# where J = -1
eansa=np.append(ABC,-1)
return (eansa)
def polyToParams3D(vec,printMe):
#gets 3D parameters of an ellipsoid. Found at http://www.juddzone.com/ALGORITHMS/least_squares_3D_ellipsoid.html
# convert the polynomial form of the 3D-ellipsoid to parameters
# center, axes, and transformation matrix
# vec is the vector whose elements are the polynomial
# coefficients A..J
# returns (center, axes, rotation matrix)
#Algebraic form: X.T * Amat * X --> polynomial form
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)
#See B.Bartoni, Preprint SMU-HEP-10-14 Multi-dimensional Ellipsoidal Fitting
# equation 20 for the following method for finding the center
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)
# Center the ellipsoid at the origin
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]
# print('normed \n',R3test)
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) #inverse is actually the transpose here
if printMe: print('\nRotation matrix\n',inve)
return (center,axes,inve)
#let us assume some definition of x, y and z
#get convex hull
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)
#fit ellipsoid on convex hull
eansa = ls_ellipsoid(hull[0],hull[1],hull[2]) #get ellipsoid polynomial coefficients
print("coefficients:" , eansa)
center,axes,inve = polyToParams3D(eansa,False) #get ellipsoid 3D parameters
print("center:" , center)
print("axes:" , axes)
print("rotationMatrix:", inve)
所以我已经尝试了上面的代码，但我看到的问题是 ls_ellipsoid() 并不总是保证通过其最小二乘法输出椭圆体。因为它只是简单地使用 二次曲面的一般公式 ，所以我经常得到一个方程，它产生一个 sheet 或两个 sheet 的双曲线（或任何其他二次曲面）。有没有办法修改这种最小二乘法，使其始终输出某种形式的椭圆？而不是任何其他二次曲面？
亲爱的 Whosebug 用户，
我遇到如下问题：我想在我的 python 脚本中将 3D 椭圆体拟合到 3D 数据点。
起始数据是一组 x、y 和 z 坐标（笛卡尔坐标）。我想得到的是3D数据点凸包最佳拟合椭球定义方程中的a和c。
方程在正确旋转和平移的坐标系中是：
所以我最想做的任务是：
求3D数据点的凸包
将best-fit椭球拟合到凸包上得到a和c
您知道 Python 中是否有一些库允许使用最少的代码行来做到这一点吗？还是我必须用我有限的数学知识（在寻找最合适的椭球体时基本上等于零）对这些步骤中的每一个进行明确编码？
在此先感谢您的帮助，祝您有愉快的一天！
好的，我通过将 scipy 的凸包算法与 this website 上的一些 python 函数相结合找到了我的解决方案。
让我们假设您得到一个 x 坐标的 numpy 向量、一个 y 坐标的 numpy 向量和一个 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):
#finds best fit ellipsoid. Found at http://www.juddzone.com/ALGORITHMS/least_squares_3D_ellipsoid.html
#least squares fit to a 3D-ellipsoid
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz = 1
#
# Note that sometimes it is expressed as a solution to
# Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1
# where the last six terms have a factor of 2 in them
# This is in anticipation of forming a matrix with the polynomial coefficients.
# Those terms with factors of 2 are all off diagonal elements. These contribute
# two terms when multiplied out (symmetric) so would need to be divided by two
# change xx from vector of length N to Nx1 matrix so we can use hstack
x = xx[:,np.newaxis]
y = yy[:,np.newaxis]
z = zz[:,np.newaxis]
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz = 1
J = np.hstack((x*x,y*y,z*z,x*y,x*z,y*z, x, y, z))
K = np.ones_like(x) #column of ones
#np.hstack performs a loop over all samples and creates
#a row in J for each x,y,z sample:
# J[ix,0] = x[ix]*x[ix]
# J[ix,1] = y[ix]*y[ix]
# etc.
JT=J.transpose()
JTJ = np.dot(JT,J)
InvJTJ=np.linalg.inv(JTJ);
ABC= np.dot(InvJTJ, np.dot(JT,K))
# Rearrange, move the 1 to the other side
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz - 1 = 0
# or
# Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz + J = 0
# where J = -1
eansa=np.append(ABC,-1)
return (eansa)
def polyToParams3D(vec,printMe):
#gets 3D parameters of an ellipsoid. Found at http://www.juddzone.com/ALGORITHMS/least_squares_3D_ellipsoid.html
# convert the polynomial form of the 3D-ellipsoid to parameters
# center, axes, and transformation matrix
# vec is the vector whose elements are the polynomial
# coefficients A..J
# returns (center, axes, rotation matrix)
#Algebraic form: X.T * Amat * X --> polynomial form
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)
#See B.Bartoni, Preprint SMU-HEP-10-14 Multi-dimensional Ellipsoidal Fitting
# equation 20 for the following method for finding the center
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)
# Center the ellipsoid at the origin
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]
# print('normed \n',R3test)
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) #inverse is actually the transpose here
if printMe: print('\nRotation matrix\n',inve)
return (center,axes,inve)
#let us assume some definition of x, y and z
#get convex hull
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)
#fit ellipsoid on convex hull
eansa = ls_ellipsoid(hull[0],hull[1],hull[2]) #get ellipsoid polynomial coefficients
print("coefficients:" , eansa)
center,axes,inve = polyToParams3D(eansa,False) #get ellipsoid 3D parameters
print("center:" , center)
print("axes:" , axes)
print("rotationMatrix:", inve)
所以我已经尝试了上面的代码，但我看到的问题是 ls_ellipsoid() 并不总是保证通过其最小二乘法输出椭圆体。因为它只是简单地使用 二次曲面的一般公式 ，所以我经常得到一个方程，它产生一个 sheet 或两个 sheet 的双曲线（或任何其他二次曲面）。有没有办法修改这种最小二乘法，使其始终输出某种形式的椭圆？而不是任何其他二次曲面？