Triangulation

In geometry, triangulation refers to subdividing planar objects into triangles, and by extension, subdividing higher-dimensional geometric objects into simplices. For a given set of points, there are many triangulations, such as:

In OI, triangulation mainly refers to perfect triangulation in two-dimensional geometry (2D Delaunay triangulation, abbreviated as DT).

Delaunay Triangulation

Definition

In mathematics and computational geometry, for a given discrete point set \(P\) in the plane, its Delaunay triangulation DT(\(P\)) satisfies:

1. Empty circle property: DT(\(P\)) is unique (no four points are cocircular). In DT(\(P\)), the circumcircle of any triangle contains no other points.
2. Maximizing the minimum angle: Among all possible triangulations of point set \(P\), DT(\(P\)) maximizes the minimum angle of the triangles formed. In this sense, DT(\(P\)) is the most regular triangulation. Specifically, for two adjacent triangles forming a convex quadrilateral, swapping the diagonal does not increase the minimum of the six interior angles of the two triangles.

Properties

1. Closest: Triangles are formed by the three closest points, and all segments (triangle edges) do not intersect.
2. Uniqueness: Regardless of where construction starts, the final result is consistent (no four points in the set are cocircular).
3. Optimality: If the diagonal of any convex quadrilateral formed by two adjacent triangles can be swapped, then the minimum angle among the six interior angles of the two triangles does not change.
4. Most regular: If the minimum angles of each triangle in the triangulation are sorted in ascending order, the Delaunay triangulation yields the largest sequence.
5. Locality: Adding, deleting, or moving a vertex only affects nearby triangles.
6. Convex hull boundary: The outermost boundary of the triangulation forms a convex polygon hull.

Divide-and-Conquer Algorithm for Constructing DT

There are many algorithms for constructing DT. Among \(O(n \log n)\) construction algorithms, the divide-and-conquer algorithm is the easiest to understand and implement.

The first step of the divide-and-conquer DT construction is to sort the given point set in ascending order by \(x\)-coordinate. The figure below shows a sorted point set of size \(10\).

Once the point set is sorted, we can recursively split it into two parts (divide-and-conquer) until each subset contains at most \(3\) points. These subsets can then be immediately triangulated into a triangle or a line segment.

Then, during the backtracking phase of the divide-and-conquer, the already triangulated left and right subsets can be merged step by step. The merged triangulation contains LL-edges (edges from the left subset), RR-edges (edges from the right subset), and LR-edges (new edges connecting the left and right triangulations), as shown: LL-edges (gray), RR-edges (red), LR-edges (blue). To maintain the DT properties after merging, we may need to delete some LL-edges and RR-edges, but we do not add new LL-edges or RR-edges during merging.

The first step in merging the left and right triangulations is to insert the base LR-edge. The base LR-edge is the lowest LR-edge that does not intersect any LL-edge or RR-edge.

Next, we need to determine the next LR-edge immediately above the base LR-edge. For example, for the right point set, the possible endpoint (right endpoint) for the next LR-edge is the other endpoint of the RR-edge connected to the right endpoint of the base LR-edge (points \(6, 7, 9\)), and the left endpoint is point \(2\).

For each candidate endpoint, we need to check against the following two criteria:

1. The angle between its corresponding RR-edge and the base LR-edge is less than \(180\) degrees.
2. The circle defined by the two endpoints of the base LR-edge and this candidate point does not contain any other candidate points.

As shown above, the green circle corresponding to candidate point \(6\) contains candidate point \(9\), while the purple circle corresponding to candidate point \(7\) contains no other candidate points. Therefore, point \(7\) is the right endpoint of the next LR-edge.

For the left point set, we perform a mirror process.

When neither the left nor right point set contains any candidate points meeting the criteria, the merge is complete. When a candidate point meets the criteria, an LR-edge is added, and any LL-edges or RR-edges intersecting this new LR-edge are deleted.

When both left and right point sets have candidate points, check if the circle defined by the left candidate point contains the right candidate point; if it does, it does not meet the criteria. The same check applies to the right candidate point. Usually, only one candidate point meets the criteria (unless four points are cocircular).

Once this LR-edge is added, it becomes the new base LR-edge, and the steps are repeated to add the next edge until the merge is complete.

Code

Implementation

#include <algorithm>
#include <cmath>
#include <cstring>
#include <list>
#include <utility>
#include <vector>

constexpr double EPS = 1e-8;
constexpr int MAXV = 10000;

struct Point {
  double x, y;
  int id;

  Point(double a = 0, double b = 0, int c = -1) : x(a), y(b), id(c) {}

  bool operator<(const Point &a) const {
    return x < a.x || (fabs(x - a.x) < EPS && y < a.y);
  }

  bool operator==(const Point &a) const {
    return fabs(x - a.x) < EPS && fabs(y - a.y) < EPS;
  }

  double dist2(const Point &b) {
    return (x - b.x) * (x - b.x) + (y - b.y) * (y - b.y);
  }
};

struct Point3D {
  double x, y, z;

  Point3D(double a = 0, double b = 0, double c = 0) : x(a), y(b), z(c) {}

  Point3D(const Point &p) { x = p.x, y = p.y, z = p.x * p.x + p.y * p.y; }

  Point3D operator-(const Point3D &a) const {
    return Point3D(x - a.x, y - a.y, z - a.z);
  }

  double dot(const Point3D &a) { return x * a.x + y * a.y + z * a.z; }
};

struct Edge {
  int id;
  std::list<Edge>::iterator c;

