计算三角形3D网格的质心

c#

我正在尝试计算 3D 三角形网格的质心。

编辑：事实证明我不是，我试图计算重心，这是不一样的

我的代码由点点滴滴组成，主要是：

• Cha Zhang和Tsuhan Chen的这篇好论文
• SO：如何计算表面由三角形组成的 3D 网格对象的体积
• SO：如何计算具有三角形面的网格的质心？

我将我的结果与 Rhino 提供的结果进行了比较。我计算质心和体积：

• 使用 Rhino 7 的参考 NURBS 体积
• 使用 Rhino 7 的 27k 三角形网格
• 使用 Rhino 7 的简化 1k 三角形网格
• 与我的代码相同的 1k 三角形网格。

如您所见，计算体积效果很好，但不适用于质心，我似乎不知道为什么。我需要误差小于 0.01。我检查了好几次，但肯定有一些明显的东西。

我不太擅长数值不稳定：

• 我应该以毫米而不是米为单位工作吗？
• 我是否应该按照第二个参考文献中 galinette 的建议，用另一个点而不是原点来计算四面体有符号体积？我试过了，并没有太大改善。

我的代码

在计算任何东西之前，我检查我的网格是否正确（未提供代码）：

• 封闭的网格，没有裸露的边缘或任何孔洞；
• 所有三角形的顶点都是一致的，即三角形正确地朝向网格的外部。

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;
    }
}