I decided to share my implementations for the basic polygon algorithms. I see almost no problems on this topic and I hope this will change in the future.

First, let's remind the definitions we will use:

1. Polygon is a plane figure that is bounded by a finite chain of straight line segments closing in a loop to form a closed chain or circuit. These segments are called its edges or sides, and the points where two edges meet are the polygon's vertices or corners (wiki).

2. Polygon is convex if a line segment connecting any two points on its boundary lies inside the polygon. Equivalently, all its interior angles are less than or equal to 180 degrees.

3. Polygon is strictly convex if in addition no three vertices lie on the same line. Equivalently, all its interior angles are less than 180 degrees.

4. Polygon is simple if its boundary doesn't cross itself.

I will present two algorithms for each problem: one for arbitrary simple polygon, and one for strictly convex polygon, that has better complexity. The priorities in implementation design were as follows:

• handle all the corner cases, except for degenerate polygons with zero area;
• perform all computations in integers;
• optimize performance;
• write as concise and clear code as possible.

Please, let me know if you find a way to improve on any of these goals in any algorithm listed. First, some references.

1. The code in this article will make no sense unless you read my previous one. I will use structures and functions defined there.

2. I highly recommend reading Computational Geometry in C (2nd Edition) by Joseph O'Rourke (google for pdf version). This is truly fundamental work on this topic! While the code in there is far from perfect, the ideas are described well in my opinion. I will indicate specific pages when I reference this book.

3. Some algorithms can be found on geomalgorithms.com, but I dislike the implementations. I took only one from there.

I define Polygon as a vector of Points. Polygon size is denoted as N in the complexity formulas. Since polygon is a chain, the next two functions help in jumping from the last vertex to the first and vice versa when necessary. After writing all the algorithms I realized it would be better to implement custom iterators to do that. Maybe one day…

template <class F> using Polygon = vector<Point<F>>;
inline int prev(int i, int n) { return i == 0 ? n-1 : i-1; }
inline int next(int i, int n) { return i == n-1 ? 0 : i+1; }
template <class T> inline int sgn(const T& x) { return (T(0) < x) - (x < T(0)); }

1. Area: O(N)

