模板-计算几何

计算几何算法相关模版, 可能有错误, 省选前持续更正中

重要的不是模版内容, 而是提供算法的实现思路.

struct Point {
	double x, y;
	Point(double x = 0, double y = 0): x(x), y(y) {}
};

typedef Point Vector;

Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
Vector operator - (Vector A, Vector B) { return Vector(A.x-B.x, A.y-B.y); }
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }

bool operator < (const Vector& a, const Vector& b) {
	return a.x < b.x || (a.x == b.x && a.y < b.y);
}

const double eps = 1e-10;

int dcmp(double x) {
	if(fabs(x) < eps) return 0;
	else return x < 0 ? -1 : 1;
}

bool operator == (const Vector& a, const Vector& b) {
	return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}

double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; }
double Area2(Point A, Point B, Point C) { return Cross(B-A, C-A); }

// 旋转公式 用 (length*cos(theta), length*sin(theta)) 表示向量即可推出
Vector Rotate(Vector A, double rad) {
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}

// 法向量
Vector Normal(Vector A) {
	double L = Length(A);
	return Vector(-A.y/L, A.x/L);
}

// P(x1, y1), v(x2, y2), Q(x3, y3), w(x4, y4), P-Q=u(x5, y5).
// P + v*t1 = Q + w*t2 // 用虚数表示如下
// x1+y1*i + t1*x2+t1*y2*i = x3+y3*i + t2*x4+t2*y4*i
// x1 + t1*x2 + (y1+t1*y2)*i = x3 + t2*x4 + (y3+t2*y4)*i

// => /x1+t1*x2 = x3+t2*x4 -> x5 + t1*x2 = t2*x4 -> x5*y4 + t1*x2*y4 = t2*x4*y4 (1)
//    \y1+t1*y2 = y3+t2*y4 -> y5 + t1*y2 = t2*y4 -> y5*x4 + t1*y2*x4 = t2*y4*x4 (2)

// (1)-(2) => (x5*y4-y5*x4) + t1(x2*y4-y2*x4) = 0
// => t1 = (x4*y5-x5*y4) / (x2*y4-y2*x4)
//       = Cross(w, u) / Cross(v, w)

Vector GetLineIntersection(Point P, Vector v, Point Q, Vector w) {
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t;
}

double DistanceToLine(Point P, Point A, Point B) {
	Vector v1 = B-A, v2 = P-A;
	return fabs(Cross(v1, v2)) / Length(v1);
}

double DistanceToSegment(Point P, Point A, Point B) {
	if(A == B) return Length(P-A); // 两点
	Vector v1 = B-A, v2 = P-A, v3 = P-B;
	if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
	if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
	return fabs(Cross(v1, v2) / Length(v1));
}

Point GetLineProjection(Point P, Point A, Point B) {
	Vector v = B-A;
	return A+v*(Dot(v, P-A) / Dot(v, v));
}

bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2) {
	double c1=Cross(a2-a1,b1-a1), c2=Cross(a2-a1,b2-a1), c3=Cross(b2-b1,a1-b1), c4=Cross(b2-b1,a2-b1);
	return (dcmp(c1)*dcmp(c2) < 0) && (dcmp(c3)*dcmp(c4) < 0);
}

bool OnSegment(Point P, Point a1, Point a2) {
	return dcmp(Cross(a1-P, a2-P)) == 0 && dcmp(Dot(a1-P, a2-P)) < 0;
}

double PolygonArea(Point* P, int n) {
	double area = 0;
	for(int i = 1; i < n-1; i++)
		area += Cross(P[i]-P[0], P[i+1]-P[0]);
	return area/2; // 有向面积
}

struct Line {
	Point p;
	Vector v;
	Point point(double a) {
		return p + v*a;
	}
};

const double PI = acos(double(-1));

struct Circle {
	Point c;
	double r;
	Circle(Point c, double r):c(c),r(r) {}
	Point point(double a) {
		return Point(c.x + cos(a)*r, c.y + sin(a)*r);
	}
};

int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector
   
   
    
    & sol) {
	double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
	double e = a*a + c*c, f = 2 * (a*b+c*d), g = b*b + d*d - C.r*C.r;
	double delta = f*f - 4*e*g;
	if(dcmp(delta) < 0) return 0;
	if(dcmp(delta) == 0) {
		t1 = t2 = -f / (2*e);
		sol.push_back(L.point(t1));
		return 1;
	}
	t1 = (-f - sqrt(delta)) / (2*e);
	sol.push_back(L.point(t1));
	t2 = (-f + sqrt(delta)) / (2*e);
	sol.push_back(L.point(t2));
	return 2;
}

