[学习笔记] 计算几何模板
存一个跑得快用得爽的板子。
所有vector不多存任何一个点或直线。
#include <bits/stdc++.h>
using namespace std;
#define lep(i, l, r) for(int i = (l); i <= (r); i ++)
#define rep(i, l, r) for(int i = (l); i >= (r); i --)
#define Lep(i, l, r) for(int i = (l); i < (r); i ++)
#define Rep(i, l, r) for(int i = (l - 1); i >= (r); i --)
#define debug(...) fprintf (stderr, __VA_ARGS__)
#define pb push_back
#define fi first
#define se second
#define gc getchar
#define pc putchar
using i64 = long long;
using uint = unsigned int;
using ui64 = unsigned long long;
using pii = pair<int, int>;
using vi = vector<int>;
using db = double;
//using db = long double;
template<typename A, typename B> inline void Min(A &x, B y) { x = x < y ? x : y; }
template<typename A, typename B> inline void Max(A &x, B y) { x = x > y ? x : y; }
namespace Geometry {
const double eps = 1e-8;
#define lt(x, y) ((x) < (y) - eps)
#define gt(x, y) ((x) > (y) + eps)
#define le(x, y) ((x) <= (y) + eps)
#define ge(x, y) ((x) >= (y) - eps)
#define eq(x, y) (le(x, y) && ge(x, y))
#define dot(x, y, z) (((y) - (x)) * ((z) - (x)))
#define cross(x, y, z) (((y) - (x)) ^ ((z) - (x)))
struct vec2 {
double x, y;
vec2 (double _x = 0, double _y = 0) : x(_x), y(_y) {}
inline vec2 operator - () const { return vec2(- x, - y); }
inline vec2 operator + (const vec2 &t) const { return vec2(x + t.x, y + t.y); }
inline vec2 operator - (const vec2 &t) const { return vec2(x - t.x, y - t.y); }
inline vec2 operator * (double k) const { return vec2(x * k, y * k); }
inline vec2 operator / (double k) const { return * this * (1.0 / k); }
inline double operator * (const vec2 & t) const { return x * t.x + y * t.y; }
inline double operator ^ (const vec2 & t) const { return x * t.y - y * t.x; }
inline double norm() { return sqrt(x * x + y * y); }
inline double norm2() { return x * x + y * y; }
inline bool operator < (const vec2 &t) const { return lt(x, t.x) || le(x, t.x) && lt(y, t.y); }
inline bool operator == (const vec2 &t) const { return eq(x, t.x) && eq(y, t.y); }
inline void rotate(db t) { * this = vec2(cos(t) * x - sin(t) * y, sin(t) * x + cos(t) * y); }
} ;
struct line {
vec2 p1, p2;
line() {}
line(vec2 _p1, vec2 _p2) : p1(_p1), p2(_p2) {}
line(double k, double b) { p1 = vec2(0, b); p2 = vec2(1000, b + k * 1000); }
inline vec2 direct() const { return p2 - p1; }
} ;
inline bool parallel(const line &a, const line &b) {
return eq(a.direct() ^ b.direct(), 0);
}
inline vec2 intersection(const line &l1, const line &l2) {
double ls = cross(l1.p1, l1.p2, l2.p1);
double rs = cross(l1.p1, l1.p2, l2.p2);
return l2.p1 + (l2.p2 - l2.p1) * ls / (ls - rs);
}
inline void graham(vector<vec2> &vec, int type) {
if(vec.size() == 0) { cerr << "Hull is empty !\n" << endl; exit(0); }
sort(vec.begin(), vec.end());
vector<vec2> rec; rec.pb(* vec.begin()); int sz = 0;
for(int i = 1; i < vec.size(); i ++) {
if(type == 0) while(sz >= 1 && le(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
else while(sz >= 1 && ge(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
rec.push_back(vec[i]); sz ++;
}
swap(vec, rec);
}
inline void graham_full(vector<vec2> &vec) {
vector<vec2> v1 = vec, v2 = vec;
graham(v1, 0); graham(v2, 1);
v1.pop_back(); for(int i = v2.size() - 1; i >= 1; i --) v1.push_back(v2[i]); swap(vec, v1);
}
inline double convDiameter(vector<vec2> &vec) {
if(vec.size() == 2) { return (vec[0] - vec[1]).norm2(); }
int j = 2, n = vec.size() - 1;
double res = 0;
for(int i = 0; i < vec.size() - 1; i ++) {
while(abs(cross(vec[i], vec[i + 1], vec[j])) < abs(cross(vec[i], vec[i + 1], vec[j + 1]))) j = (j + 1) % n;
Max(res, max((vec[i] - vec[j]).norm2(), (vec[i + 1] - vec[j]).norm2()));
}
return res;
}
inline void HPI(vector<line> &lv) {
vector<pair<line, double> > sorted(lv.size());
for(int i = 0; i < lv.size(); i ++) sorted[i].fi = lv[i], sorted[i].se = atan2(lv[i].direct().y, lv[i].direct().x);
sort(sorted.begin(), sorted.end(), [] (auto a, auto b) -> bool {
if(eq(a.se, b.se)) {
if(lt(cross(a.fi.p1, a.fi.p2, b.fi.p2), 0)) return 1;
else return 0;
}
else return a.se < b.se;
} );
for(int i = 0; i < lv.size(); i ++) lv[i] = sorted[i].fi;
deque<line> q;
q.push_back(lv[0]);
for(int i = 1; i < lv.size(); i ++) if(! parallel(lv[i], lv[i - 1])) {
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_back();
else break;
}
while(q.size() > 1) {
vec2 p = intersection(* q.begin(), * ++ q.begin());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_front();
else break;
}
q.push_back(lv[i]);
}
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(q.begin() -> p1, q.begin() -> p2, p), 0)) q.pop_back();
else break;
}
lv = vector<line> (q.size());
for(int i = 0; i < q.size(); i ++) lv[i] = q[i];
}
inline void Polarsorting(vector<vec2> &v) {
vector<pair<vec2, double> > sorted(v.size());
for(int i = 0; i < v.size(); i ++) sorted[i].fi = v[i], sorted[i].se = atan2(v[i].y, v[i].x);
sort(sorted.begin(), sorted.end());
for(int i = 0; i < v.size(); i ++) v[i] = sorted[i].fi;
}
inline vector<vec2> Minkowski(vector<vec2> v1, vector<vec2> v2) {
// v1, v2 is sorted
vector<vec2> s1(v1.size()), s2(v2.size());
for(int i = 1; i < s1.size(); i ++) s1[i - 1] = v1[i] - v1[i - 1]; s1[s1.size() - 1] = v1[0] - v1[s1.size() - 1];
for(int i = 1; i < s2.size(); i ++) s2[i - 1] = v2[i] - v2[i - 1]; s2[s2.size() - 1] = v2[0] - v2[s2.size() - 1];
vector<vec2> hull(v1.size() + v2.size() + 1);
int p1 = 0, p2 = 0, cnt = 0;
hull[cnt ++] = v1[0] + v2[0];
while(p1 < s1.size() && p2 < s2.size()) {
hull[cnt] = hull[cnt - 1] + (ge(s1[p1] ^ s2[p2], 0) ? s1[p1 ++] : s2[p2 ++]);
cnt ++;
}
while(p1 < s1.size()) hull[cnt] = hull[cnt - 1] + s1[p1 ++], cnt ++;
while(p2 < s2.size()) hull[cnt] = hull[cnt - 1] + s2[p2 ++], cnt ++;
hull.pop_back();
return hull;
}
inline bool PointInHull(const vec2 &p, const vector<vec2> &poly) {
int l = 1, r = poly.size() - 2;
while(l <= r) {
int mid = (l + r) >> 1;
double a = ( (poly[mid] - poly[0]) ^ (p - poly[0]) );
double b = ( (poly[mid + 1] - poly[0]) ^ (p - poly[0]) );
if(ge(a, 0) && le(b, 0)) {
if(ge(cross(poly[mid], poly[mid + 1], p), 0)) return 1;
else return 0;
}
if(le(a, 0)) r = mid - 1;
else l = mid + 1;
} return 0;
}
}
using Geometry :: vec2;
using Geometry :: line;
using Geometry :: graham;
using Geometry :: graham_full;
using Geometry :: convDiameter;
using Geometry :: intersection;
using Geometry :: HPI;
using Geometry :: Polarsorting;
using Geometry :: Minkowski;
using Geometry :: PointInHull;
int main() {
return 0;
}
