python - 多面体/点集中的最大内接椭圆体

Question

稍后编辑:我上传了here我的原始数据的样本。它实际上是DICOM格式的分割图像。该结构的体积约为16 mL，因此我假设内部椭球体的体积应小于该体积。从DICOM图像中提取点，我使用了以下代码:

import os
import numpy as np
import SimpleITK as sitk


def get_volume_ml(image):
    x_spacing, y_spacing, z_spacing = image.GetSpacing()
    image_nda = sitk.GetArrayFromImage(image)
    imageSegm_nda_NonZero = image_nda.nonzero()
    num_voxels = len(list(zip(imageSegm_nda_NonZero[0],
                              imageSegm_nda_NonZero[1],
                              imageSegm_nda_NonZero[2])))
    if 0 >= num_voxels:
        print('The mask image does not seem to contain an object.')
        return None
    volume_object_ml = (num_voxels * x_spacing * y_spacing * z_spacing) / 1000
    return volume_object_ml


def get_surface_points(folder_path):
    """
    :param folder_path: path to folder where DICOM images are stored
    :return: surface points of the DICOM object
    """
    # DICOM Series
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(os.path.normpath(folder_path))
    reader.SetFileNames(dicom_names)
    reader.MetaDataDictionaryArrayUpdateOn()
    reader.LoadPrivateTagsOn()
    try:
        dcm_img = reader.Execute()
    except Exception:
        print('Non-readable DICOM Data: ', folder_path)
        return None
    volume_obj = get_volume_ml(dcm_img)
    print('The volume of the object in mL:', volume_obj)
    contour = sitk.LabelContour(dcm_img, fullyConnected=False)
    contours = sitk.GetArrayFromImage(contour)
    vertices_locations = contours.nonzero()

    vertices_unravel = list(zip(vertices_locations[0], vertices_locations[1], vertices_locations[2]))
    vertices_list = [list(vertices_unravel[i]) for i in range(0, len(vertices_unravel))]
    surface_points = np.array(vertices_list)

    return surface_points

folder_path = r"C:\Users\etc\TTT [13]\20160415 114441\Series 052 [CT - Abdomen WT 1 0 I31f 3]"
points = get_surface_points(folder_path)

我在3D空间中有一组点(n> 1000)，它们描述了空心卵形形状。 我想要的是将所有点内的椭球(3D)拟合。 我正在寻找点内最大体积的椭球拟合。

我试图改写 Minimum Enclosing Ellipsoid(又称外边界椭圆形)中的代码
通过修改阈值 err > tol，从我的逻辑开始，给定椭圆方程，所有点均应小于<1。但是没有成功。

我还尝试了 mosek上的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 plot_ellipsoid(A, centroid, color, ax):
"""

:param A: matrix
:param centroid: center
:param color: color
:param ax: axis
:return:
"""
centroid = np.asarray(centroid)
A = np.asarray(A)
U, D, V = la.svd(A)
rx, ry, rz = 1. / np.sqrt(D)
u, v = np.mgrid[0:2 * np.pi:20j, -np.pi / 2:np.pi / 2:10j]
x = rx * np.cos(u) * np.cos(v)
y = ry * np.sin(u) * np.cos(v)
z = rz * np.sin(v)
E = np.dstack((x, y, z))
E = np.dot(E, V) + centroid
x, y, z = np.rollaxis(E, axis=-1)
ax.plot_wireframe(x, y, z, cstride=1, rstride=1, color=color, alpha=0.2)
ax.set_zlabel('Z-Axis')
ax.set_ylabel('Y-Axis')
ax.set_xlabel('X-Axis')

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

folder_path = r"" # path to a DICOM img folder
points = get_surface_points(folder_path) # or some random pts 

A, centroid = mvee(points)    
U, D, V = la.svd(A)    
rx_outer, ry_outer, rz_outer = 1./np.sqrt(D)
# PLOT
fig = plt.figure()
ax1 = fig.add_subplot(111, projection='3d')
ax1.scatter(points[:, 0], points[:, 1], points[:, 2], c='blue')
plot_ellipsoid(A, centroid, 'green', ax1)

这为我的采样点上的外部椭球提供了以下结果:


主要问题:如何使用Python将椭圆体(3D)放入3D点云中？

是否可以修改外部椭球的算法以获得最大的内接(内部)椭球？

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

Answer 1

问题陈述
给定许多点v₁, v₂, ..., vₙ，找到一个满足两个约束的大椭圆体:

椭球在凸包ℋ= ConvexHull(v₁，v²，...，vₙ)中。

点v₁，v²，...，v ...都不在椭球内。

我提出了一个迭代过程，以找到满足这两个约束的大椭圆体。在每次迭代中，我们需要解决一个半定编程问题。该迭代过程可以保证收敛，但是不能保证收敛到全局最大的椭球。
方法
查找单个椭球
迭代过程的核心是在每次迭代中，我们找到一个满足3个条件的椭球:

椭球包含在ConvexHull(v₁，v²，...，vₙ)= {x | Ax <= b}。

对于一组点u₁，...uₘ(其中{v₁，v²，...，vₙ}⊂{u₁，...uₘ}，即点云中的给定点属于该组点u₁ ，...uₘ)，则椭球在u₁，...uₘ中不包含任何点。我们称这个集合u₁，...uₘ为“外点”。

对于点w₁，...，wₖ的集合(其中{w₁，...，wₖ}∩{v v，v²，...，vₙ} =∅，即v₁，v²，...中没有点。 ..，vₙ属于{w₁，...，wₖ})，椭球包含所有点w₁，...，wₖ。我们称此集合w₁，...，wₖ为“内点”。

直观的想法是“内点”w₁，...，wₖ表示椭球的体积。我们将新点添加到“内部点”上，以增加椭球体的体积。
为了通过凸优化找到这样的椭球，我们将椭球参数化为

{x | xᵀPx + 2qᵀx  ≤ r}

我们将搜索 P, q, r。
将“外部点”u₁，...uₘ都在椭球体之外的条件表示为

uᵢᵀPuᵢ + 2qᵀuᵢ >= r ∀ i=1, ..., m

这是对 P, q, r的线性约束。
将“内点”w₁，...，wₖ都在椭球内的条件表示为

wᵢᵀPwᵢ + 2qᵀwᵢ <= r ∀ i=1, ..., k

这也是对 P, q, r的线性约束。
我们还施加了约束

P is positive definite

P是正定的，加上存在满足wᵢᵀPwᵢ+2qᵀwᵢ<= r的点wᵢ的约束，保证了集合{x | xᵀPx+2qᵀx≤r}是一个椭球。
我们还具有这样的约束:椭球在凸包inside = {x |内。 aᵢᵀx≤bᵢ，i = 1，...，l}(即，存在 l半空间作为H的H表示)。使用 s-lemma，我们知道包含椭球的半空间 {x|aᵢᵀx≤ bᵢ}的充要条件是

∃ λᵢ >= 0,
s.t [P            q -λᵢaᵢ/2]  is positive semidefinite.
    [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]

因此，我们可以解决以下半定编程问题，以找到包含所有“内部点”，不包含任何“外部点”并且在凸包内的椭圆体ℋ

find P, q, r, λ
s.t uᵢᵀPuᵢ + 2qᵀuᵢ >= r ∀ i=1, ..., m
    wᵢᵀPwᵢ + 2qᵀwᵢ <= r ∀ i=1, ..., k
    P is positive definite.
    λ >= 0,
    [P            q -λᵢaᵢ/2]  is positive semidefinite.
    [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]

我们称这个为 P, q, r = find_ellipsoid(outside_points, inside_points, A, b)。
该椭圆体的体积与(r +qᵀP⁻q)/功率(det(P)，1/3)成正比。
迭代过程。
我们将“外部点”初始化为点云中的所有点v 1，v 2，...，v 1，并将“内部点”初始化为凸包in中的单个点 w₁。在每次迭代中，我们在上一小节中使用 find_ellipsoid函数在within中查找包含所有“内部点”但不包含任何“外部点”的椭球。根据 find_ellipsoid中SDP的结果，我们执行以下操作

如果SDP可行。然后，我们将新发现的椭球与迄今发现的最大椭球进行比较。如果新的椭球较大，则将其接受为迄今为止找到的最大椭球。

