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

Question

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

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

我的代码都是点点滴滴，主要是：

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

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

Rhino 7 的参考 NURBS 体积
27k 三角形网格与 Rhino 7
使用 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;
    }
}

Answer 1

我认为你把它弄得有点复杂了。您需要做的就是计算出所有点的平均 x、y、z，这就是您的中心。您可以单独执行每个操作，因此伪代码为：

xSum=0;
ySum=0;
zSum=0;
pointCount=0;
for each triangle in meshArray {
 for each point in triangle {
   xSum += point.x;
   ySum += point.y;
   zSum += point.z;
   pointCount++;
 }
}
centerPointX = xSum / pointCount;
centerPointY = ySum / pointCount;
centerPointZ = zSum / pointCount;

- 编辑 -

为了平滑不平衡的网格，您可以将 meshArray 处理成大小相同的体素（3D 区域立方体），您可以根据所需的精度选择粒度。

因此对于 100x100x100 体素点云平滑：

for (x=0;x<=100;x++){
 for (y=0;y<=100;y++){
  for (z=0;z<=100;z++){
   hasPoint = false;
   for each triangle in meshArray {
    for each point in triangle {
      if point.x,y,z contained with x,y,z of this voxel {hasPoint=true;}
    }
   }
   addPoint at voxel (x,y,z) center to lowDensityPointCloud
  }
 }
}

然后你拿你的lowDensityPointCloud，用前面的方法找到它的中心点。

Answer 2

原来堆栈溢出中的和整个互联网都是错误的:

几何质心并不总是与质心重合，尽管它非常接近。

这仅适用于某些卷，如 n-dimensional simplexes and platonic solids。

我已经通过将我的结果与 Rhino 的结果进行比较，并通过计算平衡时船体的位置与现实情况进行了比较。

想了很多，现在看来还是挺明显的，不过也欢迎大家检查批评我的方法。

方法

我使用 Zhang & Chen's paper 中描述的方法的原理，使用质心。

四面体的质心与其质心重合（即顶点的等重心）。

由四面体的有符号体积加权的质心的质心是网格的质心。

代码

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 3

我将从计算体积开始，一旦建立了基本概念，我就会移动到质心。

这个答案听起来有点像黑魔法。我们可以在这里使用 Divergence_theorem 。该定理将体积积分问题转化为对物体表面的计算。

要理解它，您确实需要一些微积分。如果你知道一点积分，那么我们可以通过评估端点（区间的边界）来计算积分。

散度定理的工作原理大致相同，我们可以通过计算边界上的相关量来计算某个量的积分，但现在边界是表面本身。

这里F是一些向量值，V代表体积，S代表表面，n代表法线，∇代表发散算子。

对于我们的例子，我们希望 LHS 积分中的位仅为 1。因此我们计算 1 在体积上的积分。有几个候选者，最简单的就是 F(x,y,z) = (x,0,0)。 ( ∇ . F = df/dx + df/dy + df/dz = 1 + 0 + 0 = 0).

对于 RHS，我们现在在曲面上积分 F(x,y,z).N。由于我们的表面是三角网格，这意味着我们取总和

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

所以任务归结为找到函数在三角形上的积分。这比看起来要简单得多。对于顶点 A=(x,y,z)，F(A) = (x,0,0) 并且如果向外法向的单位是 (nx, ny, nz) 那么 F(A) 。 N = (x,0,0) 。 (nx,ny,nz) = x * nx，称这个量为 fA。还要计算 fB 和 fC 的数量。

对于三角形内的点，可以将数量计算为坐标值的线性组合。冗长的积分产生了简单的公式

  1/3 area ( fA + fB + fC)

其中面积是三角形的面积，我们就完成了。

因此将以上所有内容放入 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)

告诉你它的黑魔法。

现在对于质心，我们应用相同的技术，但对 F 使用三个不同的函数。质心的 x 坐标由

给出

F 的候选者是 F(x,y,z) = (1/2 x^2, 0, 0) 作为 ∇ 。 F = x.

事情类似，但不是我们的函数是非线性的，所以每个三角形的积分计算会更复杂。

我们可以使用Barycentric coordinate system来简化三角形的积分。

正在为

计算

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

给出三角形的积分为

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

我使用 wolfram alpha 进行了整合。计算伪代码然后是体积和质心

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

我已经将代码合并到我自己的程序中，并检查了具有已知体积的各种对象的体积和质心结果，例如椭圆体，SingSurf volume calculator。

如果我们选择 F(x,y,z) = 1/3 (x,y,z)，我们可以看到体积公式与通过对四面体的有符号体积求和得到的公式相同.

利用散度定理，r = (x,y,z)

对于线性函数 g(r)，三角形上的积分是面积乘以函数在顶点求值的平均值。

这里

现在，非标准化法线可以计算为 n = (B-A) x (C-A)。面积 = 1/2 | ñ |和单位正常 n-hat = n / |n|所以面积 * n-hat = 1/2 n.

根据需要。

Answer 4

问题出在centroid += triangleArea * (a + b + c) / 3

这给出了三角形的质心，但是对于质心和体积计算，您需要底面 (a,b,c) 的 triangular pyramid 的体积质心并倾斜原点。

但要理解，请查看为每个三角形定义的以下体积元素的属性。

上面标记为CG的是金字塔的体积质心，其加权平均值与金字塔的体积一起定义了整个形状的质心。

给定向量 A、B 和 C，您的加权平均值为

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

CG = CG/V

其中 · 是向量点积，× 是向量叉积。

这里关键是你需要

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

因为金字塔的质心距离底部的距离是 1/4。这是相关的维基百科段落

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

Calculate the centroid of a 3D mesh of triangles

c#

math

3d

mesh

方法

代码