Calculating the volume of an arbitrary mesh

Question

Given a mesh (vertices, edges, faces), I need to calculate the volume of the formed mesh. By volume, I mean the amount of space that a substance or object occupies. The resulting volume does not include any unoccupied space in or outside of the mesh.

• The given mesh may or may not be closed, that is;
  • The given mesh may or may not have holes on its surface
  • The given mesh may or may not be hollow/partially hollow
  • The given mesh may or may not be a primitive shape
  • The given mesh may or may not be composite

For example

This wireframe of a 5x5 cube has been "crushed" and has a partial hole in its center. I'm looking for the volume of this mesh without the unoccupied/hollow space.

Every vertex is defined by its 3D coordinates and normal, both of which are floating-point. Though, a 2D solution would work just as well with extension.

As described here and here, the volume of a closed mesh can be calculated fairly easily. My research hasn't brought any information on non-closed (open?) meshes.

One possible solution would be to brute-force it by checking each coordinate to a certain precision. With space partitioning, I could also determine whether a coordinate is within the solid bounds of any vertices. This provides a ratio that could be used on the rectangular bounding box volume. The main downside to this is speed, like all brute-force techniques. I'm hoping for a particularly more efficient solution.

Answer 1

If any three vertices (defined by vectors $A$, $B$, $C$) form a triangle as a part of the skin of a volume, then the volume of the triangle pyramid is formed by the three vertices, and the origin is

$${\rm volume}(A,B,C) = \frac{ A \cdot ( B \times C) }{6}$$

_{Here $\cdot$ is the vector dot product, and $\times$ the vector cross product. The above is often called the vector triple product.}

To find the total volume, split the object into triangles and add up all the volumes from above. If a face normal is away from the origin it will count as positive volume, and towards the origin, it will be negative. This way due to Green's theorem in total only the volume enclosed by the mesh will be counted.

Please also note the area of the triangle is

$${\rm area}(A,B,C)=\tfrac{1}{2}\|A\times B+B\times C+C\times A\|$$

Here is a simple example with

$$ \begin{aligned}{A} & =\begin{pmatrix}40\\ 0\\ 0 \end{pmatrix} & {B} & =\begin{pmatrix}0\\ 25\\ 0 \end{pmatrix} & {C} & =\begin{pmatrix}0\\ 0\\ 8 \end{pmatrix}\end{aligned}$$

The volume is

$$ V=\tfrac{1}{6}\begin{pmatrix}40\\ 0\\ 0 \end{pmatrix}\cdot\left(\begin{pmatrix}0\\ 25\\ 0 \end{pmatrix}\times\begin{pmatrix}0\\ 0\\ 8 \end{pmatrix}\right)=\tfrac{1}{6}\begin{pmatrix}40\\ 0\\ 0 \end{pmatrix}\cdot\begin{pmatrix}200\\ 0\\ 0 \end{pmatrix}=\frac{8000}{6}=1333.3\overline{3}\;\checkmark$$

Which is the value that matches with CAD.

---

And here is a more complex example

which I matched with CAD also. The STL file contained 5088 vertices and 1696 triangular faces. The CAD model gave me a volume of 46085.4 and a Fortran program I checked with gave me the value 46086.1

Now granted the Fortran program was intended to derive the mass properties from an STL file, such as volume, the center of mass, and mass moment of inertia tensor.

In the code there is clearly visible the line i_V = dot(A, cross(B, C))/6 which computes the volume for each triangle.

For completeness here is the processing code

subroutine calc_properties(mesh,mass,V,cg,Ic)
class(mesh3), intent(in) :: mesh
real(8), intent(in) :: mass
real(8), intent(out) :: V
type(vector3), intent(out) :: cg
type(symmatrix3), intent(out) :: Ic
real(8) :: rho, i_V
type(vector3) :: i_cg, A,B,C
type(symmatrix3) :: i_mmoi, I0
integer :: i
    V = 0
    cg = 0
    I0 = 0
    do i=1, mesh%n_faces
        A = mesh%nodes( mesh%faces(1, i) )
        B = mesh%nodes( mesh%faces(2, i) )
        C = mesh%nodes( mesh%faces(3, i) )
        i_V = dot(A, cross(B, C))/6
        V = V + i_V
        i_cg = (A+B+C)/4
        cg = cg + i_cg*i_V
        i_mmoi = (parallel(A+B) + parallel(B+C) + parallel(C+A))/20
        I0 = I0 + i_mmoi*i_V
    end do
! Find the center of mass
    cg = cg/V
! I0 now contains the MMOI at the origin (divided by density)
! We need to transfer it to the center of mass
    rho = mass/V
    Ic = rho*I0  - mass*parallel(cg)
end subroutine 

function parallel(pos) result(r)
class(vector3), intent(in) :: pos
type(symmatrix3) :: r
    r%A11 = (pos%y**2 + pos%z**2)
    r%A22 = (pos%x**2 + pos%z**2)
    r%A33 = (pos%x**2 + pos%y**2)
    r%A12 = -(pos%x*pos%y)
    r%A23 = -(pos%y*pos%z)
    r%A13 = -(pos%x*pos%z)
end function 

function dot(this,other) result(r)
class(vector3), intent(in) :: this, other
real(8) :: r
    r = this%x*other%x + this%y*other%y + this%z*other%z
end function

function cross(this,other) result(r)
class(vector3), intent(in) :: this, other
type(vector3) :: r
    r = vector3( &
        this%y*other%z - this%z*other%y, &
        this%z*other%x - this%x*other%z, &
        this%x*other%y - this%y*other%x )
end function