[学习笔记] 计算几何模板

存一个跑得快用得爽的板子。
所有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;
}
上一篇:MySQL 语句汇总


下一篇:获取 本周、本月、本季度、本年 半年的开始时间或结束时间