如果SDP不可行，则我们删除“内部点”中的最后一个点，将此点添加到“外部点”中。

在这两种情况下，我们然后在凸包take中获取一个新的采样点，将该采样点添加到“内部点”，然后再次求解SDP。
完整的算法如下

将“外部点”初始化为v₁，v2，...，vₙ，将“内部点”初始化为凸包ℋ中的单个随机点。

迭代 :

解决SDP P, q, r = find_ellipsoid(outside_points, inside_points, A, b)。

如果SDP可行并且volume(Ellipsoid(P，q，r))> maximum_volume，则设置P_best = P, q_best=q, r_best = r。

如果SDP不可行，则pt = inside_points.pop_last()，outside_points.push_back(pt)。

在Random中随机采样一个新点，将该点附加到“内部点”，迭代+ =1。转到步骤3。

代码

from scipy.spatial import ConvexHull, Delaunay
import scipy
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import dirichlet
from mpl_toolkits.mplot3d import Axes3D  # noqa


def get_hull(pts):
    dim = pts.shape[1]
    hull = ConvexHull(pts)
    A = hull.equations[:, 0:dim]
    b = hull.equations[:, dim]
    return A, -b, hull


def compute_ellipsoid_volume(P, q, r):
    """
    The volume of the ellipsoid xᵀPx + 2qᵀx ≤ r is proportional to
    r + qᵀP⁻¹q / power(det(P), 1/dim)
    We return this number.
    """
    return (r + q @ np.linalg.solve(P, q)) / \
        np.power(np.linalg.det(P), 1. / P.shape[0])


def uniform_sample_from_convex_hull(deln, dim, n):
    """
    Uniformly sample n points in the convex hull Ax<=b
    This is copied from
    https://stackoverflow.com/questions/59073952/how-to-get-uniformly-distributed-points-in-convex-hull
    @param deln Delaunay of the convex hull.
    """
    vols = np.abs(np.linalg.det(deln[:, :dim, :] - deln[:, dim:, :]))\
        / np.math.factorial(dim)
    sample = np.random.choice(len(vols), size=n, p=vols / vols.sum())

    return np.einsum('ijk, ij -> ik', deln[sample],
                     dirichlet.rvs([1]*(dim + 1), size=n))


def centered_sample_from_convex_hull(pts):
    """
    Sample a random point z that is in the convex hull of the points
    v₁, ..., vₙ. z = (w₁v₁ + ... + wₙvₙ) / (w₁ + ... + wₙ) where wᵢ are all
    uniformly sampled from [0, 1]. Notice that by central limit theorem, the
    distribution of this sample is centered around the convex hull center, and
    also with small variance when the number of points are large.
    """
    num_pts = pts.shape[0]
    pts_weights = np.random.uniform(0, 1, num_pts)
    z = (pts_weights @ pts) / np.sum(pts_weights)
    return z


