7.3.4 Convex Polygon Queries

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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199

/*

Given a convex polygon in counter-clockwise order, answer containment and tangent/intersection
queries. Point containment is logarithmic by splitting the polygon into triangles sharing `poly[0]`.
The line and tangent helpers below are written as simple linear scans; this is often shorter and
less assumption-heavy than the classic logarithmic versions, while still documenting the exact query
semantics in one place.

- `point_in_convex_polygon(poly, p, EDGE_IS_INSIDE)` returns whether point `p` lies inside the
  convex polygon `poly`. The polygon must be in counter-clockwise order, with no repeated final
  vertex. Points on the boundary count as inside unless `EDGE_IS_INSIDE` is set to `false`.
- `line_convex_polygon_intersection(poly, a, b)` describes where the infinite directed line through
  `a` $\to$ `b` hits the polygon: ${-1, -1}$ for no intersection, ${i, -1}$ for touching vertex $i$,
  ${i, i}$ for containing side $i \to i+1$, or ${i, j}$ for crossing sides $i \to i+1$ and
  $j \to j+1$.
- `convex_polygon_tangents(poly, p)` returns the two tangent vertex indices from an outside point
  `p` as (left, right). The left tangent has all polygon vertices on or to the left of the directed
  line `p` $\to$ `poly[left]`; the right tangent is analogous for the right side.

For integer-coordinate inputs, all tests are exact provided the cross products do not overflow.

Time Complexity:
- O(log n) per query, where $n$ is the number of polygon vertices.
- O(n) per call to `line_convex_polygon_intersection` and `convex_polygon_tangents`.

Space Complexity:
- O(1) auxiliary.

*/

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <vector>

template<class Pt>
auto cross(const Pt &a, const Pt &b, const Pt &o) {
  // Overflow risk for integer Pt: these products are ~O(max_coord^2); use int64_t if necessary.
  return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

template<class T>
int signum(const T &x) {
  return (T(0) < x) - (x < T(0));
}

template<class Pt>
bool on_segment(const Pt &p, const Pt &a, const Pt &b) {
  return cross(a, b, p) == 0 && std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x) &&
         std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}

template<class Pt>
bool point_in_convex_polygon(
    const std::vector<Pt> &poly, const Pt &p, const bool EDGE_IS_INSIDE = true
) {
  int n = static_cast<int>(poly.size());
  if (n == 0) {
    return false;
  }
  if (n == 1) {
    return EDGE_IS_INSIDE && p.x == poly[0].x && p.y == poly[0].y;
  }
  if (n == 2) {
    return EDGE_IS_INSIDE && on_segment(p, poly[0], poly[1]);
  }
  int left = signum(cross(poly[1], p, poly[0]));
  int right = signum(cross(poly[n - 1], p, poly[0]));
  if (left < 0 || right > 0) {
    return false;
  }
  if (left == 0) {
    return EDGE_IS_INSIDE && on_segment(p, poly[0], poly[1]);
  }
  if (right == 0) {
    return EDGE_IS_INSIDE && on_segment(p, poly[0], poly[n - 1]);
  }
  int lo = 1, hi = n - 1;
  while (hi - lo > 1) {
    int mid = (lo + hi) / 2;
    if (cross(poly[mid], p, poly[0]) >= 0) {
      lo = mid;
    } else {
      hi = mid;
    }
  }
  int s = signum(cross(poly[lo + 1], p, poly[lo]));
  return EDGE_IS_INSIDE ? s >= 0 : s > 0;
}

