Привет, codeforces!
Сегодня, ориентируясь на статью Триангуляция полигонов с сайта neerc.ifmo.ru, я реализовал триангуляцию любых многоугольников без самопересечений (в том числе не являющихся выпуклыми) на C++.
$$$\text{ }$$$ $$$\text{ }$$$
Зачем мне это понадобилось? На отбор на RuCode 7.0 в этом году дали задачу H. Территориальное размежевание. Задача ставится следующим образом: вам дана бесконечная плоскость, покрашенная зеброй (чёрная полоска высоты $$$1$$$, белая полоска высоты $$$1$$$, чёрная ..., белая ...). На плоскости дан многоугольник без самопересечений, состоящий из $$$n\text{ }\left(3\le n \le 100000\right)$$$ вершин с целочисленными координатами $$$x_i, y_i \left(-10^9 \le x_i, y_i \le 10^9\right)$$$. Вершины указаны в порядке какого-то обхода. Требуется посчитать чёрную площадь внутри многоугольника.
Сначала я решил для случая, когда фигура — треугольник, а потом подумал, почему бы не триангулировать исходный многоугольник и не вызвать решение для треугольника $$$\left(n-2\right)$$$ раза?
Итак, всё упёрлось в триангуляцию. К сожалению, в статье выше многие моменты опускаются, есть всякие ошибки (к примеру, автор говорит, что вершина типа merge соединяется только с вершиной типа merge, хотя в примере, который он приводит, это не так, а в псевдокоде заифано на merge...). Плюс псевдокод, если его реализовывать так, как там написано, не проходит стресс-тестирование.
Я всё исправил и предоставляю свою реализацию. Итак, чтобы посчитать триангуляцию, надо вызвать следующую функцию:
std::vector<Triangle> triangulate(vpii points) {
makeCounterClockwise(points);
return triangulateFast(points);
}
Она вернёт вам вектор из треугольников, с которыми вы можете делать всё, что хотите. Будьте внимательны с тем, что я вырезаю вершины, которые лежат на одной прямой с соседними вершинами, ибо они бесполезны (т.е., если $$$P_{i-1}$$$, $$$P_{i}$$$ и $$$P_{i+1}$$$ лежат на одной прямой, то $$$P_i$$$ вырезается.
Исходный код (полное решение задачи H. Территориальное размежевание)
Исходный код
#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
#define isz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
void remin(auto &x, const auto &y) { if (x > y) x = y; }
void remax(auto &x, const auto &y) { if (x < y) x = y; }
using ll = long long;
using ld = long double;
using pii = std::pair<int,int>;
using vi = std::vector<int>;
using vpii = std::vector<pii>;
using vvpii = std::vector<vpii>;
using tiii = std::tuple<int,int,int>;
using vtiii = std::vector<tiii>;
#ifndef __GEOMA_HPP__
#define __GEOMA_HPP__
namespace algos {
namespace geoma {
template<typename A, typename B>
inline std::pair<A,B> operator-(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first - b.first, a.second - b.second);
}
template<typename A, typename B>
inline std::pair<A,B> operator+(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first + b.first, a.second + b.second);
}
template<typename T, typename A, typename B>
inline T cross(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return T(a.first) * b.second - T(a.second) * b.first;
}
template<typename T, typename A, typename B>
inline T doubleSquare(const std::pair<A,B> &a, const std::pair<A,B> &b, const std::pair<A,B> &c) {
T res = cross<T>(a-b, c-b);
if (res < 0)
res = -res;
return res;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A,B> &a) {
return T(a.first)*a.first + T(a.second)*a.second;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A, B> &a, const std::pair<A, B> &b) {
return dist2<T>(a-b);
}
using Triangle = std::tuple<pii,pii,pii>;
inline bool isPointInsideTriangle(Triangle tr, pii p) {
auto [A, B, C] = tr;
auto S = doubleSquare<ll>(A,B,C);
auto S1 = doubleSquare<ll>(p,B,C);
auto S2 = doubleSquare<ll>(A,p,C);
auto S3 = doubleSquare<ll>(A,B,p);
return S == S1 + S2 + S3;
}
inline ll orientedSquare(const std::vector<pii> &points) {
ll res = cross<ll>(points.back(), points.front());
for (int i = 0; i + 1 < isz(points); i++)
res += cross<ll>(points[i], points[i+1]);
return res;
}
inline vpii readPolygon() {
int n;
vpii points;
std::cin >> n;
points.resize(n);
for (auto &[x,y] : points)
std::cin >> x >> y;
return points;
}
inline void moveToPositive(vpii &points) {
int minX = INT_MAX, minY = INT_MAX;
for (auto &[x,y] : points) {
remin(minX, x);
remin(minY, y);
}
if (minX % 2 != 0) minX--;
if (minY % 2 != 0) minY--;
for (auto &[x,y] : points)
x -= minX, y -= minY;
}
inline void transpose(vpii &points) {
for (auto &[x,y] : points)
std::swap(x, y);
}
inline void makeCounterClockwise(vpii &points) {
ll sq = orientedSquare(points);
if (sq < 0) std::reverse(all(points));
}
inline int sign(auto x) { return x < 0 ? -1 : (x > 0 ? +1 : 0); }
const int LEFT = -1, RIGHT = +1;
inline int rotation(pii first, pii mid, pii last) {
// > 0 --> right
// < 0 --> left
return sign(cross<ll>(first-mid,last-mid));
}
inline bool isPointOnSegment(pii P, pii A, pii B) {
auto [x, y] = P;
auto [xmin,xmax] = std::minmax({A.first, B.first});
auto [ymin,ymax] = std::minmax({A.second, B.second});
if (x < xmin || x > xmax || y < ymin || y > ymax)
return false;
return cross<ll>(A-P,B-P) == 0;
}
const int NONE = 0, POINT = 1, INSIDE = 2;
inline int intersects(pii A, pii B, pii C, pii D) {
if (A == C || A == D || B == C || B == D)
return POINT;
if (isPointOnSegment(A,C,D) || isPointOnSegment(B,C,D) ||
isPointOnSegment(C,A,B) || isPointOnSegment(D,A,B))
return POINT;
auto [xmin1,xmax1] = std::minmax({A.first, B.first});
auto [ymin1,ymax1] = std::minmax({A.second, B.second});
auto [xmin2,xmax2] = std::minmax({C.first, D.first});
auto [ymin2,ymax2] = std::minmax({C.second, D.second});
if (xmax1 < xmin2 || xmin1 > xmax2)
return NONE;
if (ymax1 < ymin2 || ymin1 > ymax2)
return NONE;
auto rot1 = rotation(A,B,C);
auto rot2 = rotation(A,B,D);
if (sign(rot1) * sign(rot2) < 0)
return INSIDE;
return NONE;
}
inline bool isAlmostEqual(ld x, ld y) {
const ld EPS = 1e-10;
return std::abs(x-y) / (1 + std::abs(x)) < EPS;
}
inline bool angleLessThanPi(pii L, pii M, pii R) {
return rotation(L,M,R) == LEFT;
}
inline bool angleGreaterThanPi(pii L, pii M, pii R) {
return rotation(L,M,R) == RIGHT;
}
inline auto getPrev(const auto &vec, int i) {
i--;
if (i < 0) i = isz(vec)-1;
return vec[i];
}
inline auto getNext(const auto &vec, int i) {
i++;
if (i == isz(vec))
i = 0;
return vec[i];
}
inline bool isConvex(const vpii &P) {
for (int i = 0; i < isz(P); i++)
if (rotation(getPrev(P, i), P[i], getNext(P, i)) == RIGHT)
return false;
return true;
}
inline void removeSameLine(vpii &P) {
static std::vector<bool> keep;
keep.assign(isz(P), false);
for (int i = 0; i < isz(P); i++)
keep[i] = (rotation(getPrev(P,i), P[i], getNext(P,i)) != 0);
int p = 0;
for (int i = 0; i < isz(P); i++)
if (keep[i]) P[p++] = P[i];
P.resize(p);
}
} // namespace geoma
} // namespace algos
#endif // __GEOMA_HPP__
#ifndef __TRIANGULATE_HPP__
#define __TRIANGULATE_HPP__
namespace algos {
namespace triangulate {
using namespace geoma;
inline bool liesBelow(pii P1, pii P2) {
auto [x1,y1] = P1;
auto [x2,y2] = P2;
return (y1 < y2 || (y1 == y2 && x1 > x2));
}
inline bool liesAbove(pii P1, pii P2) {
auto [x1,y1] = P1;
auto [x2,y2] = P2;
return (y1 > y2 || (y1 == y2 && x1 < x2));
}
const int START = 0, SPLIT = 1, END = 2, MERGE = 3, REGULAR = 4;
inline auto getTypeOfVertex(const auto &P, int i) {
pii L = getPrev(P, i), M = P[i], R = getNext(P, i);
if (liesBelow(L,M) && liesBelow(R,M) && angleLessThanPi(L,M,R))
return START;
if (liesBelow(L,M) && liesBelow(R,M) && angleGreaterThanPi(L,M,R))
return SPLIT;
if (liesAbove(L,M) && liesAbove(R,M) && angleLessThanPi(L,M,R))
return END;
if (liesAbove(L,M) && liesAbove(R,M) && angleGreaterThanPi(L,M,R))
return MERGE;
return REGULAR;
}
inline std::string type2str(int t) {
switch(t) {
case START: return "start";
case SPLIT: return "split";
case END: return "end";
case MERGE: return "merge";
case REGULAR: return "regular";
default: return "none";
};
}
inline bool isMonotone(const vpii &P) {
if (isz(P) <= 3)
return true;
// we should start from max vertex
int imax = -1, ymax = INT_MIN;
for (int i = 0; i < isz(P); i++)
if (P[i].second >= ymax) {
imax = i;
ymax = P[i].second;
}
int len = 0;
int i = imax;
while (len <= isz(P) && P[i].second >= getNext(P, i).second) {
len++;
i++;
if (i == isz(P)) i = 0;
}
while (len <= isz(P) && P[i].second <= getNext(P, i).second) {
len++;
i++;
if (i == isz(P)) i = 0;
}
return len >= isz(P);
}
inline vpii findDiagonals(const vpii &P) {
static vtiii pq;
pq.clear();
for (int i = 0; i < isz(P); i++) {
auto [x,y] = P[i];
pq.emplace_back(y,-x,i);
}
std::sort(all(pq), std::greater<>());
int currY, currX;
auto setCurrY = [&](int newCurrY) { currY = newCurrY; };
auto setCurrX = [&](int newCurrX) { currX = newCurrX; };
static vi helper;
helper.assign(isz(P),-1);
auto getHelper = [&](int i) {
if (i < 0) i += isz(P);
assert(0 <= i && i < isz(P));
assert(helper[i] != -1);
return helper[i];
};
auto setHelper = [&](int i, int j) {
if (i < 0) i += isz(P);
assert(0 <= i && i < isz(P));
if (helper[i] == i)
helper[i] = j;
else if (helper[i] != -1) {
if (int k = helper[i]; P[k].second > P[j].second)
helper[i] = j;
} else
helper[i] = j;
};
auto getEdge = [&](int i) {
if (i < 0) i += isz(P);
if (i >= isz(P)) i -= isz(P);
auto L = P[i], R = getNext(P, i);
auto [ymin, ymax] = std::minmax({L.second, R.second});
return tiii(ymin, ymax, i);
};
auto getEdgeId = [&](int i) {
auto [y1,y2,j] = getEdge(i);
if (y1 == y2) return -1;
return j;
};
int cachedId = -1;
ld cachedX = -1;
auto isCached = [&](int id) { return id == cachedId; };
auto setCached = [&](int id, ld x) { cachedId = id; cachedX = x; };
auto clearCached = [&](){cachedId = -1; cachedX = -1;};
auto coordOfIntersection = [&](int i, int y) -> ld {
if (i < 0) i += isz(P);
if (isCached(i))
return cachedX;
assert(0 <= i && i < isz(P));
auto [x1,y1] = P[i];
auto [x2,y2] = getNext(P, i);
//ld B = ld(x2-x1)/(y2-y1);
//ld A = x1 - y1 * B;
// x1 - y1 * (x2-x1)/(y2-y1)
// x1 - y1 * (x2-x1)/(y2-y1) + y * (x2-x1)/(y2-y1)
// (x1 * (y2 - y1) + (y - y1) * (x2 - x1)) / (y2 - y1)
// (x1 * y2 + y * x2 - y * x1 - y1 * x2) / (y2 - y1)
assert(y2 != y1);
return ld((ld)x1 * y2 + (ld)y * x2 - (ld)y * x1 - (ld)y1 * x2) / (y2 - y1);
};
auto LessByIntersection = [&](const int &j1, const int &j2) {
ld x1 = coordOfIntersection(j1, currY);
ld x2 = coordOfIntersection(j2, currY);
if (x1 > x2) return false;
if (x1 < x2) return true;
return j1 < j2;
};
// y1 + t * (y2-y1) = y => t = (y - y1)/(y2 - y1)
// x1 + t * (x2-x1) = x => x1 + (y - y1) * (x2-x1) / (y2 - y1)
// coeff = (x2 - x1) / (y2 - y1)
// x1 + (y - y1) * coeff <= x
// x1 - y1 * coeff + y * coeff <= x
// A + y * B <= x -- need to sort by these values
std::set<int, decltype(LessByIntersection)> T(LessByIntersection);
auto insertEdgeInT = [&](int i) {
int ei = getEdgeId(i);
if (ei == -1) return;
auto y = P[i].second;
setCurrY(y);
setCached(ei, coordOfIntersection(ei, currY));
T.insert(ei);
clearCached();
};
auto deleteEdgeFromT = [&](int i) {
int j = getEdgeId(i);
if (j == -1) return;
auto it = T.find(j);
assert(it != T.end());
T.erase(it);
};
auto searchEdgeInT = [&](int i) {
assert(0 <= i && i < isz(P));
const auto [x, y] = P[i];
setCurrY(y); setCurrX(x);
int cand;
setCached(isz(P), x);
auto it = T.upper_bound(isz(P));
assert(it != T.begin());
cand = *std::prev(it);
clearCached();
return cand;
};
vpii diagonals;
auto addEdgeInTriangles = [&](int u, int v) {
if (getNext(P, u) != P[v] && getPrev(P, u) != P[v])
diagonals.emplace_back(u,v);
};
for (auto [y,x,i] : pq) {
x *= -1;
setCurrY(y);
switch (getTypeOfVertex(P,i)) {
case START: {
insertEdgeInT(i);
setHelper(i, i);
break;
}
case END: {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
break;
}
case SPLIT: {
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
insertEdgeInT(i);
setHelper(i, i);
break;
}
case MERGE: {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
break;
}
case REGULAR: {
if (liesAbove(getPrev(P,i), P[i]) && liesBelow(getNext(P,i), P[i])) {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
insertEdgeInT(i);
setHelper(i, i);
} else {
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
}
break;
}
};
}
return diagonals;
}
struct Node {
const vpii * parent = nullptr;
Node * next = nullptr;
Node * prev = nullptr;
int where{}, id{};
std::list<int> diagonals{};
void setNext(Node *n) {
next = n;
n->prev = this;
}
void setPrev(Node *n) {
prev = n;
n->next = this;
}
vpii to_vector() const {
vpii result(1, (*parent)[id]);
for (Node *it = next; it != this; it = it->next)
result.push_back((*parent)[it->id]);
return result;
}
pii to_point() const {
return (*parent)[id];
}
};
template<typename Func>
inline void splitByDiag(const vpii &P, const vpii &D, const Func &f) {
if (isz(P) <= 2)
return;
if (isz(P) == 3) {
assert(D.empty());
f(P);
return;
}
static std::vector<Node *> nodes;
nodes.clear();
for (int i = 0; i < isz(P); i++)
nodes.emplace_back(new Node{&P, nullptr, nullptr, 0, i, {}});
for (int i = 0; i+1 < isz(nodes); i++)
nodes[i]->setNext(nodes[i+1]);
nodes.front()->setPrev(nodes.back());
static std::vector<std::pair<Node*, Node*>> diagonals;
diagonals.clear();
for (const auto &[u,v] : D) {
assert(0 <= u && u < isz(nodes));
assert(0 <= v && v < isz(nodes));
int currDiag = isz(diagonals);
diagonals.push_back({nodes[u], nodes[v]});
nodes[u]->diagonals.push_back(currDiag);
nodes[v]->diagonals.push_back(currDiag);
}
int cntPoly = 1;
for (int id = 0; id < isz(diagonals); id++) {
// need to split by this diagonal
auto [startFi, startSe] = diagonals[id];
auto currFi = startFi, currSe = startSe;
while (currFi != startSe && currSe != startFi) {
currFi = currFi->next;
currSe = currSe->next;
}
// need to copy ends of diagonal:
nodes.push_back(new Node{&P, nullptr, nullptr, cntPoly, startFi->id, {}});
Node *copyFi = nodes.back();
nodes.push_back(new Node{&P, nullptr, nullptr, cntPoly, startSe->id, {}});
Node *copySe = nodes.back();
// choose a list with min len:
if (currFi != startSe) {
std::swap(copyFi, copySe);
std::swap(startFi, startSe);
std::swap(currFi, currSe);
}
assert(currFi == startSe);
// iterate over all vertices and updates diagonals
for (currFi = currFi->prev; currFi != startFi; currFi = currFi->prev) {
auto &dd = currFi->diagonals;
for (auto it = dd.begin(); it != dd.end(); ) {
// this diagonal still exists?
if (*it <= id) {
it = dd.erase(it);
continue;
}
auto &[lhs, rhs] = diagonals[*it];
if (lhs == startFi) (lhs = copyFi)->diagonals.push_back(*it);
if (rhs == startFi) (rhs = copyFi)->diagonals.push_back(*it);
if (lhs == startSe) (lhs = copySe)->diagonals.push_back(*it);
if (rhs == startSe) (rhs = copySe)->diagonals.push_back(*it);
it = std::next(it);
}
// set where this vertex lies
currFi->where = cntPoly;
}
// disconnect less half:
copyFi->setNext(currFi->next);
copySe->next = copyFi;
copySe->setNext(copyFi);
startSe->prev->setNext(copySe);
startFi->setNext(startSe);
cntPoly++;
}
static std::vector<Node *> anyNode;
anyNode.resize(cntPoly);
for (auto it : nodes)
anyNode[it->where] = it;
for (int i = 0; i < cntPoly; i++) {
vpii temp = anyNode[i]->to_vector();
removeSameLine(temp);
//assert(isMonotone(temp));
f(temp);
}
for (auto node : nodes)
delete node;
}
inline vvpii splitByDiag(const vpii &P, const vpii &D) {
vvpii result;
splitByDiag(P, D, [&](const vpii &PP){ result.push_back(PP); });
return result;
}
inline vpii TriangulateMonotonePolygon(const vpii &P) {
if (isz(P) <= 3)
return {};
// we need to create chains
int ymax = INT_MIN, ymin = INT_MAX, imax = -1;
int leftChainBegin = -1, leftChainEnd = -1;
for (int i = 0; i < isz(P); i++) {
if (P[i].second >= ymax) {
imax = i;
ymax = P[i].second;
}
if (P[i].second == ymax && getNext(P, i).second != ymax) {
leftChainBegin = i;
}
remin(ymin, P[i].second);
}
const int LEFT_CHAIN = 0, RIGHT_CHAIN = 1;
assert(imax != -1);
static vi chainType; chainType.assign(isz(P), -1);
static vi chainPos; chainPos.assign(isz(P), -1);
assert(leftChainBegin != -1);
leftChainEnd = leftChainBegin;
int pos = 0;
while (!(P[leftChainEnd].second == ymin && getPrev(P, leftChainEnd).second != ymin)) {
chainType[leftChainEnd] = LEFT_CHAIN;
chainPos[leftChainEnd] = pos++;
leftChainEnd++;
if (leftChainEnd == isz(P))
leftChainEnd = 0;
}
chainType[leftChainEnd] = LEFT_CHAIN;
if (chainPos[leftChainEnd] == -1)
chainPos[leftChainEnd] = pos++;
pos = 0;
for (int i = (leftChainBegin-1+isz(P)) % isz(P); chainType[i] == -1; ) {
chainType[i] = RIGHT_CHAIN;
chainPos[i] = pos++;
i--;
if (i < 0) i = isz(P)-1;
}
// {y-coord, chainId, posInChain}
auto cmp = [&](pii a, pii b) {
auto [ya, ida] = a;
auto [yb, idb] = b;
if (ya < yb) return true;
if (ya > yb) return false;
if (chainType[ida] != chainType[idb]) {
return chainType[ida] < chainType[idb];
}
return chainPos[ida] > chainPos[idb];
};
static vpii pq;
pq.clear();
for (int i = 0; i < isz(P); i++) {
auto [x,y] = P[i];
pq.push_back({y,i});
}
std::sort(pq.rbegin(), pq.rend(), cmp);
static vi S;
S.clear();
S.push_back(std::get<1>(pq[0]));
S.push_back(std::get<1>(pq[1]));
vpii diagonals;
auto addDiagonal = [&](int i, int j) {
if (getNext(P, i) != P[j] && getPrev(P, i) != P[j]) {
diagonals.emplace_back(i, j);
return true;
}
return false;
};
auto isSameChain = [&](int i, int j) {
return chainType[i] == chainType[j];
};
for (const auto &[y, id] : pq) {
if (!isSameChain(id,S.back())) {
int last = S.back();
while (isz(S) > 1) {
addDiagonal(id, S.back());
S.pop_back();
}
S = {last, id};
} else {
int want = LEFT, skip = RIGHT;
if (chainType[id] == LEFT_CHAIN)
std::swap(want, skip);
pii first = P[id];
while (isz(S) >= 2) {
pii mid = P[S[isz(S)-1]];
pii last = P[S[isz(S)-2]];
ll rot = rotation(first, mid, last);
if (rot == skip)
break;
if (rot != want) {
S.pop_back();
continue;
}
assert(rot == want);
addDiagonal(id, S[isz(S)-2]);
S.pop_back();
}
S.push_back(id);
}
}
return diagonals;
}
template<typename Func>
inline void triangulateFast(vpii inP, const Func &f) {
ll fullSquare = orientedSquare(inP);
removeSameLine(inP);
vpii D = findDiagonals(inP);
splitByDiag(inP, D, [&](vpii PP){
removeSameLine(PP);
vpii DD = TriangulateMonotonePolygon(PP);
auto ff = [&](const vpii &t){
fullSquare -= orientedSquare(t);
f(t);
};
splitByDiag(PP, DD, ff);
});
assert(fullSquare == 0);
}
inline std::vector<Triangle> triangulateFast(const vpii &inP) {
std::vector<Triangle> answ;
triangulateFast(inP, [&](const vpii &t){
assert(isz(t) <= 3);
if (isz(t) == 3) answ.emplace_back(t[0], t[1], t[2]);
});
return answ;
}
inline std::vector<Triangle> triangulateSlow(const vpii &points)
{
// slow triangulation in time O(n^2)
std::vector<Triangle> triangles;
std::list<pii> list(all(points));
while (isz(list) > 3) {
// search a vertex which we can cut off:
for (auto it = list.begin(); it != list.end(); ) {
auto next = (std::next(it) == list.end() ? list.begin() : std::next(it));
auto prev = (it == list.begin() ? std::prev(list.end()) : std::prev(it));
ll rot = rotation(*prev, *it, *next);
if (rot == 0) {
it = list.erase(it);
} else if (rot < 0) {
// check this vertex:
bool ok = true;
for (auto p : list) {
bool bad = isPointInsideTriangle({*prev, *it, *next}, p);
bad = bad && !isPointOnSegment(p, *prev, *it);
bad = bad && !isPointOnSegment(p, *it, *next);
bad = bad && !isPointOnSegment(p, *prev, *next);
if (bad) {
ok = false;
break;
}
}
if (ok) {
triangles.emplace_back(*prev, *it, *next);
it = list.erase(it);
} else {
it++;
}
} else {
it++;
}
}
}
if (isz(list) == 3)
triangles.emplace_back(*list.begin(), *(++list.begin()), *(++++list.begin()));
else assert(isz(list) < 3);
return triangles;
}
inline std::vector<Triangle> triangulateConvex(const vpii &points) {
// points is convex polygon
// take any vertex and make triangles
std::vector<Triangle> triangles;
for (int i = 1; i + 1 < isz(points); i++)
triangles.emplace_back(points[0], points[i], points[i+1]);
return triangles;
}
inline std::vector<Triangle> triangulate(vpii points) {
makeCounterClockwise(points);
return triangulateFast(points);
}
} // namespace triangulate
} // namespace algos
#endif // __TRIANGULATE_HPP__
using namespace algos::triangulate;
using namespace algos::geoma;
pii pickPointOutside(pii A, pii B) {
pii C1(A.first, B.second);
pii C2(B.first, A.second);
pii res = cross<ll>(A-C1, B-C1) < 0 ? C1 : C2;
return res;
}
std::pair<ld,ld> rectangle(int xmin, int ymin, int xmax, int ymax) {
ll dx = xmax - xmin, dy = ymax - ymin;
ld fi = (dx + 1) / 2 * dy;
ld se = dx / 2 * dy;
if (xmin % 2 != 0) std::swap(fi,se);
return {fi, se};
}
ll sum(ll n) {
// n^2 - (n-1)^2 + (n-2)^2 - (n-3)^2 + (n-4)^2 + ...
return n * (n+1) / 2;
}
std::pair<ld, ld> solveBaseCaseLD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
ld c = y / (2.0L * x);
std::pair<ld, ld> res{c * sum(x), c * sum(x-1)};
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseRU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseRD(xmax-xmin,ymax-ymin);
if (xmin % 2 != 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> triangle(pii A, pii B, pii C) {
// B - corner
// AC - diagonal
if (doubleSquare<ll>(A,B,C) == 0)
return {0, 0};
if (A.first > C.first)
std::swap(A, C);
ll xmin = std::min({A.first, B.first, C.first});
ll xmax = std::max({A.first, B.first, C.first});
ll ymin = std::min({A.second, B.second, C.second});
ll ymax = std::max({A.second, B.second, C.second});
pii LD(xmin,ymin), LU(xmin,ymax), RD(xmax,ymin), RU(xmax,ymax);
std::pair<ld,ld> res;
if (A.second > C.second) {
if (B == LD) res = solveBaseCaseLD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseRU(xmin,ymin,xmax,ymax);
} else {
if (B == RD) res = solveBaseCaseRD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseLU(xmin,ymin,xmax,ymax);
}
return res;
}
std::pair<ld,ld> twoRectangles(pii A, pii B, pii C) {
if (A.first > B.first) std::swap(A, B);
if (A.first > C.first) std::swap(A, C);
if (B.first > C.first) std::swap(B, C);
assert(A.first <= B.first && B.first <= C.first);
auto [ymin, ymax] = std::minmax({A.second,C.second});
if (B.second >= ymax || B.second <= ymin || B.first == A.first || B.first == C.first) {
remin(ymin, B.second);
remax(ymax, B.second);
return rectangle(A.first, ymin, C.first, ymax);
}
std::pair<ld,ld> rect = rectangle(A.first, ymin, C.first, ymax);
// above or below than diagonal?
if (rotation(A,C,B) <= 0) { // above
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(A.first, B.second, B.first, ymax);
else
rect = rect - rectangle(B.first, B.second, C.first, ymax);
} else { // below
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(B.first, ymin, C.first, B.second);
else
rect = rect - rectangle(A.first, ymin, B.first, B.second);
}
return rect;
}
std::pair<ld,ld> calcTriangle(pii A, pii B, pii C) {
if (cross<ll>(A-B,C-B) > 0)
std::swap(A, C);
auto temp = twoRectangles(A,B,C);
auto P1 = pickPointOutside(A,B);
auto P2 = pickPointOutside(B,C);
auto P3 = pickPointOutside(C,A);
if (!(A.first == B.first || A.second == B.second))
temp = temp - triangle(A,P1,B);
if (!(B.first == C.first || B.second == C.second))
temp = temp - triangle(B,P2,C);
if (!(C.first == A.first || C.second == A.second))
temp = temp - triangle(C,P3,A);
return temp;
}
ld solveFast(vpii points) {
moveToPositive(points);
transpose(points);
makeCounterClockwise(points);
ld answ = 0;
triangulateFast(points, [&](const vpii &tr){
assert(isz(tr) <= 3);
if (isz(tr) == 3)
answ += calcTriangle(tr[0],tr[1],tr[2]).first;
});
return answ;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(0);
std::cout << std::fixed << std::setprecision(12);
auto points = readPolygon();
std::cout << solveFast(points) << std::endl;
}
Это решение работает от $$$280$$$ до $$$296$$$ мс в системе Яндекс.Контест при $$$n=10^5$$$.
UPD. Обновил исходный код, удалив все неиспользуемые функции и закомментированные строки.
UPD 2. В комментариях на русском языке orz предложил решение оригинальной задачи за $$$O(N)$$$. Давайте применим тот же трюк, который делается при вычислении обычной площади многоугольника. Зафиксируем точку $$$P=(0,0)$$$ и посчитаем ориентированные площади треугольников $$$P A_i A_{i+1}$$$ для каждой стороны: если при обходе точек в данном порядке был левый поворот, то прибавим к ответу чёрную площадь внутри этого треугольника, иначе вычтем её из ответа.
Linear solution of problem H
#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
#define isz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
void remin(auto &x, const auto &y) { if (x > y) x = y; }
void remax(auto &x, const auto &y) { if (x < y) x = y; }
using ll = long long;
using ld = long double;
using pii = std::pair<int,int>;
using vi = std::vector<int>;
using vpii = std::vector<pii>;
using vvpii = std::vector<vpii>;
using tiii = std::tuple<int,int,int>;
using vtiii = std::vector<tiii>;
#ifndef __GEOMA_HPP__
#define __GEOMA_HPP__
namespace algos {
namespace geoma {
template<typename A, typename B>
inline std::pair<A,B> operator-(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first - b.first, a.second - b.second);
}
template<typename A, typename B>
inline std::pair<A,B> operator+(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first + b.first, a.second + b.second);
}
template<typename T, typename A, typename B>
inline T cross(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return T(a.first) * b.second - T(a.second) * b.first;
}
template<typename T, typename A, typename B>
inline T doubleSquare(const std::pair<A,B> &a, const std::pair<A,B> &b, const std::pair<A,B> &c) {
T res = cross<T>(a-b, c-b);
if (res < 0)
res = -res;
return res;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A,B> &a) {
return T(a.first)*a.first + T(a.second)*a.second;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A, B> &a, const std::pair<A, B> &b) {
return dist2<T>(a-b);
}
inline ll orientedSquare(const std::vector<pii> &points) {
ll res = cross<ll>(points.back(), points.front());
for (int i = 0; i + 1 < isz(points); i++)
res += cross<ll>(points[i], points[i+1]);
return res;
}
inline vpii readPolygon() {
int n;
vpii points;
std::cin >> n;
points.resize(n);
for (auto &[x,y] : points)
std::cin >> x >> y;
return points;
}
inline void moveToPositive(vpii &points) {
int minX = INT_MAX, minY = INT_MAX;
for (auto &[x,y] : points) {
remin(minX, x);
remin(minY, y);
}
if (minX % 2 != 0) minX--;
if (minY % 2 != 0) minY--;
for (auto &[x,y] : points)
x -= minX, y -= minY;
}
inline void transpose(vpii &points) {
for (auto &[x,y] : points)
std::swap(x, y);
}
inline void makeCounterClockwise(vpii &points) {
ll sq = orientedSquare(points);
if (sq < 0) std::reverse(all(points));
}
inline int sign(auto x) { return x < 0 ? -1 : (x > 0 ? +1 : 0); }
const int LEFT = -1, RIGHT = +1;
inline int rotation(pii first, pii mid, pii last) {
// > 0 --> right
// < 0 --> left
return sign(cross<ll>(first-mid,last-mid));
}
} // geoma
} // algos
#endif // __GEOMA_HPP__
using namespace algos::geoma;
pii pickPointOutside(pii A, pii B) {
pii C1(A.first, B.second);
pii C2(B.first, A.second);
pii res = cross<ll>(A-C1, B-C1) < 0 ? C1 : C2;
return res;
}
std::pair<ld,ld> rectangle(int xmin, int ymin, int xmax, int ymax) {
ll dx = xmax - xmin, dy = ymax - ymin;
ld fi = (dx + 1) / 2 * dy;
ld se = dx / 2 * dy;
if (xmin % 2 != 0) std::swap(fi,se);
return {fi, se};
}
ll sum(ll n) {
// n^2 - (n-1)^2 + (n-2)^2 - (n-3)^2 + (n-4)^2 + ...
return n * (n+1) / 2;
}
std::pair<ld, ld> solveBaseCaseLD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
ld c = y / (2.0L * x);
std::pair<ld, ld> res{c * sum(x), c * sum(x-1)};
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseRU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseRD(xmax-xmin,ymax-ymin);
if (xmin % 2 != 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> triangle(pii A, pii B, pii C) {
// B - corner
// AC - diagonal
if (doubleSquare<ll>(A,B,C) == 0)
return {0, 0};
if (A.first > C.first)
std::swap(A, C);
ll xmin = std::min({A.first, B.first, C.first});
ll xmax = std::max({A.first, B.first, C.first});
ll ymin = std::min({A.second, B.second, C.second});
ll ymax = std::max({A.second, B.second, C.second});
pii LD(xmin,ymin), LU(xmin,ymax), RD(xmax,ymin), RU(xmax,ymax);
std::pair<ld,ld> res;
if (A.second > C.second) {
if (B == LD) res = solveBaseCaseLD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseRU(xmin,ymin,xmax,ymax);
} else {
if (B == RD) res = solveBaseCaseRD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseLU(xmin,ymin,xmax,ymax);
}
return res;
}
std::pair<ld,ld> twoRectangles(pii A, pii B, pii C) {
if (A.first > B.first) std::swap(A, B);
if (A.first > C.first) std::swap(A, C);
if (B.first > C.first) std::swap(B, C);
assert(A.first <= B.first && B.first <= C.first);
auto [ymin, ymax] = std::minmax({A.second,C.second});
if (B.second >= ymax || B.second <= ymin || B.first == A.first || B.first == C.first) {
remin(ymin, B.second);
remax(ymax, B.second);
return rectangle(A.first, ymin, C.first, ymax);
}
std::pair<ld,ld> rect = rectangle(A.first, ymin, C.first, ymax);
// above or below than diagonal?
if (rotation(A,C,B) <= 0) { // above
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(A.first, B.second, B.first, ymax);
else
rect = rect - rectangle(B.first, B.second, C.first, ymax);
} else { // below
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(B.first, ymin, C.first, B.second);
else
rect = rect - rectangle(A.first, ymin, B.first, B.second);
}
return rect;
}
std::pair<ld,ld> calcTriangle(pii A, pii B, pii C) {
if (cross<ll>(A-B,C-B) > 0)
std::swap(A, C);
auto temp = twoRectangles(A,B,C);
auto P1 = pickPointOutside(A,B);
auto P2 = pickPointOutside(B,C);
auto P3 = pickPointOutside(C,A);
if (!(A.first == B.first || A.second == B.second))
temp = temp - triangle(A,P1,B);
if (!(B.first == C.first || B.second == C.second))
temp = temp - triangle(B,P2,C);
if (!(C.first == A.first || C.second == A.second))
temp = temp - triangle(C,P3,A);
return temp;
}
ld solveFast(vpii points) {
moveToPositive(points);
transpose(points);
makeCounterClockwise(points);
ld answ = 0;
pii A(0,0);
for (int i = 0; i < isz(points); i++) {
pii B = points[i], C = points[(i+1)%isz(points)];
ld curr = calcTriangle(A, B, C).first;
if (rotation(A,B,C) == LEFT) answ += curr;
else answ -= curr;
}
return answ;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(0);
std::cout << std::fixed << std::setprecision(12);
auto points = readPolygon();
std::cout << solveFast(points) << std::endl;
}
UPD 3. k1r1t0 предложил использовать вместо треугольников трапеции, основания которых параллельны осям Ox или Oy. Здесь описана идея метода трапеций.
Молодец, что разобрался со сложным алгоритмом.
Отмечу, что, если ты выкладываешь вот так решение на публичное обозрение, лучше постараться причесать код. Базовое причёсывание — убрать все дефайны, инклюды, функции и переменные, которыми ты не пользуешься, а также убрать все комментарии, которые не комментируют код, а являются результатом того, что ты закомментировал какую-то строчку.
Задача выглядит так, будто в ней не нужна триангуляция. Скажи, если бы задача звучала как «найдите площадь многоугольника», ты бы тоже сводил эту задачу к задаче «найдите площадь треугольника» при помощи триангуляции? Неужели не нашлось бы чего-то менее хардкорного?
Да в том-то и дело, что я не придумал, как решать это без триангуляции. Конечно, мы можем складывать треугольники и, что может быть полезно, вычитать треугольники из какой-то объемлющей фигуры, но это ни к чему не привело.
В дополнение скажу, что там есть подгруппы. Общий случай — это 4-я подгруппа. Третья подгруппа: выпуклый многоугольник на 100000 вершин. Вторая подгруппа: произвольный многоугольник на 1000 вершин. Первая подгруппа: дан прямоугольник. За первые 3 подгруппы дают 59 баллов.
А искать площадь многоугольника без триангуляции как?
Складывать определители второго порядка в процессе обхода многоугольника. Полученная площадь будет удвоенной. Знак плюс, если против часовой, знак минус, если по часовой, и ноль, если есть пересечение причём это пересечение является одной из заданных точек.
Ну точно так же можно и чёрную площадь посчитать. Перебираешь треугольники, пусть очередной треугольник $$$A_1 A_i A_{i + 1}$$$. Считаешь в нём чёрную площадь $$$S_i$$$. Если $$$A_1 A_i A_{i + 1}$$$ — обход треугольника против часовой стрелки, то $$$S_i$$$ добавляется к ответу, а если по часовой стрелке — вычитается из ответа.
Понял, спасибо, днём попробую. Я думал, что никто не сдал эту задачу за две недели отбора из-за того, что там триангуляция, а не из-за того, что никто не знает, откуда формула площади многоугольника берётся, какая идея там заложена и как можно её применить здесь.
Я подозреваю, сложность может быть в написании функции для треугольника. Она всё-таки далеко не очевидная.
А разве нельзя вместо треугольников считать площадь трапециями? Там тогда всё совсем просто должно быть
Можно. Да, должно быть ещё проще, если брать не треугольники, а трапеции с горизонтальными или вертикальными основаниями.
Да, я долго смотрел на треугольник чтобы увидеть, что это прямоугольник, из которого вырезали три прямоугольных треугольника, и то не всегда три и не всегда из одного
Сдал, выбрав в качестве $$$A_1$$$ точку $$$(0,0)$$$ и перебрав все стороны $$$A_i A_{i+1}$$$.
С трапециями пока не понял как делать.
UPD. Похоже, понял, вот статья
Wow I got chills, literally about 10 hours ago I was looking for a valid Subquadratic Polygon Triangulation after seeing the lecture 2 days ago.
But I found this code by Takanori Maehara which does the same thing in $$$O(N)$$$
How is this $$$O(N)$$$?
Without even looking at the code, $$$O(n)$$$ triangulation is known to be notoriously complicated algorithm. Possible, but definitely not as simple
my bad. I think it's $$$N^2$$$
https://github.com/spaghetti-source/algorithm/blob/master/geometry/_geom.cc#L729
This approach can be used in $$$O(N)$$$ solution of original problem H without triangulation. I updated the blog with a source code: if order of traversal of triangle $$$O A_i A_{i+1}$$$ is clockwise, then we can take the answer for this triangle with a sign "minus", otherwise — with a sign "plus".
Also, seems that there is an easier solution with trapezoids instead of triangles.