[Laplacian Mesh Smoothing by Throwing Vertices]

日期：2025年3月9日 标签：[Programming]

在这篇博文中，我将讨论使用 Laplacian 网格平滑来平滑和模糊 3D 网格。一个很好的例子是 Adobe Substance 3D Modeler 的 [smooth tool]，我就是用 Laplacian 网格平滑来实现的。Laplacian 网格平滑的工作原理是迭代地将每个顶点移动到该顶点邻居的平均位置。其数学公式为：

vi=1NN∑j=1vj

我们可以在代码中这样实现：

for (const Vertex& vertex : mesh.vertices)
{
  vec3 avrgPosition = vec3(0);
  for (const Vertex& vertexNeighbor : neighbors(vertex))
  {
    avrgPosition += vertexNeighbor;
  }
  avrgPosition /= numNeighbors(vertex);
}

运行几次，网格就会像变魔术一样平滑。重要的是要注意，这个过程只改变网格的顶点数据，而连通性（三角形索引）保持不变。

上面的内容很简单，只是我们需要 neighbors() 和 numNeighbors() 函数。虽然一种方法是构建一个半边网格数据结构，它可以提供关于网格连通性的信息，但这篇文章将展示两种替代方法，它们可以在网格是二流形（即每条边恰好连接两个三角形）时使用。

一次性技巧：抛掷顶点

忘记邻居查找吧。相反，我们将直接迭代三角形面，利用二流形网格中每条边恰好连接两个三角形这一事实。我们不会预先计算邻居，而是在一次性遍历面的过程中动态地累积邻居位置。

这个过程通常被称为 抛掷顶点，这是像 libigl 这样的库中使用的术语。它指的是在与其相连的邻居处累积顶点位置，然后对它们进行平均。对于由顶点 A、B 和 C 组成的三角形，我们定义三个有向半边：A->B、B->C 和 C->A。当我们遍历这些半边时，我们在目标顶点位置累积源顶点位置，并计算每个顶点被更新的次数。这有效地计算了每个顶点的所有相邻顶点位置的总和。一旦我们有了这个总和，我们就除以计数以获得新的、平滑后的位置。

以下是实现代码：

std::vector<vec3> accumulated(mesh.vertices.size(), 0);
std::vector<int> numNeighbors(mesh.vertices.size(), 0);
for (const auto& triangle : mesh.triangles)
{
  int vertexIndex0 = triangle.index0;
  int vertexIndex1 = triangle.index1;
  int vertexIndex2 = triangle.index2;
  // vertexIndex0 -> vertexIndex1
  accumulated[vertexIndex1] += mesh.vertices[vertexIndex0];
  numNeighbors[vertexIndex1]++;
  // vertexIndex1 -> vertexIndex2
  accumulated[vertexIndex2] += mesh.vertices[vertexIndex1];
  numNeighbors[vertexIndex2]++;
  // vertexIndex2 -> vertexIndex0
  accumulated[vertexIndex0] += mesh.vertices[vertexIndex2];
  numNeighbors[vertexIndex0]++;
}
for (int i=0; i<mesh.vertices.size(); i++)
{
  mesh.vertices[i] = accumulated[i] / numNeighbors[i];
}

就这样：一次 Laplacian 平滑处理完成，而无需预先计算任何数据或创建半边数据结构。

如果你想让网格更平滑，只需应用更多次上述过程。如果你想减慢平滑过程，只需使用小数值：

mesh.vertices[i] = lerp(mesh.vertices[i], accumulated[i]/numNeighbors[i], fractionalAmount);

这种方法简单有效。

通过并行性榨取性能

上述方法在单线程场景中效果很好，但在多线程代码或 GPU 上却不起作用。问题在于不同线程同时写入 accumulated 和 numNeighbors 数组。这属于可变、共享数据的问题类别。例如，如果不同的线程处理单独的三角形，它们可能会同时修改这些三角形共享的同一个顶点。

原子混乱：事情变得 MESI

