快得惊人的散度定理体积计算

# 利用散度定理实现超快速体积计算 来源：https://alyssarosenzweig.ca/blog/hilariously-fast-volume-computation-with-the-divergence-theorem.html 2018年2月16日（注意：本文不含笑话） 下文提出一种针对简单闭合三角化三维网格的快速体积计算算法。该算法基于散度定理推导得出，理论上可扩展至其他网格类型，但这超出当前讨论范围。 我们首先将体积定义为区域上的常数1三重积分： $$ V = \iiint_R 1 \, \mathrm{d}V $$ 设向量场 $$ \mathbf{F} $$ 为三维空间 $$ \mathbb{R}^3 $$ 上的函数，且其散度为1。为推导方便，本文定义： $$ \mathbf{F}(x, y, z) = \begin{pmatrix} x \\ 0 \\ 0 \end{pmatrix} $$ 易证其散度： $$ \mathrm{div} \, \mathbf{F} = \frac{\partial F}{\partial x} + \frac{\partial F}{\partial y} + \frac{\partial F}{\partial z} = 1 + 0 + 0 = 1 $$ 因此体积可改写为： $$ V = \iiint_R 1 \, \mathrm{d}V = \iiint_R \mathrm{div} \, \mathbf{F}(x, y, z) \, \mathrm{d}V $$ 根据散度定理，该体积分等于闭合曲面S上的面积分： $$ V = \iint_S \mathbf{F}(x, y, z) \, \mathrm{d}\mathbf{S} $$ 此面积分定义在三维网格的表面S上，可分解为各三角形面片的贡献之和。设第i个三角形的表面为 $$ T_i $$，则： $$ V = \sum_{i=0} \iint_{T_i} \mathbf{F}(x, y, z) \, \mathrm{d}\mathbf{S} $$ 设三角形 $$ T_i $$ 的第n个顶点为 $$ T_{in} $$，定义向量： $$ \Delta_1 = T_{i1} - T_{i0}, \quad \Delta_2 = T_{i2} - T_{i0} $$ 每个三角形 $$ T_i $$ 可参数化为： $$ \mathbf{r}(u, v) = T_{i0} + u\Delta_1 + v\Delta_2 $$ 求偏导得： $$ \mathbf{r}_u = \Delta_1, \quad \mathbf{r}_v = \Delta_2 $$ 叉乘结果为： $$ \mathbf{r}_u \times \mathbf{r}_v = \Delta_1 \times \Delta_2 $$ 将参数化表达式代入面积分，结合 $$ \mathbf{F} $$ 的定义可得： $$ V = \sum_{i=0} \iint_{T_i} \mathbf{F}(x, y, z) \cdot (\mathbf{r}_u \times \mathbf{r}_v) \, \mathrm{d}A $$ $$ = \sum_{i=0} \iint_{T_i} \begin{pmatrix} x \\ 0 \\ 0 \end{pmatrix} \cdot (\Delta_{i1} \times \Delta_{i2}) \, \mathrm{d}A $$ $$ = \sum_{i=0} \iint_{T_i} (\Delta_{i1} \times \Delta_{i2})_x \cdot x \, \mathrm{d}A $$ 由于叉乘向量在三角形区域内为常量，且与 $$ \mathbf{F} $$ 的点积仅保留x分量，可提出积分号外： $$ V = \sum_{i=0} (\Delta_{i1} \times \Delta_{i2})_x \iint_{T_i} x \, \mathrm{d}A $$ 现聚焦于面积分 $$ \iint_{T_i} x \, \mathrm{d}A $$。利用参数展开： $$ \iint_{T_i} x \, \mathrm{d}A = \int_0^1 \int_0^{1-u} \left( T_{i0x} + u\Delta_{i1x} + v\Delta_{i2x} \right) dv \, du $$ 直接计算二重积分（视顶点坐标为常数）： $$ = T_{i0x} \int_0^1 \int_0^{1-u} dv \, du + \Delta_{i1x} \int_0^1 \int_0^{1-u} u \, dv \, du + \Delta_{i2x} \int_0^1 \int_0^{1-u} v \, dv \, du $$ $$ = T_{i0x} \cdot \frac{1}{2} + \Delta_{i1x} \cdot \frac{1}{6} + \Delta_{i2x} \cdot \frac{1}{6} $$ 代回向量定义 $$ \Delta_{ix} = T_{i,x} - T_{i0x} $$ 并化简： $$ = \frac{1}{6} \left( T_{i0x} + T_{i1x} + T_{i2x} \right) $$ 最终得到体积的紧凑计算公式： $$ V = \frac{1}{6} \sum_{i=0} (\Delta_{i1} \times \Delta_{i2})_x \cdot \left( T_{i0x} + T_{i1x} + T_{i2x} \right) $$ ## 性能分析 该算法完全规避了数值积分与微分。相较于传统网格体积计算（通常需渲染网格后进行采样），本算法仅需单次三角形遍历，时间复杂度为 $$ O(n) $$。 单次三角形计算仅需7次加法和3次乘法，循环外还有1次乘法。对于包含 $$ n $$ 个三角形的网格，总运算量为 $$ 8n-1 $$ 次加法、$$ 3n+1 $$ 次乘法，即 $$ 11n $$ 次浮点运算。 此效率极高。以基准估算：若需在60fps高性能应用中实时计算体积（仅依靠CPU，不依赖GPU），即便是35美元的树莓派（性能参考：https://raspberrypi.stackexchange.com/questions/55862/what-is-the-performance-and-the-performance-per-watt-of-raspberry-pi-3-in-gflops），每帧仍可处理约3000万三角形。 ## 动机 向量微积分考试临近需要复习，加之谁不爱3D图形呢？ 若该算法具有创新性，那将是（令人惊喜的）。但后续查阅发现Cha Zheng与Tsuhan Chen的论文《Efficient Feature Extraction for 2D/3D Objects in Mesh Representation》（http://chenlab.ece.cornell.edu/Publication/Cha/icip01_Cha.pdf）已描述类似算法，只是推导路径不同。 推导过程很有趣，但文献先已存在。 [返回主页](https://alyssarosenzweig.ca/)