一个二维计算几何的模板, 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;
}
}