If you find this article contains errors or problems rendering it unreadable (missing images or files, mangled code, improper text formatting, etc) please contact the editor so corrections can be made. Thank you for helping us improve this resource |
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
where the path l encloses the area a. The dell vector is defined by:
dell = (∂/∂x, ∂/∂y)
When F is set to:
F(x, y) = (x, 0)
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
which equates to the area of a because:
dell ∙ F = x ∙ ∂/∂x + 0 ∙ ∂/∂y = ∂x/∂x = 1
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 {v_{1}, v_{2}} is given by:
e = v_{2} – v_{1}
l(t) = v_{1} + te
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 = (-e_{y}, e_{x})dt
The flux of F through the edge {v_{1}, v_{2}} now follows:
Φ = ∫ F(l(t)) ∙ (-e_{y}, e_{x})dt
Φ = ∫ -e_{y}(v_{1}_{x} + te_{x})dt
where the integral is evaluated through t on [0, 1]. The analytical solution to equation 1.10 can be written as:
Φ = ½ ∙ (v_{1}_{y} – v_{2}_{y}) ∙ (v_{1}_{x} +v_{2}_{x})
The right hand side of equation 1.11 yields the flux of F through the edge {v_{1}, v_{2}}. Given a closed and clockwise wound edge set {v_{i1},
v_{i2}} for i = {0, 1, 2,…,n} the enclosed area is computed by evaluating the following sum:
area = ½ ∑ (v_{i1}_{y} – v_{i2}_{y}) ∙ (v_{i1}_{x} +v_{i2}_{x})
C++ Implementation (2D)
<span class="codecomment">// 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.</span> <span class="codekeyword">float</span> CalculateArea(<span class="codekeyword">const</span> VECTOR2 *v, <span class="codekeyword">const unsigned int</span> *i, <span class="codekeyword">const unsigned int</span> n) { <span class="codekeyword">unsigned int</span> j; VECTOR2 v1; VECTOR2 v2; <span class="codekeyword">float</span> area = 0.0f; <span class="codekeyword">for</span>(j = 0; j < n; j+=2) { v1 = v[i[j ]]; v2 = v[i[j+1]]; area += (v1.y-v2.y) * (v1.x+v2.x); } <span class="codekeyword">return</span> 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
where the surface s encloses the volume v. The dell vector is defined by:
dell = (∂/∂x, ∂/∂y, ∂/∂z)
When F is set to:
F(x, y, z) = (x, 0, 0)
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
which equates to the volume of v because:
dell ∙ F = x ∙ ∂/∂x + 0 ∙ ∂/∂y + 0 ∙ ∂/∂z = ∂x/∂x = 1
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 {v_{1}, v_{2},
v_{3}} is given by:
e_{1} = v_{2} – v_{1}
e_{2} = v_{3} – v_{1}
s(u, v)= v_{1} + ue_{1} + ve_{2}
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 = (e_{1} cross e_{2})dvdu
The flux of F through the triangle {v_{1}, v_{2}, v_{3}} now follows:
Φ = ∫∫F(s(u, v)) ∙ (e_{1} cross e_{2})dvdu
Φ = ∫∫ (v_{1x} + ue_{1x} + ve_{2x})(e_{1y}e_{2z} – e_{1z}e_{2y})dvdu
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 ((v_{2}_{y}–v_{1}_{y})(v_{3}_{z}–v_{1}_{z}) –
(v_{2}_{z}–v_{1}_{z})(v_{3}_{y}–v_{1}_{y}))(v_{1}_{x} +
v_{2}_{x} + v_{3}_{x})
The right hand side of equation 2.12 yields the flux of F through the triangle {v_{1}, v_{2}, v_{3}}. Given a closed and clockwise wound
triangle set {v_{i1}, v_{i2}, v_{i3}} for i = {0, 1, 2,…,n} the enclosed volume is computed by evaluating the following sum:
volume = 1/6 ∑ ((v_{i2}_{y}–v_{i1}_{y})(v_{i3}_{z}–v_{i1}_{z}) –
(v_{i2}_{z}–v_{i1}_{z})(v_{i3}_{y}–v_{i1}_{y}))(v_{i1}_{x} +
v_{i2}_{x} + v_{i3}_{x})
C++ Implementation (3D)
<span class="codecomment">// 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).</span> <span class="codekeyword">float</span> CalculateVolume(<span class="codekeyword">const</span> VECTOR3 *v, <span class="codekeyword">const unsigned int</span> *i, <span class="codekeyword">const unsigned int</span> n) { <span class="codekeyword">unsigned int</span> j; VECTOR3 v1; VECTOR3 v2; VECTOR3 v3; <span class="codekeyword">float</span> volume = 0.0f; <span class="codekeyword">for</span>(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); } <span class="codekeyword">return</span> 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.