回答问题

稍后编辑:我在这里上传了我原始数据的样本。它实际上是 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 Enclosure Ellipsoid(又名外边界椭球)调整代码

通过修改阈值err > tol,我的逻辑开始认为所有点都应该小于 < 1 给定椭球方程。但没有成功。

我还在mosek上尝试了 Loewner-John 改编,但我不知道如何描述超平面与 3D 多面体(Ax <u003d 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中的代码。

Answers

问题陈述

给定许多点v₁, v₂, ..., vₙ,找到一个满足两个约束的大椭球:

1. 椭球在凸包中 H u003d ConvexHull(v1, v2, ..., vn)。

2. 点 v1、v2、...、vn 均不在椭球内。

我提出了一个迭代过程来找到一个满足这两个约束的大椭球。在每次迭代中,我们需要解决一个半定规划问题。此迭代过程保证收敛,但不保证收敛到_globally_ 最大的椭球。

方法

求单个椭球

我们迭代过程的核心是,在每次迭代中,我们找到一个满足 3 个条件的椭球:

• 椭球包含在 ConvexHull(v1, v2, ..., vn) u003d {x |斧头<u003db}。

• 对于一组点u1, ... um(其中{v1, v2, ..., vn} ⊂ {u1, ... um},即点云中的给定点属于这组点u1 , ... um),椭球不包含 u1, ... um 中的任何点。我们称这个集合 u1, ... um 为“外点”。

• 对于一组点 w1,..., wk (其中 {w1,..., wk} ∩ {v1, v2, ..., vn} u003d ∅, 即 v1, v2, . .., vn 属于 {w1,..., wk}),椭球包含所有点 w1,..., wk。我们称这个集合 w1,..., wk 为“内点”。

直观的想法是“内点” w1,..., wk 表示椭圆体的体积。我们将新点添加到“内部点”以增加椭球体积。

为了通过凸优化找到这样的椭球,我们将椭球参数化为

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

我们将搜索P, q, r。

“外点” u1, ... um 都在椭球之外的条件被表述为

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

这是对P, q, r的线性约束。

“内点” w1,..., wk 都在椭球内的条件被表述为

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

这也是对P, q, r的线性约束。

我们还施加了约束

P is positive definite

P是正定的,加上存在满足 wiTPwi + 2qTwi <u003d r 的点 wi 的约束保证集合 {x | xTPx + 2qTx ≤ r} 是一个椭球体。

我们也有椭球在凸包内的约束Hu003d{x | aiTx≤ bi, iu003d1,...,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]

因此,我们可以解决以下半定规划问题,以找到包含所有“内点”、不包含任何“外点”且在凸包内的椭球H

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 + qTP−1q)/power(det(P), 1/3) 成正比。

迭代过程。

我们将“外部点”初始化为点云中的所有点 v1、v2、...、vn,将“内部点”初始化为凸包 H 中的单个点w₁。在每次迭代中,我们使用前面小节中的find_ellipsoid函数来查找 H 内包含所有“内点”但不包含任何“外点”的椭球。根据find_ellipsoid中 SDP 的结果,我们执行以下操作

• 如果SDP可行。然后,我们将新发现的椭球与迄今为止发现的最大椭球进行比较。如果这个新的椭球体更大,那么接受它作为迄今为止发现的最大的椭球体。

• 如果SDP不可行,那么我们将“inside points”中的最后一个点去掉,将这个点添加到“outside point”中。

在这两种情况下,我们都在凸包 H 中取一个新的样本点,将该样本点添加到“内部点”,然后再次求解 SDP。

完整算法如下

1. 将“外点”初始化为 v1, v2, ..., vn,将“内点”初始化为凸包 H 中的单个随机点。

2. 当 iter < max_iterations:

3.解决SDPP, q, r = find_ellipsoid(outside_points, inside_points, A, b)。

4.如果SDP可行且volume(Ellipsoid(P, q, r)) >最大_volume,设置P_best = P, q_best=q, r_best = r。

5.如果SDP不可行,pt u003d inside\points.pop\last(), outside\points.push_back(pt)。

6.在H中随机采样一个新点,将点追加到“inside points”中,iter +u003d 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
    power(r + qᵀP⁻¹q, dim/2)/sqrt(det(P))
    We return this number.
    """
    dim = P.shape[0]
    return np.power(r + q @ np.linalg.solve(P, q)), dim/2) / \
        np.sqrt(np.linalg.det(P))


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 = hull_vertices[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 示例的结果,假设我们要找到一个不包含下图中的五个红点的大椭球。这是第一次迭代后的结果。红点是“外点”(初始外点是 v1, v2, ..., vn),绿点是初始“内点”。

在第二次迭代中,椭球变为

。通过再增加一个“内点”(绿点),椭圆体变得更大。

这个 gif 显示了 10 个迭代的动画。