def find_ellipsoid(outside_pts, inside_pts, A, b):
    """
    For a given sets of points v₁, ..., vₙ, find the ellipsoid satisfying
    three constraints:
    1. The ellipsoid is within the convex hull of these points.
    2. The ellipsoid doesn't contain any of the points.
    3. The ellipsoid contains all the points in @p inside_pts
    This ellipsoid is parameterized as {x | xᵀPx + 2qᵀx ≤ r }.
    We find this ellipsoid by solving a semidefinite programming problem.
    @param outside_pts outside_pts[i, :] is the i'th point vᵢ. The point vᵢ
    must be outside of the ellipsoid.
    @param inside_pts inside_pts[i, :] is the i'th point that must be inside
    the ellipsoid.
    @param A, b The convex hull of v₁, ..., vₙ is Ax<=b
    @return (P, q, r, λ) P, q, r are the parameterization of this ellipsoid. λ
    is the slack variable used in constraining the ellipsoid inside the convex
    hull Ax <= b. If the problem is infeasible, then returns
    None, None, None, None
    """
    assert(isinstance(outside_pts, np.ndarray))
    (num_outside_pts, dim) = outside_pts.shape
    assert(isinstance(inside_pts, np.ndarray))
    assert(inside_pts.shape[1] == dim)
    num_inside_pts = inside_pts.shape[0]

    constraints = []
    P = cp.Variable((dim, dim), symmetric=True)
    q = cp.Variable(dim)
    r = cp.Variable()

    # Impose the constraint that v₁, ..., vₙ are all outside of the ellipsoid.
    for i in range(num_outside_pts):
        constraints.append(
            outside_pts[i, :] @ (P @ outside_pts[i, :]) +
            2 * q @ outside_pts[i, :] >= r)
    # P is strictly positive definite.
    epsilon = 1e-6
    constraints.append(P - epsilon * np.eye(dim) >> 0)

    # Add the constraint that the ellipsoid contains @p inside_pts.
    for i in range(num_inside_pts):
        constraints.append(
            inside_pts[i, :] @ (P @ inside_pts[i, :]) +
            2 * q @ inside_pts[i, :] <= r)

    # Now add the constraint that the ellipsoid is in the convex hull Ax<=b.
    # Using s-lemma, we know that the constraint is
    # ∃ λᵢ > 0,
    # s.t [P            q -λᵢaᵢ/2]  is positive semidefinite.
    #     [(q-λᵢaᵢ/2)ᵀ     λᵢbᵢ-r]
    num_faces = A.shape[0]
    lambda_var = cp.Variable(num_faces)
    constraints.append(lambda_var >= 0)
    Q = [None] * num_faces
    for i in range(num_faces):
        Q[i] = cp.Variable((dim+1, dim+1), PSD=True)
        constraints.append(Q[i][:dim, :dim] == P)
        constraints.append(Q[i][:dim, dim] == q - lambda_var[i] * A[i, :]/2)
        constraints.append(Q[i][-1, -1] == lambda_var[i] * b[i] - r)

    prob = cp.Problem(cp.Minimize(0), constraints)
    try:
        prob.solve(verbose=False)
    except cp.error.SolverError:
        return None, None, None, None

    if prob.status == 'optimal':
        P_val = P.value
        q_val = q.value
        r_val = r.value
        lambda_val = lambda_var.value
        return P_val, q_val, r_val, lambda_val
    else:
        return None, None, None, None


def draw_ellipsoid(P, q, r, outside_pts, inside_pts):
    """
    Draw an ellipsoid defined as {x | xᵀPx + 2qᵀx ≤ r }
    This ellipsoid is equivalent to
    |Lx + L⁻¹q| ≤ √(r + qᵀP⁻¹q)
    where L is the symmetric matrix satisfying L * L = P
    """
    fig = plt.figure()
    dim = P.shape[0]
    L = scipy.linalg.sqrtm(P)
    radius = np.sqrt(r + q@(np.linalg.solve(P, q)))
    if dim == 2:
        # first compute the points on the unit sphere
        theta = np.linspace(0, 2 * np.pi, 200)
        sphere_pts = np.vstack((np.cos(theta), np.sin(theta)))
        ellipsoid_pts = np.linalg.solve(
            L, radius * sphere_pts - (np.linalg.solve(L, q)).reshape((2, -1)))
        ax = fig.add_subplot(111)
        ax.plot(ellipsoid_pts[0, :], ellipsoid_pts[1, :], c='blue')

        ax.scatter(outside_pts[:, 0], outside_pts[:, 1], c='red')
        ax.scatter(inside_pts[:, 0], inside_pts[:, 1], s=20, c='green')
        ax.axis('equal')
        plt.show()
    if dim == 3:
        u = np.linspace(0, np.pi, 30)
        v = np.linspace(0, 2*np.pi, 30)

        sphere_pts_x = np.outer(np.sin(u), np.sin(v))
        sphere_pts_y = np.outer(np.sin(u), np.cos(v))
        sphere_pts_z = np.outer(np.cos(u), np.ones_like(v))
        sphere_pts = np.vstack((
            sphere_pts_x.reshape((1, -1)), sphere_pts_y.reshape((1, -1)),
            sphere_pts_z.reshape((1, -1))))
        ellipsoid_pts = np.linalg.solve(
            L, radius * sphere_pts - (np.linalg.solve(L, q)).reshape((3, -1)))
        ax = plt.axes(projection='3d')
        ellipsoid_pts_x = ellipsoid_pts[0, :].reshape(sphere_pts_x.shape)
        ellipsoid_pts_y = ellipsoid_pts[1, :].reshape(sphere_pts_y.shape)
        ellipsoid_pts_z = ellipsoid_pts[2, :].reshape(sphere_pts_z.shape)
        ax.plot_wireframe(ellipsoid_pts_x, ellipsoid_pts_y, ellipsoid_pts_z)
        ax.scatter(outside_pts[:, 0], outside_pts[:, 1], outside_pts[:, 2],
                   c='red')
        ax.scatter(inside_pts[:, 0], inside_pts[:, 1], inside_pts[:, 2], s=20,
                   c='green')
        ax.axis('equal')
        plt.show()


