polytope/set 个点中的最大体积内接椭球体
Maximum volume inscribed ellipsoid in a polytope/set of points
稍后编辑：我上传了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]160415 114441\Series 052 [CT - Abdomen WT 1 0 I31f 3]"
points = get_surface_points(folder_path)
我有一组 3D 点 (n > 1000) space 描述了一个类似空心卵形的形状。 我想要拟合所有点内的椭圆体 (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 中理想地寻找代码。
表示椭圆体（的表面）的一种也许是数学标准的方法是它是集合
{ X | (X-a)'*inv(C)*(X-a) = 1}
the solid ellipsoid is then
{ X | (X-a)'*inv(C)*(X-a) <= 1}
这里C是3x3对称正定矩阵，a是椭球的'centre'
我们可以通过使用 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 进行变换。
最简单的做法就是用常数 s 缩放 M，并保持 a 不变。这相当于缩放所有 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'
你也可以换一个。我留给感兴趣的细节 reader...
问题陈述
给定多个点 v₁, v₂, ..., vₙ，找到满足两个约束的大椭圆体：
- 椭球在凸包内 ℋ = ConvexHull(v₁, v₂, ..., vₙ).
- None 点 v₁, v₂, ..., vₙ 在椭球内。
我提出了一个迭代过程来找到满足这两个约束的大椭圆体。在每次迭代中，我们需要解决一个半定规划问题。此迭代过程保证收敛，但不能保证收敛到 全局 最大的椭圆体。
方法
找到一个椭圆体
我们迭代过程的核心是在每次迭代中，我们找到一个满足3个条件的椭球：
- 椭圆体包含在 ConvexHull(v₁, v₂, ..., vₙ) = {x | ax<=b}.
- 对于一组点u₁, ... uₘ (其中{v₁, v2, ..., vₙ} ⊂ {u₁, ... uₘ}，即点云中的给定点属于此点集 u₁, ... uₘ)，椭球不包含 u₁, ... uₘ 中的任何点。我们称这个集合 u₁, ... uₘ 为“外点”。
- 对于一组点w₁,...,wₖ（其中{w₁,...,wₖ}∩{v₁,v₂,...,vₙ}=∅，即none的v₁, v2, ..., 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}是一个椭球体。
我们也有椭圆体在凸包内的约束ℋ={x | aᵢᵀx≤ bᵢ, i=1,...,l}（即有l个半空间作为ℋ的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)/power(det(P), 1/3).
迭代过程。
我们将“外部点”初始化为点云中的所有点 v₁、v2、...、vₙ，并将“内部点”初始化为凸包 ℋ 中的单个点 w₁。在每次迭代中，我们使用上一小节中的 find_ellipsoid 函数来查找 ℋ 内包含所有“内部点”但不包含任何“外部点”的椭圆体。根据 find_ellipsoid 中 SDP 的结果，我们执行以下操作
- 如果SDP可行。然后我们将新发现的椭圆体与迄今为止发现的最大椭圆体进行比较。如果这个新的椭球体更大，那么接受它作为迄今为止发现的最大的椭球体。
- 如果SDP不可行，那么我们去掉“inside points”中的最后一个点，将这个点添加到“outside points”中。
在这两种情况下，我们然后在凸包ℋ中取一个新的样本点，将该样本点添加到“内部点”，然后再次求解SDP。
完整算法如下
- 将“外部点”初始化为v₁, v2, ..., vₙ，将“内部点”初始化为凸包ℋ中的单个随机点。
- while iter < max_iterations:
- 解决 SDP
P, q, r = find_ellipsoid(outside_points, inside_points, A, b)。
- 如果 SDP 可行且 volume(Ellipsoid(P, q, r)) > largest_volume，设置
P_best = P, q_best=q, r_best = r.
- 如果SDP不可行，pt = inside_points.pop_last(), outside_points.push_back(pt).
- 在ℋ中随机抽取一个新的点，将该点追加到“内点”，iter += 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
@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
结果
这是一个简单的二维示例的结果，假设我们要找到一个不包含下图中五个红点的大椭圆体。这是第一次迭代后的结果
。红色的点是“外点”（初始的外点是v₁, v2, ..., vₙ），绿色的点是初始的“内点”。
在第二次迭代中，椭圆体变为
。通过再添加一个“内部点”（绿点），椭圆体会变大。
此 gif 显示了 10 次迭代的动画。
这个答案是否有效取决于你的数据中有多少噪音。这个想法是首先找到点云的凸包，然后找到适合这个凸包的最大椭圆体。如果你的大部分点都靠近它们描述的椭球表面，那么这个近似值就不会是 "too bad".
为此，请注意凸包可以用一组线性不等式来描述 Ax<=b。
请注意，边界椭球可以用 E={Bx+d for ||x||_2<=1} 来描述，其中 B 是一个半正定矩阵，描述了椭球的拉伸方式和方向，d 是一个描述其偏移量的向量。
请注意，椭圆体的体积由 det(B^-1) 给出。如果我们试图最大化或最小化这个行列式，我们会失败，因为这会产生一个非凸问题。但是，应用对数变换 log(det(B^-1)) 会使问题再次凸出。我们将要使用的优化程序不允许矩阵求逆，但很容易证明上述等价于 -log(det(B)).
最后，一些支撑代数操作给了我们优化问题：
minimize -log(det(B))
s.t. ||B*A_i||_2 + a_i^T * d <= b_i, i = 1, ..., m
B is PSD
我们可以使用 CVXPY 在 Python 中解决这个问题，如下所示：
#!/usr/bin/env python3
from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import sklearn.datasets
#From:
def random_point_ellipsoid(a,b,c,x0,y0,z0):
"""Generate a random point on an ellipsoid defined by a,b,c"""
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
def random_point_ellipse(W,d):
# random angle
alpha = 2 * np.pi * np.random.random()
# vector on that angle
pt = np.array([np.cos(alpha),np.sin(alpha)])
# Ellipsoidize it
return W@pt+d
def GetRandom(dims, Npts):
if dims==2:
W = sklearn.datasets.make_spd_matrix(2)
d = np.array([2,3])
points = np.array([random_point_ellipse(W,d) for i in range(Npts)])
elif dims==3:
points = np.array([random_point_ellipsoid(3,5,7,2,3,3) for i in range(Npts)])
else:
raise Exception("dims must be 2 or 3!")
noise = np.random.multivariate_normal(mean=[0]*dims, cov=0.2*np.eye(dims), size=Npts)
return points+noise
def GetHull(points):
dim = points.shape[1]
hull = ConvexHull(points)
A = hull.equations[:,0:dim]
b = hull.equations[:,dim]
return A, -b, hull #Negative moves b to the RHS of the inequality
def Plot(points, hull, B, d):
fig = plt.figure()
if points.shape[1]==2:
ax = fig.add_subplot(111)
ax.scatter(points[:,0], points[:,1])
for simplex in hull.simplices:
plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
display_points = np.array([random_point_ellipse([[1,0],[0,1]],[0,0]) for i in range(100)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1])
elif points.shape[1]==3:
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:,1], points[:,2])
display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
plt.show()
def FindMaximumVolumeInscribedEllipsoid(points):
"""Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
dim = points.shape[1]
A,b,hull = GetHull(points)
B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
d = cp.Variable(dim) #Center
constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
optval = prob.solve()
if optval==np.inf:
raise Exception("No solution possible!")
print(f"Optimal value: {optval}")
Plot(points, hull, B.value, d.value)
return B.value, d.value
FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=2, Npts=100))
FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=3, Npts=100))
解的计算速度很快。
在视觉上，这给出（对于 2D）：
请注意，我添加了很多噪音来强调正在发生的事情。
对于 3D：
虽然上面的代码是为两个或三个维度编写的，但您可以轻松地将其调整为任意数量的维度，尽管可视化会变得更加困难。
如果凸包不好，而你想要某种 "interior convex hull"，那会更难：这个包没有明确定义。但是，您可以use alpha shapes尝试找到这样的船体，然后使用上面的算法来解决它。
另请注意，由于我们使用凸多胞形来限制椭圆，而不是点本身，即使这些点完美地描述了一个椭圆体，我们最终也会低估体积。我们可以想象一下，如下所示：
如果正方形的顶点是点，那么正方形就是它们的凸包。以船体为界的圆明显小于仅以点为界的圆。
编辑：要获得音量，您需要将像素索引转换为 DICOM 图像的坐标系，就像这样（注意：我不确定我是否我们已经用正确的值缩放了正确的坐标，但是根据您对数据的了解，您将能够解决这个问题):
from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import os
import sklearn.datasets
import SimpleITK as sitk
import code
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(dcm_img):
x_spacing, y_spacing, z_spacing = dcm_img.GetSpacing()
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)
surface_points = surface_points.astype(np.float64)
surface_points[:,0] *= x_spacing/10
surface_points[:,1] *= y_spacing/10
surface_points[:,2] *= z_spacing/10
return surface_points
def get_dcm_image(folder_path):
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:
raise Exception('Non-readable DICOM Data: ', folder_path)
return dcm_img
def GetHull(points):
dim = points.shape[1]
hull = ConvexHull(points)
A = hull.equations[:,0:dim]
b = hull.equations[:,dim]
return A, -b, hull #Negative moves b to the RHS of the inequality
def FindMaximumVolumeInscribedEllipsoid(points):
"""Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
dim = points.shape[1]
A,b,hull = GetHull(points)
B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
d = cp.Variable(dim) #Center
constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
optval = prob.solve()
if optval==np.inf:
raise Exception("No solution possible!")
print(f"Optimal value: {optval}")
return B.value, d.value
#From:
def random_point_ellipsoid(a,b,c,x0,y0,z0):
"""Generate a random point on an ellipsoid defined by a,b,c"""
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
def Plot(points, B, d):
hull = ConvexHull(points)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:,1], points[:,2], marker=".")
display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
plt.show()
folder_path = r"data"
dcm_img = get_dcm_image(folder_path)
points = get_surface_points(dcm_img)
B, d = FindMaximumVolumeInscribedEllipsoid(points)
Plot(points, B, d)
ball_vol = 4/3.0*np.pi*(1.0**3)
print("DCM vol: ", get_volume_ml(dcm_img))
print("Ellipsoid Volume: ", np.linalg.det(B) * ball_vol)
这给
DCM vol: 16.2786318359375
Ellipsoid Volume: 11.947614772444393
我认为如果你可以假设椭圆体的质心和你的点是相同的，你可以解决椭圆体通过最近或最远的 n 点的方程质心。
我不确定我是否有时间加强这个答案，但这种方法应该很容易用标准 Python 工具实现，例如：
https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.center_of_mass.html
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html
也许还有 SymPy 来求解解析方程。
稍后编辑：我上传了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]160415 114441\Series 052 [CT - Abdomen WT 1 0 I31f 3]"
points = get_surface_points(folder_path)
我有一组 3D 点 (n > 1000) space 描述了一个类似空心卵形的形状。 我想要拟合所有点内的椭圆体 (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 中理想地寻找代码。
表示椭圆体（的表面）的一种也许是数学标准的方法是它是集合
{ X | (X-a)'*inv(C)*(X-a) = 1}
the solid ellipsoid is then
{ X | (X-a)'*inv(C)*(X-a) <= 1}
这里C是3x3对称正定矩阵，a是椭球的'centre'
我们可以通过使用 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 进行变换。
最简单的做法就是用常数 s 缩放 M，并保持 a 不变。这相当于缩放所有 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'
你也可以换一个。我留给感兴趣的细节 reader...
问题陈述
给定多个点 v₁, v₂, ..., vₙ，找到满足两个约束的大椭圆体：
- 椭球在凸包内 ℋ = ConvexHull(v₁, v₂, ..., vₙ).
- None 点 v₁, v₂, ..., vₙ 在椭球内。
我提出了一个迭代过程来找到满足这两个约束的大椭圆体。在每次迭代中，我们需要解决一个半定规划问题。此迭代过程保证收敛，但不能保证收敛到 全局 最大的椭圆体。
方法
找到一个椭圆体
我们迭代过程的核心是在每次迭代中，我们找到一个满足3个条件的椭球：
- 椭圆体包含在 ConvexHull(v₁, v₂, ..., vₙ) = {x | ax<=b}.
- 对于一组点u₁, ... uₘ (其中{v₁, v2, ..., vₙ} ⊂ {u₁, ... uₘ}，即点云中的给定点属于此点集 u₁, ... uₘ)，椭球不包含 u₁, ... uₘ 中的任何点。我们称这个集合 u₁, ... uₘ 为“外点”。
- 对于一组点w₁,...,wₖ（其中{w₁,...,wₖ}∩{v₁,v₂,...,vₙ}=∅，即none的v₁, v2, ..., 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}是一个椭球体。
我们也有椭圆体在凸包内的约束ℋ={x | aᵢᵀx≤ bᵢ, i=1,...,l}（即有l个半空间作为ℋ的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)/power(det(P), 1/3).
迭代过程。
我们将“外部点”初始化为点云中的所有点 v₁、v2、...、vₙ，并将“内部点”初始化为凸包 ℋ 中的单个点 w₁。在每次迭代中，我们使用上一小节中的 find_ellipsoid 函数来查找 ℋ 内包含所有“内部点”但不包含任何“外部点”的椭圆体。根据 find_ellipsoid 中 SDP 的结果，我们执行以下操作
- 如果SDP可行。然后我们将新发现的椭圆体与迄今为止发现的最大椭圆体进行比较。如果这个新的椭球体更大，那么接受它作为迄今为止发现的最大的椭球体。
- 如果SDP不可行，那么我们去掉“inside points”中的最后一个点，将这个点添加到“outside points”中。
在这两种情况下，我们然后在凸包ℋ中取一个新的样本点，将该样本点添加到“内部点”，然后再次求解SDP。
完整算法如下
- 将“外部点”初始化为v₁, v2, ..., vₙ，将“内部点”初始化为凸包ℋ中的单个随机点。
- while iter < max_iterations:
- 解决 SDP
P, q, r = find_ellipsoid(outside_points, inside_points, A, b)。 - 如果 SDP 可行且 volume(Ellipsoid(P, q, r)) > largest_volume，设置
P_best = P, q_best=q, r_best = r. - 如果SDP不可行，pt = inside_points.pop_last(), outside_points.push_back(pt).
- 在ℋ中随机抽取一个新的点，将该点追加到“内点”，iter += 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
@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
结果
这是一个简单的二维示例的结果，假设我们要找到一个不包含下图中五个红点的大椭圆体。这是第一次迭代后的结果
在第二次迭代中，椭圆体变为
此 gif 显示了 10 次迭代的动画。
这个答案是否有效取决于你的数据中有多少噪音。这个想法是首先找到点云的凸包，然后找到适合这个凸包的最大椭圆体。如果你的大部分点都靠近它们描述的椭球表面，那么这个近似值就不会是 "too bad".
为此，请注意凸包可以用一组线性不等式来描述 Ax<=b。
请注意，边界椭球可以用 E={Bx+d for ||x||_2<=1} 来描述，其中 B 是一个半正定矩阵，描述了椭球的拉伸方式和方向，d 是一个描述其偏移量的向量。
请注意，椭圆体的体积由 det(B^-1) 给出。如果我们试图最大化或最小化这个行列式，我们会失败，因为这会产生一个非凸问题。但是，应用对数变换 log(det(B^-1)) 会使问题再次凸出。我们将要使用的优化程序不允许矩阵求逆，但很容易证明上述等价于 -log(det(B)).
最后，一些支撑代数操作给了我们优化问题：
minimize -log(det(B))
s.t. ||B*A_i||_2 + a_i^T * d <= b_i, i = 1, ..., m
B is PSD
我们可以使用 CVXPY 在 Python 中解决这个问题，如下所示：
#!/usr/bin/env python3
from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import sklearn.datasets
#From:
def random_point_ellipsoid(a,b,c,x0,y0,z0):
"""Generate a random point on an ellipsoid defined by a,b,c"""
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
def random_point_ellipse(W,d):
# random angle
alpha = 2 * np.pi * np.random.random()
# vector on that angle
pt = np.array([np.cos(alpha),np.sin(alpha)])
# Ellipsoidize it
return W@pt+d
def GetRandom(dims, Npts):
if dims==2:
W = sklearn.datasets.make_spd_matrix(2)
d = np.array([2,3])
points = np.array([random_point_ellipse(W,d) for i in range(Npts)])
elif dims==3:
points = np.array([random_point_ellipsoid(3,5,7,2,3,3) for i in range(Npts)])
else:
raise Exception("dims must be 2 or 3!")
noise = np.random.multivariate_normal(mean=[0]*dims, cov=0.2*np.eye(dims), size=Npts)
return points+noise
def GetHull(points):
dim = points.shape[1]
hull = ConvexHull(points)
A = hull.equations[:,0:dim]
b = hull.equations[:,dim]
return A, -b, hull #Negative moves b to the RHS of the inequality
def Plot(points, hull, B, d):
fig = plt.figure()
if points.shape[1]==2:
ax = fig.add_subplot(111)
ax.scatter(points[:,0], points[:,1])
for simplex in hull.simplices:
plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
display_points = np.array([random_point_ellipse([[1,0],[0,1]],[0,0]) for i in range(100)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1])
elif points.shape[1]==3:
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:,1], points[:,2])
display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
plt.show()
def FindMaximumVolumeInscribedEllipsoid(points):
"""Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
dim = points.shape[1]
A,b,hull = GetHull(points)
B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
d = cp.Variable(dim) #Center
constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
optval = prob.solve()
if optval==np.inf:
raise Exception("No solution possible!")
print(f"Optimal value: {optval}")
Plot(points, hull, B.value, d.value)
return B.value, d.value
FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=2, Npts=100))
FindMaximumVolumeInscribedEllipsoid(GetRandom(dims=3, Npts=100))
解的计算速度很快。
在视觉上，这给出（对于 2D）：
请注意，我添加了很多噪音来强调正在发生的事情。
对于 3D：
虽然上面的代码是为两个或三个维度编写的，但您可以轻松地将其调整为任意数量的维度，尽管可视化会变得更加困难。
如果凸包不好，而你想要某种 "interior convex hull"，那会更难：这个包没有明确定义。但是，您可以use alpha shapes尝试找到这样的船体，然后使用上面的算法来解决它。
另请注意，由于我们使用凸多胞形来限制椭圆，而不是点本身，即使这些点完美地描述了一个椭圆体，我们最终也会低估体积。我们可以想象一下，如下所示：
如果正方形的顶点是点，那么正方形就是它们的凸包。以船体为界的圆明显小于仅以点为界的圆。
编辑：要获得音量，您需要将像素索引转换为 DICOM 图像的坐标系，就像这样（注意：我不确定我是否我们已经用正确的值缩放了正确的坐标，但是根据您对数据的了解，您将能够解决这个问题):
from mpl_toolkits.mplot3d import axes3d
from scipy.spatial import ConvexHull
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
import os
import sklearn.datasets
import SimpleITK as sitk
import code
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(dcm_img):
x_spacing, y_spacing, z_spacing = dcm_img.GetSpacing()
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)
surface_points = surface_points.astype(np.float64)
surface_points[:,0] *= x_spacing/10
surface_points[:,1] *= y_spacing/10
surface_points[:,2] *= z_spacing/10
return surface_points
def get_dcm_image(folder_path):
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:
raise Exception('Non-readable DICOM Data: ', folder_path)
return dcm_img
def GetHull(points):
dim = points.shape[1]
hull = ConvexHull(points)
A = hull.equations[:,0:dim]
b = hull.equations[:,dim]
return A, -b, hull #Negative moves b to the RHS of the inequality
def FindMaximumVolumeInscribedEllipsoid(points):
"""Find the inscribed ellipsoid of maximum volume. Return its matrix-offset form."""
dim = points.shape[1]
A,b,hull = GetHull(points)
B = cp.Variable((dim,dim), PSD=True) #Ellipsoid
d = cp.Variable(dim) #Center
constraints = [cp.norm(B@A[i],2)+A[i]@d<=b[i] for i in range(len(A))]
prob = cp.Problem(cp.Minimize(-cp.log_det(B)), constraints)
optval = prob.solve()
if optval==np.inf:
raise Exception("No solution possible!")
print(f"Optimal value: {optval}")
return B.value, d.value
#From:
def random_point_ellipsoid(a,b,c,x0,y0,z0):
"""Generate a random point on an ellipsoid defined by a,b,c"""
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
def Plot(points, B, d):
hull = ConvexHull(points)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:,1], points[:,2], marker=".")
display_points = np.array([random_point_ellipsoid(1,1,1,0,0,0) for i in range(1000)])
display_points = display_points@B+d
ax.scatter(display_points[:,0], display_points[:,1], display_points[:,2])
plt.show()
folder_path = r"data"
dcm_img = get_dcm_image(folder_path)
points = get_surface_points(dcm_img)
B, d = FindMaximumVolumeInscribedEllipsoid(points)
Plot(points, B, d)
ball_vol = 4/3.0*np.pi*(1.0**3)
print("DCM vol: ", get_volume_ml(dcm_img))
print("Ellipsoid Volume: ", np.linalg.det(B) * ball_vol)
这给
DCM vol: 16.2786318359375
Ellipsoid Volume: 11.947614772444393
我认为如果你可以假设椭圆体的质心和你的点是相同的，你可以解决椭圆体通过最近或最远的 n 点的方程质心。
我不确定我是否有时间加强这个答案，但这种方法应该很容易用标准 Python 工具实现，例如：
https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.center_of_mass.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html
也许还有 SymPy 来求解解析方程。