7.4.2 Half-Plane Intersection

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 directed lines, compute the convex polygon contained in every left half-plane.
Half-plane intersection sorts the lines by direction, keeps only the tightest representative among
parallel lines, and maintains a deque of lines whose pairwise intersections form the current
boundary. Directions are ordered by half-plane and cross-product signs, without epsilon comparisons,
to satisfy `std::sort`'s strict ordering requirement. `EPS` is then used outside the comparator to
group parallel directions and retain only the tightest half-plane in each group.

- `half_plane_intersection(planes)` returns the polygon cut out by the half-planes in
  counter-clockwise order. Each half-plane is represented by `HalfPlane(p, q)`, meaning the closed
  region to the left of the directed line `p` $\to$ `q`. This implementation is intended for
  intersections that form a bounded polygon with positive area; empty, unbounded, or degenerate
  inputs may return an empty vector. Add explicit bounding-box half-planes when a bounded polygon is
  required. The two points defining each half-plane must differ.

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

Space Complexity:
- O(n) auxiliary.

*/

#include <algorithm>
#include <cmath>
#include <deque>
#include <type_traits>
#include <vector>

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

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); }
  Point operator+(const Point &p) const { return {x + p.x, y + p.y}; }
  Point operator-(const Point &p) const { return {x - p.x, y - p.y}; }
  Point operator*(double k) const { return {x * k, y * k}; }
  double cross(const Point &p) const { return x * p.y - y * p.x; }
};

struct HalfPlane {
  Point p, dir;

  HalfPlane() = default;
  HalfPlane(Point p, Point q) : p(p), dir(q - p) {}

  bool out(const Point &q) const { return LT(dir.cross(q - p), 0); }

  Point intersect(const HalfPlane &h) const {
    double t = h.dir.cross(p - h.p) / dir.cross(h.dir);
    return p + dir * t;
  }

  bool operator<(const HalfPlane &h) const {
    auto half = [](const Point &d) { return d.y < 0 || (d.y == 0 && d.x < 0); };
    if (half(dir) != half(h.dir)) {
      return half(dir) < half(h.dir);
    }
    return dir.cross(h.dir) > 0;
  }
};

std::vector<Point> half_plane_intersection(std::vector<HalfPlane> planes) {
  std::sort(planes.begin(), planes.end());
  std::vector<HalfPlane> unique;
  for (const HalfPlane &h : planes) {
    if (!unique.empty() && EQ(unique.back().dir.cross(h.dir), 0)) {
      if (unique.back().dir.cross(h.p - unique.back().p) > 0) {
        unique.back() = h;
      }
      continue;
    }
    unique.push_back(h);
  }
  std::deque<HalfPlane> dq;
  auto bad_back = [&](const HalfPlane &h) {
    return dq.size() >= 2 && h.out(dq[dq.size() - 2].intersect(dq.back()));
  };
  auto bad_front = [&](const HalfPlane &h) {
    return dq.size() >= 2 && h.out(dq[0].intersect(dq[1]));
  };
  for (const HalfPlane &h : unique) {
    while (bad_back(h)) {
      dq.pop_back();
    }
    while (bad_front(h)) {
      dq.pop_front();
    }
    dq.push_back(h);
  }
  while (dq.size() >= 3 && dq.front().out(dq[dq.size() - 2].intersect(dq.back()))) {
    dq.pop_back();
  }
  while (dq.size() >= 3 && dq.back().out(dq[0].intersect(dq[1]))) {
    dq.pop_front();
  }
  if (dq.size() < 3) {
    return {};
  }
  std::vector<Point> res;
  for (int i = 0; i < static_cast<int>(dq.size()); i++) {
    if (EQ(dq[i].dir.cross(dq[(i + 1) % dq.size()].dir), 0)) {
      return {};
    }
    res.push_back(dq[i].intersect(dq[(i + 1) % dq.size()]));
  }
  return res;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

double polygon_area(const vector<Point> &p) {
  double area = 0;
  for (int i = 0, j = static_cast<int>(p.size()) - 1; i < static_cast<int>(p.size()); j = i++) {
    area += p[j].cross(p[i]);
  }
  return fabs(area) / 2.0;
}

int main() {
  vector<HalfPlane> box{
      HalfPlane(Point(0, 0), Point(4, 0)),
      HalfPlane(Point(4, 0), Point(4, 4)),
      HalfPlane(Point(4, 4), Point(0, 4)),
      HalfPlane(Point(0, 4), Point(0, 0)),
  };
  auto square = half_plane_intersection(box);
  assert(square.size() == 4);
  assert(EQ(polygon_area(square), 16));

  auto tighter = box;
  tighter.push_back(HalfPlane(Point(0, 1), Point(4, 1)));  // y >= 1
  assert(EQ(polygon_area(half_plane_intersection(tighter)), 12));

  box.push_back(HalfPlane(Point(2, 0), Point(2, 4)));  // x <= 2
  auto rect = half_plane_intersection(box);
  assert(rect.size() == 4);
  assert(EQ(polygon_area(rect), 8));

  vector<HalfPlane> empty{
      HalfPlane(Point(0, 0), Point(1, 0)),  // y >= 0
      HalfPlane(Point(1, 1), Point(0, 1)),  // y <= 1
      HalfPlane(Point(0, 2), Point(1, 2)),  // y >= 2
  };
  assert(half_plane_intersection(empty).empty());
  return 0;
}