Motivation

In recent months, the indie game programming scene has witnessed the rising popularity of a pressure volume model described in a paper written by Maciej Matyka and Mark Ollila. The paper can be found here: http://panoramix.ift.uni.wroc.pl/~maq/eng/index.php. In the model's first presentation the authors used an axis aligned bounding box to approximate the volume of a body. Such a rough approximation is easily computed but highly inaccurate. A second paper followed sometime later describing a two dimensional implementation that employed Gauss's Theorem to calculate the exact area of a simulated body. The author mentioned the possibility of a three dimensional version using the same technique but restricted his presentation to the two dimensional case. It would seem that a three dimensional volume calculator should have surfaced shortly afterward; however, no such algorithm has been detailed by Matyka or Ollila at the time of this writing. The goal of this article is to fill in the details left out of Matyka's second paper.

Problem Space

The two dimensional case: Given a set of edges that form a closed loop calculate the enclosed area.

The three dimensional case: Given a set of triangles that form a closed surface calculate the enclosed volume.

Derivation (2D)

Gauss's Theorem relates the divergence of a vector field within an area to the flux of a vector field through a closed path by the following:

??a dell ?F da = ?l F ? dn
equation 1.1

where the path l encloses the area a. The dell vector is defined by:

dell = (?/?x, ?/?y)
equation 1.2

When F is set to:

F(x, y) = (x, 0)
equation 1.3

the innards of the double integral in equation 1.1 can be reduced to one. This reduction lets us write:

??a 1 da = ?l F ? dn
equation 1.4

which equates to the area of a because:

dell ? F = x ? ?/?x + 0 ? ?/?y = ?x/?x = 1
equation 1.5

Equation 1.4 tells us that the flux of the vector field F through the closed path l is equal to the area enclosed by l. When l is defined by a set of edges the enclosed area is equal to the sum of the flux of F through every edge. A straight forward parameterization of the edge {v1, v2} is given by:

e = v2 - v1
equation 1.6

l(t) = v1 + te
equation 1.7

where the evaluation of l(t) yields any point on the edge for t on [0, 1]. Under this parameterization dn from equation 2.4 can be written as:

dn = cross dl = cross e dt = (-ey, ex)dt
equation 1.8

The flux of F through the edge {v1, v2} now follows:

? = ? F(l(t)) ? (-ey, ex)dt
equation 1.9

? = ? -ey(v1x + tex)dt
equation 1.10

where the integral is evaluated through t on [0, 1]. The analytical solution to equation 1.10 can be written as:

? = 1/2 ? (v1y - v2y) ? (v1x +v2x)
equation 1.11

The right hand side of equation 1.11 yields the flux of F through the edge {v1, v2}. Given a closed and clockwise wound edge set {vi1, vi2} for i = {0, 1, 2,...,n} the enclosed area is computed by evaluating the following sum:

area = 1/2 ? (vi1y - vi2y) ? (vi1x +vi2x)
equation 1.12

C++ Implementation (2D)

Listing 1: Area Calculation:

// v: A pointer to the array of vertices
// i: A pointer to the array of indices
// n: The number of indices (multiple of 2)
// This function uses Gauss's Theorem to calculate the area (or field)
// of a body enclosed by a set of edges. The edge set must form a
// closed path in R2 and be outward facing. An outward facing edge set
// wraps clockwise around the body where the x-axis points left and the
// y-axis points up.
float CalculateArea(const VECTOR2      *v,
                    const unsigned int *i,
                    const unsigned int  n)
{
   unsigned int j;
   VECTOR2      v1;
   VECTOR2      v2;
   float        area = 0.0f;

   for(j = 0; j < n; j+=2)
   {
      v1 = v[i[j ]];
      v2 = v[i[j+1]];

      area += (v1.y-v2.y) * (v1.x+v2.x);
   }

   return area / 2.0f;
}

Derivation (3D)

Gauss's Theorem relates the divergence of a vector field within a volume to the flux of a vector field through a closed surface by the following:

???v dell ? F dv = ??s F ? da
equation 2.1

where the surface s encloses the volume v. The dell vector is defined by:

dell = (?/?x, ?/?y, ?/?z)
equation 2.2

When F is set to:

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

the innards of the triple integral in equation 1.1 can be reduced to one. This reduction lets us write:

???v 1 dv = ??s F ? da
equation 2.4

which equates to the volume of v because:

dell ? F = x ? ?/?x + 0 ? ?/?y + 0 ? ?/?z = ?x/?x = 1
equation 2.5

Equation 2.4 tells us that the flux of the vector field F through the closed surface s is equal to the volume enclosed by s. When s is defined by a set of triangles the enclosed volume is equal to the sum of the flux of F through every triangle. A straightforward parameterization of the triangle {v1, v2, v3} is given by:

e1 = v2 - v1
equation 2.6

e2 = v3 - v1
equation 2.7

s(u, v)= v1 + ue1 + ve2
equation 2.8

where the evaluation of s(u, v) yields any point on the triangle for u on [0, 1] and v on [0, 1-u]. Under this parameterization da from equation 2.4 can be written as:

dn = (e1 cross e2)dvdu
equation 2.9

The flux of F through the triangle {v1, v2, v3} now follows:

? = ??F(s(u, v)) ? (e1 cross e2)dvdu
equation 2.10

? = ?? (v1x + ue1x + ve2x)(e1ye2z - e1ze2y)dvdu
equation 2.11

where the integral is evaluated through u on [0, 1] and v on [0, 1-u]. The analytical solution to equation 2.11 can be written as:

? = 1/6 ((v2y-v1y)(v3z-v1z) - (v2z-v1z)(v3y-v1y))(v1x + v2x + v3x)
equation 2.12

The right hand side of equation 2.12 yields the flux of F through the triangle {v1, v2, v3}. Given a closed and clockwise wound triangle set {vi1, vi2, vi3} for i = {0, 1, 2,...,n} the enclosed volume is computed by evaluating the following sum:

volume = 1/6 ? ((vi2y-vi1y)(vi3z-vi1z) - (vi2z-vi1z)(vi3y-vi1y))(vi1x + vi2x + vi3x)

equation 2.13

C++ Implementation (3D)

Listing 2: Volume Calculation

// v: A pointer to the array of vertices
// i: A pointer to the array of indices
// n: The number of indices (multiple of 3)
// This function uses Gauss's Theorem to calculate the volume of a body
// enclosed by a set of triangles. The triangle set must form a closed
// surface in R3 and be outward facing. Outward facing triangles index
// their vertices in a counterclockwise order where the x-axis points
// left, the y-axis point up and the z-axis points toward you (rhs).
float CalculateVolume(const VECTOR3      *v,
                      const unsigned int *i,
                      const unsigned int  n)
{
   unsigned int j;
   VECTOR3      v1;
   VECTOR3      v2;
   VECTOR3      v3;
   float        volume = 0.0f;

   for(j = 0; j < n; j+=3)
   {
      v1 = v[i[j ]];
      v2 = v[i[j+1]];
      v3 = v[i[j+2]];

      volume += ((v2.y-v1.y)*(v3.z-v1.z)-
                 (v2.z-v1.z)*(v3.y-v1.y) )*(v1.x+v2.x+v3.x);

   }

   return volume / 6.0f;
}

Conclusion

And there you have it, a fast and precise way to find the area or volume of a body given its edge or triangle set. I have tested the code on a variety of basic primitives and am 99% certain it is bug free. You can contact me with question, comments and bug reports at kwoxrf@umr.edu.