Mesh is Art (4): Subdivision

MeshIsArt

2020-07-19

前言

上文我们介绍了半边网格的底层架构，介绍了点、半边和面所携带的信息。本文我们来探讨一个半边网格的应用——网格细分。网格细分是一种高效地表达几何的方法。

网格细分简述

在 Mesh is Art（1) 中，我们讲解了网格数据大多数由三角形或四边形组成。网格细分技术为分割曲面提供了解决方案。这种技术的核心是关注如何通过细分算法（Subdivision）计算来用大量的较小的网格面来替代原来的曲面，从而细分并优化输入的基础网格面。并且由此产生的精密网格还可以使用相同的算法来产生更精密的网格，如此可以反复迭代。随着迭代步的增加，网格的数量也不断地翻倍，从而更加逼近于精确、光滑的初始曲面。

网格细分技术已经被广泛地应用于游戏和动画产业当中，在项目中常采用非常粗糙的基础网格来表达一个物体（大概只需要几百个多边形），因此存储时所占用的计算机资源非常少。我们知道动画和游戏中常常需要渲染模型，对于模型的渲染时，并不是所有模型都要达到几十万网格面级别再去渲染，这样不仅增加了模型场景中的负担，还浪费了许多时间。在一个场景中，通常只有距离摄像机近、需要突出表达细节的模型才会被渲染地非常清晰，其他背景、场景往往不需要那么高精度的表达。因此，网格细分技术变得十分重要，设计师只需要建立一些粗糙网格模型，当该模型需要被表达时，将其细分，不需要表达细节时，就保持不变。

Figure 1：随着网格细分迭代次数增加，网格模型逐渐变得精细（从右至左）图片来源：Pixie3D | Animation Software
Figure 2：随着网格细分迭代次数增加，网格模型逐渐变得精细（从左至右）图片来源：《Subdivision for Modeling and Animation》

在建筑设计中，许多建筑在设计阶段是由NURBS曲面这样的光滑连续曲面建模技术来设计的。然而，这些模型在几何学方面很容易操作（很容易在计算机中建模），但它们很难被精确地制造出来，由于加工尺度问题，在实际项目中不可避免地要将它们分割成小构件，从而制作出能安装在结构上的表面嵌板（如图3）。在建造中，即使最后完成时是一个单体结构，通常也需要采用离散单元来拼合的传统方法。因此，这些项目在建模时经常使用离散网格来表示曲面（如图4、5），来代替一张参数化连续的NURBS曲面（至少在工程和制造阶段上是这样）。在 Mesh is Art（2）中，也简述了一个由细分算法主导的建筑设计实例。
Figure 3：细分网格在建筑中的应用图片来源：Design and Panelization of Architectural Freeform-Surfaces by Planar Quadrilateral Meshes
Figure 4：用细分网格表达光滑曲面（一）：从左至右依次为：初始光滑曲面、手动创建的初始网格、通过投影算法匹配到初始光滑曲面后的网格、细分之后的网格图片来源：Design and Panelization of Architectural Freeform-Surfaces by Planar Quadrilateral Meshes
Figure 5：用细分网格表达光滑曲面（二）：初始网格（左）和细分后的网格（右）图片来源：Design and Panelization of Architectural Freeform-Surfaces by Planar Quadrilateral Meshes
网格细分除了在建造上表达曲面的优势外，还常被用于分析和优化。

在幕墙优化中，细分曲面常用来控制嵌板的尺寸以及辅助其他优化算法。比如当设计师设计了一张曲面，需要将该曲面铺上嵌板，设计师可以使用细分算法对该曲面不断地迭代细分，直至嵌板的尺寸符合可加工的尺寸，接着可进行更进一步地平板优化和规格优化，当然这两个问题先放到后面了。

在对于需要执行多目标优化的曲面来说，网格细分同样非常有用。在设计和施工的过程中，通常要对建筑曲面进行优化来提高其性能。结构优化通常要求在最大限度节省材料的前提下尽可能小地改变形状，优化的实现可以通过将结构构件仅放置在所需要的地方来尽可能减少结构构件的数量，或者通过调整结构构件以便有效地将载荷传递到支撑来减少结构构件的大小。改变建筑物的形状可以减小风和雪带来的荷载，当然也可以通过考虑热环境、通风和采光来优化建筑的场地、形状和遮阳，从而提高建筑的环境性能。

这些优化目标都需要对建筑模型采用一种特殊表式方法来分析其性能。这些分析需要大量离散的类似网格的表示方法，并且所需要的网格的密度会取决于正在执行的分析类型。例如有限元结构分析可能会要求网格的离散方式能够反映给出的支撑结构，计算风流分析可能需要更精密的分析网格，而太阳能生成计算会相对精确，可能只需要比较粗糙的网格来表示即可。连续壳体的有限元结构分析需要形状函数来在两节点间插入几何结构，而细分曲面也可用于此。

