7.3.2 Point-in-Polygon

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

/*

Given a point `p` and a polygon, determines whether `p` lies inside using ray casting. A ray cast
from `p` must cross the polygon's boundary an odd number of times exactly when `p` is inside, so the
function counts crossings between a horizontal ray and each polygon edge. The function is templated
on the point type `Pt` and accepts `Point`, `PointI`, or `PointL` from 7.1.1, or any compatible type
with numeric `.x` and `.y` fields. All comparisons are exact (no epsilon): the ray-crossing parity
needs a consistent comparison to be correct, and the result is exact for integer-coordinate points.

- `point_in_polygon(p, lo, hi, include_boundary = true)` returns whether `p` lies within the polygon
  with vertices specified by the range $[`lo`, `hi`)$ of points in either clockwise or
  counter-clockwise order. The `include_boundary` flag controls whether points exactly on an edge or
  vertex count as inside.

Overflow warning: the edge orientation test forms a cross product that grows like the squared
coordinate magnitude. For integer point types use a 64-bit coordinate type (e.g. `PointL` from
7.1.1) once coordinates exceed roughly a few tens of thousands.

Time Complexity:
- O(n) per call, where $n$ is the distance between `lo` and `hi`.

Space Complexity:
- O(1) auxiliary.

*/

template<typename 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);
}

// Detection is exact for integer-coordinate Pt.
template<typename Pt, typename It>
bool point_in_polygon(const Pt &p, It lo, It hi, bool include_boundary = true) {
  if (lo == hi) {
    return false;
  }
  bool res = false;
  for (It i = lo, j = hi - 1; i != hi; j = i++) {
    if (i->y == p.y && (i->x == p.x || (j->y == p.y && (i->x <= p.x || j->x <= p.x)))) {
      return include_boundary;
    }
    if ((p.y < i->y) != (p.y < j->y)) {
      auto det = cross(*i, *j, p);
      if (det == 0) {
        return include_boundary;
      }
      if ((0 < det) != (i->y < j->y)) {
        res = !res;
      }
    }
  }
  return res;
}

/*** Example Usage ***/

#include <cassert>
#include <vector>
using namespace std;

struct Point {
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
};

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

int main() {
  vector<Point> poly{Point(-1, 3), Point(1, 3), Point(2, 1), Point(0, 0)};
  assert(point_in_polygon(Point(1, 2), poly.begin(), poly.end()));
  assert(point_in_polygon(Point(0, 3), poly.begin(), poly.end()));
  assert(!point_in_polygon(Point(0, 3), poly.begin(), poly.end(), false));
  assert(!point_in_polygon(Point(0, 3.01), poly.begin(), poly.end()));
  assert(!point_in_polygon(Point(2, 2), poly.begin(), poly.end()));

  // Integer-coordinate points: exact detection.
  vector<PointI> ipoly{{0, 0}, {4, 0}, {4, 4}, {0, 4}};
  assert(point_in_polygon(PointI(2, 2), ipoly.begin(), ipoly.end()));
  assert(!point_in_polygon(PointI(4, 2), ipoly.begin(), ipoly.end(), false));
  assert(!point_in_polygon(PointI(5, 2), ipoly.begin(), ipoly.end()));
  return 0;
}