当给定多面体的顶点时，计算重心和体积

Question

给定一个凸多面体（3D）的顶点位置，我需要计算该多面体的重心和体积。以下代码可在Mathworks网站上找到。

function C = centroid(P)
k=convhulln(P);
if length(unique(k(:)))<size(P,1)
    error('Polyhedron is not convex.');
end
T = delaunayn(P);
n = size(T,1);
W = zeros(n,1);
C=0;
for m = 1:n
    sp = P(T(m,:),:);
    [null,W(m)]=convhulln(sp);
    C = C + W(m) * mean(sp);
end
C=C./sum(W);
return
end

代码优雅但速度极慢。我需要计算数千个多面体的体积和重心，需要进行数百次计算。在当前状态下使用此代码不可行。有人知道更好的方法或者能否使这段代码更快？我可以想到一些小的更改，例如用表达式来代替mean。

- John Smith

Answer 1

有一种更简单的方法可以用最少的工作量计算体积。第一种方法使用多面体的3个本地拓扑信息集，即边缘的切向单位矢量、在此切向上的平面法向量的单位矢量以及多面体本身的单位矢量（这些都很容易从顶点中提取出来）。请参考Polyhedron的体积以获取更多详细信息。

第二种方法使用面积、法向量和面重心来计算多面体的体积，可参考Wikipedia文章。这两种算法都相当简单，非常易于实现，并且通过简单的求和结构易于向量化。我认为这两种方法都比对多面体进行充分的镶嵌要快得多。

然后，可以通过应用散度定理，将整个多面体体积的积分转换为多面体表面的积分来计算多面体的重心。详细描述可以在在三维空间中计算多面体的体积和重心中找到。我没有检查过将多面体镶嵌成三角形是否真的有必要，或者可以使用多面体更复杂的多边形表面进行计算，但无论如何，面积的镶嵌都比体积的镶嵌简单得多。总的来说，这种组合方法应该比体积法更快。

- Rainer

Answer 2

如果您认为quickhull不够好，那么您唯一的选择就是cudahull，如果您想要精确的解决方案。尽管如此，即使这样，您最多只能获得约40倍的增长（看起来是这样）。

我假设您制作的凸包至少有10个顶点（如果少于这个数量，您无法做太多）。如果您不介意“足够接近”的解决方案。您可以创建一个版本的quickhull，限制每个多边形的顶点数。您限制计算的顶点数也将允许在需要时计算最大误差。

问题在于，随着凸包上的顶点数量趋近于无穷大，您最终会得到一个球体。这意味着由于quick hull的工作方式，您添加到凸包中的每个额外顶点对其影响都比之前的顶点小*。

*根据quickhull的编码方式，这可能只是一般意义上的真实情况。要使其在实践中成立，需要修改quickhull的递归算法，因此，虽然“下一个顶点”总是被计算出来（除了在添加最后一个顶点之后或该部分没有剩余点时），但实际上顶点是按照最大化多面体体积增加的顺序添加到凸包中的（可能按照最远到最近的顺序）。跟踪添加顶点的顺序会带来一些性能成本，但只要待处理的凸包点与待处理的点的比率足够高，就应该值得这样做。至于误差，最好的选择可能是在达到实际凸包或最大体积增加小于当前总体积的某个分数时停止算法。如果性能更重要，则只需限制每个多边形的凸包点数即可。

您还可以查看各种近似凸包算法，但我上面概述的方法应该适用于体积/质心近似，并具有确定误差的能力。

- Nuclearman

Answer 3

你能加速代码的程度取决于你计算质心的方式。关于质心计算的选项，请参考这个回答。事实证明，如果你需要计算立体多面体的质心，那么你基本上是没有办法的。然而，如果只有多面体的顶点具有权重，那么你可以简单地写成：

[k,volume] = convhulln(P);
centroid = mean(P(k,:));

- Jonas