Hilariously Fast Volume Computation with the Divergence Theorem (2018)

Author: Rosenzweig
Date: February 16, 2018

Warning: This article contains high-level humor. Correction: No, there actually won't be any jokes.

This paper outlines a highly efficient algorithm designed to calculate the volume of a simple, closed, 3D mesh composed of triangles. While the current focus is on triangulated meshes, the logic could potentially be extended to other mesh types in the future.

---

1. Mathematical Foundation

The starting point is the standard definition of volume $V$ as a triple integral over a region $R$ where the integrand is a constant one:

$V = \iiint_R 1 \mathrm{d}V$

To simplify this, we introduce a vector field $\mathbf{F}$ in $\mathbb{R}^3$ such that its divergence is exactly one. For this specific implementation, we define:

F(x, y, z) = (x, 0, 0)

We can verify the divergence $\text{div} \mathbf{F}$ as follows:

$\mathrm{div} \mathbf{F} = \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z} = 1 + 0 + 0 = 1$

By substituting this back into our volume equation, we get:

$V = \iiint_R \mathrm{div} \mathbf{F}(x, y, z) \mathrm{d}V$

Applying the Divergence Theorem, we can transform this volume integral into a surface integral over the boundary $S$:

$V = \iint_S \mathbf{F}(x, y, z) \cdot \mathrm{d}\mathbf{S}$

Decomposing the Mesh

Since the 3D mesh is a collection of piecewise triangles, the total volume is simply the sum of the integrals over each individual triangle $T_i$:

$V = \sum_{i = 0} \iint_{T_i} \mathbf{F}(x, y, z) \cdot \mathrm{d}\mathbf{S}$

---

2. Triangle Parametrization

To solve the integral for a single triangle $T_i$, we define its vertices as $T_{i0}, T_{i1},$ and $T_{i2}$. We then establish two edge vectors:

• $\Delta_1 = T_{i1} - T_{i0}$
• $\Delta_2 = T_{i2} - T_{i0}$

The surface of the triangle can be parametrized by $\mathbf{r}(u, v)$:

$\mathbf{r}(u, v) = T_{i0} + u\Delta_1 + v\Delta_2$

From this, we derive the partial derivatives:

• $\mathbf{r}_u = \Delta_1$
• $\mathbf{r}_v = \Delta_2$

The cross product $\mathbf{r}_u \times \mathbf{r}_v = \Delta_1 \times \Delta_2$ represents the normal vector, which remains constant across the triangle's surface.

Simplifying the Surface Integral

Substituting our parametrization and the definition of $\mathbf{F}$ into the integral:

$V = \sum_{i = 0} \iint_{T_i} (x, 0, 0) \cdot (\Delta_{i1} \times \Delta_{i2}) dA$

Because $\mathbf{F}$ only has an $x$-component, only the $x$-component of the cross product $(\Delta_{i1} \times \Delta_{i2})_x$ contributes to the result. The formula simplifies to:

$V = \sum_{i = 0} (\Delta_{i1} \times \Delta_{i2})_x \iint_{T_i} x dA$

---

3. Solving the Integral of $x$

We now evaluate the remaining integral $\iint_{T_i} x dA$ using the parametrization:

$\iint_{T_i} x dA = \int_{0}^{1} \int_{0}^{1-u} (T_{i0x} + u \Delta_{i1x} + v \Delta_{i2x}) dv du$

By treating the vertex coordinates as constants and integrating, we find:

$= T_{i0x} \left(\frac{1}{2}\right) + \Delta_{i1x} \left(\frac{1}{6}\right) + \Delta_{i2x} \left(\frac{1}{6}\right)$

Substituting $\Delta_{i1x} = (T_{i1x} - T_{i0x})$ and $\Delta_{i2x} = (T_{i2x} - T_{i0x})$:

$= T_{i0x} \left(\frac{1}{2}\right) + (T_{i1x} - T_{i0x})\left(\frac{1}{6}\right) + (T_{i2x} - T_{i0x})\left(\frac{1}{6}\right)$ $= \frac{1}{6}(T_{i0x} + T_{i1x} + T_{i2x})$

Final Volume Formula

To optimize performance and avoid division inside the loop, we factor out $1/6$:

$V = \frac{1}{6} \sum_{i = 0} (\Delta_{i1} \times \Delta_{i2})_x (T_{i0x} + T_{i1x} + T_{i2x})$

---

4. Implementation and Performance

Workflow Visualization

Complexity Analysis

Unlike naive sampling methods (which involve rendering the mesh and sampling pixels—a very costly process), this algorithm is strictly linear.

Hardware Example: On a budget-friendly $35 Raspberry Pi (CPU only), this algorithm could process roughly 30 million triangles per frame while maintaining a 60 FPS target.

Implementation Checklist

• Ensure mesh is closed (watertight).
• Compute $\Delta_{i1}$ and $\Delta_{i2}$ for every triangle.
• Extract the $x$-component of the cross product.
• Sum the $x$-coordinates of the three vertices.
• Apply the final $1/6$ scalar.

---

5. Conclusion & Motivation

Why was this written?

"The vector calculus exam is soon, and I need to study."

Further Reading: After publication, it was noted that similar logic is described in the paper "Efficient Feature Extraction for 2D/3D Objects in Mesh Representation" by Cha Zheng and Tsuhan Chen.

Summary Code Snippet

double calculateVolume(Mesh mesh) {
    double totalSum = 0;
    for (auto& tri : mesh.triangles) {
        Vector3 d1 = tri.v1 - tri.v0;
        Vector3 d2 = tri.v2 - tri.v0;
        double crossX = (d1.y * d2.z) - (d1.z * d2.y);
        double sumX = tri.v0.x + tri.v1.x + tri.v2.x;
        totalSum += crossX * sumX;
    }
    return totalSum / 6.0;
}

Figure 1: Conceptual representation of converting a volume integral to a surface sum.