平面剖分三角形Delaunay算法实现，以及扩展到立体空间中的实现（附源码）

平面剖分三角形Delaunay算法实现

前端版实现：https://github.com/ccinos/TriangleMesh

Delaunay2D.cs下载

  void Triangulate() {
    /* 基于Bowyer-Watson算法原理
     * 通过所有顶点和现有三角形进行判断, 对包含该点的三角形进行拆分
     * 首先构建初始的三角形,包含全部点在内(最后会删除)
     * 
     * 遍历所有顶点 {
     *    遍历所有现有三角形 {
     *        首先判断是否 [顶点] 在 [三角形] 内, 如果不在则跳过
     *        [顶点] 和 [三角形] 的3条边,生成新的3个三角形,加入现有三角形集合中
     *    }
     * }
     */
    #region 构建初始三角形,可以包含住所有的点
    float minX = Vertices[0].Position.x;
    float minY = Vertices[0].Position.y;
    float maxX = minX;
    float maxY = minY;

    foreach (var vertex in Vertices) {
      if (vertex.Position.x < minX) minX = vertex.Position.x;
      if (vertex.Position.x > maxX) maxX = vertex.Position.x;
      if (vertex.Position.y < minY) minY = vertex.Position.y;
      if (vertex.Position.y > maxY) maxY = vertex.Position.y;
    }

    float dx = maxX - minX;
    float dy = maxY - minY;
    float deltaMax = Mathf.Max(dx, dy) * 2;

    Vertex p1 = new Vertex(new Vector2(minX - 1, minY - 1));
    Vertex p2 = new Vertex(new Vector2(minX - 1, maxY + deltaMax));
    Vertex p3 = new Vertex(new Vector2(maxX + deltaMax, minY - 1));

    Triangles.Add(new Triangle(p1, p2, p3));
    #endregion

    #region 通过顶点对三角形进行拆分
    //遍历所有顶点
    foreach (var vertex in Vertices) {
      List<Edge> polygon = new List<Edge>();
      //遍历所有现有三角形
      foreach (var t in Triangles) {
        if (t.CircumCircleContains(vertex.Position)) {//如果包含该点
          //标记删除
          t.IsBad = true;
          //记录下三条边
          polygon.Add(new Edge(t.A, t.B));
          polygon.Add(new Edge(t.B, t.C));
          polygon.Add(new Edge(t.C, t.A));
        }
      }
      //删除标记的三角形
      Triangles.RemoveAll((Triangle t) => t.IsBad);
      //删除近似的边
      for (int i = 0; i < polygon.Count; i++) {
        for (int j = i + 1; j < polygon.Count; j++) {
          if (Edge.AlmostEqual(polygon[i], polygon[j])) {
            polygon[i].IsBad = true;
            polygon[j].IsBad = true;
          }
        }
      }
      polygon.RemoveAll((Edge e) => e.IsBad);
      //通过点和边构建新的三角形,加入到现有三角形集合中
      foreach (var edge in polygon) {
        Triangles.Add(new Triangle(edge.U, edge.V, vertex));
      }
    }
    //删除所有包含初始点的三角形
    Triangles.RemoveAll((Triangle t) => t.ContainsVertex(p1.Position) || t.ContainsVertex(p2.Position) || t.ContainsVertex(p3.Position));
    #endregion

    #region 提取边信息
    HashSet<Edge> edgeSet = new HashSet<Edge>();

    foreach (var t in Triangles) {
      var ab = new Edge(t.A, t.B);
      var bc = new Edge(t.B, t.C);
      var ca = new Edge(t.C, t.A);

      if (edgeSet.Add(ab)) {
        Edges.Add(ab);
      }

      if (edgeSet.Add(bc)) {
        Edges.Add(bc);
      }

      if (edgeSet.Add(ca)) {
        Edges.Add(ca);
      }
    }
    #endregion
  }

扩展立体空间中，三角形剖分变为四面体剖分：

Delaunay3D.cs下载

  void Triangulate() {
    /* 基于Bowyer-Watson算法原理
     * 通过所有顶点和现有四面体进行判断, 对包含该点的四面体进行拆分
     * 首先构建初始的四面体,包含全部点在内,并在最后删除
     * 
     * 遍历所有顶点 {
     *    遍历所有现有四面体 {
     *        首先判断是否 [顶点] 在 [四面体] 内, 如果不在则跳过
     *        [顶点] 和 [四面体] 的四个三角形,生成新的四个四面体,加入现有四面体集合中
     *    }
     * }
     * 
     * 通过以上步骤可以得到所有点相关的四面体
     * 将这些四面体的三角形信息拿出来即可
     */
    #region 构建初始四面体,可以包含住所有点的四面体
    float minX = Vertices[0].Position.x;
    float minY = Vertices[0].Position.y;
    float minZ = Vertices[0].Position.z;
    float maxX = minX;
    float maxY = minY;
    float maxZ = minZ;

    //1.计算顶点中最大和最小的坐标 x,y,z => minX,minY,minZ,maxX,maxY,maxZ
    foreach (var vertex in Vertices) {
      if (vertex.Position.x < minX) minX = vertex.Position.x;
      if (vertex.Position.x > maxX) maxX = vertex.Position.x;
      if (vertex.Position.y < minY) minY = vertex.Position.y;
      if (vertex.Position.y > maxY) maxY = vertex.Position.y;
      if (vertex.Position.z < minZ) minZ = vertex.Position.z;
      if (vertex.Position.z > maxZ) maxZ = vertex.Position.z;
    }
    //2.计算坐标差 dx,dy,dz
    float dx = maxX - minX;
    float dy = maxY - minY;
    float dz = maxZ - minZ;
    float deltaMax = Mathf.Max(dx, dy, dz) * 2;

    Vertex p1 = new Vertex(new Vector3(minX - 1, minY - 1, minZ - 1));
    Vertex p2 = new Vertex(new Vector3(maxX + deltaMax, minY - 1, minZ - 1));
    Vertex p3 = new Vertex(new Vector3(minX - 1, maxY + deltaMax, minZ - 1));
    Vertex p4 = new Vertex(new Vector3(minX - 1, minY - 1, maxZ + deltaMax));

    this.Tetrahedra.Add(new Tetrahedron(p1, p2, p3, p4));
    #endregion
    
    #region 通过顶点对四面体进行拆分
    //遍历所有顶点
    foreach (var vertex in Vertices) {
      List<Triangle> triangles = new List<Triangle>();
      //遍历所有现有四面体
      foreach (var t in this.Tetrahedra) {
        if (t.CircumCircleContains(vertex.Position)) {
          //如果包含则标记该四面体需要删除
          t.IsBad = true;
          //记录四面体的4个三角形
          triangles.Add(new Triangle(t.A, t.B, t.C));
          triangles.Add(new Triangle(t.A, t.B, t.D));
          triangles.Add(new Triangle(t.A, t.C, t.D));
          triangles.Add(new Triangle(t.B, t.C, t.D));
        }
      }
      //判断四面体4个三角形是否近似相同,相同则都删除
      for (int i = 0; i < triangles.Count; i++) {
        for (int j = i + 1; j < triangles.Count; j++) {
          if (Triangle.AlmostEqual(triangles[i], triangles[j])) {
            triangles[i].IsBad = true;
            triangles[j].IsBad = true;
          }
        }
      }
      //删除带标记的四面体和三角形
      this.Tetrahedra.RemoveAll((Tetrahedron t) => t.IsBad);
      triangles.RemoveAll((Triangle t) => t.IsBad);
      //将剩余三角形和该顶点重新构成四面体
      foreach (var triangle in triangles) {
        this.Tetrahedra.Add(new Tetrahedron(triangle.U, triangle.V, triangle.W, vertex));
      }
    }

    //四面体中删除所有含有初始四面体点的四面体
    this.Tetrahedra.RemoveAll((Tetrahedron t) => t.ContainsVertex(p1) || t.ContainsVertex(p2) || t.ContainsVertex(p3) || t.ContainsVertex(p4));
    #endregion

    #region 提取三角形和边的信息
    HashSet<Triangle> triangleSet = new HashSet<Triangle>();
    HashSet<Edge> edgeSet = new HashSet<Edge>();

    foreach (var t in this.Tetrahedra) {
      var abc = new Triangle(t.A, t.B, t.C);
      var abd = new Triangle(t.A, t.B, t.D);
      var acd = new Triangle(t.A, t.C, t.D);
      var bcd = new Triangle(t.B, t.C, t.D);

      if (triangleSet.Add(abc)) {
        this.Triangles.Add(abc);
      }

      if (triangleSet.Add(abd)) {
        this.Triangles.Add(abd);
      }

      if (triangleSet.Add(acd)) {
        this.Triangles.Add(acd);
      }

      if (triangleSet.Add(bcd)) {
        this.Triangles.Add(bcd);
      }

      var ab = new Edge(t.A, t.B);
      var bc = new Edge(t.B, t.C);
      var ca = new Edge(t.C, t.A);
      var da = new Edge(t.D, t.A);
      var db = new Edge(t.D, t.B);
      var dc = new Edge(t.D, t.C);

      if (edgeSet.Add(ab)) {
        Edges.Add(ab);
      }

      if (edgeSet.Add(bc)) {
        Edges.Add(bc);
      }

      if (edgeSet.Add(ca)) {
        Edges.Add(ca);
      }

      if (edgeSet.Add(da)) {
        Edges.Add(da);
      }

      if (edgeSet.Add(db)) {
        Edges.Add(db);
      }

      if (edgeSet.Add(dc)) {
        Edges.Add(dc);
      }
    }
    #endregion
  }

源码：

Delaunay2D.cs

Delaunay3D.cs