Probably, everyone knows how to compute polygon area as a sum of triangle or trapezoid areas. If not, you can read [geomalgorithms] or [O'Rourke, p.16] for more detailed explanation. However, as it's shown in area2D_Polygon on geomalgorithms, this computation can be optimized, resulting in the following code:

template <class F>
F area(const Polygon<F>& poly) {
  int n = static_cast<int>(poly.size());
  F area = F(0);
  for (int i = 0; i < n; ++i)
    area += poly[i].x * (poly[next(i, n)].y - poly[prev(i, n)].y);
  return area;
}

Note that this function returns doubled oriented area. Doubled, because for any polygon with integer coordinated the area may be not integer but half-integer. Oriented means that it's positive if polygon vertices are listed in counter-clockwise (ccw) order and negative otherwise.

2. Orientation

Polygon vertices are given in either counter-clockwise or clockwise order. All the algorithms below assume that orientation is counter-clockwise. Because of that, we need a way to verify this condition. Common solution is to check the sign of area return value. However, it's not necessary to compute the area.

2.1 Simple polygon: O(N)

As noted in [O'Rourke, p.12], it's enough to check orientation at one of extreme polygon vertices, for example lower-left. Complexity remains the same, but operations are faster: comparisons instead of multiplication and summation.

// True if orientation of a simple polygon is counter-clockwise.
template <class F> 
bool orientation(const Polygon<F>& poly) {
  int n = static_cast<int>(poly.size());
  int i = static_cast<int>(min_element(begin(poly), end(poly)) - begin(poly));
  return ccw(poly[prev(i, n)], poly[next(i, n)], poly[i]) > 0;
}

2.2 Convex polygon: O(1)

All vertices of the convex polygon are extreme, so we can check orientation at an arbitrary vertex.

3. Convex hull: O(NlogN)

We will be using convex polygons, so first let's learn to build one from a set of points. There are various algorithms with the same complexity, so I chose the one with the smallest constant: Andrew's monotone chain algorithm [geomalgorithms]. It's faster, because points are sorted lexicographically, not using more complex angle comparison as in sortByAngle in my previous article. The implementation was taken from here. The difference is that I build right and left chains, not upper and lower. Also, I handle the degenerate case when all the points coincide. In such case the original code outputs 2 points instead of one.

template <class F>
Polygon<F> convexHull(Polygon<F> points) {
  sort(begin(points), end(points));
  Polygon<F> hull;
  hull.reserve(points.size() + 1);
  for (int phase = 0; phase < 2; ++phase) {
    auto start = hull.size();
    for (auto& point : points) {
      while (hull.size() >= start+2 &&
             ccw(point, hull.back(), hull[hull.size()-2]) <= 0)
        hull.pop_back();
      hull.push_back(point);
    }
    hull.pop_back();
    reverse(begin(points), end(points));
  }
  if (hull.size() == 2 && hull[0] == hull[1]) hull.pop_back();
  return hull;
}

4. Inclusion tests

Given a point, sometimes we need to check whether it lies inside or outside the given polygon. All the functions below return

• -1 if the point is strictly inside the polygon;
• 0 if the point is on the boundary;
• 1 if the point is outside the polygon.

4.1 Triangle: O(1)

The easiest example of a polygon is triangle. It's an important object in computer graphics and deserves special treating :) All the cases are visualized in [O'Rourke, p.235], from where I derived the following conditions:

• point is inside if all the signs are equal;
• point is outside if there are two opposite signs;
• otherwise point is on the boundary (two signs are equal and the third one is 0).

template <class F1, class F2>
int pointVsTriangle(const Point<F1>& point, const Polygon<F2>& triangle) {
  assert(triangle.size() == 3);
  int signs[3];
  for (int i = 0; i < 3; ++i)
    signs[i] = sgn(ccw(point, triangle[next(i, 3)], triangle[i]));
  if (signs[0] == signs[1] && signs[1] == signs[2]) return -1;
  for (int i = 0; i < 3; ++i) if (signs[i] * signs[next(i, 3)] == -1) return 1;
  return 0;
}

4.2 Convex polygon: O(logN)

The popular approach here is to triangulate the polygon by drawing diagonals from one vertex to all the others, find the angle where the given point lies using binary search, and then check if it's inside the triangle or not. However, I got the best time on the problem 166B - Polygons using the same approach as in convex hull construction (and I'm not the only one using fast I/O there). Also, I have less special cases to check (even if we don't need to distinguish the border from the interior).

If the point is below (or on the same horizontal line but to the right of) the lower-left vertex, it's definitely outside, similarly for upper-right. Now, let's connect the lower-left and the upper-right polygon vertices. If the given point lies on this line, it's either on the polygon boundary (if it coincides with one of the endpoints or this line is a polygon edge) or inside. If the point is to the right (analogously to the left), we need to check if it's to the left of the right chain built by the convex hull algorithm. The corresponding edge is found using binary search, again comparing points lexicographically.

Current implementation assumes that lower-left index has index 0, and we know the index of upper-right vertex (called top, it can be found once in linear time). This assumption is perfectly reasonable, because if fits the output of convex hull algorithm. Picture shows points locations that will be compared with each polygon edge. For example, point P with be compared with edge 2-3, and it's inside the polygon because it's to the left of this edge.

template <class F1, class F2>
int pointVsConvexPolygon(const Point<F1>& point, const Polygon<F2>& poly, int top) {
  if (point < poly[0] || point > poly[top]) return 1;
  auto orientation = ccw(point, poly[top], poly[0]);
  if (orientation == 0) {
    if (point == poly[0] || point == poly[top]) return 0;
    return top == 1 || top + 1 == poly.size() ? 0 : -1;
  } else if (orientation < 0) {
    auto itRight = lower_bound(begin(poly) + 1, begin(poly) + top, point);
    return sgn(ccw(itRight[0], point, itRight[-1]));
  } else {
    auto itLeft = upper_bound(poly.rbegin(), poly.rend() - top-1, point);
    return sgn(ccw(itLeft == poly.rbegin() ? poly[0] : itLeft[-1], point, itLeft[0]));
  }
}

Exercise: verify that this code works correctly when horizontal edges are present.

4.3 Simple polygon: O(N)

We will use the winding number described for example on [geomalgorithms], the same number is the answer for the problem Timus 1599. It's defined as a number of full counter-clockwise turns around the given point we make if we follow the polygon edges. In case of a simple polygon it's either 1 when the point is inside or 0 if outside.

Winding number can be computed in the following way. We draw an arbitrary ray from the given point and each time the polygon edge crosses it we add either +1 if the crossing was from the right to the left side and -1 otherwise. Such edges are marked with arrows on the picture. For optimization purposes I chose the horizontal ray going in increasing x direction. The ray itself is included in the upper half-plane. For a horizontal edge, I check whether it contains the given point. If the edge endpoints are on the different sides of the line containing the ray, I check if the edge crosses this line to the right side of the given point.

template <class F1, class F2>
int pointVsPolygon(const Point<F1>& point, const Polygon<F2>& poly) {
  int n = static_cast<int>(poly.size()), windingNumber = 0;
  for (int i = 0; i < n; ++i) {
    if (point == poly[i]) return 0;
    int j = next(i, n);
    if (poly[i].y == point.y && poly[j].y == point.y) {
      if (min(poly[i].x, poly[j].x) <= point.x &&
          point.x <= max(poly[i].x, poly[j].x)) return 0;
    } else {
      bool below = poly[i].y < point.y;
      if (below != (poly[j].y < point.y)) {
        auto orientation = ccw(point, poly[j], poly[i]);
        if (orientation == 0) return 0;
        if (below == (orientation > 0)) windingNumber += below ? 1 : -1;
      }
    }
  }
  return windingNumber == 0 ? 1 : -1;
}

Exercise: think of an example when the point == poly[i] check is necessary.

5. Extreme points

I define an extreme vertex with respect to the given direction as following: if we draw a line through this vertex in the given direction, there are no polygon vertices to the right of it. This direction can be fixed for all vertices, or it can be for example direction from a given point to the vertex. In the latter case extreme vertex is the right tangent vertex (for example vertices 1 and 2 for point P in the picture).

5.1 Convex polygon: O(logN)

The binary search algorithm was described in [O'Rourke, p.270] and implemented in functions [polyMax_2D] and [tangent_PointPolyC]. My contribution here is writing one general function instead of three similar ones. I packed all the cases handling into one if statement. At least for me, it's easier to understand and remember. Also, I believe the cases when direction coincides with polygon side (as in the picture) aren't handled correctly. These cases are (in O'Rourke notation):

if C points neither up nor down
	if A points up
		then [a, b] <- [a, c]
		else [a, b] <- [c, b]
if A points neither up nor down
	if C points up
		then [a, b] <- [c, b]
		else [a, b] <- [a, c]

There are at most 2 vertices where this can happen (1 and 5 in the picture) and one of them is always extreme (number 1). This implies A and C can't both point neither up nor down at the same time.

Now I'll explain my implementation a bit. vertexCmp returns -1 if vertex i is to the right of vertex j with respect to the given direction. isExtreme just follows my definition, but in case when there are two extreme vertices on the same line, it returns true only for the first in counter-clockwise order (true for 1 and 5, false for 2 and 0 in the picture). Once extreme vertex is encountered, it's immediately returned. For the algorithm to work correctly, the binary search isn't executed in case vertex 0 is extreme.

template <class F, class Function>
int extremeVertex(const Polygon<F>& poly, Function direction) {
  int n = static_cast<int>(poly.size()), left = 0, leftSgn;
  auto vertexCmp = [&poly, direction](int i, int j) {
    return sgn(ccw(direction(poly[j]), poly[j] - poly[i])); };
  auto isExtreme = [n, vertexCmp](int i, int& iSgn) {
    return (iSgn = vertexCmp(next(i, n), i)) >= 0 && vertexCmp(i, prev(i, n)) < 0; };
  for (int right = isExtreme(0, leftSgn) ? 1 : n; left + 1 < right;) {
    int middle = (left + right) / 2, middleSgn;
    if (isExtreme(middle, middleSgn)) return middle;
    if (leftSgn != middleSgn ? leftSgn < middleSgn
        : leftSgn == vertexCmp(left, middle)) right = middle;
    else left = middle, leftSgn = middleSgn;
  }
  return left;
}

Probably, you have a question regarding direction function. Below is the code showing how this function looks like to find the right and left tangents from the given point. Fixed direction will be demonstrated in Section 6.1.

template <class F1, class F2>
pair<int, int> tangentsConvex(const Point<F1>& point, const Polygon<F2>& poly) {
  return {
    extremeVertex(poly, [&point](const Point<F2>& q) { return q - point; }),
    extremeVertex(poly, [&point](const Point<F2>& q) { return point - q; })};
}

5.2 Simple polygon

Extreme points always lie on the convex hull. If you're processing multiple queries for such points, it's worthwhile to build the convex hull for the given points beforehand and apply the algorithm from Section 5.1. Convex hull can be built using algorithm from Section 3, but there's more efficient solution for simple polygons I have yet to analyze.

6. Polygon stabbing

The problem of finding the intersection of a geometric object with a line is often called the "stabbing" problem.

6.1 Convex polygon: O(logN)

OK, here's the place where I got lazy. I just implemented an idea expressed in [O'Rourke, p.271] and didn't test it well. The result is either empty or a segment. First, I find two extreme points with respect to the line direction, and then I find intersection segment endpoints on two chains between them using binary search. Care must be taken, because one of these chains passes through 0 vertex. After the correct edge is found, I call the intersectLines function from my previous article. I believe extremeVertex can be modified to find each intersection point in one binary search, but I'll leave this idea for the future.

