🎦 👨🏿‍🍳 👩🏻‍🏭 搜索多面体的体积和质心的算法〽️ 🤙🏻 👩🏾‍🚒

也许每个人都知道这种算法，但是“当局对我隐瞒了”。我在英语论坛的自动翻译档案中的搜索引擎第三页上找到了他的口头描述。在我看来，它的详细说明（以及代码）值得habrosta使用。

因此，例如，您需要为玩具和过程中的某个位置生成生物，以清除那些没有站起来的人。为此，您需要找到生物的质心（这几乎与找到其体积相同），并确保它位于生物腿上方的某处。

暴民是多面体，为简单起见，我们认为多面体仅由三角形组成（在内部算法中为高斯面积公式，因此您可以将其扩展为任何多面体，但是为什么...）。此外，多面体不应具有自相交的形状，并应限制封闭的体积，因为它适合于体面的多面体。

（好吧，那样）

一个小的UPD解释了为什么在KDPV上右侧暴民不好，但是左侧暴民好：
正确的图片不好，因为暴民会往前掉，因为它的质心超出了支撑区域。立于表面上的多边形的支撑区域定义为表面上所有点都位于其中的最小多边形。在左图的情况下，支架的区域移到质心中心（因为恐龙的爪子较大）而变大了，在右图的情况下，支架的区域本身变小并且更靠近尾巴。
参考面积与质心的比率将如下所示：

我将立即开始使用体积搜索代码（Python，输入-点列表和过渡矩阵）：

一些代码

def RecSetDirsTriangles(para, Connects, TR): """ ,           """ #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 #    break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def FindV(dots, Connects): """ .   - dots     [x, y, z], Connects -  , Connects[i][j]=1      i, j,  =0 """ #1.      TR = [] for i in range(1,len(Connects)):#        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.        V = 0 for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... return math.fabs(V)

该算法的本质是考虑形成“下降”到xy平面上的多面体面的图形的体积。为此，您需要知道三角形的投影区域以及用于添加图形体积的符号（截短的棱镜）。实际上，如果预先订购了三角形，则体积和正负号都将减少为一次计算。

因此，输入递归函数收集三角形的第一件事。组装的方式是，从多面体“向外”看时，绕三角形的方向是相同的（理想情况是逆时针；如果顺时针方向，则结果是正确的，但为负-因此，将体积模数赋予返回值）。

这很容易实现-取一个三角形（点a1，a2，a3），寻找其邻居，并以相反的顺序列出两个匹配的顶点（例如，如下所示：a2，a1，b1）。
原来是这样的：

现在，如果将这样的三角形投影到xy平面上，则“上”三角形的投影的遍历顺序将与最初选择的三角形重合，而“下”三角形的投影的遍历顺序将更改其方向。结果，它将改变由高斯公式计算的三角形的符号和面积。在这里，“下部”三角形（一个有条件的概念）意味着紧接其下方的体积不包含在多面体的体积中。非凸多面体的“下”三角形可以高于“上”的三角形。

在完成这些初步步骤之后，为了计算多面体的总体积，您只需要添加（考虑符号，即“通过其自身获得”）从平面和这些平面在xy平面上的投影收集的所有截顶棱柱体的体积。棱镜的体积被认为是面积（高斯，带一个符号）与三角形顶点的算术平均z坐标的乘积。

如果多面体与xy平面相交，则在计算体积时，所有符号将相互抵消，结果仍然正确（您只需要在没有模块的情况下获取棱镜高度）。

（以某种方式，“上”截顶棱镜看起来像）

通过寻找质心，一切都差不多。同样，我们需要找到每个截顶棱柱的质心并乘以棱柱体的体积，以坐标方式进行汇总（假定质量在整个体积中均匀分布，并且一个可以被另一个代替）。为了找到截顶棱柱的质心，有必要计算两个四面体（+1函数）和一个普通棱柱的质心。如果多面体穿过xy平面，算法也不会“恶化”（这里可能是Magritte的复制品）。

（这两个分别标记为红色和红色的四面体与一个三角棱镜（在红色四面体下方）一起形成了所需的截顶棱镜。我们需要找到所有三个图形的质心和体积中心。这些名称大致对应于代码中的名称）

可以计算以下内容的代码：

多一点代码

 def RecSetDirsTriangles(para, Connects, TR): #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def TetrV(mas):#dot1, dot2, dot3, dot4): """   """ M = np.zeros((3,3),float) for i in range(1,4): for j in range(0,3): M[i-1][j] = mas[i][j] - mas[0][j] #print(M) return math.fabs(np.linalg.det(M)/6) def FindVandCM(dots, Connects): """     """ #1.      TR = [] for i in range(1,len(Connects)): #        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.   ,          V = 0 CM = [0, 0, 0] for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... #c        c1 = ((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0])/3, (dots[T[0]][1]+ dots[T[1]][1]+ dots[T[2]][1])/3) #    hm = min([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) hM = max([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) indM = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hM) indm = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hm) V3 = S * hm if indM == indm: # ! CM[0] += V3*c1[0] CM[1] += V3*c1[1] CM[2] += V3*hm/2 continue L = [0,1,2] L.remove(indM) L.remove(indm) indmidle = L[0] dots1 = [dots[T[0]], dots[T[1]], dots[T[2]], (dots[T[indM]][0], dots[T[indM]][1] , hm)] #  V1 = TetrV(dots1) if S < 0: V1 = -V1 V2 = S * ( dots[T[indmidle]][2] - hm)/3 #V3 = S * hm CM[0] += V1*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indM]][0])/4) + V2*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indmidle]][0])/4) + V3*c1[0] CM[1] += V1*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indM]][1])/4) + V2*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indmidle]][1])/4) + V3*c1[1] CM[2] += V1*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V2*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V3*hm/2 CM[0] = CM[0]/V CM[1] = CM[1]/V CM[2] = CM[2]/V return (math.fabs(V), CM)

考虑到三角形的方向并用来理解外部和内部体积的算法的一部分是非常有力的一步，在处理多面体时可以使用很多。例如，如果您需要计算法线“向外”的方向-只需知道一张脸的“逆时针”方向-瞧！

（猜电影！）

搜索多面体的体积和质心的算法

More articles: