How can I calculate the volume of a mesh?

General Usage
,
1

How can I calculate the volume of a mesh?

The mesh is a closed triangles mesh, with 3-D coordinates; it is a surface mesh.

I have tried the following, but I am getting an error:

using GLMakie

using FileIO

using MeshIO

file_path = “my_file.stl”

my_mesh = load(assetpath(file_path))

volume(my_mesh)

ERROR: Makie.convert_arguments for the plot type Volume{Tuple{GeometryBasics.Mesh{3, Float32, GeometryBasics.TriangleP{3, Float32, GeometryBasics.PointMeta{3, Float32, Point{3, Float32}, (:normals,), Tuple{Vec{3, Float32}}}}, GeometryBasics.FaceView{GeometryBasics.TriangleP{3, Float32, GeometryBasics.PointMeta{3, Float32, Point{3, Float32}, (:normals,), Tuple{Vec{3, Float32}}}}, GeometryBasics.PointMeta{3, Float32, Point{3, Float32}, (:normals,), Tuple{Vec{3, Float32}}}, GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}, StructArrays.StructVector{GeometryBasics.PointMeta{3, Float32, Point{3, Float32}, (:normals,), Tuple{Vec{3, Float32}}}, NamedTuple{(:position, :normals), Tuple{Vector{Point{3, Float32}}, Vector{Vec{3, Float32}}}}, Int64}, Vector{GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}}}}}} and its conversion trait VolumeLike() was unsuccessful.

2

I don’t know whether there is a Julia package for it; others might.

Anyhow, is the mesh guaranteed to be convex (no “indents”)?

If yes, a simple way would be to calculate the sum of the volumes of all triangular pyramids with their apex being a (the same) point inside the mesh. If the mesh is not convex, things are more complex.

1 Like
3

The mesh is not necessarily convex.
I thought about creating a 3D-mesh (a mesh of tetrahedron) from the 2-D mesh (mesh of triangles), and then calculating the sum of the tetrahedron volumes. But I also do not know how to get the 3D-mesh from the triangle-mesh.

4

Isn’t there some theorem that proves the volume can be calculated by taking an arbitrary center, computing the signed volumes of tetrahedra constructed by connecting the center with each triangular facet, and adding them together? The facets must represent a consistent cover of the boundary.

5

Do you need convexity for this? Suppose the region is two spheres intersecting only at one boundayr point, does “singed volume” fix this up?

6

In this case the boundary is empty, isn’t it?

7

I’m thinking of two separate spheres, so the boundary is the union of the two boundaries.

8

In that case, wouldn’t that be a union of two convex shapes sharing a single commong point. Two convex shapes period, really?

9

Yes. Does the volume computation work for such a geometry?

10

Or let them overlap a bit and call it a peanut.

11

Solution described here seems simple to implement

2 Likes
12

I still think it would work. I believe the principle is the same as that of the polar planimeter.

13

It appears to be a variant of the Green’s theorem.

1 Like
14

I found a Matlab implementation:

It seems to be the same idea as the link you mention.
I have ported the code to Julia, and it gives me the same volume as Meshlab for the two meshes I tested on.

using LinearAlgebra

# ------------------------------------------------
# original code in Matlab:
# https://ch.mathworks.com/matlabcentral/fileexchange/93875-volume-surface-area-and-centroid-of-polyhedron
# ------------------------------------------------
# This function find the volume, surface area, and centroid of a polyhedron. 
# The polyhedron faces should be divided into triangles. The numbering of all triangles faces should be consistent. 
# Either all of them in clock wise or counter clock wise direction as seen from the outside.
# 
# Input:
# V: (Nv,3) vertices matrix
# F: (Nf,3) faces matrix
# Output:
# Vol: volume
# A: surface area
# Cx, Cy, and Cz: x,y, and z coordinates of centroid, respectively.
# 
# ------------------------------------------------
function Polyhedron_VAC(V,F)
  Vol=0;
  A=0;
  Mx=0;
  My=0;
  Mz=0;
  
  for i=1:1:size(F,1)
      
    x1=V[F[i,1],1];
    y1=V[F[i,1],2];
    z1=V[F[i,1],3];

    x2=V[F[i,2],1];
    y2=V[F[i,2],2];
    z2=V[F[i,2],3];
      
    x3=V[F[i,3],1];
    y3=V[F[i,3],2];
    z3=V[F[i,3],3];
     
    vi=1/6*det([x1 x2 x3 0; y1 y2 y3 0; z1 z2 z3 0; 1 1 1 1])
    
    # ai=1/2*norm( cross(V(F(i,2),:)-V(F(i,1),:),V(F(i,3),:)-V(F(i,1),:)));
    ai=1/2*norm( cross(V[F[i,2],:] - V[F[i,1],:], V[F[i,3],:]-V[F[i,1],:]) )
    
    xi=1/4*(x1+x2+x3);
    yi=1/4*(y1+y2+y3);
    zi=1/4*(z1+z2+z3);
    Vol=Vol+vi;
    A=A+ai;
    Mx=Mx+xi*vi;
    My=My+yi*vi;
    Mz=Mz+zi*vi;
  end
  
  Cx=Mx/Vol
  Cy=My/Vol
  Cz=Mz/Vol
  S=sign(Vol)
  Vol=Vol*S
end
```
1 Like