网格细分可以看作是一种圆滑的过程，细分得到的网格子顶点坐标是由周围相邻父顶点的平滑平均值产生。用这种方式，网格的坐标被圆滑，形成更接近初始的近似曲面。这一点已经被证明（Zorin et al.，2000），这种生成的近似曲面在非边界和非奇异顶点处都是C2连续的，这意味着形态表面（没有缝隙）、表面切线（没有折痕）以及切线的变化率（在视觉反射中没有变形）都没有突变（这三点分别是C0，C1和C2连续性的判断标准）。这也是细分网格看上去特别美观的原因。

圆滑过程的另一个优点是，细分不仅能可以圆滑顶点坐标，任何一组与顶点关联的数据都会被圆滑，如果一个初始的粗糙网格中每一个顶点附带着一个颜色的信息，所有连续细分的网格都能够自动匹配圆滑曲面上的颜色。可以类似地应用于更多实际的参数，例如顶点携带的贴图信息、应力值、立面渗透率或透明度值，百叶窗角度或者覆层的偏移距离。在基础网格上，这些数据可以被定义在所需要的顶点上，细分网格后，这些信息将完全以与顶点位置相同的C2连续的方式均匀地分布在整个曲面上。

细分算法原理

实现网格细分有多种方案，有基于三角形网格的，基于四边形网格，也有基于高阶多边形及混合网格的。细分算法的基本过程首先是从一张网格面开始，用某种方式在网格上添加点，并对拓扑结构进行优化。本文中着重讨论两个细分算法的原理及C#实现。

Loop细分算法

原理

Loop细分算法（Loop，1992）是一种应用在三角形网格的算法，如图6所示，这是从一个球形网格中提取的一张网格面，蓝色顶点为父顶点，即网格面原有的顶点，红色顶点为新创建的子顶点，每个三角形上的每条边都被新引入的子顶点分割成了两段，然后每一个三角形都被四个新的小三角形代替。

Figure 6：Loop细分算法的实现步骤：（a）基础网格；（b）拓扑划分；（c）圆滑子顶点；（d）圆滑父顶点图片来源：《Shell Structures for Architecture: Form Finding and Optimization》
内部子顶点：设内部边的两端点 $v_0$ ， $v_1$ ，其相对的两个顶点为 $v_2$ ， $v_3$ ，要创建的子顶点为 $v_{new}$ ，则

v_{new} = \frac{3}{8}(v_0+v_1)+\frac{1}{8}(v_2+v_3)

内部父顶点：设内部父顶点 $v_0$ 的相邻点为 $v_1$ ， $v_2$ ，...， $v_n$ ，更新后的父顶点为 $v_{new}$ ，w为权重值， $n$ 为顶点的价（该顶点的相邻点的数量），则

w=\frac{1}{n}(\frac{5}{8}-(\frac{3}{8}+\frac{1}{4}cos\frac{2\pi}{n})^2)

v_{new} = (1-nw)v_0+w\sum_{i=1}^{n} {v_i}

在细分算法里曲面的边缘通常需要被特殊处理，比如一条只通过一个面的边（面的边缘）。如图8-a所示，沿着边界创建的子顶点会被放置在边界上的边的中点。如图8-b所示，边界上的边的父顶点会被放置在以原始位置为 $\frac{3}{4}$ 权重值和相邻两个顶点位置各占 $\frac{1}{8}$ 权重值的位置，而不会考虑内部邻近顶点的位置。

Figure 8：边界顶点的权重值图片来源：《Shell Structures for Architecture: Form Finding and Optimization》
边界子顶点：设边界边的两个端点为 $v_0$ 和 $v_1$ ，要创建的子顶点为 $v_{new}$ ，则

v_{new} = \frac{v_0+v_1}{2}

边界父顶点：设边界父顶点 $v_0$ 的相邻点为 $v_1$ 和 $v_2$ ，更新后的父顶点为 $v_{new}$ ，w为权重值， $n$ 为顶点的价（该顶点的相邻点的数量），则

v_{new} = \frac{3}{4}v_{0}+\frac{1}{8}(v_1+v_2)

代码实现

在Mesh is Art（3）中，我们讲解了半边网格的特性，即半边网格在处理邻域查找问题上尤为高效。因此在实现网格细分算法中，半边网格是最适合的网格数据结构，这也是为什么我一直在强调使用半边数据结构，类似细分一样大量使用网格邻域查找的算法还有很多。简略版的Loop细分的C#源码如下：

 public Mesh Loop(Mesh m)
  {
    PlanktonMesh P = m.ToPlanktonMesh();
    PlanktonMesh NP = new PlanktonMesh();
    PlanktonMesh TP = P;

    for (int v = 0; v < P.Vertices.Count; v++)
    {
      // Loop
      int[] Neighbours = P.Vertices.GetVertexNeighbours(v);

      int Valence = Neighbours.Length;
      double Beta = 0.625 - (Math.Pow((3 + 2 * Math.Cos((2 * Math.PI) / Valence)), 2) * 0.015625);
      double Mult = Beta / Valence;

      Point3d NewPosition = P.Vertices[v].ToPoint3d() * (1 - Beta);

      foreach (int Neighbour in Neighbours)
      {
        NewPosition += P.Vertices[Neighbour].ToPoint3d() * Mult;
      }
      NP.Vertices.Add(NewPosition);
    }
    for (int HE = 0; HE < P.Halfedges.Count; HE += 2)
    {
      int Pair = P.Halfedges.GetPairHalfedge(HE);
      int ThisOpp = P.Halfedges[HE].PrevHalfedge;
      int PairOpp = P.Halfedges[Pair].PrevHalfedge;

      int HeSV = P.Halfedges[HE].StartVertex;
      int PrSV = P.Halfedges[Pair].StartVertex;

      // split the halfedges in the topology lookup mesh
      TP.Halfedges.SplitEdge(HE);
      // ensure that each face has a starting halfedge on an even vertex
      if (P.Halfedges[Pair].AdjacentFace > -1) 
      	TP.Faces[P.Halfedges[Pair].AdjacentFace].FirstHalfedge = TP.Halfedges.Count - 1;

      // Loop
      Point3d NewPosition = P.Vertices[HeSV].ToPoint3d() * 0.375 +
      P.Vertices[PrSV].ToPoint3d() * 0.375 +
      P.Vertices[P.Halfedges[ThisOpp].StartVertex].ToPoint3d() * 0.125 +
      P.Vertices[P.Halfedges[PairOpp].StartVertex].ToPoint3d() * 0.125;
      NP.Vertices.Add(NewPosition);
    }
    for (int f = 0; f < P.Faces.Count; f++)
    {
      // add new triangulated faces
      int[] FaceVerts = TP.Faces.GetFaceVertices(f);
      NP.Faces.AddFace(FaceVerts[0], FaceVerts[1], FaceVerts[5]);
      NP.Faces.AddFace(FaceVerts[2], FaceVerts[3], FaceVerts[1]);
      NP.Faces.AddFace(FaceVerts[4], FaceVerts[5], FaceVerts[3]);
      NP.Faces.AddFace(FaceVerts[1], FaceVerts[3], FaceVerts[5]);
    }
    return NP.ToRhinoMesh();
  }

Catmull-Clark细分算法

原理

Catmull-Clark细分是一种基于四边形的网格细分算法，它所产生的网格能够达到 $C^{2}$ 连续，奇异点处为 $C^1$ 连续。和Loop类似，其实就是父的更新法则和子顶点的创建法则有不同。对于内部顶点，位于面中的子顶点创建法则为该面顶点的加权平均，位于边界的子顶点的创建法则是内部边两端点的 $\frac{3}{8}$ ，过该边的两张网格面的其他顶点占 $\frac{1}{16}$ 。父顶点的更新法则分两种，与父顶点相邻的顶点的权重为β/k（β=3/2k，k为顶点的价），与父顶点共面却不相邻的顶点的权重为γ/k（γ=1/4k），父顶点挪动前的权重为1-β-γ。边界顶点的处理和Loop算法完全相同，这里不再多做解释了。
Figure 9：Catmull-Clark细分法则图片来源：《Subdivision for Modeling and Animation》
位于面中的内部子顶点：设四边形的四个顶点分别为 $v_0$ ， $v_1$ ， $v_2$ ， $v_3$ ，要创建的子顶点为 $v_{new}$ ，则

v_{new} = \frac{v_0+v_1+v_2+v_3}{4}

位于内部边上的子顶点：设内部网格边的两端为 $v_0$ 和 $v_1$ ，与该边相邻的两个四边形顶点分别为 $v_0$ - $v_1$ - $v_2$ - $v_3$ 和 $v_0$ - $v_1$ - $v_4$ - $v_5$ 要创建的子顶点为 $v_{new}$ ，则

v_{new} = \frac{3}{8}(v_0+v_1)+\frac{1}{16}(v_2+v_3+v_4+v_5)

内部父顶点：设内部父顶点 $v_0$ 的相邻点为 $v_1$ ， $v_2$ ，...， $v_2n$ ，更新后的父顶点为 $v_{new}$ ，β和γ均为权重值， $n$ 为顶点的价（该顶点的相邻点的数量），则

β=\frac{3}{2n}$$$$γ=\frac{1}{4n}$$$$v_{new} = (1-β-γ)v_0+\frac{β}{n}\sum_{i=1}^{n} {v_{2i}}+\frac{γ}{n}\sum_{i=1}^{n} {v_{2i-1}}

边界子顶点：设边界边的两个端点为 $v_0$ 和 $v_1$ ，要创建的子顶点为 $v_{new}$ ，则

v_{new} = \frac{v_0+v_1}{2}

边界父顶点：设边界父顶点 $v_0$ 的相邻点为 $v_1$ 和 $v_2$ ，更新后的父顶点为 $v_{new}$ ，w为权重值， $n$ 为顶点的价（该顶点的相邻点的数量），则

v_{new} = \frac{3}{4}v_{0}+\frac{1}{8}(v_1+v_2)

代码实现

简略版的Catmull-Clark细分的C#源码如下：

public Mesh CatmullClark(Mesh m)
  {
    PlanktonMesh P = m.ToPlanktonMesh();
    PlanktonMesh NP = new PlanktonMesh();
    PlanktonMesh TP = m.ToPlanktonMesh();

    for (int v = 0; v < P.Vertices.Count; v++)
    {
      // Catmull-Clark
      int[] InHEs = P.Vertices.GetIncomingHalfedges(v);

      double Beta = 3.0 / (2 * InHEs.Length);
      double Delta = 1.0 / (4 * InHEs.Length);
      double BetaK = Beta / InHEs.Length;
      double DeltaK = Delta / InHEs.Length;

      Point3d NewPosition = P.Vertices[v].ToPoint3d() * (1 - Beta - Delta);

      foreach (int InHE in InHEs)
      {
        NewPosition += P.Vertices[P.Halfedges[InHE].StartVertex].ToPoint3d() * BetaK;
        NewPosition += P.Vertices[P.Halfedges[P.Halfedges[InHE].PrevHalfedge].StartVertex].ToPoint3d() * DeltaK;
      }
      NP.Vertices.Add(NewPosition);
    }
    for (int HE = 0; HE < P.Halfedges.Count; HE += 2)
    {
      int Pair = P.Halfedges.GetPairHalfedge(HE);
      int ThisOpp = P.Halfedges[HE].PrevHalfedge;
      int PairOpp = P.Halfedges[Pair].PrevHalfedge;

      int HeSV = P.Halfedges[HE].StartVertex;
      int PrSV = P.Halfedges[Pair].StartVertex;

      // split the halfedges in the topology lookup mesh
      TP.Halfedges.SplitEdge(HE);
      // ensure that each face has a starting halfedge on an even vertex
      if (P.Halfedges[Pair].AdjacentFace > -1) TP.Faces[P.Halfedges[Pair].AdjacentFace].FirstHalfedge = TP.Halfedges.Count - 1;

      // Catmull-Clark
      Point3d NewPosition = P.Vertices[HeSV].ToPoint3d() * 0.375 +
        P.Vertices[PrSV].ToPoint3d() * 0.375 +
        P.Vertices[P.Halfedges[P.Halfedges.GetPairHalfedge(P.Halfedges[HE].NextHalfedge)].StartVertex].ToPoint3d() * 0.0625 +
        P.Vertices[P.Halfedges[P.Halfedges[HE].PrevHalfedge].StartVertex].ToPoint3d() * 0.0625 +
        P.Vertices[P.Halfedges[P.Halfedges.GetPairHalfedge(P.Halfedges[Pair].NextHalfedge)].StartVertex].ToPoint3d() * 0.0625 +
        P.Vertices[P.Halfedges[P.Halfedges[Pair].PrevHalfedge].StartVertex].ToPoint3d() * 0.0625;
      NP.Vertices.Add(NewPosition);
    }
    for (int f = 0; f < P.Faces.Count; f++)
    {
      int CenterIdx = NP.Vertices.Count;
      NP.Vertices.Add(P.Faces.GetFaceCenter(f).ToPoint3d());
      int[] FaceVerts = TP.Faces.GetFaceVertices(f);
      for (int nf = 0; nf < FaceVerts.Length; nf += 2)
      {
        int LastVert = nf - 1;
        if (nf == 0) LastVert = FaceVerts.Length - 1;
        NP.Faces.AddFace(FaceVerts[nf], FaceVerts[nf + 1], CenterIdx, FaceVerts[LastVert]);
      }
    }
    return NP.ToRhinoMesh();
  }

Mesh is Art (4): Subdivision

前言

网格细分简述

细分算法原理

Loop细分算法

原理

代码实现

Catmull-Clark细分算法

原理

代码实现

【Mesh is Art】第一章常见几何表示方法

Mesh is Art (3): Half-Edge Data Structure

Mesh is Art (4): Subdivision

前言

网格细分简述

细分算法原理

Loop细分算法

原理

代码实现

Catmull-Clark细分算法

原理

代码实现

【Mesh is Art】第一章 常见几何表示方法

Mesh is Art (3): Half-Edge Data Structure

【Mesh is Art】第一章常见几何表示方法