template <class F1, class F2, class F3>
bool stabConvexPolygon(const Line<F1>& line, const Polygon<F2>& poly, Line<F3>& res) {
  assert(line);
  int right = extremeVertex(poly, [&line](const Point<F2>&) { return line.ab; });
  int left = extremeVertex(poly, [&line](const Point<F2>&) { return -line.ab; });
  auto vertexCmp = [&line](const Point<F2>& vertex) {
    return sgn(ccw(line.ab, vertex - line.a)); };
  int rightSgn = vertexCmp(poly[right]), leftSgn = vertexCmp(poly[left]);
  if (rightSgn < 0 || leftSgn > 0) return false;
  auto intersectChain = [&line, &poly, vertexCmp](int first, int last,
                                                  int firstSgn, Point<F3>& res) {
    int n = static_cast<int>(poly.size());
    while (next(first, n) != last) {
      int middle = (first + last + (first < last ? 0 : n)) / 2;
      if (middle >= n) middle -= n;
      if (vertexCmp(poly[middle]) == firstSgn) first = middle;
      else last = middle;
    }
    intersectLines<0, 0, 0, 0>(line, makeLine(poly[first], poly[last]), res);
  };
  intersectChain(left, right, leftSgn, res.a);
  intersectChain(right, left, rightSgn, res.ab);
  res.ab -= res.a;
  return true;
}

