多面体/点集合中最大内接椭圆形的体积

Question

我在3D空间中有一组描述椭球的点（> 1000）。 我想在所有点内拟合一个椭圆体（3D）。这是描述点的最大椭圆体，而在点所描述的多面体内部（所有点都在外部）拟合的椭球）。

[从最小封闭椭球体（又称外边界椭球体）中查看代码，我试图对其进行调整通过修改阈值Port MATLAB bounding ellipsoid code to Python来确定err > tol，根据我的逻辑，给定椭圆方程，所有点都应小于<1。

[我还发现了Moek（https://docs.mosek.com/9.2/pythonfusion/case-studies-ellipsoids.html）上的Loewner-John适应，但是我没有弄清楚如何描述超平面与3D多面体（Ax <= b表示）的交点，因此我可以将其用于3D（不是我有任何保证可以在我的示例中使用）。

来自外部椭球的代码：

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

pi = np.pi
sin = np.sin
cos = np.cos

def mvee(points, tol = 0.001):
    """
    Finds the ellipse equation in "center form"
    (x-c).T * A * (x-c) = 1
    """
    N, d = points.shape
    Q = np.column_stack((points, np.ones(N))).T
    err = tol+1.0
    u = np.ones(N)/N
    while err > tol:
        # assert u.sum() == 1 # invariant
        X = np.dot(np.dot(Q, np.diag(u)), Q.T)
        M = np.diag(np.dot(np.dot(Q.T, la.inv(X)), Q))
        jdx = np.argmax(M)
        step_size = (M[jdx]-d-1.0)/((d+1)*(M[jdx]-1.0))
        new_u = (1-step_size)*u
        new_u[jdx] += step_size
        err = la.norm(new_u-u)
        u = new_u
    c = np.dot(u,points)        
    A = la.inv(np.dot(np.dot(points.T, np.diag(u)), points)
               - np.multiply.outer(c,c))/d
    return A, c

#some random points
points = np.array([[ 0.53135758, -0.25818091, -0.32382715], 
                   [ 0.58368177, -0.3286576,  -0.23854156,], 
                   [ 0.18741533,  0.03066228, -0.94294771], 
                   [ 0.65685862, -0.09220681, -0.60347573],
                   [ 0.63137604, -0.22978685, -0.27479238],
                   [ 0.59683195, -0.15111101, -0.40536606],
                   [ 0.68646128,  0.0046802,  -0.68407367],
                   [ 0.62311759,  0.0101013,  -0.75863324]])

# Singular matrix error!
# points = np.eye(3)

A, centroid = mvee(points)    
U, D, V = la.svd(A)    
rx, ry, rz = 1./np.sqrt(D)
u, v = np.mgrid[0:2*pi:20j, -pi/2:pi/2:10j]

def ellipse(u,v):
    x = rx*cos(u)*cos(v)
    y = ry*sin(u)*cos(v)
    z = rz*sin(v)
    return x,y,z

E = np.dstack(ellipse(u,v))
E = np.dot(E,V) + centroid
x, y, z = np.rollaxis(E, axis = -1)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(x, y, z, cstride = 1, rstride = 1, alpha = 0.05)
ax.scatter(points[:,0],points[:,1],points[:,2])

plt.show()

哪个给我这个结果：

问题：如何使用Python将椭圆体（3D）放入3D点云中？

是否可以修改外椭球的算法以获得最大的内接（内）椭球？

我理想地在Python中寻找代码。

Answer 1

表示椭圆的（表面）的一种可能是数学上标准的方法是它是集合

{ X | (X-a)'*inv(C)*(X-a) = 1}
the solid ellipsoid is then 
{ X | (X-a)'*inv(C)*(X-a) <= 1}

这里C是一个3x3对称正定矩阵，a是椭球的“中心”。

我们可以通过使用cholesky分解使它更容易处理，即找到一个较低的三角矩阵L，这样

C = L*L'

并且使用L的倒数M（L是较低的三角形，这很容易计算）。我们比实心椭球是

{ X | (M*(X-a))'*(M*(X-a)) <= 1}
= { | ||M*(X-a))|| <= 1} where ||.|| is the euclidean 

规范

我们有一堆点X []和一个包含它们的椭球（C，a），即

for all i ||M*(X[i]-a)|| <= 1
i.e. for all i ||Y[i]|| <= 1 where Y[i] = M*(X[i]-a)

现在，我们要变换椭球体（即更改C和a），以便所有点都在变换的椭球体之外。我们也可以变换M和a。

最简单的方法就是将M缩放为常数s，然后不做任何处理。这相当于缩放所有Y []，在这种情况下，很容易看出要使用的缩放因子将是|| Y [i] ||最小值中的一个。这样，所有点都将在变换的椭球体的外部或上方，并且某些点将在其上，因此变换的椭圆形将尽可能大。

就D而言，新的椭圆就是

D = (1/(s*s))*C

如果这种简单的方法给出可接受的结果，那就是我要使用的。

[我认为，不移动中心，最普遍的事情就是改变

M to N*M

约束是N是上三角并且对角线上有正数。我们需要N的那个

N*Y[i] >= 1 for all i

我们需要选择N的标准。一个应该是尽可能减少体积，即行列式（对于较低的三角形矩阵，它只是对角线元素的乘积）应受约束而尽可能小。

[也许有很多软件包可以做这种事情，但是我不知道哪个软件包（应该把它更多地看作是我的无知，而不是表明没有这样的软件包）。] >

一旦找到N，则变换后的C矩阵为

D = L*inv(N)*inv(N')*L'
您也可以更改一个。我把细节留给感兴趣的读者...