def find_large_ellipsoid(pts, max_iterations):
    """
    We find a large ellipsoid within the convex hull of @p pts but not
    containing any point in @p pts.
    The algorithm proceeds iteratively
    1. Start with outside_pts = pts, inside_pts = z where z is a random point
       in the convex hull of @p outside_pts.
    2. while num_iter < max_iterations
    3.   Solve an SDP to find an ellipsoid that is within the convex hull of
         @p pts, not containing any outside_pts, but contains all inside_pts.
    4.   If the SDP in the previous step is infeasible, then remove z from
         inside_pts, and append it to the outside_pts.
    5.   Randomly sample a point in the convex hull of @p pts, if this point is
         outside of the current ellipsoid, then append it to inside_pts.
    6.   num_iter += 1
    When the iterations limit is reached, we report the ellipsoid with the
    maximal volume.
    @param pts pts[i, :] is the i'th points that has to be outside of the
    ellipsoid.
    @param max_iterations The iterations limit.
    @return (P, q, r) The largest ellipsoid is parameterized as
    {x | xᵀPx + 2qᵀx ≤ r }
    """
    dim = pts.shape[1]
    A, b, hull = get_hull(pts)
    hull_vertices = pts[hull.vertices]
    deln = pts[Delaunay(hull_vertices).simplices]

    outside_pts = pts
    z = centered_sample_from_convex_hull(pts)
    inside_pts = z.reshape((1, -1))

    num_iter = 0
    max_ellipsoid_volume = -np.inf
    while num_iter < max_iterations:
        (P, q, r, lambda_val) = find_ellipsoid(outside_pts, inside_pts, A, b)
        if P is not None:
            volume = compute_ellipsoid_volume(P, q, r)
            if volume > max_ellipsoid_volume:
                max_ellipsoid_volume = volume
                P_best = P
                q_best = q
                r_best = r
            else:
                # Adding the last inside_pts doesn't increase the ellipsoid
                # volume, so remove it.
                inside_pts = inside_pts[:-1, :]
        else:
            outside_pts = np.vstack((outside_pts, inside_pts[-1, :]))
            inside_pts = inside_pts[:-1, :]

        # Now take a new sample that is outside of the ellipsoid.
        sample_pts = uniform_sample_from_convex_hull(deln, dim, 20)
        is_in_ellipsoid = np.sum(sample_pts.T*(P_best @ sample_pts.T), axis=0)\
            + 2 * sample_pts @ q_best <= r_best
        if np.all(is_in_ellipsoid):
            # All the sampled points are in the ellipsoid, the ellipsoid is
            # already large enough.
            return P_best, q_best, r_best
        else:
            inside_pts = np.vstack((
                inside_pts, sample_pts[np.where(~is_in_ellipsoid)[0][0], :]))
            num_iter += 1

    return P_best, q_best, r_best

if __name__ == "__main__":
    pts = np.array([[0., 0.], [0., 1.], [1., 1.], [1., 0.], [0.2, 0.4]])
    max_iterations = 10
    P, q, r = find_large_ellipsoid(pts, max_iterations)

我也将代码放在 github repo中
结果
这是一个简单的2D示例的结果，假设我们要找到一个不包含下图中五个红色点的大椭圆体。这是第一次迭代后的结果
。红色点是“外部点”(初始外部点是v are，v 2，...，v，)，绿色点是初始“内部点”。
在第二次迭代中，椭球变为
。通过再添加一个“内部点”(绿点)，椭球会变大。
此gif显示了10个迭代的动画。