【模板】 計算幾何大模板
阿新 • • 發佈:2019-01-01
一個二維計算幾何的模板, 13kb
#include <cstdio> #include <cstring> #include <cmath> #include <string> #include <algorithm> namespace g2{ using std :: swap; using std :: pair; using std :: max; using std :: min; //const double M_PI = 3.1415926535898; // global variables const int MAXT = 50000; const double EPS = 0.000001; bool tmp[MAXT]; const int CS_SINT = (1<<3) - 1; // strictly intersect const int CS_INT = (1<<4) - 1; // intersect const int CS_MOV = 1<<15; enum CASE {EQUAL = 1, CROSS = 1<<1, COVER = 1<<2, VERTX = 1<<3, SEPAR = 1<<4, NONE = CS_MOV - 1, ZERO = CS_MOV, ONE = CS_MOV+1, TWO = CS_MOV+2, THREE = CS_MOV+3, MULTP = CS_MOV+CS_MOV, INFI = CS_MOV*256, SUCCESS, FAIL}; // end // announcements of functions ans structures struct Vector ; struct Line ; inline double abs(const double & a) {return a<0 ? - a : a;} inline bool equal(const double & a, const double & b) {return abs(a - b) < EPS;} inline bool equal(const Vector & a, const Vector & b); inline bool equal(const Line & a, const Line & b); Vector operator + (const Vector a, const Vector b) ; Vector operator - (const Vector a, const Vector b) ; Vector operator * (const Vector a, const double b) ; Vector operator / (const Vector a, const double b) ; Vector operator += (Vector & a, const Vector b) ; Vector operator -= (Vector & a, const Vector b) ; Vector operator *= (Vector & a, const Vector b) ; Vector operator /= (Vector & a, const Vector b) ; double cross(const Vector a, const Vector b) ; double dot (const Vector a, const Vector b) ; Vector * sortPolygonVec(Vector * arr, int len) ; Vector * getPolygonEdges(Vector * ans, const Vector * arr, int len) ; int getVolumes(Vector * ans, const Vector * arr, int sz, const Vector * polyvec, int psz) ; double getPolygonSize (const Vector * arr, int sz) ; inline bool operator / (Line a, Line b) ; inline CASE haveIntc(Line a, Line b) ; //end struct Vector { double x,y; void * ad; Vector(): x(0.0), y(0.0), ad(NULL){} Vector(const double a, const double b, void * const c = NULL): x(a), y(b), ad(c){} Vector init(const double & a, const double & b, void * const c = NULL) {x = a, y = b; ad = c; return *this;} Vector getRotateR() const {return Vector(y, -x, ad);} double getLengthSquared() const {return x*x + y*y;} double getLength() const {return sqrt(x*x + y*y);} Vector getUnit() const {return (*this) / getLength();} Vector getRotate(double v) { return Vector(x*cos(v) - y*sin(v), y*cos(v) + x*sin(v), ad); } void print() { printf("(%.2lf, %.2lf)", x, y); } }; const Vector VUNITX(1, 0), VUNITY(0, 1); const Vector VORIGIN(0, 0); inline Vector operator - (const Vector a) {return Vector(-a.x, -a.y);} inline Vector operator + (const Vector a, const Vector b) {return Vector(a.x + b.x, a.y + b.y);} inline Vector operator - (const Vector a, const Vector b) {return Vector(a.x - b.x, a.y - b.y);} inline Vector operator * (const Vector a, const double b) {return Vector(a.x * b, a.y * b);} inline Vector operator / (const Vector a, const double b) {return Vector(a.x / b, a.y / b);} inline Vector operator += (Vector & a, const Vector b) {return a = a+b;} inline Vector operator -= (Vector & a, const Vector b) {return a = a-b;} inline Vector operator *= (Vector & a, const double b) {return a = a*b;} inline Vector operator /= (Vector & a, const double b) {return a = a/b;} inline bool equal(const Vector & a, const Vector & b) {return equal(a.x, b.x) && equal(a.y, b.y);} double cross(const Vector a, const Vector b) {return a.x*b.y - a.y*b.x;} double dot (const Vector a, const Vector b) {return a.x*b.x + a.y*b.y;} double getDistance(const Vector & a, const Vector & b) {return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y) * (a.y - b.y));} Vector _stv; double _dtmp; inline bool _cmp(const Vector & a, const Vector & b) { return cross(a - _stv, b - _stv) < 0; } Vector * sortPolygonVec(Vector * arr, int len) { _stv = arr[0]; std :: sort(arr, arr + len, _cmp); return arr; } Vector * getPolygonEdges(Vector * ans, const Vector * arr, int len) { for(int i = 0; i < len - 1; i++) ans[i] = arr[i+1] - arr[i]; ans[len - 1] = arr[0] - arr[len - 1]; return ans; } Vector mst[MAXT]; Vector tmpv[MAXT]; inline bool _hcmp(const Vector & a, const Vector & b) { if(equal(a.x, b.x)) return a.y < b.y; return a.x < b.x; } int getHull(Vector * ans, Vector * arr, int len) { for(int i = 0; i<len; i++) tmpv[i] = arr[i]; std :: sort(tmpv, tmpv + len, _hcmp); //for(int i = 0; i<len; i++) // tmpv[i].print(); int top = -1; mst[++top] = tmpv[0]; mst[++top] = tmpv[1]; for(int i = 2; i<len; i++){ _dtmp = cross(tmpv[i] - mst[top - 1], mst[top] - mst[top - 1]); if(!equal(_dtmp, 0) && _dtmp > 0){ if(top == 1) mst[top] = tmpv[i]; else -- top, i --; } else mst[++top] = tmpv[i]; } for(int i = len - 2; i>=0; i--){ _dtmp = cross(tmpv[i] - mst[top - 1], mst[top] - mst[top - 1]); if(!equal(_dtmp, 0) && _dtmp > 0){ if(top == 1) mst[top] = tmpv[i]; else -- top, i++; } else mst[++top] = tmpv[i]; } for(int i = 0; i<=top; i++) ans[i] = mst[i]; return top; } int getVolumes (Vector * ans, const Vector * arr, int sz, const Vector * polyvec, int psz) { // O(N^2) memset(tmp, 0, sz * (sizeof tmp[0])); double t; for(int i = 1; i<psz; i++) for(int j = 0; j<sz; j++){ if(tmp[j]) continue; t = cross(polyvec[i] - polyvec[i - 1], arr[j] - polyvec[i-1]); if(t <= 0 && !equal(t, 0.0)) tmp[j] = true; } for(int j = 0; j<sz; j++){ if(tmp[j]) continue; t = cross(polyvec[0] - polyvec[psz - 1], arr[j] - polyvec[psz - 1]); if(t <= 0 && !equal(t, 0.0)) tmp[j] = true; } int tot = 0; for(int i = 0; i<sz; i++) if(!tmp[i]) ans[tot++] = arr[i]; return tot; } double getPolygonSize (const Vector * arr, int sz) { double ans = 0; for(int i = 0; i<sz - 1; i++) ans += cross(arr[i+1], arr[i]); ans += cross(arr[0], arr[sz - 1]); return abs(ans); } struct Line { Vector a,b; bool isSeg; Line(){} Line(const Vector & aa, const Vector & bb, const bool var = true): a(aa), b(bb), isSeg(var) {} Vector getVector() const {return b - a;} Line getLine() const {return Line(a, b, false);} Line getSeg() const {return Line(a, b, true);} double getLength() const {return (b - a).getLength();} double getLengthSquared() const {return (b - a).getLengthSquared();} Line init(const Vector & aa, const Vector & bb, const bool var = false) {a = aa; b = bb; isSeg = var; return * this;} void print() { printf("[Line] "); a.print(); putchar(','); b.print(); printf(",[%d]\n", isSeg); } }; inline bool operator / (Vector a, Vector b) { return equal(cross(a, b), 0); } inline bool operator / (Line a, Line b) { return equal(cross(a.getVector(), b.getVector()), 0); } inline bool equal(const Line & a, const Line & b) {return equal(a.a, b.a) && equal(a.b, b.b);} inline CASE haveIntc(Line a, Line b){ if(!a.isSeg && !b.isSeg) { if(a / b) { if(equal(a.a.getUnit(), b.a.getUnit()) || equal(a.a.getUnit(), b.a.getUnit())) return EQUAL; return SEPAR; } else return CROSS; } if(a.isSeg && b.isSeg) // only returns SEPAR or CROSS return ((((int)haveIntc(a.getLine(), b)) & CS_INT ) &&(((int)haveIntc(b.getLine(), a)) & CS_INT )) ? CROSS : SEPAR; if(a.isSeg && !b.isSeg) swap(a, b); if(a / b) { if(equal(cross(a.a, b.a), 0)) return COVER; else return SEPAR; } //a.print(); //b.print(); double v1 = cross(b.b - a.a, a.getVector()), v2 = cross(b.a - a.a, a.getVector()); if(equal(v1, 0) || equal(v2, 0)) return VERTX; return ((v1 < 0) != (v2 < 0)) ? CROSS : SEPAR; } inline double getTrangArea(Vector a, Vector b, Vector c) { return abs(cross(a-c, b-c)) / 2; } inline Vector getIntc(Line a, Line b, bool flag = true) { if(flag && !((int)haveIntc(a,b) & CS_INT)) return *((Vector*)NULL) ; double S1 = getTrangArea(a.a, b.a, a.b); double S2 = getTrangArea(a.a, b.b, a.b); if(equal(S1, 0.0)) return b.a; if(equal(S2, 0.0)) return b.b; return b.a + (b.b - b.a) * (S1 / (S1+S2)); } struct Circle { Vector c; double r; Circle(){r = 0;} Circle(const Vector & vc, double rr): c(vc), r(rr){} }; inline bool equal(const Circle & a, const Circle & b) {return equal(a.r, b.r) && equal(a.c, b.c);} CASE getIntc(Circle a, Line b, pair<Vector, Vector> & p) { double dist = dot(a.c - b.a, b.b - b.a) / (b.b - b.a).getLength(); Vector pmid = b.a + (b.getVector().getUnit() * dist); double dish = cross(pmid - b.a, a.c - b.a) / dist; if(abs(dish) > a.r && !equal(dish, a.r)){ p = pair<Vector, Vector>(a.c, a.c); return NONE; } double row = sqrt(a.r*a.r - dish*dish); Vector inca = pmid - (pmid - b.a).getUnit() * row, incb = pmid + (pmid - b.a).getUnit() * row; if(equal(inca, incb)) { p = pair<Vector, Vector>(inca, a.c); return ONE; } p = pair<Vector, Vector> (inca, incb); return TWO; } Line getPerpend(Vector v, Line l){ if(l.isSeg) return *(Line*)NULL; return Line(v, v + l.getVector().getRotate(0.5*M_PI), false); } int getPublcTangent(Circle a, Circle b, Line * ans) { int tot = -1; double dist = getDistance(a.c, b.c); if(dist < abs(a.r - b.r) && !equal(dist, a.r - b.r)) return 0; if(abs(a.r - b.r) > dist || equal(dist, abs(a.r - b.r))){ double sinv = (b.r - a.r) / dist; double rlen = acos(sinv); Vector p2 = (b.c - a.c).getUnit().getRotate(rlen)*b.r; //TODO Vector p1 = (a.c - b.c).getUnit().getRotate(-M_PI-asin(sinv))*a.r; ans[++tot] = Line(a.c + p1, b.c + p2, false); ans[++tot] = Line(a.c - p1, b.c - p2, false); } if(equal(dist, abs(a.r - b.r))) -- tot; if(equal(a.r + b.r, dist)){ Vector pmid = a.c + (b.c - a.c) * (a.r / (a.r + b.r)); ans[++ tot] = getPerpend(pmid, Line(a.c, b.c, false)); } if(!equal(a.r + b.r, dist) && dist > a.r+b.r) { double sinv = (a.r + b.r) / dist; double rlen = acos(sinv); Vector p1 = (b.c - a.c).getUnit().getRotate(-rlen)*b.r; Vector p2 = (a.c - b.c).getUnit().getRotate(-rlen)*a.r; ans[++ tot] = Line(p1, p2, false); p1 = (b.c - a.c).getUnit().getRotate(rlen)*b.r; p2 = (a.c - b.c).getUnit().getRotate(rlen)*a.r;; ans[++ tot] = Line(p1, p2, false);; } return tot + 1; } }