最小凸多边形的质心

Question

我有一个物种在一个地区的发生点。我想计算包含该物种所有出现点的最小凸多边形的质心。谁能告诉我如何在 R 中做到这一点？

Answer 1

首先请注意，凸多边形的质心不等于其顶点的质心，除非在特殊情况下，例如当多边形是三角形时。 （如果您实际上想要顶点的质心而不是多边形的质心，则使用下面计算的ch 形式的sapply(ch[c("x", "y")], mean) 或上一节中计算的ix 形式的colMeans(X[ix, ])。）

要使质心将凸包多边形的面积分成三角形（Delaunay 三角化），请计算每个三角形的质心和面积，然后使用面积作为权重对质心进行加权平均。

tripack 包中的voronoi.mosaic 将提供三角形顶点的索引和三角形的面积。在其输出中p1、p2 和p3 是三角形顶点输入X 的索引，area 是相应的区域。由此我们计算x 和y 作为三角形质心的坐标。然后在他们计算之后的行中，我们采用他们的加权平均值来获得整体质心cen。最后，我们绘制所有内容。

library(tripack)

set.seed(43)

# test data
n <- 25
X <- matrix(rnorm(2 * n), ncol = 2)

vm <- voronoi.mosaic(xy.coords(X))
x <- with(vm, rowMeans(cbind(X[p1, 1], X[p2, 1], X[p3, 1])))
y <- with(vm, rowMeans(cbind(X[p1, 2], X[p2, 2], X[p3, 2])))
cen <- apply(cbind(x, y), 2, weighted.mean, vm$area)  # centroid of conv hull

tri <- tri.mesh(xy.coords(X)) # triangularization
ch <- convex.hull(tri) # ch$i gives indexes of vertices of conv hull

# plot points & Delauney triangularization, conv hull (red) and centroid (red)
plot(tri)
lines(X[c(ch$i, ch$i[1]), ], col = "red", lwd = 2)
points(cen[1], cen[2], col = "red", pch = 20, cex = 2)

（剧情后续）

删除目录

tripack 的替代方法是使用 deldir 包来获取质心。 deldir 函数在其输出中提供del.wts，这是对每个输入点进行加权的数量，使得它们的加权平均值是质心。 cen.dd 等于上面cen 的浮点近似值，情节也相似。 chull 来自 R 的基础。

library(deldir)

dd <- deldir(xy.coords(X))
cen.dd <- with(dd$summary, sapply(list(x, y), weighted.mean, del.wts))

# plot points & triangularization, compute & plot conv hull (red) and 
#  centroid (red)
plot(dd, wlines = "triang")
ix <- chull(X)
lines(X[c(ix, ix[1]), ], col = "red", lwd = 2)
points(cen.dd[1], cen.dd[2], col = "red", pch = 20, cex = 2)

更新

已经完全修改了答案。

Answer 2

正如G. Grothendieck 指出的那样，多边形的质心不是顶点的平均值。凸多边形质心的公式并不难，我还没有看到它们发布：

set.seed(42)
pts <- matrix(rnorm(40, 5, 1), 20, 2)
plot(pts)
# The centroid of the data points
points(t(pts.mn), pch=8, col="darkgreen", cex=2)
verts <- chull(pts)
poly <- pts[verts,]
polygon(poly)
cent <- colMeans(poly)
# The centroid of the vertices
points(t(cent), pch=8, col="blue", cex=2)

现在是面积和质心的公式。这些假设多边形是闭合的（第一行与最后一行相同）并且点逆时针遍历多边形：

poly <- rbind(poly, poly[1, ])
n <- nrow(poly)
x <- rev(poly[, 1])
y <- rev(poly[, 2])
i <- 1:(n-1)
# Area of the polygon
A <- sum(c(x[i] * y[i+1] - x[i+1] * y[i])) / 2
# Coordinates of the centroid
Cx <- sum((x[i] + x[i+1]) * (x[i] * y[i+1] - x[i+1] * y[i])) / (6 * A)
Cy <- sum((y[i] + y[i+1]) * (x[i] * y[i+1] - x[i+1] * y[i])) / (6 * A)
# The centroid of the polygon
points(Cx, Cy, pch=8, col="red", cex=2)