std::vector<std::atomic<vec3>> accumulated(mesh.vertices.size());
std::vector<std::atomic<int>> numNeighbors(mesh.vertices.size());

这可以允许多个线程并行运行平滑代码。但是，原子操作并非没有代价，缓存争用会降低 CPU 和 GPU 的速度。在我自己的 CPU 测试中，我发现使用多线程的这种方法比单线程非原子版本慢。这可能不适用于你的情况，所以在尝试之前/之后请对你的情况进行分析。

一种更简洁的方法：预先计算邻居

一个更好的解决方案是预先计算邻居列表一次，并在多次平滑处理中重用它。 如下所示：

// Find the number of neighbors of each vertex
std::vector<int> numNeighbors(mesh.vertices.size(), 0);
int totalNeighbors = 0;
for (const auto& triangle : mesh.triangles)
{
  numNeighbors[triangle.index0]++;
  numNeighbors[triangle.index1]++;
  numNeighbors[triangle.index2]++;
  totalNeighbors += 3;
}
// Prefix sum
std::vector<int> neighborsStartIndex(mesh.vertices.size(), 0);
int currentIndex = 0;
for (int i=0; i<numNeighbors.size(); i++)
{
  neighborsStartIndex[i] = currentIndex;
  currentIndex += numNeighbors[i];
}
// Populate the neighbors array
std::vector<int> neighbors(totalNeighbors, 0);
std::fill(numNeighbors.begin(), numNeighbors.end(), 0);
for (const auto& triangle : mesh.triangles)
{
  // vertexIndex0 -> vertexIndex1
  neighbors[neighborsStartIndex[triangle.index1] + numNeighbors[triangle.index1]] = triangle.index0;
  numNeighbors[triangle.index1]++;
  // vertexIndex1 -> vertexIndex2
  neighbors[neighborsStartIndex[triangle.index2] + numNeighbors[triangle.index2]] = triangle.index1;
  numNeighbors[triangle.index2]++;
  // vertexIndex2 -> vertexIndex0
  neighbors[neighborsStartIndex[triangle.index0] + numNeighbors[triangle.index0]] = triangle.index2;
  numNeighbors[triangle.index0]++;
}

有了这个，平滑处理可以并行处理每个顶点。

std::vector<vec3> result(mesh.vertices.size(), 0);
// this for() loop can now operate on each vertex in parallel;
// the result[] array holds the new mesh vertices
for (int i=0; i<mesh.vertices.size(); i++)
{
  for (int j=0; j<numNeighbors[i]; j++)
  {
    int neighborIndex = neighbors[neighborStartIndex[i] + j];
    result[i] += mesh.vertices[neighborIndex] / numNeighbors[i];
  }
}

一旦构建了 neighbors 和 numNeighbors 列表，我们就再也不需要碰它了。

顶点法线

顶点平滑后，需要重新计算顶点法线。可以使用相同的技巧来完成：迭代面，计算面法线，并将其抛掷到该面的顶点。

std::vector<vec3> accumulated(mesh.vertices.size(), 0);
for (const auto& triangle : mesh.triangles)
{
  vec3 v0 = mesh.vertices[triangles.index0];
  vec3 v1 = mesh.vertices[triangles.index1];
  vec3 v2 = mesh.vertices[triangles.index2];
  vec3 faceNormal = normalize(cross(v1-v0, v2-v0));
  accumulated[triangles.index0] += faceNormal;
  numNeighbors[triangles.index0]++;
  accumulated[triangles.index1] += faceNormal;
  numNeighbors[triangles.index1]++;
  accumulated[triangles.index2] += faceNormal;
  numNeighbors[triangles.index2]++;
}
for (int i=0; i<mesh.vertices.size(); i++)
{
  mesh.vertices[i] = normalize(accumulated[i]);
}

如果想要更准确的结果，可以使用角度或面积加权的法线。那是另一个深坑，但一个好的参考是 [Weighted Vertex Normals]。