Python:将3D椭圆体(扁圆/扁长)拟合到3D点 [英] Python: fit 3D ellipsoid (oblate/prolate) to 3D points

问题描述

尊敬的stackoverflow用户，

Dear fellow stackoverflow users,

我面临如下问题:我想在我的python脚本中将3D椭球体适合3D数据点.

I face a problem as follows: I would like to fit a 3D ellipsoid to 3D data points within my python script.

起始数据是一组x，y和z坐标(笛卡尔坐标).我想得到的是3D数据点的凸包的最佳拟合椭球的定义方程中的a和c.

The starting data are a set of x, y and z coordinates (cartesian coordinates). What I would like to get are a and c in the defining equation of the best-fit ellipsoid of the convex hull of the 3D data points.

等式在正确旋转和平移的坐标系中:

The equation is, in the properly rotated and translated coordinate system:

所以我理想地要做的任务是:

So the tasks I would ideally like to do are:

1. 查找3D数据点的凸包

1. Find convex hull of 3D data points

将最佳椭球体拟合到凸包，并得到a和c

Fit best-fit ellipsoid to the convex hull and get a and c

您知道是否有一些库允许使用最少的代码行在Python中执行此操作?还是我必须用有限的数学知识(在寻找最合适的椭球体时实际上为零)显式地编码所有这些步骤?

Do you know whether there is some library allowing to do this in Python with minimal lines of code? Or do I have to explicitly code every of these steps with my limited math knowledge (which essentially amounts to zero when it comes to find best fit ellipsoid)?

提前感谢您的帮助，祝您愉快！

Thanks in advance for your help and have a nice day!

推荐答案

好吧，我找到了我的解决方案，方法是将scipy的凸包算法与

All right, I found my solution by combining the convex hull algorithm of scipy with some python function found on this website.

让我们假设您得到一个x坐标的numpy向量，一个y坐标的numpy向量以及一个z坐标的numpy向量，分别命名为x，y和z.这对我有用:

Let us assume that you get a numpy vector of x coordinates, a numpy vector of y coordinates, and a numpy vector of z coordinates, named x, y and z. This worked for me:

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)

这篇关于Python:将3D椭圆体(扁圆/扁长)拟合到3D点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！