template<class Pt>
std::array<int, 2> line_convex_polygon_intersection(
    const std::vector<Pt> &poly, const Pt &a, const Pt &b
) {
  int n = static_cast<int>(poly.size());
  std::vector<int> on, cross_edges;
  bool has_pos = false, has_neg = false;
  for (int i = 0; i < n; i++) {
    int s = signum(cross(b, poly[i], a));
    has_pos |= s > 0;
    has_neg |= s < 0;
    if (s == 0) {
      on.push_back(i);
    }
  }
  if (!has_pos && !has_neg) {
    return n >= 2 ? std::array<int, 2>{0, 0} : std::array<int, 2>{-1, -1};
  }
  if (!(has_pos && has_neg)) {
    for (int i = 0; i < static_cast<int>(on.size()); i++) {
      int j = (i + 1) % static_cast<int>(on.size());
      if ((on[i] + 1) % n == on[j]) {
        return {on[i], on[i]};
      }
    }
    return on.empty() ? std::array<int, 2>{-1, -1} : std::array<int, 2>{on[0], -1};
  }
  for (int i = 0; i < n; i++) {
    int j = (i + 1) % n;
    int si = signum(cross(b, poly[i], a)), sj = signum(cross(b, poly[j], a));
    if (si == 0 || si * sj < 0) {
      cross_edges.push_back(i);
    }
  }
  if (cross_edges.size() == 1) {
    return {cross_edges[0], -1};
  }
  auto line_parameter = [&](int i) {
    int j = (i + 1) % n;
    auto ex = poly[j].x - poly[i].x, ey = poly[j].y - poly[i].y;
    auto ax = a.x, ay = a.y, dx = b.x - a.x, dy = b.y - a.y;
    auto det = dx * ey - dy * ex;
    return static_cast<double>((poly[i].x - ax) * ey - (poly[i].y - ay) * ex) / det;
  };
  if (line_parameter(cross_edges[1]) < line_parameter(cross_edges[0])) {
    std::swap(cross_edges[0], cross_edges[1]);
  }
  return {cross_edges[0], cross_edges[1]};
}

template<class Pt>
std::array<int, 2> convex_polygon_tangents(const std::vector<Pt> &poly, const Pt &p) {
  int n = static_cast<int>(poly.size());
  std::array<int, 2> res{-1, -1};
  for (int i = 0; i < n; i++) {
    bool left = true, right = true;
    for (int j = 0; j < n; j++) {
      int s = signum(cross(poly[i], poly[j], p));
      left &= s >= 0;
      right &= s <= 0;
    }
    if (left) {
      res[0] = i;
    }
    if (right) {
      res[1] = i;
    }
  }
  return res;
}

/*** Example Usage ***/

#include <algorithm>
using namespace std;

struct PointI {
  int64_t x, y;
  PointI(int64_t x = 0, int64_t y = 0) : x(x), y(y) {}
};

int main() {
  vector<PointI> square{{0, 0}, {4, 0}, {4, 4}, {0, 4}};
  assert(point_in_convex_polygon(square, PointI(2, 2)));
  assert(point_in_convex_polygon(square, PointI(4, 2)));
  assert(!point_in_convex_polygon(square, PointI(4, 2), false));
  assert(!point_in_convex_polygon(square, PointI(5, 2)));

  vector<PointI> pentagon{{0, 0}, {3, 0}, {5, 2}, {2, 5}, {-1, 3}};
  assert(point_in_convex_polygon(pentagon, PointI(2, 3)));
  assert(point_in_convex_polygon(pentagon, PointI(0, 0)));
  assert(!point_in_convex_polygon(pentagon, PointI(0, 0), false));
  assert(!point_in_convex_polygon(pentagon, PointI(5, 4)));

  assert(
      (line_convex_polygon_intersection(square, PointI(-1, 2), PointI(5, 2)) == array<int, 2>{3, 1})
  );
  assert(
      (line_convex_polygon_intersection(square, PointI(0, 0), PointI(4, 0)) == array<int, 2>{0, 0})
  );
  assert((
      line_convex_polygon_intersection(square, PointI(5, 0), PointI(5, 4)) == array<int, 2>{-1, -1}
  ));

  auto tangents = convex_polygon_tangents(square, PointI(6, 2));
  assert((tangents == array<int, 2>{2, 1}));
  return 0;
}