计算几何 综合模板
1、基本函数
1.1 Point 定义
const double eps = 1e-8; const double PI = acos(-1.0); int sgn(double x) { if(fabs(x) < eps)return 0; if(x < 0)return -1; else return 1; } struct point { double x,y; point() {} point(double _x,double _y) { x = _x; y = _y; } point operator -(const point &b)const { return point(x - b.x,y - b.y); } //叉积 double operator ^(const point &b)const { return x*b.y - y*b.x; } //点积 double operator *(const point &b)const { return x*b.x + y*b.y; } //绕原点旋转角度B（弧度值），后x,y的变化 void transXY(double B) { double tx = x,ty = y; x = tx*cos(B) - ty*sin(B); y = tx*sin(B) + ty*cos(B); } };
1.2 Line 定义
struct Line { point s,e; Line() {} Line(point _s,point _e) { s = _s; e = _e; } //两直线相交求交点 //第一个值为0表示直线重合，为1表示平行，为0表示相交,为2是相交 //只有第一个值为2时，交点才有意义 pair<int,point> operator &(const Line &b)const { point res = s; if(sgn((s-e)^(b.s-b.e)) == 0) { if(sgn((s-b.e)^(b.s-b.e)) == 0) return make_pair(0,res);//重合 else return make_pair(1,res);//平行 } double t = ((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e)); res.x += (e.x-s.x)*t; res.y += (e.y-s.y)*t; return make_pair(2,res); } };
1.4 判断：线段相交
bool inter(Line l1,Line l2) { return max(l1.s.x,l1.e.x) >= min(l2.s.x,l2.e.x) && max(l2.s.x,l2.e.x) >= min(l1.s.x,l1.e.x) && max(l1.s.y,l1.e.y) >= min(l2.s.y,l2.e.y) && max(l2.s.y,l2.e.y) >= min(l1.s.y,l1.e.y) && sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0 && sgn((l1.s-l2.e)^(l2.s-l2.e))*sgn((l1.e-l2.e)^(l2.s-l2.e)) <= 0; }
1.5 判断：直线和线段相交
//判断直线和线段相交 bool Seg_inter_line(Line l1,Line l2) //判断直线l1和线段l2是否相交 { return sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0; }
1.6 点到直线距离
//点到直线距离 //返回为result,是点到直线最近的点 Point PointToLine(Point P,Line L) { Point result; double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s)); result.x = L.s.x + (L.e.x-L.s.x)*t; result.y = L.s.y + (L.e.y-L.s.y)*t; return result; }
1.7 点到线段距离
Point NearestPointToLineSeg(Point P,Line L) { Point result; double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s)); if(t >= 0 && t <= 1) { result.x = L.s.x + (L.e.x - L.s.x)*t; result.y = L.s.y + (L.e.y - L.s.y)*t; } else { if(dist(P,L.s) < dist(P,L.e)) result = L.s; else result = L.e; } return result; }
7.1.7 求 两 线 段 间 最 短 距离
double dis_segments(Segment seg1,Segment seg2) { double m1=dis_point_segment(seg1.s,seg2); //这个函数是求点到线段的距离 double m2=dis_point_segment(seg1.e,seg2); double m3=dis_point_segment(seg2.s,seg1); double m4=dis_point_segment(seg2.e,seg1); return min(min(m1,m2),min(m3,m4)); }
1.8 计算多边形面积
//计算多边形面积 //点的编号从0~n-1 double CalcArea(Point p[],int n) { double res = 0; for(int i = 0; i < n; i++) res += (p[i]^p[(i+1)%n])/2; return fabs(res); }
1.9 判断点在线段上
//*判断点在线段上 bool OnSeg(point P,Line L) { return sgn((L.s-P)^(L.e-P)) == 0 && sgn((P.x - L.s.x) * (P.x - L.e.x)) <= 0 && sgn((P.y - L.s.y) * (P.y - L.e.y)) <= 0; }
1.11 判断点在任意多边形内
//*判断点在任意多边形内 //射线法，poly[]的顶点数要大于等于3,点的编号0~n-1 //返回值 //-1:点在凸多边形外 //0:点在凸多边形边界上 //1:点在凸多边形内 int inPoly(point p,point poly[],int n) { int cnt; Line ray,side; cnt = 0; ray.s = p; ray.e.y = p.y; ray.e.x = -100000000000.0;//-INF,注意取值防止越界 for(int i = 0; i < n; i++) { side.s = poly[i]; side.e = poly[(i+1)%n]; if(OnSeg(p,side))return 0; //如果平行轴则不考虑 if(sgn(side.s.y - side.e.y) == 0) continue; if(OnSeg(side.s,ray)) { if(sgn(side.s.y - side.e.y) > 0)cnt++; } else if(OnSeg(side.e,ray)) { if(sgn(side.e.y - side.s.y) > 0)cnt++; } else if(inter(ray,side)) cnt++; } if(cnt % 2 == 1)return 1; else return -1; }
7.2.1 判 断 线 段 是 否 与 矩 形 相 交 （ 包 括 线 段 在 矩 形 内 部 ）
struct Rectangle//这好像是矩形 { Point lt;//lefttop Point rb;//rightbottom }; int segement_rectangle_intersect(Segment l,Rectangle r) { Segment d1,d2;//retangle ’s diagonal d1.s=r.lt; d1.e=r.rb; d2.s.x=d1.e.x; d2.s.y=d1.s.y; d2.e.x=d1.s.x; d2.e.y=d1.e.y; if (l.s.x>=r.lt.x&&l.s.x<=r.rb.x&& l.s.y<=r.lt.y&&l.s.y>=r.rb.y|| l.e.x>=r.lt.x&&l.e.x<=r.rb.x&& l.e.y<=r.lt.y&&l.e.y>=r.rb.y) return 1; if (segment_intersect(l,d1)|| segment_intersect(l,d2)) //这个函数是判断线段是否相交 //segment_intersect(endpoint inclusive) return 1; return 0; }
1.12 判断凸多边形
//判断凸多边形 //允许共线边 //点可以是顺时针给出也可以是逆时针给出 //点的编号1~n-1 bool isconvex(Point poly[],int n) { bool s[3]; memset(s,false,sizeof(s)); for(int i = 0; i < n; i++) { s[sgn( (poly[(i+1)%n]-poly[i])^(poly[(i+2)%n]-poly[i]) )+1] = true; if(s[0] && s[2])return false; } return true; }
7.2.7 判 断 两 凸 多 边 形 是 否 相 交 （graham scan 求 凸 包+ 枚 举 边 、 点)
#include <iostream > #include <cstdio > #include <algorithm > #define N 110 using namespace std; struct Point { int x,y; }; struct Polygon { Point p[N]; int n; }; Point pt[N]; int stack[N]; int cross(Point a,Point b,Point s) { return (a.x-s.x)*(b.y-s.y)-(b.x-s.x)*(a.y-s.y); } int dist(Point a,Point b) { return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y); } int cmp(Point a,Point b) { if (cross(a,b,pt[0]) >0||cross(a,b,pt[0])==0&&dist(a,pt[0])<dist(b,pt[0])) return 1; else return 0; } int on_segment(Point s,Point e,Point o) { if(cross(s,e,o)==0&& o.x>=min(s.x,e.x)&& o.x<=max(s.x,e.x)&& o.y>=min(s.y,e.y)&& o.y<=max(s.y,e.y)) return 1; else return 0; } int graham_scan(int n) { int i,top,t; if (n<=1) { stack[0]=0; return n; } for (t=0,i=1; i<n; i++) if (pt[i].y<pt[t].y||pt[i].y==pt[t].y&&pt[i].x<pt[t].x) t=i; swap(pt[0],pt[t]); sort(pt+1,pt+n,cmp); top=2; for (i=0; i<2; i++) stack[i]=i; for (i=2; i<n; i++) { while (top >1&&cross(pt[stack[top -1]],pt[i],pt[stack[top -2]]) <=0) top--; stack[top++]=i; } return top; } int segment_intersect(Point s1,Point e1,Point s2,Point e2) { if (max(s1.x,e1.x)>=min(s2.x,e2.x)&& max(s2.x,e2.x)>=min(s1.x,e1.x)&& max(s1.y,e1.y)>=min(s2.y,e2.y)&& max(s2.y,e2.y)>=min(s1.y,e2.y)&& (double)cross(s2,e1,s1)*(double)cross(e2,e1,s1)<=0&& (double)cross(s1,e2,s2)*(double)cross(e1,e2,s2)<=0) return 1; else return 0; } int point_inside(Point o,Polygon pln) { int i,a1=0,a2=0,n=pln.n; pln.p[n]=pln.p[0]; for (i=0; i<n; i++) a1+=abs(cross(pln.p[i],pln.p[i+1],o)); for (i=1; i<n; i++) a2+=abs(cross(pln.p[i],pln.p[i+1],pln.p[0])); if (a1==a2) return 1; else return 0; } int convex_polygon_intersect(Polygon pln1,Polygon pln2) { int i,j; pln1.p[pln1.n]=pln1.p[0]; pln2.p[pln2.n]=pln2.p[0]; for (i=0; i<pln1.n; i++) for (j=0; j<pln2.n; j++) if (segment_intersect(pln1.p[i],pln1.p[i+1],pln2.p[j],pln2.p[j+1])) return 1; if (point_inside(pln1.p[0],pln2)||point_inside(pln2.p[0],pln1)) return 1; return 0; } int main() { int n,m,i,vertexnum,ans; Polygon pln1,pln2; while (cin>>n>>m,n||m) { for (i=0; i<n; i++) cin>>pt[i].x>>pt[i].y; vertexnum=graham_scan(n); pln1.n=vertexnum; for (i=0; i<vertexnum; i++) pln1.p[i]=pt[stack[i]]; for (i=0; i<m; i++) cin>>pt[i].x>>pt[i].y; vertexnum=graham_scan(m); pln2.n=vertexnum; for (i=0; i<vertexnum; i++) pln2.p[i]=pt[stack[i]]; if (pln1.n==1&&pln2.n==1) ans=1; else if (pln1.n==1&&pln2.n==2) { if (on_segment(pln2.p[0],pln2.p[1],pln1.p[0])) ans=0; else ans=1; } else if (pln1.n==2&&pln2.n==1) { if (on_segment(pln1.p[0],pln1.p[1],pln2.p[0])) ans=0; else ans=1; } else if (pln1.n==2&&pln2.n==2) { if (segment_intersect(pln1.p[0],pln1.p[1],pln2.p[0],pln2.p[1])) ans=0; else ans=1; } else if (pln1.n==1) { if (point_inside(pln1.p[0],pln2)) ans=0; else ans=1; } else if (pln2.n==1) { if (point_inside(pln2.p[0],pln1)) ans=0; else ans=1; } else { if (convex_polygon_intersect(pln1,pln2)==0) ans=1; else ans=0; } if (ans==0) cout <<"NO"<<endl; else cout <<"YES"<<endl; } return 0; }
7.3.3 两 个 简 单 多 边 形 求 面 积 并 、交
#include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #define eps 1e-8 #define N 550 using namespace std; struct Point { double x,y; Point () {} Point (double a,double b):x(a),y(b) {} Point operator - (const Point a) const { return Point(x-a.x,y-a.y); } }; Point zero(0,0); int dlcmp(double x) { return x<-eps?-1:x>eps; } double cross(Point v1,Point v2) { return v1.x*v2.y-v2.x*v1.y; } struct Polygon { Point p[N]; int n; Polygon():n(0) {} void clear() { n=0; } void add(Point a) { p[n++]=a; } double area() { double res=0; for (int i=1; i<n-1; i++) res+=cross(p[i]-p[0],p[i+1]-p[0]); return fabs(res/2); } }; Polygon A,B,rec; Point line_intersect(Point s1,Point e1,Point s2,Point e2) //两直线交点 { Point v1,v2,res; double r; v1=s1-e1; v2=s2-e2; r=((s1.x-s2.x)*v2.y-(s1.y-s2.y)*v2.x)/(v1.x*v2.y-v1.y*v2.x); res.x=-r*v1.x+s1.x; res.y=-r*v1.y+s1.y; return res; } void cut(Point s,Point e) //半平面交 { int i,j,d1,d2; Polygon ker; for (i=0; i<rec.n; i++) { j=(i+1)%rec.n; d1=dlcmp(cross(e-s,rec.p[i]-s)); d2=dlcmp(cross(e-s,rec.p[j]-s)); if (d1>=0) ker.add(rec.p[i]); if (d1*d2<0) ker.add(line_intersect(s,e,rec.p[i],rec.p[j])); } rec=ker; } double calc(Point p1,Point p2,Point q1,Point q2) { int dp=dlcmp(cross(p1,p2)),dq=dlcmp(cross(q1,q2)); int sgn=dp*dq; if (sgn==0) return 0; rec.clear(); rec.add(zero); rec.add(p1); rec.add(p2); if (dp<0) swap(rec.p[1],rec.p[2]); if (dq>0) { cut(zero,q1); cut(q1,q2); cut(q2,zero); } else { cut(zero,q2); cut(q2,q1); cut(q1,zero); } return sgn*rec.area(); } double solve() { double res=A.area()+B.area(); double sum=0; int i,j; //对两个多边形三角剖分，分别求两个三角形的面积交 for (i=0; i<A.n; i++) for (j=0; j<B.n; j++) sum+=calc(A.p[i],A.p[(i+1)%A.n],B.p[j],B.p[(j+1)%B.n]); res-=fabs(sum); //fabs(sum)为两个多边形的面积交 return res; //面积并 } //hdu3060(题目数据有误) int main() { int n,m,i; Point pt; double ans; while (scanf("%d%d",&n,&m)!=EOF) { A.clear(); B.clear(); for (i=0; i<n; i++) { scanf("%lf%lf",&pt.x,&pt.y); A.add(pt); } for (i=0; i<m; i++) { scanf("%lf%lf",&pt.x,&pt.y); B.add(pt); } ans=solve(); printf("%.2f\n",ans+eps); } return 0; }
简单的极角排序（以第一个点为基准点）
const int maxn=55; point List[maxn]; double dist(point p0,point p1) { return (double)sqrt((p1.x-p0.x)*(p1.x-p0.x)+(p1.y-p0.y)*(p1.y-p0.y)); } bool _cmp(point p1,point p2) { double tmp = (p1-List[0])^(p2-List[0]); if(sgn(tmp) > 0)return true; else if(sgn(tmp) == 0 && sgn(dist(p1,List[0]) - dist(p2,List[0])) <= 0) return true; else return false; } sort(List+1,List+n,_cmp);
2、凸包（ps：这个算法前面的的排序是先找到最左下角的一个点，算法复杂度O（nlogn）
/* * 求凸包，Graham算法 * 点的编号0~n-1 * 返回凸包结果Stack[0~top-1]为凸包的编号 */ const int MAXN = 1010; point List[MAXN]; int Stack[MAXN],top; //相对于List[0]的极角排序 bool _cmp(point p1,point p2) { double tmp = (p1-List[0])^(p2-List[0]); if(sgn(tmp) > 0)return true; else if(sgn(tmp) == 0 && sgn(dist(p1,List[0]) - dist(p2,List[0])) <= 0) return true; else return false; } void Graham(int n) { point p0; int k = 0; p0 = List[0]; //找最下边的一个点 for(int i = 1; i < n; i++) { if( (p0.y > List[i].y) || (p0.y == List[i].y && p0.x > List[i].x) ) { p0 = List[i]; k = i; } } swap(List[k],List[0]); sort(List+1,List+n,_cmp); if(n == 1) { top = 1; Stack[0] = 0; return; } if(n == 2) { top = 2; Stack[0] = 0; Stack[1] = 1; return ; } Stack[0] = 0; Stack[1] = 1; top = 2; for(int i = 2; i < n; i++) { while(top > 1 && sgn((List[Stack[top-1]]-List[Stack[top-2]])^(List[i]-List[Stack[top-2]])) <= 0)top--; Stack[top++] = i; } }
3、平面最近点对（HDU 1007）
#include <stdio.h> #include <string.h> #include <algorithm> #include <iostream> #include <math.h> using namespace std; const double eps = 1e-6; const int MAXN = 100010; const double INF = 1e20; struct Point { double x,y; }; double dist(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } Point p[MAXN]; Point tmpt[MAXN]; bool cmpxy(Point a,Point b) { if(a.x != b.x)return a.x < b.x; else return a.y < b.y; } bool cmpy(Point a,Point b) { return a.y < b.y; } double Closest_Pair(int left,int right) { double d = INF; if(left == right)return d; if(left + 1 == right) return dist(p[left],p[right]); int mid = (left+right)/2; double d1 = Closest_Pair(left,mid); double d2 = Closest_Pair(mid+1,right); d = min(d1,d2); int k = 0; for(int i = left; i <= right; i++) { if(fabs(p[mid].x - p[i].x) <= d) tmpt[k++] = p[i]; } sort(tmpt,tmpt+k,cmpy); for(int i = 0; i <k; i++) { for(int j = i+1; j < k && tmpt[j].y - tmpt[i].y < d; j++) { d = min(d,dist(tmpt[i],tmpt[j])); } } return d; } int main() { int n; while(scanf("%d",&n)==1 && n) { for(int i = 0; i < n; i++) scanf("%lf%lf",&p[i].x,&p[i].y); sort(p,p+n,cmpxy); printf("%.2lf\n",Closest_Pair(0,n-1)/2); } return 0; }
4、旋转卡壳
4.1 求解平面最远点对（POJ 2187 Beauty Contest）
struct Point { int x,y; Point(int _x = 0,int _y = 0) { x = _x; y = _y; } Point operator -(const Point &b)const { return Point(x - b.x, y - b.y); } int operator ^(const Point &b)const { return x*b.y - y*b.x; } int operator *(const Point &b)const { return x*b.x + y*b.y; } void input() { scanf("%d%d",&x,&y); } }; //距离的平方 int dist2(Point a,Point b) { return (a-b)*(a-b); } //******二维凸包，int*********** const int MAXN = 50010; Point list[MAXN]; int Stack[MAXN],top; bool _cmp(Point p1,Point p2) { int tmp = (p1-list[0])^(p2-list[0]); if(tmp > 0)return true; else if(tmp == 0 && dist2(p1,list[0]) <= dist2(p2,list[0])) return true; else return false; } void Graham(int n) { Point p0; int k = 0; p0 = list[0]; for(int i = 1; i < n; i++) if(p0.y > list[i].y || (p0.y == list[i].y && p0.x > list[i].x)) { p0 = list[i]; k = i; } swap(list[k],list[0]); sort(list+1,list+n,_cmp); if(n == 1) { top = 1; Stack[0] = 0; return; } if(n == 2) { top = 2; Stack[0] = 0; Stack[1] = 1; return; } Stack[0] = 0; Stack[1] = 1; top = 2; for(int i = 2; i < n; i++) { while(top > 1 && ((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <= 0) top--; Stack[top++] = i; } } //旋转卡壳，求两点间距离平方的最大值 int rotating_calipers(Point p[],int n) { int ans = 0; Point v; int cur = 1; for(int i = 0; i < n; i++) { v = p[i]-p[(i+1)%n]; while((v^(p[(cur+1)%n]-p[cur])) < 0) cur = (cur+1)%n; ans = max(ans,max(dist2(p[i],p[cur]),dist2(p[(i+1)%n],p[(cur+1)%n]))); } return ans; } Point p[MAXN]; int main() { int n; while(scanf("%d",&n) == 1) { for(int i = 0; i < n; i++)list[i].input(); Graham(n); for(int i = 0; i < top; i++)p[i] = list[Stack[i]]; printf("%d\n",rotating_calipers(p,top)); } return 0; }
4.2 求解平面点集最大三角形
//旋转卡壳计算平面点集最大三角形面积 int rotating_calipers(Point p[],int n) { int ans = 0; Point v; for(int i = 0; i < n; i++) { int j = (i+1)%n; int k = (j+1)%n; while(j != i && k != i) { ans = max(ans,abs((p[j]-p[i])^(p[k]-p[i]))); while( ((p[i]-p[j])^(p[(k+1)%n]-p[k])) < 0 ) k = (k+1)%n; j = (j+1)%n; } } return ans; } Point p[MAXN]; int main() { int n; while(scanf("%d",&n) == 1) { if(n == -1)break; for(int i = 0; i < n; i++)list[i].input(); Graham(n); for(int i = 0; i < top; i++)p[i] = list[Stack[i]]; printf("%.2f\n",(double)rotating_calipers(p,top)/2); } return 0; }
4.3 求解两凸包最小距离（POJ 3608）
const double eps = 1e-8; int sgn(double x) { if(fabs(x) < eps)return 0; if(x < 0)return -1; else return 1; } struct Point { double x,y; Point(double _x = 0,double _y = 0) { x = _x; y = _y; } Point operator -(const Point &b)const { return Point(x - b.x, y - b.y); } double operator ^(const Point &b)const { return x*b.y - y*b.x; } double operator *(const Point &b)const { return x*b.x + y*b.y; } void input() { scanf("%lf%lf",&x,&y); } }; struct Line { Point s,e; Line() {} Line(Point _s,Point _e) { s = _s; e = _e; } }; //两点间距离 double dist(Point a,Point b) { return sqrt((a-b)*(a-b)); } //点到线段的距离，返回点到线段最近的点 Point NearestPointToLineSeg(Point P,Line L) { Point result; double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s)); if(t >=0 && t <= 1) { result.x = L.s.x + (L.e.x - L.s.x)*t; result.y = L.s.y + (L.e.y - L.s.y)*t; } else { if(dist(P,L.s) < dist(P,L.e)) result = L.s; else result = L.e; } return result; } /* * 求凸包，Graham算法 * 点的编号0~n-1 * 返回凸包结果Stack[0~top-1]为凸包的编号 */const int MAXN = 10010; Point list[MAXN]; int Stack[MAXN],top; //相对于list[0]的极角排序 bool _cmp(Point p1,Point p2) { double tmp = (p1-list[0])^(p2-list[0]); if(sgn(tmp) > 0)return true; else if(sgn(tmp) == 0 && sgn(dist(p1,list[0]) - dist(p2,list[0])) <= 0) return true; else return false; } void Graham(int n) { Point p0; int k = 0; p0 = list[0]; //找最下边的一个点 for(int i = 1; i < n; i++) { if( (p0.y > list[i].y) || (p0.y == list[i].y && p0.x > list[i].x) ) { p0 = list[i]; k = i; } } swap(list[k],list[0]); sort(list+1,list+n,_cmp); if(n == 1) { top = 1; Stack[0] = 0; return; } if(n == 2) { top = 2; Stack[0] = 0; Stack[1] = 1; return ; } Stack[0] = 0; Stack[1] = 1; top = 2; for(int i = 2; i < n; i++) { while(top > 1 && sgn((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <= 0) top--; Stack[top++] = i; } } //点p0到线段p1p2的距离 double pointtoseg(Point p0,Point p1,Point p2) { return dist(p0,NearestPointToLineSeg(p0,Line(p1,p2))); }//平行线段p0p1和p2p3的距离 double dispallseg(Point p0,Point p1,Point p2,Point p3) { double ans1 = min(pointtoseg(p0,p2,p3),pointtoseg(p1,p2,p3)); double ans2 = min(pointtoseg(p2,p0,p1),pointtoseg(p3,p0,p1)); return min(ans1,ans2); } //得到向量a1a2和b1b2的位置关系 double Get_angle(Point a1,Point a2,Point b1,Point b2) { return (a2-a1)^(b1-b2); } double rotating_calipers(Point p[],int np,Point q[],int nq) { int sp = 0, sq = 0; for(int i = 0; i < np; i++) if(sgn(p[i].y - p[sp].y) < 0) sp = i; for(int i = 0; i < nq; i++) if(sgn(q[i].y - q[sq].y) > 0) sq = i; double tmp; double ans = dist(p[sp],q[sq]); for(int i = 0; i < np; i++) { while(sgn(tmp = Get_angle(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])) < 0) sq = (sq+1)%nq; if(sgn(tmp) == 0) ans = min(ans,dispallseg(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])); else ans = min(ans,pointtoseg(q[sq],p[sp],p[(sp+1)%np])); sp = (sp+1)%np; } return ans; } double solve(Point p[],int n,Point q[],int m) { return min(rotating_calipers(p,n,q,m),rotating_calipers(q,m,p,n)); } Point p[MAXN],q[MAXN]; int main() { int n,m; while(scanf("%d%d",&n,&m) == 2) { if(n == 0 && m == 0)break; for(int i = 0; i < n; i++) list[i].input(); Graham(n); for(int i = 0; i < top; i++) p[i] = list[i]; n = top; for(int i = 0; i < m; i++) list[i].input(); Graham(m); for(int i = 0; i < top; i++) q[i] = list[i]; m = top; printf("%.4f\n",solve(p,n,q,m)); } return 0; }
5、半平面交
5.1 半平面交模板(from UESTC)
const double eps = 1e-8; const double PI = acos(-1.0); int sgn(double x) { if(fabs(x) < eps) return 0; if(x < 0) return -1; else return 1; } struct point { double x,y; point() {} point(double _x,double _y) { x = _x; y = _y; } point operator -(const point &b)const { return point(x - b.x, y - b.y); } double operator ^(const point &b)const { return x*b.y - y*b.x; } double operator *(const point &b)const { return x*b.x + y*b.y; } }; struct Line { point s,e; double k; Line() {} Line(point _s,point _e) { s = _s; e = _e; k = atan2(e.y - s.y,e.x - s.x); } point operator &(const Line &b)const { point res = s; double t = ((s - b.s)^(b.s - b.e))/((s - e)^(b.s - b.e)); res.x += (e.x - s.x)*t; res.y += (e.y - s.y)*t; return res; } }; //半平面交，直线的左边代表有效区域 //这个好像和给出点的顺序有关 bool HPIcmp(Line a,Line b) { if(fabs(a.k - b.k) > eps)return a.k < b.k; return ((a.s - b.s)^(b.e - b.s)) < 0; } Line Q[110]; //第一个位代表半平面交的直线，第二个参数代表直线条数，第三个参数是相交以后把 //所得点压栈，第四个参数是栈有多少个元素 void HPI(Line line[], int n, point res[], int &resn) { int tot = n; sort(line,line+n,HPIcmp); tot = 1; for(int i = 1; i < n; i++) if(fabs(line[i].k - line[i-1].k) > eps) line[tot++] = line[i]; int head = 0, tail = 1; Q[0] = line[0]; Q[1] = line[1]; resn = 0; for(int i = 2; i < tot; i++) { if(fabs((Q[tail].e-Q[tail].s)^(Q[tail-1].e-Q[tail-1].s)) < eps || fabs((Q[head].e-Q[head].s)^(Q[head+1].e-Q[head+1].s)) < eps) return; while(head < tail && (((Q[tail]&Q[tail-1]) - line[i].s)^(line[i].e-line[i].s)) > eps) tail--; while(head < tail && (((Q[head]&Q[head+1]) - line[i].s)^(line[i].e-line[i].s)) > eps) head++; Q[++tail] = line[i]; } while(head < tail && (((Q[tail]&Q[tail-1]) - Q[head].s)^(Q[head].e-Q[head].s)) > eps) tail--; while(head < tail && (((Q[head]&Q[head-1]) - Q[tail].s)^(Q[tail].e-Q[tail].e)) > eps) head++; if(tail <= head + 1)return; for(int i = head; i < tail; i++) res[resn++] = Q[i]&Q[i+1]; if(head < tail - 1) res[resn++] = Q[head]&Q[tail]; }
5.2 普通半平面交写法
POJ 1750 const double eps = 1e-18; int sgn(double x) { if(fabs(x) < eps)return 0; if(x < 0)return -1; else return 1; } struct Point { double x,y; Point() {} Point(double _x,double _y) { x = _x; y = _y; } Point operator -(const Point &b)const { return Point(x - b.x, y - b.y); } double operator ^(const Point &b)const { return x*b.y - y*b.x; } double operator *(const Point &b)const { return x*b.x + y*b.y; } }; //计算多边形面积 double CalcArea(Point p[],int n) { double res = 0; for(int i = 0; i < n; i++) res += (p[i]^p[(i+1)%n]); return fabs(res/2); } //通过两点，确定直线方程 void Get_equation(Point p1,Point p2,double &a,double &b,double &c) { a = p2.y - p1.y; b = p1.x - p2.x; c = p2.x*p1.y - p1.x*p2.y; } //求交点 Point Intersection(Point p1,Point p2,double a,double b,double c) { double u = fabs(a*p1.x + b*p1.y + c); double v = fabs(a*p2.x + b*p2.y + c); Point t; t.x = (p1.x*v + p2.x*u)/(u+v); t.y = (p1.y*v + p2.y*u)/(u+v); return t; } Point tp[110]; void Cut(double a,double b,double c,Point p[],int &cnt) { int tmp = 0; for(int i = 1; i <= cnt; i++) { //当前点在左侧，逆时针的点 if(a*p[i].x + b*p[i].y + c < eps)tp[++tmp] = p[i]; else { if(a*p[i-1].x + b*p[i-1].y + c < -eps) tp[++tmp] = Intersection(p[i-1],p[i],a,b,c); if(a*p[i+1].x + b*p[i+1].y + c < -eps) tp[++tmp] = Intersection(p[i],p[i+1],a,b,c); } } for(int i = 1; i <= tmp; i++) p[i] = tp[i]; p[0] = p[tmp]; p[tmp+1] = p[1]; cnt = tmp; } double V[110],U[110],W[110]; int n; const double INF = 100000000000.0; Point p[110]; bool solve(int id) { p[1] = Point(0,0); p[2] = Point(INF,0); p[3] = Point(INF,INF); p[4] = Point(0,INF); p[0] = p[4]; p[5] = p[1]; int cnt = 4; for(int i = 0; i < n; i++) if(i != id) { double a = (V[i] - V[id])/(V[i]*V[id]); double b = (U[i] - U[id])/(U[i]*U[id]); double c = (W[i] - W[id])/(W[i]*W[id]); if(sgn(a) == 0 && sgn(b) == 0) { if(sgn(c) >= 0)return false; else continue; } Cut(a,b,c,p,cnt); } if(sgn(CalcArea(p,cnt)) == 0)return false; else return true; } int main() { while(scanf("%d",&n) == 1) { for(int i = 0; i < n; i++) scanf("%lf%lf%lf",&V[i],&U[i],&W[i]); for(int i = 0; i < n; i++) { if(solve(i))printf("Yes\n"); else printf("No\n"); } } return 0; }
6、三点求圆心坐标（三角形外心）
//过三点求圆心坐标 Point waixin(Point a,Point b,Point c) { double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2; double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2; double d = a1*b2 - a2*b1; return Point(a.x + (c1*b2 - c2*b1)/d, a.y + (a1*c2 -a2*c1)/d); }
7.2.3 求 简 单 多 边 形 重 心
Point get_center(Point pt[],int n) { double sum,area; Point res(0,0),o(0,0); int i; sum=0; for (i=0; i<n; i++) { area=cross(pt[i]-o,pt[(i+1)%n]-o); res=res+(pt[i]+pt[(i+1)%n])/3*area; sum+=area; } res=res/sum; return res; }
7、求两圆相交的面积
//两个圆的公共部分面积 double Area_of_overlap(Point c1,double r1,Point c2,double r2) { double d = dist(c1,c2); if(r1 + r2 < d + eps)return 0; if(d < fabs(r1 - r2) + eps) { double r = min(r1,r2); return PI*r*r; } double x = (d*d + r1*r1 - r2*r2)/(2*d); double t1 = acos(x / r1); double t2 = acos((d - x)/r2); return r1*r1*t1 + r2*r2*t2 - d*r1*sin(t1); }
8、Pick 公式
顶点坐标均是整点的简单多边形：面积=内部格点数目+边上格点数目/2-1
7.4 圆
7.4.1 点类
struct Point { double x,y; Point() {} Point(double a,double b):x(a),y(b) {} Point operator + (const Point a) const { return Point(x+a.x,y+a.y); } Point operator - (const Point a) const { return Point(x-a.x,y-a.y); } Point operator * (const double a) const { return Point(x*a,y*a); } Point operator / (const double a) const { return Point(x/a,y/a); } bool operator < (const Point a) const { if (dlcmp(x-a.x)==0) return dlcmp(x-a.y)<0; else return dlcmp(x-a.x)<0; } bool operator == (const Point a) const { return !dlcmp(x-a.x)&&!dlcmp(y-a.y); } //向量长度定为d Point trunc(double d) { double dis(Point,Point); double len=dis(*this,Point(0,0)); return Point(x*d/len,y*d/len); } //坐标逆时针旋转a度 Point rotate(double a) { return Point(x*cos(a)-y*sin(a),y*cos(a)+x*sin(a)); } }; double dis(Point a,Point b) { return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)); } double cross(Point a,Point b,Point s) { double x1=a.x-s.x,y1=a.y-s.y; double x2=b.x-s.x,y2=b.y-s.y; return x1*y2-x2*y1; } double cross(Point a,Point b) { return a.x*b.y-b.x*a.y; } double dot(Point a,Point b,Point s) { double x1=a.x-s.x,y1=a.y-s.y; double x2=b.x-s.x,y2=b.y-s.y; return x1*x2+y1*y2; } double dot(Point a,Point b) { return a.x*b.x+a.y*b.y; }
7.4.2 圆 类
struct Circle { Point o; double r; Circle() {} Circle(Point a,double l):o(a),r(l) {} double area() { return sqr(r)*PI; } }; //判断圆a是否含于圆b int inner_circle(Circle a,Circle b) { if (dlcmp(a.r-b.r)>0) return 0; return dlcmp(dis(a.o,b.o)+a.r-b.r)<=0; } //以base点为基点，极角排序，排序前base需赋初值 Point base; int cmp(const Point a,const Point b) { return atan2(a.y-base.y,a.x-base.x)<atan2(b.y-base.y,b.x-base.x); } //向量a,b的夹角 double vec_angle(Point a,Point b) { double tmp=dot(a,b)/(dis(a,Point(0,0))*dis(b,Point(0,0))); if (dlcmp(tmp -1) >=0) tmp=1; if (dlcmp(tmp+1) <=0) tmp=-1; return acos(tmp); } //计算由a到b逆时针方向的弓形面积 double arc_area(Point a,Point b,Circle c) { double theta=vec_angle(a-c.o,b-c.o); double sf=sqr(c.r)*theta/2.0; double st=sqr(c.r)*sin(theta)/2.0; if (dlcmp(cross(a,b,c.o))>0) return sf-st; else return c.area()-sf+st; } double arc_area(double th,double r) { return 0.5*sqr(r)*(th-sin(th)); }
7.4.7 最 小 圆 覆盖（平面上n个点，求一个半径最小的圆，能够覆盖所有的点。）
#include <Bits/stdc++.h> using namespace std; #define eps 1e-8 #define MAX_P 2000 struct Point { double x,y; Point operator -(Point &a) { Point t; t.x=x-a.x; t.y=y-a.y; return t; } }; struct Circle { double r; Point center; }; struct Triangle { Point t[3]; }; Point pt[MAX_P]; //点集 Circle c; //最小圆 double distance(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } double cross(Point a,Point b) { return a.x*b.y-b.x*a.y; } double triangle_area(Triangle tri) //三角形距离 { Point v1=tri.t[1]-tri.t[0]; Point v2=tri.t[2]-tri.t[0]; return fabs(cross(v1,v2))/2; } Circle circumcircle_triangle(Triangle tri) //三角形外接圆 { Circle res; double a,b,c,c1,c2; double xA,yA,xB,yB,xC,yC; a=distance(tri.t[0],tri.t[1]); b=distance(tri.t[1],tri.t[2]); c=distance(tri.t[2],tri.t[0]); res.r=a*b*c/triangle_area(tri)/4; xA=tri.t[0].x; yA=tri.t[0].y; xB=tri.t[1].x; yB=tri.t[1].y; xC=tri.t[2].x; yC=tri.t[2].y; c1=(xA*xA+yA*yA-xB*xB-yB*yB)/2; c2 = (xA*xA+yA*yA-xC*xC-yC*yC)/2; res.center.x=(c1*(yA-yC)-c2*(yA-yB))/ ((xA-xB)*(yA-yC)-(xA-xC)*(yA-yB)); res.center.y = (c1*(xA-xC)-c2*(xA-xB))/ ((yA-yB)*(xA-xC)-(yA-yC)*(xA-xB)); return res; } Circle mincircle_triangle(int trinum,Triangle tri) { Circle res; if (trinum==0) res.r=-2; else if (trinum==1) { res.center=tri.t[0]; res.r=0; } else if (trinum==2) { res.center.x=(tri.t[0].x+tri.t[1].x)/2; res.center.y=(tri.t[0].y+tri.t[1].y)/2; res.r=distance(tri.t[0],tri.t[1])/2; } else if (trinum==3) res=circumcircle_triangle(tri); return res; } void mincircle_pointset(int m,int trinum,Triangle tri) //求点集的最小覆盖圆 { int i,j; Point tmp; c=mincircle_triangle(trinum,tri); if (trinum==3) return; for (i=0; i<m; i++) if (distance(pt[i],c.center)>c.r) { tri.t[trinum]=pt[i]; mincircle_pointset(i,trinum+1,tri); tmp=pt[i]; for (j=i; j>=1; j--) pt[j]=pt[j-1]; pt[0]=tmp; } } int main() { int n,i,f1,f2; Triangle tri; while (scanf("%d%d%d",&f1,&f2,&n)!=EOF) { for (i=0; i<n; i++) scanf("%lf%lf",&pt[i].x,&pt[i].y); mincircle_pointset(n,0,tri); printf("%lf %lf %lf\n",c.center.x,c.center.y,c.r); } return 0; }
7.4.8 单 位 圆 覆盖（一个单位圆最多能覆盖平面上多少点）
#include <math.h> #include <bits/stdc++.h> using namespace std; #define eps 1e-8 #define MAX_P 505 const double r=1.0;//单位圆半径 struct Point { double x,y; }; Point pt[MAX_P]; double distance(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //sqrt函数速度较慢，应尽量避免出现，此处可优化为距离的平方和的形式 } Point find_center(Point a,Point b) { Point v,mid,center; double d,s,ang; v.x=a.x-b.x; v.y=a.y-b.y; mid.x=(a.x+b.x)/2; mid.y=(a.y+b.y)/2; d=distance(a,mid); s=sqrt(r*r-d*d); //优化为s=sqrt(r*r-d); if (fabs(v.y)<eps) { center.x=mid.x; center.y=mid.y+s; } else { ang=atan(-v.x/v.y); center.x=mid.x+s*cos(ang); center.y=mid.y+s*sin(ang); } return center; } int main() { int n,i,j,k,ans,cnt; double tmp; Point center; while (~scanf("%d",&n),n) { for (i=0; i<n; i++) scanf("%lf%lf",&pt[i].x,&pt[i].y); ans=1; for (i=0; i<n; i++) for (j=i+1; j<n; j++) { if (distance(pt[i],pt[j])>2*r) //优化为distance(pt[i],pt[j])>2*2*r*r continue; cnt=0; center=find_center(pt[i],pt[j]); for (k=0; k<n; k++) if (distance(pt[k],center)<=r+eps) cnt++; if (ans<cnt) ans=cnt; } printf("%d\n",ans); } return 0; }
7.5 模 拟 退 火
7.5.1 求 多 边 形 费 马点
#include <iostream > #include <cstdio > #include <cmath > #define eps 1e-6 #define N 105 using namespace std; struct Point { double x,y; }; double point_dis(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } double sum_dis(Point pt[],int n,Point o) { double res=0; int i; for (i=0; i<n; i++) res+=point_dis(pt[i],o); return res; } double polygon_Fermatpoint(Point pln[],int n) { Point cp,np,tmp; double min,step,d; int flag; cp=pln[0]; //cp保存当前更新后最优的费马点 min=sum_dis(pln,n,cp); step=10000; //选取坐标范围的最大值 while (step>eps) { flag=1; while (flag) { flag=0; np=cp; tmp=cp,tmp.x+=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp, flag=1; tmp=cp,tmp.x-=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; tmp=cp,tmp.y+=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; tmp=cp,tmp.y-=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; cp=np; } step*=0.98; //系数根据精度要求修改 } return min; } int main() { int n,i; double min; Point pt[N]; cin>>n; for (i=0; i<n; i++) cin>>pt[i].x>>pt[i].y; min=polygon_Fermatpoint(pt,n); printf("%.0f\n",min); return 0; }
7.5.2 最 小 球 覆盖（求能覆盖所有点的最小球的半径。）
#include <iostream > #include <cstdio > #include <cmath > #define oo 1e20 #define eps 1e-10 #define N 105 using namespace std; struct Point { double x,y,z; }; double dis(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z)); } int max_dis(Point pt[],int n,Point o) { int i,res; double max,tmp; max=0; res=0; for (i=0; i<n; i++) { tmp=dis(pt[i],o); if (max<tmp) { max=tmp; res=i; } } return res; } int main() { Point pt[N],o; int n,i,t; double dx,dy,dz,step,r,tmp; cin>>n; for (i=0; i<n; i++) cin>>pt[i].x>>pt[i].y>>pt[i].z; step=10000; //step选取最大的坐标范围 r=oo; if (n==1) { o.x=pt[0].x; o.y=pt[0].y; o.z=pt[0].z; } else { o.x=o.y=o.z=0; while (step>eps) { t=max_dis(pt,n,o); tmp=dis(pt[t],o); if (r>tmp) r=tmp; dx=(pt[t].x-o.x)/tmp; dy=(pt[t].y-o.y)/tmp; dz=(pt[t].z-o.z)/tmp; o.x+=step*dx; o.y+=step*dy; o.z+=step*dz; step*=0.9993; //系数的选取根据具体精度调整 } } printf("%.6f %.6f %.6f\n",o.x,o.y,o.z); return 0; }
8.求四面体体积//hdu 1411
double solve(double ab,double ac,double ad,double bc,double bd,double cd) { double x=(ad*ad+bd*bd-ab*ab)/(2*ad*bd); double y=(bd*bd+cd*cd-bc*bc)/(2*bd*cd); double z=(ad*ad+cd*cd-ac*ac)/(2*ad*cd); double v=ad*bd*cd*sqrt(1.0+2.0*x*y*z-x*x-y*y-z*z)/6.0; return v; }
