Rosenzweig – Hilariously Fast Volume Computation with the Divergence Theorem

(No,
there won’t be jokes.)

The following presents a fast algorithm for volume computation of a
simple, closed, triangulated 3D mesh. This assumption is a consequence
of the divergence theorem. Further extensions may generalise to other
meshes as well, although that is presently out of scope.

We begin with the definition of volume as the triple integral over a
region of the constant one:

$$V = \iiint_R 1 \mathrm⚡V$$

Let
$\mathbf🔥$
be a function in
$\mathbb🔥^3$
such that its divergence is equal to one. For the purposes of this
paper, we choose:

$$\mathbf💬(x, y, z) =$$

It can easily be verified that

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

Therefore,

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

By the Divergence Theorem, this is equal to the surface integral:

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

This surface integral, defined over the surface S of the 3D mesh, is
equal to the sum of its piecewise triangle parts. Let
$T_i$
denote the surface of the
$i$’th
triangle in the mesh. Then,

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

Let
$T_{in}$
represent the
$n$’th
vertex of the
$i$’th
triangle. Let
$\Delta_1$
equal the vector difference between
$T_{i1}$
and
$T_{i0}$,
and
$\Delta_2$
likewise equal to
$T_{i2} – T{i0}$.
Each individual triangle
$T_i$
may thus be parametrised as:

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

Then, simple differentiation yields:

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

Therefore,

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

Thus, the surface integral can be rewritten in terms of this
parametrisation, substituting in the definition of
$\mathbf{F}$
as needed:

$$V = \sum_{i = 0} \iint_{T_i} \mathbf{F}(x, y, z) (\mathbf{r}_u \times \mathbf{r}_v) dA$$
$$= \sum_{i = 0} \iint_{T_i} \mathbf{F}(x, y, z) \dot (\Delta_{i1} \times \Delta_{i2}) dA$$
$$= \sum_{i = 0} \iint_{T_i} \dot (\Delta_{i1} \times \Delta_{i2}) dA$$

This cross product is constant throughout the triangle and easy to
calculate from the vertex data. Only the X component of the cross
product should be calculated; the others are equal to zero due to the
dot product with the zero components of
$\mathbf{F}$.
$V$
can be thus be rewritten as:

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

We now focus on the surface integral
$\iint_{T_i} x dA$.
Expanding with the parametrisation yields:

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

This integral can be directly evaluated, treating vertex data as
constants:

$$\int_{0}^{1} \int_{0}^{1-u} (T_{i0x} + u \Delta_{i1x} + v \Delta_{i2x}) 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} (\frac{1}{2}) + \Delta_{i1x} (\frac{1}{6}) + \Delta_{i2x} (\frac{1}{6})$$
$$= T_{i0x} (\frac{1}{2}) + (T_{i1x} – T_{i0x})(\frac{1}{6}) + (T_{i2x} – T_{i0x})(\frac{1}{6})$$
$$= T_{i0x} (\frac{1}{6}) + (T_{i1x})(\frac{1}{6}) + (T_{i2x})(\frac{1}{6})$$
$$= \frac{1}{6}(T_{i0x} + T_{i1x} + T_{i2x})$$

Substituting into the original sum and pulling out a constant factor
of
$\frac{1}{6}$
to avoid the inner loop division, this yields the following compact
formula for the volume:

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

Performance analysis

The final algorithm contains no numerical integration nor
differentiation. In contrast to common naive algorithms for volume,
which are equivalent to rendering the mesh and then sampling the render,
an expensive operation, there is only a single loop in this algorithm,
over the triangles. Thus, this algorithm for volume computation is O(n)
to the number of the triangles. Furthermore, the per-triangle
calculation is similarly efficient: given the natural expansion of the
cross product, the inner part contains seven additions and three
multiplications. On the outside of the loop is only a single
multiplication. Thus, for a mesh of
$n$
triangles, the algorithm requires
$8n – 1$
additions and
$3n + 1$
multiplications, or
$11n$
floating point operations. This is very fast.

For a ballpark number, if volume needs to be calculated every frame
in a high-performance 60 frames per second application, without the aid
of a GPU, only using the CPU capabilities of a $35
Raspberry Pi, around 30 million triangles could be measured every
frame.

Motivation

The vector calculus exam is soon, and I need to study. Plus, who
doesn’t love 3D graphics?!

I would be (pleasantly) surprised if the algorithm is
novel. Further research after posting reveals the paper
Efficient
Feature Extraction for 2D/3D Objects in Mesh Representation by Cha
Zheng and Tsuhan Chen, which appears to describe the same algorithm,
although the derivation is different. It was fun while it lasted!

Retourner à l’acceuil