7.3.3 Convex Hull and Diametral Pair

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

/*

The convex hull of points in two dimensions is the smallest convex set containing them. A diametral
pair is a pair of input points at maximum distance. Monotone chain computes the hull by sorting the
points lexicographically and building the lower and upper boundaries in one pass each, popping any
point that would fail to create a counter-clockwise turn. Rotating calipers then walks two antipodal
pointers around the hull, advancing whichever increases the separation and visiting every candidate
diametral pair in linear time. Both functions accept either floating-point or integral coordinates,
and use exact comparisons for integral points.

- `convex_hull(p)` returns the convex hull of points `p` in counter-clockwise order. Duplicate input
  points and collinear points along hull edges are omitted. To instead return the hull points in
  clockwise order, replace the cross product comparisons `<= 0` with `>= 0`.
- `diametral_pair(p)` returns the maximum-distance pair of points in `p`.

Overflow warning: `cross()` and the squared distance in `diametral_pair()` grow like the squared
coordinate magnitude. With 32-bit `int` coordinates they overflow once coordinates exceed a few tens
of thousands; use a 64-bit (`int64_t`) coordinate type for larger integer inputs.

Time Complexity:
- O(n log n) per call, where $n$ is the number of points.

Space Complexity:
- O(n) auxiliary for storage of the convex hull.

*/

#include <algorithm>
#include <cmath>
#include <utility>
#include <vector>

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

// Convex hull: exact for integer-coordinate points.
template<typename Pt>
std::vector<Pt> convex_hull(std::vector<Pt> p) {
  std::sort(p.begin(), p.end());
  auto same_point = [](const Pt &a, const Pt &b) { return !(a < b) && !(b < a); };
  p.erase(std::unique(p.begin(), p.end(), same_point), p.end());
  int n = static_cast<int>(p.size());
  if (n <= 1) {
    return p;
  }
  int k = 0;
  std::vector<Pt> res(2 * n);
  for (const Pt &q : p) {
    while (k >= 2 && cross(res[k - 1], q, res[k - 2]) <= 0) {
      k--;
    }
    res[k++] = q;
  }
  int t = k + 1;
  for (int i = n - 2; i >= 0; i--) {
    while (k >= t && cross(res[k - 1], p[i], res[k - 2]) <= 0) {
      k--;
    }
    res[k++] = p[i];
  }
  res.resize(k - 1);
  return res;
}

// Diametral pair: squared-distance comparisons are exact for integer-coordinate points.
template<typename Pt>
std::pair<Pt, Pt> diametral_pair(const std::vector<Pt> &p) {
  auto h = convex_hull(p);
  int m = static_cast<int>(h.size());
  if (m == 0) {
    return std::pair<Pt, Pt>{};
  }
  if (m == 1) {
    return std::pair<Pt, Pt>{h[0], h[0]};
  }
  if (m == 2) {
    return std::pair<Pt, Pt>{h[0], h[1]};
  }
  int k = 1;
  while (std::abs(cross(h[0], h[(k + 1) % m], h[m - 1])) > std::abs(cross(h[0], h[k], h[m - 1]))) {
    k++;
  }
  auto sqdist = [](const Pt &a, const Pt &b) {
    // Overflow risk for integer Pt: ~O(max_coord^2); use int64_t if necessary.
    auto dx = a.x - b.x, dy = a.y - b.y;
    return dx * dx + dy * dy;
  };
  auto maxsq = sqdist(h[0], h[0]);
  std::pair<Pt, Pt> res{h[0], h[0]};
  for (int i = 0, j = k; i <= k && j < m; i++) {
    auto d = sqdist(h[i], h[j]);
    if (d > maxsq) {
      maxsq = d;
      res = {h[i], h[j]};
    }
    while (j < m && std::abs(cross(h[(i + 1) % m], h[(j + 1) % m], h[i])) >
                        std::abs(cross(h[(i + 1) % m], h[j], h[i]))) {
      d = sqdist(h[i], h[(j + 1) % m]);
      if (d > maxsq) {
        maxsq = d;
        res = {h[i], h[(j + 1) % m]};
      }
      j++;
    }
  }
  return res;
}

/*** Example Usage ***/

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

struct PointI {
  int x, y;
  PointI(int x = 0, int y = 0) : x(x), y(y) {}
  bool operator==(const PointI &o) const { return x == o.x && y == o.y; }
  bool operator!=(const PointI &o) const { return !(*this == o); }
  bool operator<(const PointI &o) const { return x != o.x ? x < o.x : y < o.y; }
  bool operator>(const PointI &o) const { return o < *this; }
};

int main() {
  {
    vector<PointI> v{{1, 3}, {1, 2}, {2, 1}, {0, 0}, {-1, 3}};
    mt19937 rng(1234567);  // Fixed seed for reproducibility.
    shuffle(v.begin(), v.end(), rng);
    vector<PointI> h{{-1, 3}, {0, 0}, {2, 1}, {1, 3}};
    assert(convex_hull(v) == h);
  }
  {
    auto [p1, p2] = diametral_pair(vector<PointI>{{0, 0}, {3, 0}, {0, 3}, {1, 1}, {4, 4}});
    assert(p1 == PointI(0, 0) && p2 == PointI(4, 4));
  }
  {
    vector<PointI> v{{0, 0}, {4, 0}, {4, 4}, {0, 4}, {2, 2}};
    auto h = convex_hull(v);
    assert(h.size() == 4);  // interior point (2,2) excluded
    auto [p1, p2] = diametral_pair(v);
    // diametral pair is exact
    assert(
        (p1 == PointI(0, 0) && p2 == PointI(4, 4)) || (p1 == PointI(4, 4) && p2 == PointI(0, 0)) ||
        (p1 == PointI(4, 0) && p2 == PointI(0, 4)) || (p1 == PointI(0, 4) && p2 == PointI(4, 0))
    );
  }
  {
    assert((convex_hull(vector<PointI>{{2, 3}, {2, 3}, {2, 3}}) == vector<PointI>{{2, 3}}));
  }
  return 0;
}