double angle(Vector v) {
	return atan2(v.y, v.x);
}

int getCircleCircleIntersection(Circle C1, Circle C2, vector
    
    
     
     & sol) {
	double d = Length(C1.c - C2.c);
	if(dcmp(d) == 0) {
		if(dcmp(C1.r-C2.r) == 0) return -1;
		return 0;
	}
	if(dcmp(C1.r+C2.r-d) < 0) return 0;
	if(dcmp(fabs(C1.r-C2.r)-d > 0)) return 0;

	double a = angle(C2.c-C1.c);
	double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));

	Point p1 = C1.point(a-da), p2 = C1.point(a+da);
	sol.push_back(p1);
	if(p1 == p2) return 1;
	sol.push_back(p2);
	return 2;
}

int getTangents(Point p, Circle C, Vector* v) {
	Vector u = C.c-p;
	double dist = Length(u);
	if(dist < C.r) return 0;
	else if(dcmp(dist-C.r) == 0) {
		v[0] = Rotate(u, PI/2);
		return 1;
	} else {
		double ang = asin(C.r / dist);
		v[0] = Rotate(u, -ang);
		v[1] = Rotate(u, +ang);
		return 2;
	}
}

// 可用 dcmp 比较
// 如果不希望在凸包边上有输入点, 可以将 <= 换为 <
int ConvexHull(Point* p, Point* ch, int n) {
	sort(p, p+n); // x 为第一关键字, y 为第二关键字
	int m = 0;
	for(int i = 0; i < n; i++) {
		while(m > 1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
		ch[m++] = p[i];
	}
	int k = m;
	for(int i = n-2; i >= 0; i--) {
		while(m > k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
		ch[m++] = p[i];
	}
	if(n > 1) m--;
	return m;
}

bool OnLeft(const Line& L, const Point& p) {
	return Cross(L.v, p-L.P) > 0;
}

vector
     
     
      
       HalfplaneIntersection(vector
      
      
       
        L) {
	int n = L.size();
	sort(L.begin(), L.end()); // 按极角排序

	int first, last;          // 双端队列的第一个元素和最后一个元素的下标
	vector
       
       
         p(n); // p[i]为q[i]和q[i+1]的交点 vector 
        
          q(n); // 双端队列 vector 
         
           ans; // 结果 q[first=last=0] = L[0]; // 双端队列初始化为只有一个半平面L[0] for(int i = 1; i < n; i++) { while(first < last && !OnLeft(L[i], p[last-1])) last--; while(first < last && !OnLeft(L[i], p[first])) first++; q[++last] = L[i]; if(fabs(Cross(q[last].v, q[last-1].v)) < eps) { // 两向量平行且同向，取内侧的一个 last--; if(OnLeft(q[last], L[i].P)) q[last] = L[i]; } if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]); } while(first < last && !OnLeft(q[first], p[last-1])) last--; // 删除无用平面 if(last - first <= 1) return ans; // 空集 p[last] = GetLineIntersection(q[last], q[first]); // 计算首尾两个半平面的交点 // 从deque复制到输出中 for(int i = first; i <= last; i++) ans.push_back(p[i]); return ans; } int diameter2(vector 
          
            & points) { vector 
           
             p = ConvexHull(points); int n = p.size(); if(n == 1) return 0; if(n == 2) return Dist2(p[0], p[1]); p.push_back(p[0]); // 免得取模 int ans = 0; for(int u = 0, v = 1; u < n; u++) { // 一条直线贴住边p[u]-p[u+1] while(1) { // 当 Area(p[u], p[u+1], p[v+1]) <= Area(p[u], p[u+1], p[v])时停止旋转 // 即 Cross(p[u+1]-p[u], p[v+1]-p[u]) - Cross(p[u+1]-p[u], p[v]-p[u]) <= 0 // 根据 Cross(A,B) - Cross(A,C) = Cross(A,B-C) // 化简得 Cross(p[u+1]-p[u], p[v+1]-p[v]) <= 0 int diff = Cross(p[u+1]-p[u], p[v+1]-p[v]); if(diff <= 0) { ans = max(ans, Dist2(p[u], p[v])); // u和v是对踵点 if(diff == 0) ans = max(ans, Dist2(p[u], p[v+1])); // diff == 0时u和v+1也是对踵点 break; } v = (v + 1) % n; } } return ans; }