稍后编辑：我在这里上传了我的原始数据样本。它实际上是 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）。 我正在寻找点内拟合的最大体积椭球。
我试图
通过修改阈值来调整来自最小封闭椭圆体（又名外边界椭圆体）的代码，我的逻辑开始是所有点都应该小于 < 1 给定椭圆体方程。但没有成功。err > tol
我还尝试了对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理想的代码。
