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
173
174
175
176
177
178
179
180
181
182
183

/*

Given a list of distinct points in two-dimensions, order them into a valid polygon and determine the
area. The sorting comparators order two points by the sign of their cross product around a chosen
center, that is, by angle. The area functions use the shoelace formula, summing the cross products
of consecutive vertex pairs. The functions are templated on the point type. The local `Point` struct
(`double` coordinates) is the default; replace it with `Point`/`PointD`/ `PointI` from 7.1.1 or any
struct with numeric `.x` and `.y` fields and `<` / `==` operators.

- `cw_comp(a, b, c)` returns whether point `a` compares clockwise before `b` about `c`.
- `CWComparator(c)` / `CCWComparator(c)` return comparators for `std::sort`.
- `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`, where $n$ is the number of points.
- O(n) per call to `polygon_centroid`.
- O(1) per call to `cw_comp` and the comparators.

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

*/

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

const double EPS = 1e-9;

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

// 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<class 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;
  if (det == 0) {
    auto acnorm = acx * acx + acy * acy;
    auto bcnorm = bcx * bcx + bcy * bcy;
    return acnorm > bcnorm;
  }
  return det < 0;
}

// Returns 2 * area. Result is exact (integer) for integer-coordinate points.
template<class 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 += (T)(j->x - i->x) * (T)(j->y + i->y);
  }
  return area < T(0) ? -area : area;
}

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

template<class It>
std::pair<double, double> polygon_centroid(It lo, It hi) {
  if (hi - lo < 3) {
    throw std::runtime_error("Cannot compute centroid of a degenerate polygon.");
  }
  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 += (j->x + i->x) * cross;
    cy += (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 and Output:

Clockwise vertices: (1, 3) (1, 2) (2, 1) (0, 0) (-1, 3) 
Counterclockwise vertices: (-1, 3) (0, 0) (2, 1) (1, 2) (1, 3) 

***/

#include <cassert>
#include <iostream>
#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 EQ(x, p.x) && EQ(y, p.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<class It>
Point mean_center(It lo, It hi) {
  if (lo == hi) {
    throw std::runtime_error("Cannot get center of an empty range.");
  }
  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);
  std::mt19937 rng(1234567);
  std::shuffle(v.begin(), v.end(), rng);
  Point c = mean_center(v.begin(), v.end());
  assert(EQ(c.x, 0.6) && EQ(c.y, 1.8));
  sort(v.begin(), v.end(), [c](const Point &a, const Point &b) { return cw_comp(a, b, c); });
  cout << "Clockwise vertices: ";
  for (const auto &p : v) {
    cout << "(" << p.x << ", " << p.y << ") ";
  }
  cout << endl;
  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); });
  cout << "Counterclockwise vertices: ";
  for (const auto &p : v) {
    cout << "(" << p.x << ", " << p.y << ") ";
  }
  cout << endl;
  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;
}