三角剖分

In geometry, triangulation refers to subdividing planar objects into triangles, and subdividing high-dimensional geometric objects into simplex using extension. For a given point set, there are many forms of triangulation, such as:

3 types of triangulation

The triangulation in OI mainly refers to the perfect triangulation in two-dimensional geometry (two-dimensional Delaunay triangulation, DT for abbreviation).

Delaunay triangulation

Definition

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

  1. Empty circumcircle: DT( P ) is unique (any four points cannot be in a circle), and there exists no other points in the circumcircle of any triangles in DT( P ).
  2. Maximize the smallest angle: in the possible triangulations of the point set P , the smallest angle of the triangle formed by DT( P ) is the largest. In this sense, DT( P ) is the triangulation that is closest to regularization. Specifically, a convex quadrilateral diagonals formed by two adjacent triangles. After exchange, the smallest angle of the two internal angles no longer increases.

Properties

  1. Closest: a triangle is formed by the three closest points, and all line segments (sides of the triangle) do not intersect.
  2. Uniqueness: no matter which point in the area is used as the start point of the construction, the final result will be consistent (any four points in the point set cannot be a circle).
  3. Optimality: if the diagonals of a convex quadrilateral formed by any two adjacent triangles are interchangeable, the smallest angle among the six internal angles of the two triangles will not change.
  4. The extreme rule: If the smallest angle of each triangle in the triangulation is sorted in ascending order, the Delaunay triangulation will get the largest value.
  5. Regionality: adding, deleting, or moving any vertex will only affect the adjacent triangle.
  6. Convex shell: the outermost boundary of the triangulation forms a convex polygonal shell.

Construct a divide and conquer algorithm for DT

DT has many construction algorithms. Among the construction algorithms with time complexity of O(n \log n) , the divide and conquer algorithm is the easiest to understand and implement.

The first step of constructing DT by divide and conquer is to sort the given point set according to the coordinates of x in ascending order. The following figure shows the sorted point set of size 10 .

Sorted point set of size 1010

Once the point set is ordered, we can continuously divide it into two parts (divide and conquer) until the size of the sub-point set does not exceed 3 . Then these sub-point sets can be immediately divided into a triangle or line segments.

Divide and conquer as a point set containing 2 or 3 points

Then in the process of divide and conquer backtracking, the left and right sub-point sets that have been divided can be merged sequentially. The merged division contains LL-edge (the edge of the left sub-point set), RR-edge (the edge of the right sub-point set), and LR-edge (the new edge generated by connecting the left and right divisions), as shown in the figure LL-edge (gray), RR-edge (red), LR-edge (blue). For the merged subdivision, in order to maintain the properties of DT, we may need to delete some LL-edge and RR-edge, but we will not add LL-edge and RR-edge when merging.

Edge

The first step to merge the left and right splits is to insert the base LR-edge, which is the bottom LR-edge that does not intersects any LL-edge and RR-edge.

Merge left and right splits

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

Next LR-edge

For possible endpoints, we need to check according to the following two standards:

  1. The angle between the corresponding RR-edge and base LR-edge is less than 180 degrees.
  2. The circle formed by the two end points of the base LR-edge and the three possible points does not contain any other possible points.

Check possible points

As shown in the figure above, the green circle corresponding to the possible point of 6 contains the possible point of 9 , and the purple circle corresponding to the possible point of 7 does not contain any other possible points, so the point of 7 is the right end point of the next LR-edge.

For the left point set, we can do mirror processing.

As shown in the figure above, the green circle corresponding to the possible point of 6 contains the possible point of 9 , and the purple circle corresponding to the possible point of 7 does not contain any other possible points, so the point of 7 is the next LR-edge The right end point.

For the left point set, we can do mirror processing.

Check possible points on the left side

When the left and right point sets no longer contain possible points that meet the standard, the merge is completed. When a possible point meets the standard, an LR-edge needs to be added, and the LL-edge and RR-edge that intersect with the LR-edge that needs to be added are deleted.

When there are possible points in the left and right point sets, check whether the circle corresponding to the left point contains the right point, and if it does, there is no match; This works the same for the right point. Generally, only one possible point meets the standard (unless the four points are in a circle).

Next LR-edge

After this LR-edge is added, use it as the base LR-edge and repeat the above steps, and continue to add the next one until the merge is completed.

Merge

Code

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

const double EPS = 1e-8;
const 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

The Voronoi diagram is composed of a set of continuous polygons composed of vertical bisectors connecting two adjacent points. According to n seed points that do not overlap on the plane, the plane is divided into n regions, so that the points in each region The distance to the seed point in its area is shorter than the distance to the seed points in other areas.

The Voronoi diagram is the dual graph of Delaunay triangulation. You can use the divide and conquer algorithm to construct the Delaunay triangulation to find the triangulation, and then use the leftmost line algorithm to find the dual graph. It is implemented in O(n \log n) for constructing the Voronoi diagram.

Problems

SGU 383 Caravans Triangulation + multiplication

CodeVS 2819. Endless destroy Triangulation to find dual graph for building the Voronoi diagram

References

buildLast update and/or translate time of this article，Check the history
editFound smelly bugs? Translation outdated? Wanna contribute with us? Edit this Page on Github
peopleContributor of this article xehoth
translateTranslator of this article Visit the original article!
copyrightThe article is available under CC BY-SA 4.0 & SATA ; additional terms may apply.

Comments