6.2 Simple polygon: O(N)

The problem of finding the intersection length was given at the first educational round 598F - Cut Length. Surprisingly, I didn't notice a single comment mentioning that it can be solved in linear time (per line). Recall that the line is represented in parametric form , then the intersection is the set of segments for the value of t. It's enough to find the total length of these segments and multiply it by |AB|. The solution is based on the simple fact: when traversing the polygon in counter-clockwise direction, every time we cross the line from right to left we get an end of some t segment, so this value is added to the answer, and from left to right we get a beginning, so it's subtracted from the answer. We only sort this values if we're asked to produce all the intersection segments, not just their total length.

The problem isn't easy because some vertices and edges may lie on the line, for a total of 6 cases illustrated in the picture. Points we're interested at are marked blue. The general idea to process these cases is the following. If we're crossing from right to left, we need the furthest point in the line direction. Suppose we're processing a vertex that lies on the line, should we include it in the answer? It's not the furthest only in case one of its adjacent edges lies on the line, and the second vertex is further. This condition is checked using dot product. To determine the direction in which we're moving, it's enough to know on which side of the line previous and next vertices lie (variables prevSgn and nextSgn take values 1 on the left, 0 on the line, -1 on the right). For example, if prevSgn == 0 && nextSgn == 1, we're moving from right to left.

I'm not going to go through all the cases, but I'll explain how the values of t are computed. If an edge crosses the line not at a vertex, the formula from intersectLines applies. If vertex V lies on the line, we have , and taking dot product of each side with results in . This is the only place in this article where I needed non-integer computations.

template <class F1, class F2, class F = distF<F1, F2>>
F stabPolygonLength(const Line<F1>& line, const Polygon<F2>& poly) {
  assert(line);
  F tSum = F(0);
  int n = static_cast<int>(poly.size());
  auto vertexSgn = [&line, &poly](int i) { return sgn(line.ab ^ (poly[i] - line.a)); };
  int prevSgn = vertexSgn(n - 1), iSgn = vertexSgn(0), nextSgn;
  for (int i = 0; i < n; ++i, prevSgn = iSgn, iSgn = nextSgn) {
    nextSgn = vertexSgn(next(i, n));
    if (iSgn == 0) {
      if (prevSgn == 0) {
        if (nextSgn != 0 && nextSgn == sgn((poly[i] - poly[prev(i, n)]) * line.ab))
          tSum += nextSgn * static_cast<F>((poly[i] - line.a) * line.ab) / abs(line.ab);
      } else {
        if ((nextSgn != 0 && nextSgn != prevSgn) ||
            (nextSgn == 0 && prevSgn == sgn((poly[next(i, n)] - poly[i]) * line.ab)))
          tSum -= prevSgn * static_cast<F>((poly[i] - line.a) * line.ab) / abs(line.ab);
      }
    } else if (nextSgn == -iSgn) {
      auto vect = poly[next(i, n)] - poly[i];
      tSum += nextSgn * static_cast<F>((poly[i] - line.a) ^ vect) / (line.ab ^ vect);
    }
  }
  return tSum * norm(line.ab);
}

7. Maximum distance (convex polygon): O(N)

As a final chord, I decided to include my implementation for the maximum distance (squared) between two vertices of the convex polygon. It uses the well known rotating calipers technique. If the input is not a convex polygon, call convexHull from Section 3 first.

template <class F>
F maxDist2(const Polygon<F>& poly) {
  int n = static_cast<int>(poly.size());
  F res = F(0);
  for (int i = 0, j = n < 2 ? 0 : 1; i < j; ++i)
    for (;; j = next(j, n)) {
      res = max(res, dist2(poly[i], poly[j]));
      if (ccw(poly[i+1] - poly[i], poly[next(j, n)] - poly[j]) >= 0) break;
    }
  return res;
}