从多变形面积到多面体体积：鞋带公式的3D版本_鞋带公式推导

怎样根据多边形顶点坐标计算多边形面积，以及3D情况下怎样根据多面体的顶点坐标计算多面体的体积？

多边形的面积计算方法中，有一条像变戏法一样的计算方式：鞋带定理，可以直接根据多边形的各个顶点坐标得到多变形的面积。本文将会介绍鞋带公式的原理，以及从2D扩展到3D的方式，即使用鞋带公式计算3D多面体的体积。

2D情况下的鞋带公式可以写作：（其中 i 是关于多边形的各个顶点求和）

A = ∑ i v e r t e x e s 1 2 ( x i y i + 1 − x i + 1 y i ) = ∑ i v e r t e x e s 1 2 ∣ x i x i + 1 y i y i + 1 ∣ = ∑ i v e r t e x e s 1 2 d e t ( x i ⃗ ,   x i + 1 ⃗ ) A = \sum_{i}^{vertexes}{\frac{1}{2}(x_{i}y_{i+1}-x_{i+1}y_{i})} = \sum_{i}^{vertexes}{ \frac{1}{2}\left|

$$\begin{matrix}x_{i} & x_{i+1}\\y_{i} & y_{i+1} \end{matrix}$$ \right| } = \sum_{i}^{vertexes}{\frac{1}{2} det(\vec{x_{i}},\ \vec{x_{i+1}})} A=i∑vertexes​21​(xi​yi+1​−xi+1​yi​)=i∑vertexes​21​∣ ∣​xi​yi​​xi+1​yi+1​​∣ ∣​=i∑vertexes​21​det(xi​ ​, xi+1​ ​)

图1：任意多边形

其中 $\overrightarrow{x_i}$ 表示 $i$ 节点的位置矢量。以上的求和过程从 $i$ 到 $i+1$ 必须保持逆时针或者顺时针的顺序，如上图所示，否则会计算出错误的结果。

鞋带公式可以看有两种直观的理解方式：

• 格林公式的特例（积分外围边界是多边形的情况下，可以从格林公式推导出鞋带公式）
• 三角形面积的加减 （没错！这是一种神奇的理解方式！鞋带公式可以表示为三角形的加加减减，而且几何上很直观，不需要多少推导操作，而且后面可以更快推导到三维体积的情况。下面来看看为什么）

图2：鞋带公式计算多边形面积的三角形解释（从坐标原点出发，指向各个顶点的矢量，与原来的边长矢量构成三角形，三角形的面积的正负号由叉乘顺序决定，正负三角面积的叠加得到多边形面积）

如上图所示，表示了鞋带公式为什么相当于用三角形面积的加减得到多边形的面积（从坐标原点出发，指向各个顶点的矢量，与原来的边长矢量构成三角形），三角形的面积可以用叉乘计算，分别是$\frac{1}{2}O{P}_1\times O{P}_2,\frac{1}{2}O{P}_2\times O{P}_3,...$，这些三角形的面积有正有负，累加之后就是多边形的面积了！

多面体的体积：二维的鞋带公式推广到三维的情况

从二维的鞋带公式推广到三维，可以在给定多面体的顶点的坐标的情况下计算多面体的体积。这里把问题简化为一个完全由三角形面片构成物体外表面的多面体（如果一个面不是三角形而是任意多边形，那么把这个面分割成三角形就行了），于是问题变成：怎样计算一个由一片片三角形构成外表面的物体的体积？（给定各三角形的顶点坐标的情况下）。

类比到二维的情况，从原点出发指向各个顶点构成三角形，通过三角形面积的加加减减可以得到多边形的面积；于是在三维的时候，可以从原点出发指向各个顶点构成四面体，通过四面体体积的加加减减就可以得到多面体的体积。（类比图2，把原来连线构成的三角形想象成四面体）

首先给出四面体的体积公式：

图3：四面体的体积公式

如上图所示，对于四面体OABC, 其体积可以表示为$$V_{OABC}=\frac{1}{6}(\overrightarrow{OA}\times \overrightarrow{OB})\cdot \overrightarrow{OC}=\frac{1}{6}\left|\overrightarrow{OA}\overrightarrow{OB}\overrightarrow{OC}\right|=\frac{1}{6}det(\overrightarrow{OA},\overrightarrow{OB},\overrightarrow{OC})$$ 其中，$(\overrightarrow{OA}\times \overrightarrow{OB})\cdot \overrightarrow{OC}$表示由这三个矢量构成的平行六面体的体积，而四面体则是这个平行六面体体积的1/6（乘以1/6的原因：锥体体积可以看作底乘以高，而底部三角形面积是平行六面体底面积的1/2, 乘以高之后，由于是锥体，体积还得再乘以1/3, 所以 1/2 x 1/3 = 1/6)

所以，三角形面片构成的多面体的体积，可以通过一系列的从原点出发构造的四面体的体积进行加加减减得到（即3D情况下的鞋带公式，原理类似于图2）：

$$\begin{gathered} V=\frac{1}{6}\sum _i^{facets}det(\overrightarrow{r_{i1}},\overrightarrow{r_{i2}},\overrightarrow{r_{i3}})=\frac{1}{6}\sum _i^{facets}\left|\overrightarrow{r_{i1}}\overrightarrow{r_{i2}}\overrightarrow{r_{i3}}\right| \\ facetnormalpointstotheoutsideofthesurface \end{gathered}$$

其中，上式的 $i$ 表示不同的三角形面片的编号，1， 2， 3 表示这个三角形三个顶点的编号，三个顶点的编号顺序必须按照右手螺旋定则指向物体表面朝外的方向，必须所有所有面片都指向朝外的方向（或者所有面片朝里指向物体内部，不可以某些面朝外某些面朝里），这样算出来的四面体体积才可以正确地把正体积加上负体积之后得到多面体的体积。

经过测试，这个公式可以处理 stl 格式的文件得到物体体积，即根据 stl 格式的文件中各个三角形面片的顶点坐标,以及右手定则指向物体外部的三角形节点排列方式，计算出这个物体的体积。（而对于由多边形构成的多面体，把每个面的多边形分割成三角形之后就可以用同样的方法计算）

下面再贴上导入 stl 格式的文件之后计算物体体积的python代码（stl 格式是其中一种表示3d物体的标准文件格式）

import numpy as np

def stlVolume(dataFile):
    facets = []
    with open (dataFile, "r") as file:
        readIn = False
        for line in file:
            if "endloop" in line:
                readIn = False
                continue
            if readIn:
                facets.extend(
                    list(map(float, line.split()[1:]))
                )
            if "outer loop" in line:
                readIn = True
    facets = np.array(facets).reshape((-1, 3, 3))

    ### 3D version of shoelace formular
    return sum(np.linalg.det(x) for x in facets) / 6.  


if __name__ == "__main__":
    computeVolume = "y"
    while computeVolume in ["y", ""]:
        dataFile = input("\033[32;1m please input the file name (include the path): \033[0m")
        volume = stlVolume(dataFile)
        print("volume = \033[40;33;1m {} \033[0m".format(volume))
        computeVolume = input("\033[35;1m continue to next file? (y/n): \033[0m")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

参考资料

https://en.wikipedia.org/wiki/Shoelace_formula