  Edge(int id = 0) { this->id = id; }
};

int cmp(double v) { return fabs(v) > EPS ? (v > 0 ? 1 : -1) : 0; }

double cross(const Point &o, const Point &a, const Point &b) {
  return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

Point3D cross(const Point3D &a, const Point3D &b) {
  return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
                 a.x * b.y - a.y * b.x);
}

int inCircle(const Point &a, Point b, Point c, const Point &p) {
  if (cross(a, b, c) < 0) std::swap(b, c);
  Point3D a3(a), b3(b), c3(c), p3(p);
  b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
  Point3D f = cross(b3, c3);
  return cmp(p3.dot(f));  // check same direction, in: < 0, on: = 0, out: > 0
}

int intersection(const Point &a, const Point &b, const Point &c,
                 const Point &d) {  // seg(a, b) and seg(c, d)
  return cmp(cross(a, c, b)) * cmp(cross(a, b, d)) > 0 &&
         cmp(cross(c, a, d)) * cmp(cross(c, d, b)) > 0;
}

class Delaunay {
 public:
  std::list<Edge> head[MAXV];  // graph
  Point p[MAXV];
  int n, rename[MAXV];

  void init(int n, Point p[]) {
    memcpy(this->p, p, sizeof(Point) * n);
    std::sort(this->p, this->p + n);
    for (int i = 0; i < n; i++) rename[p[i].id] = i;
    this->n = n;
    divide(0, n - 1);
  }

  void addEdge(int u, int v) {
    head[u].push_front(Edge(v));
    head[v].push_front(Edge(u));
    head[u].begin()->c = head[v].begin();
    head[v].begin()->c = head[u].begin();
  }

  void divide(int l, int r) {
    if (r - l <= 2) {  // #point <= 3
      for (int i = l; i <= r; i++)
        for (int j = i + 1; j <= r; j++) addEdge(i, j);
      return;
    }
    int mid = (l + r) / 2;
    divide(l, mid);
    divide(mid + 1, r);

    std::list<Edge>::iterator it;
    int nowl = l, nowr = r;

    for (int update = 1; update;) {
      // find left and right convex, lower common tangent
      update = 0;
      Point ptL = p[nowl], ptR = p[nowr];
      for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
        Point t = p[it->id];
        double v = cross(ptR, ptL, t);
        if (cmp(v) > 0 || (cmp(v) == 0 && ptR.dist2(t) < ptR.dist2(ptL))) {
          nowl = it->id, update = 1;
          break;
        }
      }
      if (update) continue;
      for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
        Point t = p[it->id];
        double v = cross(ptL, ptR, t);
        if (cmp(v) < 0 || (cmp(v) == 0 && ptL.dist2(t) < ptL.dist2(ptR))) {
          nowr = it->id, update = 1;
          break;
        }
      }
    }

    addEdge(nowl, nowr);  // add tangent

    for (int update = 1; true;) {
      update = 0;
      Point ptL = p[nowl], ptR = p[nowr];
      int ch = -1, side = 0;
      for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
        if (cmp(cross(ptL, ptR, p[it->id])) > 0 &&
            (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
          ch = it->id, side = -1;
        }
      }
      for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
        if (cmp(cross(ptR, p[it->id], ptL)) > 0 &&
            (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
          ch = it->id, side = 1;
        }
      }
      if (ch == -1) break;  // upper common tangent
      if (side == -1) {
        for (it = head[nowl].begin(); it != head[nowl].end();) {
          if (intersection(ptL, p[it->id], ptR, p[ch])) {
            head[it->id].erase(it->c);
            head[nowl].erase(it++);
          } else {
            it++;
          }
        }
        nowl = ch;
        addEdge(nowl, nowr);
      } else {
        for (it = head[nowr].begin(); it != head[nowr].end();) {
          if (intersection(ptR, p[it->id], ptL, p[ch])) {
            head[it->id].erase(it->c);
            head[nowr].erase(it++);
          } else {
            it++;
          }
        }
        nowr = ch;
        addEdge(nowl, nowr);
      }
    }
  }

  std::vector<std::pair<int, int>> getEdge() {
    std::vector<std::pair<int, int>> ret;
    ret.reserve(n);
    std::list<Edge>::iterator it;
    for (int i = 0; i < n; i++) {
      for (it = head[i].begin(); it != head[i].end(); it++) {
        if (it->id < i) continue;
        ret.push_back(std::make_pair(p[i].id, p[it->id].id));
      }
    }
    return ret;
  }
};

Voronoi Diagram

A Voronoi diagram consists of a set of continuous polygons formed by the perpendicular bisectors of lines connecting neighboring points. Given \(n\) distinct seed points in the plane, the plane is divided into \(n\) regions such that every point in a region is closer to its seed point than to any other seed point.

The Voronoi diagram is the dual graph of the Delaunay triangulation. It can be constructed in \(O(n \log n)\) time by first constructing the Delaunay triangulation using a divide-and-conquer algorithm and then applying the leftmost turn algorithm to obtain its dual graph.

Problems

SGU 383 Caravans Triangulation + Binary Lifting

ContestHunter. Endless Destruction Construct Voronoi diagram via dual graph of triangulation

Codeforces Gym 103485M. Constellation collection Graph construction after triangulation and Floodfill

References and Further Reading

1. Wikipedia - Triangulation (geometry)
2. Wikipedia - Delaunay triangulation
3. Samuel Peterson -Computing Constrained Delaunay Triangulations in 2-D (1997-98)

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用