计算几何模板
既然咱负责计算几何，就多刷一些题吧ww
模板改编自刘汝佳
点，线基础部分：
const double pi=acos(-1.0); int dcmp(double x) {if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } struct Vector { double x, y; Vector (double x=0, double y=0) :x(x),y(y) {} Vector operator + (const Vector &B) const { return Vector (x+B.x,y+B.y); } Vector operator - (const Vector &B) const { return Vector(x - B.x, y - B.y); } Vector operator * (const double &p) const { return Vector(x*p, y*p); } Vector operator / (const double &p) const { return Vector(x/p, y/p); } double operator * (const Vector &B) const { return x*B.x + y*B.y;}//点积 double operator ^ (const Vector &B) const { return x*B.y - y*B.x;}//叉积 bool operator < (const Vector &b) const { return x < b.x || (x == b.x && y < b.y); } bool operator ==(const Vector &b) const { return dcmp(x-b.x) == 0 && dcmp(y-b.y) == 0; } }; typedef Vector Point; Point Read(){double x, y;scanf("%lf%lf", &x, &y);return Point(x, y);} double Length(Vector A){ return sqrt(A*A); }//向量的模 double Angle(Vector A, Vector B){return acos(A*B / Length(A) / Length(B)); }//向量的夹角，返回值为弧度 double Area2(Point A, Point B, Point C){ return (B-A)^(C-A); }//向量AB叉乘AC的有向面积 Vector VRotate(Vector A, double rad){return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));}//向量A旋转rad弧度 Point PRotate(Point A, Point B, double rad){return A + VRotate(B-A, rad);}//将B点绕A点旋转rad弧度 Vector Normal(Vector A){double l = Length(A);return Vector(-A.y/l, A.x/l);}//求向量A向左旋转90°的单位法向量,调用前确保A不是零向量 Point GetLineIntersection/*求直线交点，调用前要确保两条直线有唯一交点*/(Point P, Vector v, Point Q, Vector w){double t = (w^(P - Q)) / (v^w);return P + v*t;}//在精度要求极高的情况下，可以自定义分数类 double DistanceToLine/*P点到直线AB的距离*/(Point P, Point A, Point B){Vector v1 = B - A, v2 = P - A;return fabs(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(v1*v2)<0) return Length(v2);else if (dcmp(v1*v3)>0) return Length(v3);else return fabs(v1^v2)/Length(v1); } Point GetLineProjection/*点在直线上的射影*/(Point P, Point A, Point B) { Vector v=B-A; return A+v*((v*(P-A))/(v*v)); } bool OnSegment/*判断点是否在线段上(含端点)*/(Point P,Point a1,Point a2) { Vector v1=a1-P,v2=a2-P; if (dcmp(v1^v2)==0 && min(a1.x,a2.x)<=P.x && P.x<=max(a1.x,a2.x) && min(a1.y,a2.y)<=P.y && P.y<=max(a1.y,a2.y)) return true; return false; } bool SegmentInter/*线段相交判定*/(Point a1, Point a2, Point b1, Point b2) { //if (OnSegment(a1,b1,b2) || OnSegment(a2,b1,b2) || OnSegment(b1,a1,a2) || OnSegment(b2,a1,a2)) return 1; //如果只判断线段规范相交（不算交点）,上面那句可以删掉 double c1=(a2-a1)^(b1-a1),c2=(a2-a1)^(b2-a1); double c3=(b2-b1)^(a1-b1),c4=(b2-b1)^(a2-b1); return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0; } bool InTri/*判断点是否在三角形内*/(Point P, Point a,Point b,Point c) { if (dcmp(fabs((c-a)^(c-b))-fabs((P-a)^(P-b))-fabs((P-b)^(P-c))-fabs((P-a)^(P-c)))==0) return true; return false; } double PolygonArea/*求多边形面积，注意凸包P序号从0开始*/(Point *P ,int n) { double ans = 0.0; for(int i=1;i<n-1;i++) ans+=(P[i]-P[0])^(P[i+1]-P[0]); return ans/2; } bool CrossOfSegAndLine/*判断线段是否与直线相交*/(Point a1,Point a2,Point b1,Vector b2) { if (OnSegment(b1,a1,a2) || OnSegment(b1+b2,a1,a2)) return true; return dcmp(b2^(a1-b1))*dcmp(b2^(a2-b1))<0; }
半平面交
struct Line//有向直线 { Point p; Vector v; double ang; Line() { } Line(Point p, Vector v): p(p), v(v) { ang = atan2(v.y, v.x); } Point point(double t) { return p + v*t; } bool operator < (const Line& L) const { if (fabs(ang-L.ang)<eps) { return ((v)^(L.p-p))<0; } return ang < L.ang; } }; bool OnRight(Line L,Point p)//判断是否p点在有向直线右边 { return ((L.v^(p-L.p))<0); } Point GetIntersection (Line a,Line b)//求两个有向直线的交点 { Vector u=a.p-b.p; double t=(b.v^u)/(a.v^b.v); return a.p+a.v*t; } bool IsParallel(Line a,Line b)//判断两条有向直线是否平行 { if (dcmp(a.v ^ b.v)==0) return 1; return 0; } Line q[MAXN]; int HalfPlane(Line *L,int n,Point *poly)//半平面交，输入直线数组L(需要保证逆时针顺序)，n是直线的个数，poly用于输出，返回值为poly中点的个数 { sort(L,L+n); int m,i; for (m=i=1;i<n;i++) { if (dcmp(L[i].ang-L[i-1].ang)!=0) L[m++]=L[i]; } n=m; int l=0,r=1; q[0]=L[0];q[1]=L[1]; for (int i=2;i<n;i++) { //if (IsParallel(q[r],q[r-1]) || IsParallel(q[l],q[l+1])) return 0; while (l < r && OnRight(L[i],GetIntersection(q[r-1],q[r]))) r--; while (l < r && OnRight(L[i],GetIntersection(q[l],q[l+1]))) l++; q[++r]=L[i]; } while (l<r && OnRight(q[l],GetIntersection(q[r-1],q[r]))) r--; while (l<r && OnRight(q[r],GetIntersection(q[l],q[l+1]))) l++; if (r-l<=1 ) return 0; q[r+1]=q[l]; m=0; for (int i=l;i<=r;i++) poly[m++]=GetIntersection(q[i],q[i+1]); return m; }
凸包
int ConvexHull(Point* p,int n,Point* ch) { sort(p,p+n); int m=0; for(int i=0;i<n;++i) { while(m>1&&((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&&((ch[m-1]-ch[m-2])^(p[i]-ch[m-2]))<=0) m--; ch[m++]=p[i]; } if(n>1) m--; return m; }
基础模板补充：
double Cross/*B-A和C-A的叉积*/(Point A, Point B,Point C) { return (B-A)^(C-A); } double dis_pair_seg/*两条线段间的最短距离*/(Point p1, Point p2, Point p3, Point p4) { return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)), min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2))); }
旋转卡壳
double rotating_calipers(Point *ch,int n) { int q=1; double ans=0; ch[n]=ch[0]; for(int p=0; p<n; p++) { while(((ch[q+1]-ch[p+1])^(ch[p]-ch[p+1]))>((ch[q]-ch[p+1])^(ch[p]-ch[p+1]))) q=(q+1)%n; ans=max(ans,max(Length(ch[p]-ch[q]),Length(ch[p+1]-ch[q+1]))); } return ans; }
为什么是大于号而不是大于等于号？因为ch默认是凸多边形
转角法判定点P是否在多边形内部（多边形不一定是凸多边形）
int isPointInPolygon(Point P, Point* Poly, int n)//转角法判定点P是否在多边形内部 { int wn=0; for(int i = 0; i < n; ++i) { if(OnSegment(P, Poly[i], Poly[(i+1)%n])) return -1; //在边界上 int k = dcmp((Poly[(i+1)%n] - Poly[i])^( P - Poly[i])); int d1 = dcmp(Poly[i].y - P.y); int d2 = dcmp(Poly[(i+1)%n].y - P.y); if(k > 0 && d1 <= 0 && d2 > 0) wn++; if(k < 0 && d2 <= 0 && d1 > 0) wn--; } if(wn != 0) return 1; //内部 return 0; //外部 }
关于圆的一些模板
圆的定义：
struct Circle{ Point c; double r; Circle() {} 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); } void read(){ c.read(); scanf("%lf", &r); } };
