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
156
157
158
159
160
161
162
163
164
165

/*

Given a list of points in two dimensions, computes the convex hull using the monotone chain
algorithm and the diametral pair using rotating calipers. Monotone chain sorts the points
lexicographically and builds the lower and upper hulls 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, visiting every candidate diametral
pair in linear time. The functions are templated on the iterator type; the value type of the
iterator is used as the point type. Both functions accept either floating-point or integral
coordinates, and use exact comparisons for integral points.

- `convex_hull(lo, hi)` returns the convex hull in counter-clockwise order for a range
  $[`lo`, `hi`)$ of points. 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(lo, hi)` returns the maximum-distance pair of points.

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 distance between `lo` and `hi`.

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

*/

#include <algorithm>
#include <cmath>
#include <iterator>
#include <random>
#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 It>
auto convex_hull(It lo, It hi) {
  using Pt = typename std::iterator_traits<It>::value_type;
  std::vector<Pt> p(lo, hi);
  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 It>
auto diametral_pair(It lo, It hi) {
  using Pt = typename std::iterator_traits<It>::value_type;
  auto h = convex_hull(lo, hi);
  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 <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.begin(), v.end()) == h);
  }
  {
    vector<PointI> v{{0, 0}, {3, 0}, {0, 3}, {1, 1}, {4, 4}};
    auto [p1, p2] = diametral_pair(v.begin(), v.end());
    assert(p1 == PointI(0, 0));
    assert(p2 == PointI(4, 4));
  }
  {
    vector<PointI> v{{0, 0}, {4, 0}, {4, 4}, {0, 4}, {2, 2}};
    auto h = convex_hull(v.begin(), v.end());
    assert(h.size() == 4);  // interior point (2,2) excluded
    auto [p1, p2] = diametral_pair(v.begin(), v.end());
    // 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))
    );
  }
  {
    vector<PointI> v{{2, 3}, {2, 3}, {2, 3}};
    vector<PointI> expected{{2, 3}};
    assert(convex_hull(v.begin(), v.end()) == expected);
  }
  return 0;
}