Calculate the centroid of a 3D mesh of triangles

Question

I'm trying to calculate the centroid of a 3D mesh of triangles.

EDIT : it turns out I was not, I was trying to calculate the center of gravity, which is not the same

My code is made of bits and pieces, mainly :

• This nice paper by Cha Zhang and Tsuhan Chen
• SO : How to calculate the volume of a 3D mesh object the surface of which is made up triangles
• SO : How to compute the centroid of a mesh with triangular faces?

I compared my results to those provided by Rhino. I calculate the centroid and volume :

• of the reference NURBS volume with Rhino 7
• of a 27k triangle mesh with Rhino 7
• of a simplified 1k triangle mesh with Rhino 7
• of the same 1k triangle mesh with my code.

As you can see, it works great to calculate the volume, but not for the centroid, and i can't seem to know why. I need the error to be less than 0.01. I checked everything several times, but there must be something obvious.

I'm not great with numerical instability :

• should I work in milimeters instead of meters ?
• should I calculate the tetrahedrons signed volume with another point than the origin, as suggested by galinette in the second reference ? I tried and it didn't improve much.

MY CODE

Before calculationg anything, I check that my mesh is correct (code not provided) :

• closed mesh, no naked edges or any holes ;
• the vertices of all triangles are ordered consistently, i.e. triangles are correctly oriented towards the outside of the mesh.

using HelixToolkit.Wpf;
using System.Collections.Generic;
using System.Windows.Media.Media3D;

internal static class CentroidHelper
{
    public static Point3D Centroid(this List<MeshGeometry3D> meshes, out double volume)
    {
        Vector3D centroid = new Vector3D();
        volume = 0;

        foreach (var mesh in meshes)
        {
            var c = mesh.Centroid(out double v);
            volume += v;
            centroid += v *c ;
        }

        return (centroid / volume).ToPoint3D();
    }

    public static Vector3D Centroid(this MeshGeometry3D mesh, out double volume)
    {
        Vector3D centroid = new Vector3D();
        double totalArea = 0;
        volume = 0;

        for (int i = 0; i < mesh.TriangleIndices.Count; i += 3)
        {
            var a = mesh.Positions[mesh.TriangleIndices[i + 0]].ToVector3D();
            var b = mesh.Positions[mesh.TriangleIndices[i + 1]].ToVector3D();
            var c = mesh.Positions[mesh.TriangleIndices[i + 2]].ToVector3D();
            var triangleArea = AreaOfTriangle(a, b, c);
            totalArea += triangleArea;
            centroid += triangleArea * (a + b + c) / 3;
                
            volume += SignedVolumeOfTetrahedron(a, b, c);
        }
        return centroid / totalArea;
    }

    private static double SignedVolumeOfTetrahedron(Vector3D a, Vector3D b, Vector3D c)
    {
        return Vector3D.DotProduct(a, Vector3D.CrossProduct(b, c)) / 6.0d;
    }

    private static double AreaOfTriangle(Vector3D a, Vector3D b, Vector3D c)
    {
        return 0.5d * Vector3D.CrossProduct(b - a, c - a).Length;
    }
}    

Answer 1

It turns out the commonly accepted answer in stack overflow and all over the internet is wrong :

The geometric centroid DOES NOT always coincide with the center of mass, although it is very close.

It is only true for some volumes, as n-dimensional simplexes and platonic solids.

I have checked this by comparing my results to Rhino's and to reality by calculating a ship's hull position at equilibrium.

After thinking about it a lot, it seems quite obvious to me now, but you are very welcome to check and criticise my method.

METHOD

I use the principles of the method described in Zhang & Chen's paper, with the center of mass.

The center of mass of a tetrahedron does coincide with its centroid (i.e. the isobarycenter of its vertices).

The barycenter of the centroids weighted by the signed volume of the tetrahedrons is the centroid of the mesh.

CODE

internal static class CentroidHelper
{
    public static Point3D Centroid(this List<MeshGeometry3D> meshes, out double volume)
    {
        Vector3D centroid = new Vector3D();
        volume = 0;

        foreach (var mesh in meshes)
        {
            var c = mesh.Centroid(out double v);
            volume += v;
            centroid += v * c;
        }
        return (centroid / volume).ToPoint3D();
    }

    public static Vector3D Centroid(this MeshGeometry3D mesh, out double volume)
    {
        Vector3D centroid = new Vector3D();
        volume = 0;

        for (int i = 0; i < mesh.TriangleIndices.Count; i += 3)
        {
            var a = mesh.Positions[mesh.TriangleIndices[i + 0]].ToVector3D();
            var b = mesh.Positions[mesh.TriangleIndices[i + 1]].ToVector3D();
            var c = mesh.Positions[mesh.TriangleIndices[i + 2]].ToVector3D();
            var tetrahedronVolume = SignedVolumeOfTetrahedron(a, b, c);
            centroid += tetrahedronVolume * (a + b + c) / 4.0d;
            volume += tetrahedronVolume;
        }
        return centroid / volume;
    }

