计算一个由三角形组成的三维网格的质心

Question

计算一个由三角形组成的三维网格的质心

c#math3dmesh

8

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

编辑：事实证明我并没有这样做，我是在尝试计算重心，这不是一回事

我的代码由各种部分组成，主要包括：

我将我的结果与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;
    }
}

- geriwald

我可能已经找到了那个问题的答案：似乎我使用的方法不起作用。我仍然需要检查，但我认为它只适用于凸固体。我应该在数学论坛上发布这个问题。我还需要进行更多的检查，但似乎通过计算按其代数质量加权的单位四面体重心的平均值可以解决问题。一旦我检查完我的结果，我会详细解释它。 - geriwald

请参考我在Stack Overflow上的这篇文章，了解如何计算STL文件（三角网格）的质量属性。链接：https://stackoverflow.com/a/58286695/380384 - John Alexiou

4个回答

2

问题在于 centroid += triangleArea * (a + b + c) / 3 这给出了一个三角形的重心，但对于质心和体积计算，您需要具有基础为(a，b，c)，顶点为原点的三角形金字塔的体积重心。

我有更详细的答案

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

trig

上面标记为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处。以下是相关的维基百科段落。

- John Alexiou

是的，这就是我在被接受的答案中所做的。 - geriwald

1

虽然你的答案是正确的，但你的演示不正确。你应该看一下金字塔的定义。你画的是四面体。这是另一个维基百科页面。 - geriwald

@geriwald 这是一个三角锥。 - John Alexiou

哦，我读你的答案太快了。我被你的第一行搞糊涂了：在我的第一个版本中，我并不是试图计算四面体的重心，而是三角形的重心，正如我在问题开头引用的已接受答案所建议的那样。我没有解释清楚，这是真的。 - geriwald

1

一旦建立了基本思想，我将开始计算体积并移至质心。

这个答案听起来有点像黑魔法。我们可以在这里使用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)

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

所以，为了将上述内容都转化为伪代码：

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。

类似地进行计算，但我们的函数是非线性的，因此每个三角形的积分计算会更加复杂。

我们可以使用重心坐标系来简化三角形上的积分。

对于以下函数的计算：

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体积计算器）检查了体积和重心结果。

如果我们选择F(x,y,z) = 1/3 (x,y,z)，则我们可以看到体积的公式与求四面体有符号体积的总和所得的公式完全相同。

使用散度定理，其中 r = (x,y,z)。

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

这里

现在可以计算非规范化的法向量，公式为 n = (B-A) x (C-A)。面积 = 1/2 | n |，单位法向量 n-hat = n / |n|，因此面积 * n-hat = 1/2 n。

按要求。

- Salix alba

这里的1/3因子是不正确的。金字塔体积元的正确因子是1/4。 - John Alexiou

我不明白为什么你的结果偏向于x坐标。这看起来很可疑。 - geriwald

@JohnAlexiou 我不是在计算金字塔元素的体积。我正在三角形上积分一个量。这个算法的好处是我们不需要将体积分成3个简单形式，只需将表面三角化即可。 - Salix alba

请注意，ABC 三角形的（非单位）法向量为 n = (a × b + b × c + c × d)，三角形的面积为 A = 1/2 |(a × b + b × c + c × d)|。 - John Alexiou

1

我已经有几天没上 Stack Overflow 了，因为我在测试我的软件（顺便说一下，我的解决方案非常好，与 Rhino 相比给出的相对差异小于 1e-9）。 - geriwald

显示剩余4条评论

0

我认为你把它想得比必要的复杂了。你只需要计算所有点的平均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处理成大小相等的体素（三维区域立方体），其粒度取决于所需的精度。

因此，对于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
  }
 }
}

然后您拿低密度点云，并使用先前的方法找到其中心点。

- James

1

嗯，我认为这不起作用。想象一下一个立方体，其中一个面有4个顶点，而另一面则是非常密集的小三角形网格。甚至比第一个面多一个顶点。 - geriwald

你说得对。那我建议将体积处理成体素，以获取每个体素的平滑中心值，然后您可以处理体素结果并取其平均值。有意义吗？ - James

我找到了解决方案，明天我会发布，这有点复杂，无法通过评论来解释。 - geriwald

我需要非常高的准确性。这可能不会在计算上非常高效。 - geriwald

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- geriwald · Accepted Answer

事实证明，Stack Overflow和整个互联网上通常被接受的答案是错误的：

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

这仅适用于某些体积，如n维单形体和柏拉图立体。

我通过将我的结果与Rhino和计算船体在平衡时的位置的真实结果进行比较来验证了这一点。

经过深思熟虑，现在对我来说似乎很显然，但您可以随时检查和批评我的方法。

方法

我使用Zhang＆Chen的论文中描述的方法原则，并使用重心。

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

以四面体的有向体积加权的质心为重心，可以得到网格的质心。

代码

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