7.3.1 Polygon Sorting and Area

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

/*

Given a list of distinct points in two dimensions, order them into a valid polygon and determine the
area. The sorting comparator orders two points by the sign of their cross product around a chosen
center, that is, by angle. To form a simple polygon, choose a center inside the points' convex hull
that differs from every input point; the arithmetic mean used in the example is usually convenient.
The area functions use the shoelace formula, summing the cross products of consecutive vertex pairs.
All functions accept `Point`, `PointI`, or `PointL` from 7.1.1, or any compatible point type with
numeric `.x` and `.y` fields and the comparison operators used by the requested operation.

- `cw_comp(a, b, c)` returns whether point `a` compares clockwise before `b` about `c`.
- `polygon_area_2x(lo, hi)` returns exactly double the area of the polygon with vertices specified
  by the range $[`lo`, `hi`)$ of points in either clockwise or counter-clockwise order. The return
  value is integral or floating-point, depending on the input point type. For integer vertices,
  divide by $2$ in the caller if the exact area is needed.
- `polygon_area(lo, hi)` returns the area as `double`.
- `polygon_centroid(lo, hi)` returns the centroid (center of mass) of a non-degenerate simple
  polygon with vertices in boundary order as an `std::pair<double, double>`. It uses signed shoelace
  weights, so either clockwise or counter-clockwise inputs will work.

Overflow warning: `cw_comp` and `polygon_area_2x` form cross products that grow like the squared
coordinate magnitude (and the shoelace sum accumulates over all vertices). For integer point types
use a 64-bit coordinate type (e.g. `PointL` from 7.1.1) for large or numerous coordinates.

Time Complexity:
- O(n) per call to `polygon_area_2x()`, `polygon_area()`, and `polygon_centroid()`, where $n$ is the
  number of points.
- O(1) per call to `cw_comp()` and the comparators.

Space Complexity:
- O(1) auxiliary for all operations.

*/

#include <algorithm>
#include <cassert>
#include <cmath>
#include <random>
#include <stdexcept>
#include <type_traits>
#include <utility>

const double EPS = 1e-9;

template<typename T, typename U, typename C = std::common_type_t<T, U>>
bool EQ(T a, U b) {
  if constexpr (std::is_floating_point_v<C>) return C(a) == C(b) || std::fabs(C(a) - C(b)) <= EPS;
  return C(a) == C(b);
}

// Comparator (clockwise angular order about c). Exact: comparisons are pure sign tests on
// coordinate differences and the cross product, so it is a valid strict weak ordering and is
// exact for integer-coordinate points.
template<typename Pt>
bool cw_comp(const Pt &a, const Pt &b, const Pt &c) {
  if (a.x - c.x >= 0 && b.x - c.x < 0) {
    return true;
  }
  if (a.x - c.x < 0 && b.x - c.x >= 0) {
    return false;
  }
  if (a.x - c.x == 0 && b.x - c.x == 0) {
    if (a.y - c.y >= 0 || b.y - c.y >= 0) {
      return a.y > b.y;
    }
    return b.y > a.y;
  }
  auto acx = a.x - c.x, acy = a.y - c.y;
  auto bcx = b.x - c.x, bcy = b.y - c.y;
  auto det = acx * bcy - acy * bcx;  // Overflow warning.
  if (det == 0) {
    auto acnorm = acx * acx + acy * acy;  // Overflow warning.
    auto bcnorm = bcx * bcx + bcy * bcy;
    return acnorm > bcnorm;
  }
  return det < 0;
}

// Returns 2 * area. Result is exact (integer) for integer-coordinate points.
template<typename It>
auto polygon_area_2x(It lo, It hi) {
  using T = decltype(lo->x * lo->y);
  if (lo == hi) {
    return T{0};
  }
  T area = 0;
  for (It i = lo, j = hi - 1; i != hi; j = i++) {
    area += static_cast<T>(j->x - i->x) * static_cast<T>(j->y + i->y);  // Overflow warning.
  }
  return area < T{0} ? -area : area;
}

template<typename It>
double polygon_area(It lo, It hi) {
  return static_cast<double>(polygon_area_2x(lo, hi)) / 2.0;
}

template<typename It>
std::pair<double, double> polygon_centroid(It lo, It hi) {
  assert(hi - lo >= 3);
  double cx = 0, cy = 0, area2 = 0;
  for (It i = lo, j = hi - 1; i != hi; j = i++) {
    double cross = static_cast<double>(j->x) * i->y - static_cast<double>(i->x) * j->y;
    cx += (static_cast<double>(j->x) + i->x) * cross;
    cy += (static_cast<double>(j->y) + i->y) * cross;
    area2 += cross;
  }
  if (EQ(area2, 0)) {
    throw std::runtime_error("Cannot compute centroid of zero-area polygon.");
  }
  return {cx / (3 * area2), cy / (3 * area2)};
}

/*** 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) {}
  bool operator==(const Point &p) const { return x == p.x && y == p.y; }
  friend bool EQ(const Point &a, const Point &b) { return EQ(a.x, b.x) && EQ(a.y, b.y); }
  bool operator!=(const Point &p) const { return !(*this == p); }
  bool operator<(const Point &p) const { return x != p.x ? x < p.x : y < p.y; }
  bool operator>(const Point &p) const { return p < *this; }
};

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 !(*this == p); }
  bool operator<(const PointI &p) const { return x != p.x ? x < p.x : y < p.y; }
  bool operator>(const PointI &p) const { return p < *this; }
};

template<typename It>
Point mean_center(It lo, It hi) {
  assert(lo != hi);
  double x_sum = 0, y_sum = 0, n = hi - lo;
  for (It it = lo; it != hi; ++it) {
    x_sum += it->x;
    y_sum += it->y;
  }
  return Point(x_sum / n, y_sum / n);
}

int main() {
  vector<Point> points{Point(1, 3), Point(1, 2), Point(2, 1), Point(0, 0), Point(-1, 3)};
  vector<Point> v(points);
  mt19937 rng(1234567);  // Fixed seed for reproducibility.
  shuffle(v.begin(), v.end(), rng);
  Point c = mean_center(v.begin(), v.end());
  assert(EQ(c, Point(0.6, 1.8)));
  sort(v.begin(), v.end(), [c](const Point &a, const Point &b) { return cw_comp(a, b, c); });
  assert((v == vector<Point>{{1, 3}, {1, 2}, {2, 1}, {0, 0}, {-1, 3}}));
  assert(EQ(polygon_area(v.begin(), v.end()), 5));

  sort(v.begin(), v.end(), [c](const Point &a, const Point &b) { return cw_comp(b, a, c); });
  assert((v == vector<Point>{{-1, 3}, {0, 0}, {2, 1}, {1, 2}, {1, 3}}));
  assert(EQ(polygon_area(v.begin(), v.end()), 5));

  // Integer points: polygon_area_2x is exact (no float arithmetic).
  vector<PointI> iv{{0, 0}, {4, 0}, {0, 3}};            // right triangle, area = 6
  assert(polygon_area_2x(iv.begin(), iv.end()) == 12);  // exact int
  assert(EQ(polygon_area(iv.begin(), iv.end()), 6.0));
  auto centroid = polygon_centroid(iv.begin(), iv.end());
  assert(EQ(centroid.first, 4.0 / 3.0) && EQ(centroid.second, 1.0));
  return 0;
}