    private static double SignedVolumeOfTetrahedron(Vector3D a, Vector3D b, Vector3D c)
    {
        return Vector3D.DotProduct(a, Vector3D.CrossProduct(b, c)) / 6.0d;
    }
}

Answer 2

I'll start by calculating the volume and move onto the centroid once the basic idea is established.

This answer is going to sound a little like black magic. We can use the Divergence_theorem here. This theorem translates problems of integration over a volume with calculation over the surface of the object.

To understand it you really need some calculus. If you know a little integration then we can calculate the integral just by evaluating the end points (the boundary of the interval).

The divergence theorem works much the same, and we can calculate the integral of some quantity by evaluating a related quantity over the boundary, but now the boundary is the surface itself.

Here F is some vector valued quantity, V stands for the volume, S the surface, n the normal and ∇ the divergence operator.

For our case we would want the bit inside the LHS integral to just be 1. So we calculate the integral of 1 over the volume. There are a few candidates for this and the simplest is just F(x,y,z) = (x,0,0). ( ∇ . F = df/dx + df/dy + df/dz = 1 + 0 + 0 = 0).

For the RHS we now integrate F(x,y,z).N over the surface. As our surface is a triangulated mesh this means that we take the sum

sum_over_all_triangles (integral of F.N for that triangle).

So the task boils down to finding the integral of the function over the triangle. This is a lot simpler than it seems. For a vertex A=(x,y,z), F(A) = (x,0,0) and if the unit outwards normal is (nx, ny, nz) then F(A) . N = (x,0,0) . (nx,ny,nz) = x * nx, call this quantity fA. Calculate the quantities fB and fC as well.

For points inside the triangle the quantity can be calculated as a linear conbination of the values at the coordinates. A long-winded integrations yields the simple formula

  1/3 area ( fA + fB + fC)

where area is the area of the triangle, and we are done.

So to put all the above in psudocode

total = 0
for each triangle:
    let A = (ax,ay,az), B, C be three verticies
    let N = (nx,ny,nz) be the unit outward normal
    let area be the triangle area
    fA = ax * nx
    fB = bx * nx
    fC = cx * nx
    total += 1/3 * area * (fA + fB + fC)

Told you its black magic.

---

Now for the centroid we apply the same technique but use three different functions for F. The x coordinate of the centroid is given by

A candidate for F is then F(x,y,z) = (1/2 x^2, 0, 0) as ∇ . F = x.

Things work similarly but not our function is non-linear, so the calculation for the integral of each triangle will be more complex.

We can use Barycentric coordinate system to make the integral over the triangle easier.

Calculating this for

f(x,y,z) = F(x,y,z) . (nx,ny,nz) = 1/2 x^2 * nx

gives the integration over a triangle as

1/12 area * nx (ax^2 + bx^2 + cx^2 + ax bx + ax cx + bx cx)

I used wolfram alpha to do the integration. The pseudocode for calculation volume and centroid is then

volume = 0
Cx = 0
Cy = 0
Cz = 0
for each triangle:
    let A = (ax,ay,az), B, C be three verticies
    let N = (nx,ny,nz) be the unit outward normal
    let area be the triangle area

    volume += 1/3 * area * (ax + bx + cx) * nx
    Cx += 1/12 area * nx (ax^2 + bx^2 + cx^2 + ax bx + ax cx + bx cx)
    Cy += 1/12 area * ny (ay^2 + by^2 + cy^2 + ay by + ay cy + by cy)
    Cz += 1/12 area * nz (az^2 + bz^2 + cz^2 + az bz + az cz + bz cz)

Cx /= volume
Cy /= volume
Cz /= volume

I've incorporated to code in my own program and checked volumes and centroids results with various objects with known volumes, such as elipsoids, SingSurf volume calculator.

---

If we choose F(x,y,z) = 1/3 (x,y,z), we can see that the formula for volumes works out to be identical to that found by summing signed volumes of tetrahedrons.

Using the divergence theorem, with r = (x,y,z)

For a linear function, g(r), the integral over the triangle is the area times the mean of the function evaluated at the vertices.

Here

Now a non-normalised normal can be calculated as n = (B-A) x (C-A). The area = 1/2 | n | and the unit normal n-hat = n / |n| so area * n-hat = 1/2 n.

As required.

Answer 3

The problem is in centroid += triangleArea * (a + b + c) / 3

This gives the centroid of a triangle, but for center of mass and volume calculations, you need the volume centroid of the triangular pyramid with base (a,b,c) and tip the origin.

• I have this answer with more details

but to understand look at properties of the following volume element defined for each triangle.

What is marked as CG above is the volume centroid of the pyramid and together with the volume of the pyramid its weighted average defines the center of mass of the entire shape.

Given vectors A, B, and C you have the weighted average as

for all triangles
  ΔV = A · (B × C)/6
  V += ΔV
  CG += ΔV * (A + B + C)/4
next

CG = CG/V

where · is the vector dot product, and × the vector cross product.

The key here is that you need

centroid += triangleVolume * (a + b + c) / 4

since the centroid of a pyramid is 1/4 the distance from the base. Here is the relevant Wikipedia paragraph