7.3.5 Minimum Enclosing Circle

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

/*

Given a set of points in two dimensions, finds the unique circle of minimum radius (equivalently,
minimum area) containing all points.

This implementation uses Welzl's randomized incremental algorithm. The points are shuffled and
processed one at a time while maintaining the current minimum enclosing circle. Whenever a point
lies outside the current circle, the circle is rebuilt with that point constrained to lie on the
boundary. The final circle is determined by at most three boundary points.

- `minimum_enclosing_circle(lo, hi)` returns the minimum enclosing circle for the range
  $[`lo`, `hi`)$ of points, where `lo` and `hi` must be random-access iterators. The input range is
  shuffled in place. The point type may be any type exposing numeric `.x` and `.y` members. The
  returned circle always uses `double` coordinates, so integer-coordinate inputs are accepted.

Time Complexity:
- O(n) expected time per call, where $n$ is the distance between `lo` and `hi`, or O(n^3) in the
  worst-case for a particular shuffle order.

Space Complexity:
- O(1) auxiliary.

*/

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

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

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

template<typename T, typename U>
bool LE(T a, U b) {
  return !LT(b, a);
}

double sqnorm(double x, double y) {
  return x * x + y * y;
}

// SFINAE guard: valid only when Pt exposes numeric .x/.y members. This keeps the templated point
// constructors from hijacking calls like Circle(double, double, double).
template<typename Pt>
using if_point = decltype(std::declval<const Pt &>().x, std::declval<const Pt &>().y, void());

struct Circle {
  double h, k, r;

  Circle() : h(0), k(0), r(0) {}
  Circle(double h, double k, double r) : h(h), k(k), r(fabs(r)) {}

  template<typename Pt, typename = if_point<Pt>>
  Circle(const Pt &a, const Pt &b) {
    h = (a.x + b.x) / 2.0;
    k = (a.y + b.y) / 2.0;
    r = std::hypot(a.x - h, a.y - k);
  }

  template<typename Pt, typename = if_point<Pt>>
  Circle(const Pt &a, const Pt &b, const Pt &c) {
    double an = sqnorm(b.x - c.x, b.y - c.y);
    double bn = sqnorm(a.x - c.x, a.y - c.y);
    double cn = sqnorm(a.x - b.x, a.y - b.y);
    double wa = an * (bn + cn - an), wb = bn * (an + cn - bn), wc = cn * (an + bn - cn);
    double w = wa + wb + wc;
    if (EQ(w, 0)) {
      throw std::runtime_error("No circumcircle from collinear points.");
    }
    h = (wa * a.x + wb * b.x + wc * c.x) / w;
    k = (wa * a.y + wb * b.y + wc * c.y) / w;
    r = std::hypot(a.x - h, a.y - k);
  }

  template<typename Pt, typename = if_point<Pt>>
  bool contains(const Pt &p) const {
    return LE(sqnorm(p.x - h, p.y - k), r * r);
  }
};

// Input points can be any type with .x and .y fields; Circle output is always double.
template<typename It>
Circle minimum_enclosing_circle(It lo, It hi) {
  if (lo == hi) {
    return Circle(0, 0, 0);
  }
  if (lo + 1 == hi) {
    return Circle(lo->x, lo->y, 0);
  }
  static std::mt19937 rng(std::random_device{}());
  std::shuffle(lo, hi, rng);
  Circle res(*lo, *(lo + 1));
  for (It i = lo + 2; i != hi; ++i) {
    if (res.contains(*i)) {
      continue;
    }
    res = Circle(*lo, *i);
    for (It j = lo + 1; j != i; ++j) {
      if (res.contains(*j)) {
        continue;
      }
      res = Circle(*i, *j);
      for (It k = lo; k != j; ++k) {
        if (!res.contains(*k)) {
          res = Circle(*i, *j, *k);
        }
      }
    }
  }
  return res;
}

/*** 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; }
};

struct PointI {
  int x, y;
  PointI(int x = 0, int y = 0) : x(x), y(y) {}
};

int main() {
  vector<Point> empty;
  Circle none = minimum_enclosing_circle(empty.begin(), empty.end());
  assert(EQ(none.h, 0) && EQ(none.k, 0) && EQ(none.r, 0));

  vector<Point> one{{2, 3}};
  Circle single = minimum_enclosing_circle(one.begin(), one.end());
  assert(EQ(single.h, 2) && EQ(single.k, 3) && EQ(single.r, 0));

  vector<Point> two{{0, 0}, {4, 0}};
  Circle diameter = minimum_enclosing_circle(two.begin(), two.end());
  assert(EQ(diameter.h, 2) && EQ(diameter.k, 0) && EQ(diameter.r, 2));

  vector<Point> v{{0, 0}, {0, 1}, {1, 0}, {1, 1}};
  Circle res = minimum_enclosing_circle(v.begin(), v.end());
  assert(EQ(res.h, 0.5) && EQ(res.k, 0.5) && EQ(res.r, 1 / sqrt(2)));

  // Integer-coordinate input: Circle output is always double.
  vector<PointI> iv{{0, 0}, {0, 2}, {2, 0}, {2, 2}};
  Circle ic = minimum_enclosing_circle(iv.begin(), iv.end());
  assert(EQ(ic.h, 1.0) && EQ(ic.k, 1.0));
  return 0;
}