【楕円の近似(Python)】英語質問サイトより，コードを引用したが，中身がわかりません．

Question

前提・実現したいこと

使用言語：python
やりたいこと
x, y, zの3座標から楕円の近似をする．
プログラムコードは英語の質問サイトの方より参考にした．
コードないの「conf.」が何を表しているかわかりません．

わかる方，教えていただきたいです．

よろしくお願いいたします．

該当のソースコード

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)

コード引用：https://stackoverflow.com/questions/58501545/python-fit-3d-ellipsoid-oblate-prolate-to-3d-points

Answer 1

https://stackoverflow.com/questions/58501545/python-fit-3d-ellipsoid-oblate-prolate-to-3d-pointsから引用したとのことですが、

質問者の"The starting data are a set of x, y and z coordinates (cartesian coordinates). "
から始まる質問文を素直に読めば
「与えられた点群にフィットする楕円体を算出するための係数を求めたいんだけど」っていう質問なわけです。

ですから、、confっていうのは、その「所与の点群」の座標に該当するんじゃないでしょうか。

なので

# let us assume some definition of x, y and z

の下に下記のコードを入れてみればよい。Configurationの中身はテキトー。

class Configuration():
    def __init__(self, px ,py, pz):
        self.x = px
        self.y = py
        self.z = pz

conf = Configuration([0,1,1,0.8],[1,0,-2,0.9],[-1,0,1.5,0.4])

出力結果：

coefficients: [-9.      1.5     0.6875  1.75    0.      0.4375  3.75    0.1875 -0.25
 -1.    ]
center: [ 0.1880728  -0.20839408  0.24812539]
axes: [0.27735711 0.65571365 1.04693673]
rotationMatrix: [[-0.99658785  0.08251814 -0.00184948]
 [-0.07994754 -0.97062807 -0.22691307]
 [-0.0205196  -0.22599095  0.97391326]]

Comments

[y08]

理解しました．
あとは自力で頑張って見ます．丁寧な回答，ありがとうございました．