7.4.5 Delaunay Triangulation (Simple)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

/*

Given a set $P$ of distinct two-dimensional points, a Delaunay triangulation of $P$ is a
triangulation of the convex hull of P such that no point of P lies strictly inside the circumcircle
of any triangle in the triangulation. The Delaunay triangulation is not necessarily unique when four
or more points are cocircular.

This implementation produces one valid triangulation using a simple brute-force algorithm. Each
candidate triangle is tested against the empty-circumcircle condition, and candidates whose edges
would properly cross already accepted triangles are rejected.

- `delaunay_triangulation(lo, hi)` returns a Delaunay triangulation for a range [`lo`, `hi`) of
  distinct points, where `lo` and `hi` must be random-access iterators, or an empty vector if a
  triangulation does not exist.

Time Complexity:
- O(n^6) per call to `delaunay_triangulation(lo, hi)`, where $n$ is the distance between `lo` and
  `hi`. The empty-circumcircle test contributes O(n^4); the remaining factor comes from checking
  candidate triangles against previously accepted triangles.

Space Complexity:
- O(n) auxiliary heap space, excluding the returned triangulation.

*/

#include <algorithm>
#include <array>
#include <cmath>
#include <utility>
#include <vector>

const double EPS = 1e-9;

#define EQ(a, b) (fabs((a) - (b)) <= EPS)
#define LT(a, b) ((a) < (b) - EPS)
#define LE(a, b) ((a) <= (b) + EPS)

// Specialized version of seg_intersection() from 7.2.3, simplified for TOUCH_IS_INTERSECT = false,
// i.e. we're detecting proper crossings of two non-parallel segments whose interiors intersect.
template<class Pt>
bool proper_seg_intersection(const Pt &a, const Pt &b, const Pt &c, const Pt &d) {
  auto cross = [](const Pt &o, const Pt &p, const Pt &q) {
    return (p.x - o.x) * (q.y - o.y) - (p.y - o.y) * (q.x - o.x);  // Overflow warning!
  };
  auto c1 = cross(a, b, c);
  auto c2 = cross(a, b, d);
  auto c3 = cross(c, d, a);
  auto c4 = cross(c, d, b);
  return ((LT(c1, 0) && LT(0, c2)) || (LT(c2, 0) && LT(0, c1))) &&
         ((LT(c3, 0) && LT(0, c4)) || (LT(c4, 0) && LT(0, c3)));
}

struct Triangle {
  double ax, ay, bx, by, cx, cy;

  Triangle(double ax, double ay, double bx, double by, double cx, double cy)
      : ax(ax), ay(ay), bx(bx), by(by), cx(cx), cy(cy) {}

  bool operator==(const Triangle &t) const {
    return EQ(ax, t.ax) && EQ(ay, t.ay) && EQ(bx, t.bx) && EQ(by, t.by) && EQ(cx, t.cx) &&
           EQ(cy, t.cy);
  }
};

// Accepts any point type with numeric .x/.y; coordinates are copied to double arrays.
template<class It>
std::vector<Triangle> delaunay_triangulation(It lo, It hi) {
  using Pt = typename std::iterator_traits<It>::value_type;
  int n = hi - lo;
  std::vector<Pt> pts;
  std::vector<double> x, y, z;
  for (It it = lo; it != hi; ++it) {
    pts.emplace_back(it->x, it->y);
    x.emplace_back(it->x);
    y.emplace_back(it->y);
    z.emplace_back(static_cast<double>(it->x) * it->x + static_cast<double>(it->y) * it->y);
  }
  std::vector<Triangle> res;
  std::vector<std::array<int, 3>> res_idx;
  for (int i = 0; i < n - 2; i++) {
    for (int j = i + 1; j < n; j++) {
      for (int k = i + 1; k < n; k++) {
        if (j == k) {
          continue;
        }
        double nx = (y[j] - y[i]) * (z[k] - z[i]) - (y[k] - y[i]) * (z[j] - z[i]);
        double ny = (x[k] - x[i]) * (z[j] - z[i]) - (x[j] - x[i]) * (z[k] - z[i]);
        double nz = (x[j] - x[i]) * (y[k] - y[i]) - (x[k] - x[i]) * (y[j] - y[i]);
        if (LE(0, nz)) {
          continue;
        }
        std::vector<Pt> s1{pts[i], pts[j], pts[k], pts[i]};
        for (int m = 0; m < n; m++) {
          if (LT(0, nx * (x[m] - x[i]) + ny * (y[m] - y[i]) + nz * (z[m] - z[i]))) {
            goto skip;
          }
        }
        // Reject candidates whose edges properly cross already accepted triangles.
        for (const auto &tri : res_idx) {
          std::vector<Pt> s2{pts[tri[0]], pts[tri[1]], pts[tri[2]], pts[tri[0]]};
          for (int u = 0; u < 3; u++) {
            for (int v = 0; v < 3; v++) {
              if (proper_seg_intersection(s1[u], s1[u + 1], s2[v], s2[v + 1])) {
                goto skip;
              }
            }
          }
        }
        res.emplace_back(pts[i].x, pts[i].y, pts[j].x, pts[j].y, pts[k].x, pts[k].y);
        res_idx.push_back({i, j, k});
      skip:;
      }
    }
  }
  return res;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

struct Point {
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
  bool operator==(const Point &p) const { return EQ(x, p.x) && EQ(y, p.y); }
  bool operator<(const Point &p) const { return x != p.x ? x < p.x : y < p.y; }
};

struct PointI {
  int x, y;
  PointI(int x = 0, int y = 0) : x(x), y(y) {}
  bool operator==(const PointI &p) const { return x == p.x && y == p.y; }
  bool operator<(const PointI &p) const { return x != p.x ? x < p.x : y < p.y; }
};

int main() {
  vector<Point> v{{1, 3}, {1, 2}, {2, 1}, {0, 0}, {-1, 3}};
  vector<Triangle> t{
      Triangle(1, 3, 1, 2, -1, 3), Triangle(1, 3, 2, 1, 1, 2), Triangle(1, 2, 2, 1, 0, 0),
      Triangle(1, 2, 0, 0, -1, 3)
  };
  assert(delaunay_triangulation(v.begin(), v.end()) == t);

  // Integer-coordinate inputs are accepted (triangulation computed in double).
  vector<PointI> iv{{1, 3}, {1, 2}, {2, 1}, {0, 0}, {-1, 3}};
  assert(delaunay_triangulation(iv.begin(), iv.end()) == t);
  return 0;
}