Scipy: 凸包的中心点

Question

我该如何使用Python和Scipy计算凸包的重心？我找到的所有方法都是计算面积和体积的。

祝好，Frank。

- OD IUM

---

如果你正在寻找凸壳所包围的体积的实际重心，假设密度均匀，请不要使用下面的答案。它们只会给出描述凸壳的点坐标的平均值。 - Raketenolli

Answer 1

假设您使用scipy.spatial.ConvexHull构建了凸壳，返回的对象应该含有点的位置信息，因此质心可以简单地表示为：

import numpy as np
from scipy.spatial import ConvexHull

points = np.random.rand(30, 2)   # 30 random points in 2-D
hull = ConvexHull(points)

#Get centoid
cx = np.mean(hull.points[hull.vertices,0])
cy = np.mean(hull.points[hull.vertices,1])

您可以按照以下方式绘制：

import matplotlib.pyplot as plt
#Plot convex hull
for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')

#Plot centroid
plt.plot(cx, cy,'x',ms=20)
plt.show()

Scipy凸包基于Qhull，Qhull的文档中应该有centrum方法，文档中说明：

Centrum是面的超平面上的一个点。Centrum是面顶点的平均值。如果每个centrum都在相邻面的超平面下，则相邻面是凸面。

对于简单面，centrum与重心相同，

对于具有d个顶点的单形面，centrum等效于重心或重心.

由于Scipy似乎没有提供这个功能，您可以在子类到凸壳中定义自己的centrum方法。

class CHull(ConvexHull):

    def __init__(self, points):
        ConvexHull.__init__(self, points)

    def centrum(self):

        c = []
        for i in range(self.points.shape[1]):
            c.append(np.mean(self.points[self.vertices,i]))

        return c

 hull = CHull(points)
 c = hull.centrum()

- Ed Smith

---

谢谢帮忙。我将检查计算出的质心是否始终位于外壳内... 另外，你知道如何在3D中绘制外壳吗？我正在使用Axes3D和scatter。只需散布简单形状即可吗？ - OD IUM

我刚刚添加了一个带有附加功能的方法，因此应该可以在N维中工作。我猜你可以使用matplotlib将顶点绘制为散点图，如http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#scatter-plots或者使用上面提到的简单形式绘制线图，但是在3D中... - Ed Smith

Answer 2

如果凸壳的点在周长上不是均匀分布的，简单的“平均”方法是不正确的，或者至少会给出一个倾斜的答案。我使用的最佳方法是计算凸壳点的德劳内三角形的重心。这将会以计算形状的COM（质心）作为加权计算质心，而不仅仅是顶点的平均值:

from scipy.spatial import ConvexHull, Delaunay

def _centroid_poly(poly):
    
    T = Delaunay(poly).simplices
    n = T.shape[0]
    W = np.zeros(n)
    C = 0
    
    for m in range(n):
        sp = poly[T[m, :], :]
        W[m] = ConvexHull(sp).volume
        C += W[m] * np.mean(sp, axis=0)
    
    return C / np.sum(W)

poly = np.random.rand(30, 2)
# will calculate the centroid of the convex hull of poly
centroid_hull = _centroid_poly(poly)

这样的东西应该可以正常工作。

- Alec Day

Answer 3

要找到船体顶点的几何中心，只需使用以下方法：

# Calculate geometric centroid of convex hull 
hull = ConvexHull(points)   
centroid = np.mean(points[hull.vertices, :], axis=0)

要绘制船体，请尝试以下操作：

import numpy as np
import pylab as pl
import scipy as sp
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3

#  Plot polyhedra    
ax = a3.Axes3D(pl.figure())
facetCol = sp.rand(3)
for simplex in hull.simplices:
    vtx = [points[simplex[0],:], points[simplex[1],:], points[simplex[2],:]]
    tri = a3.art3d.Poly3DCollection([vtx], linewidths = 2, alpha = 0.8)
    tri.set_color(facetCol)
    tri.set_edgecolor('k')
    ax.add_collection3d(tri)

# Plot centroid      
ax.scatter(centroid0], centroid[1], centroid[2]) 
plt.axis('equal')
plt.axis('off')
plt.show()

- DrBwts

Answer 4

即使外壳周长上的点是不规则分布的，此解决方案也是正确的。

import numpy as np
from scipy.spatial import ConvexHull

def areaPoly(points):
    area = 0
    nPoints = len(points)
    j = nPoints - 1
    i = 0
    for point in points:
        p1 = points[i]
        p2 = points[j]
        area += (p1[0]*p2[1])
        area -= (p1[1]*p2[0])
        j = i
        i += 1

    area /= 2
    return area

def centroidPoly(points):
    nPoints = len(points)
    x = 0
    y = 0
    j = nPoints - 1
    i = 0

    for point in points:
        p1 = points[i]
        p2 = points[j]
        f = p1[0]*p2[1] - p2[0]*p1[1]
        x += (p1[0] + p2[0])*f
        y += (p1[1] + p2[1])*f
        j = i
        i += 1

    area = areaPoly(hull_points)
    f = area*6
    return [x/f, y/f]

# 'points' is an array of tuples (x, y)
points = np.array(points)
hull = ConvexHull(points)
hull_points = points[hull.vertices]
centroid = centroidPoly(hull_points)

- exsurge-domine

---

感谢您抽出时间回答问题。由于问题通常会收到多个答案，因此鼓励您为您的答案提供一些上下文，例如解释为什么您的解决方案可能比其他解决方案更好。 - DaveL17