7.4.1 Convex Polygon Cut

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

/*

Given a convex polygon and a directed line from `p` to `q`, clips the polygon to the closed left
half-plane of that line. The polygon is walked edge by edge: vertices on the left side of the line
are kept, and whenever an edge crosses the line, the intersection point is appended to the output.

- `convex_cut(lo, hi, p, q)` returns the portion of the polygon lying on or to the left of the
  directed line from `p` to `q`. The input range $[`lo`, `hi`)$ must contain the vertices of a
  convex polygon in boundary order, either clockwise or counterclockwise. The returned polygon
  preserves that boundary order. The points `p` and `q` must differ.

The function is templated on the input point type. Side classification is done with cross products.
Edge-line intersection points are computed in floating point, and the returned polygon uses `Point`
with `double` coordinates.

Overflow warning: For integer-coordinate inputs, classification is exact only if the intermediate
products do not overflow.

Time Complexity:
- O(n) per call, where $n$ is the distance between `lo` and `hi`.

Space Complexity:
- O(n) for the returned polygon and O(1) auxiliary.

*/

#include <cassert>
#include <cmath>
#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); }
};

template<typename PtA, typename PtB, typename PtO>
auto cross(const PtA &a, const PtB &b, const PtO &o) {
  return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

template<typename PtA, typename PtO, typename PtB>
int turn(const PtA &a, const PtO &o, const PtB &b) {
  auto c = cross(a, b, o);
  return LT(c, 0) ? 1 : (LT(0, c) ? -1 : 0);
}

template<typename PtA, typename PtB>
int line_intersection(
    const PtA &p1, const PtA &p2, const PtB &p3, const PtB &p4, Point *p = nullptr
) {
  double a1 = static_cast<double>(p2.y) - p1.y, b1 = static_cast<double>(p1.x) - p2.x;
  double c1 = -(static_cast<double>(p1.x) * p2.y - static_cast<double>(p2.x) * p1.y);
  double a2 = static_cast<double>(p4.y) - p3.y, b2 = static_cast<double>(p3.x) - p4.x;
  double c2 = -(static_cast<double>(p3.x) * p4.y - static_cast<double>(p4.x) * p3.y);
  double x = -(c1 * b2 - c2 * b1), y = -(a1 * c2 - a2 * c1);
  double det = a1 * b2 - a2 * b1;
  if (EQ(det, 0)) {
    return (EQ(x, 0) && EQ(y, 0)) ? 1 : -1;
  }
  if (p != nullptr) {
    *p = Point(x / det, y / det);
  }
  return 0;
}

template<typename It, typename Pt>
std::vector<Point> convex_cut(It lo, It hi, const Pt &p, const Pt &q) {
  assert(!EQ(p.x, q.x) || !EQ(p.y, q.y));
  if (lo == hi) {
    return {};
  }
  std::vector<Point> res;
  for (It i = lo, j = hi - 1; i != hi; j = i++) {
    Point pj(static_cast<double>(j->x), static_cast<double>(j->y));
    Point pi(static_cast<double>(i->x), static_cast<double>(i->y));
    int d1 = turn(q, p, *j), d2 = turn(q, p, *i);
    if (d1 <= 0) {
      res.push_back(pj);
    }
    if (d1 * d2 < 0) {
      Point r;
      line_intersection(p, q, *j, *i, &r);
      res.push_back(r);
    }
  }
  return res;
}

/*** Example Usage ***/

using namespace std;

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

bool EQ(const vector<Point> &a, const vector<Point> &b) {
  return a.size() == b.size() &&
         equal(a.begin(), a.end(), b.begin(), [](const Point &p, const Point &q) {
           return EQ(p, q);
         });
}

int main() {
  {
    vector<Point> v{{1, 3}, {2, 2}, {2, 1}, {0, 0}, {-1, 3}};
    // Cut using the vertical line through (0, 0).
    vector<Point> c{{-1, 3}, {0, 3}, {0, 0}};
    assert(EQ(convex_cut(v.begin(), v.end(), Point(0, 0), Point(0, 1)), c));
  }
  {  // On a non-convex input, the result may be multiple disjoint polygons!
    vector<Point> v{{0, 0}, {2, 2}, {0, 4}, {3, 4}, {3, 0}};
    vector<Point> c{{1, 0}, {0, 0}, {1, 1}, {1, 3}, {0, 4}, {1, 4}};
    assert(EQ(convex_cut(v.begin(), v.end(), Point(1, 0), Point(1, 4)), c));
  }
  {
    vector<PointI> v{{0, 0}, {4, 0}, {4, 4}, {0, 4}};
    vector<Point> c{{0, 4}, {0, 0}, {2, 0}, {2, 4}};
    assert(EQ(convex_cut(v.begin(), v.end(), PointI(2, 0), PointI(2, 4)), c));
  }
  return 0;
}