计算几何细节梳理&模板 - Flash_Hu - 博客园 (cnblogs.com) 必看,还有一些知识点我漏掉了(细节)
1、二维几何
点和向量 point Vector(直线、线段、射线)
#include<bits/stdc++.h> using namespace std; const double pi=acos(-1.0); const double eps=1e-8; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } int dcmp(double x,double y){ //比较两个浮点数,0为相等,-1为小于,1为大于 if(fabs(x-y)<eps) return 0; else return x<y? -1:1; } //点,点用坐标(x,y)表示 struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} point operator / (double k){return point(x/k,y/k);} bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;} bool operator < (point B){ return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);} }; //两点之间的距离 double Distance(point a,point b){ return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //或者用库函数hypot() return hypot(a.x-b.x,a.y-b.y); }
点积和叉积
点积dot(a,b) 叉积cross(a,b)
//向量:有大小有方向,如果把其中一个点移到源点,那么就可以用点来表示了 //向量的运算:运算符重载,加减乘除 typedef point Vector; //点积和叉积 //点积是A*B=|A||B|cosC C是A,B之间的夹角 //A*B=A.x*B.x+A.y*B.y double dot(Vector A,Vector B) { return A.x*B.x+A.y*B.y; } //应用 //1.判断A,B的夹角是钝角还是锐角 //如果dot(A,B)>0 ,那么夹角为锐角 如果 dot(A,B)<0 ,那么夹角为钝角 如果dot(A,B)=0 ,那么夹角为直角 //2.求向量A的长度 double len(Vector A){ return sqrt(dot(A,A)); //或者求长度的平方 }
double len2(Vector A){
return dot(A,A);
}
//3.求向量A,B的夹角大小 double angle(Vector A,Vector B){ return acos(dot(A,B)/len(A)/len(B)); } //叉积是A*B=|A||B|sinC C表示A旋转到B的夹角 //A*B=A.x*B.y-A.y*B.x 有正负号,几何意义是A和B形成的平行四边形的“有向”面积 //计算过程中A,B是有顺序的,A*B与B*A是相反的 ------小窍门,左边的边向逆时针转到第二条边,就是正的 double cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; } //应用: //1.判断方向关系:如果A*B>0 那么B在A的逆时针方向 // 如果A*B<0 那么B在A的顺时针方向 // 如果A*B=0 那么B与A共线,可能同向,也可能反向 //2。计算两向量构成的平行四边形的有向面积 //3个点A,B,C以A为公共点,得到两个向量B-A,C-A,构成的平行四边形面积为: double area2(Vector A,Vector B,Vector C){ return cross(B-A,C-A); } //3.三个点构成的三角形面积=平行四边形的一半 //4.向量旋转 //使向量(x,y)绕起点逆时针旋转,角度为C,那么旋转后的向量为 //x`=xcosC-ysinC //y`=xsinC+ycosC Vector rotate(Vector A,double rad){ //rad是旋转角度 return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad)); } //如果逆时针旋转90,rotate(a,pi/2),返回vector(-a.y,a.x) //如果顺时针旋转90,rotate(a,-pi/2),返回vector(a.y,-a.x) //如果求单位法向量:逆时针转90然后取单位值 Vector normal(Vector A){ return Vector(-A.y/len(A),A.x/len(A)); } //5.用叉积检查两个向量是不是平行或者重合 bool parallel(Vector A,Vector B){ return sgn(cross(A,B))==0; }
点和线:
1.直线的表示、线段的表示
2、点和直线(位置关系、距离、投影、对称点) 3、点和线段(位置关系、距离)
4、两条直线(位置关系、交点) 5、两条线段(是否相交、交点)
//直线的表示: //1.用直线上的两个点来表示 //2.ax+by+c=0 //3.y=kx+b //4.P=P0+vt 点向式:比较适合计算机处理,表示直线(t无限制)、线段(t=0~1)、射线(t>0) struct Line{ point p1,p2; Line(){} //两个点 Line(point p1,point p2):p1(p1),p2(p2){} //一个点,一个倾斜角angle在0,pi之间 Line(point p,double angle){ p1=p; if(sgn(angle-pi/2)==0) {p2=(p1+point(0,1));} else{p2=(p1+point(1,tan(angle)));} } Line(double a,double b,double c){ if(sgn(a)==0){ p1=point(0,-c/b); p2=point(1,-c/b); } if(sgn(b)==0){ p1=point(-c/a,0); p2=point(-c/a,1); } else{ p1=point(0,-c/b); p2=point(1,((-c-a)/b)); } } }; //点和直线的关系 int point_line_relation(point p,Line v){ int c=sgn(cross(p-v.p1,v.p2-v.p1)); if(c<0) return 1; //1 p在v的左边 if(c>0) return 2; //2 p在v的右边 return 0; //p在v上 } //点和直线的距离 p v(p1,p2) //用叉积求p,p1,p2构成的平行四边形面积,然后除以平行四边形的底边长,就得到了高 double dis_point_line(point p,Line v){ return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2); } //点在直线上的投影 point point_line_proj(point p,Line v){ double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1); return v.p1+(v.p2-v.p1)*k; } //点关于直线的对称点 point point_line_symmetry(point p,Line v){ point q=point_line_proj(p,v); return point(2*q.x-p.x,2*q.y-p.y); } //两条直线的关系 int line_relation(Line v1,Line v2){ if(sgn(cross(v1.p2-v1.p1,v2.p2-v2.p1))==0){ if(point_line_relation(v1.p1,v2)==0) return 1; //重合 else return 0; //平行 } return 2; //相交 } x //两条直线的交点 point cross_point(point a,point b,point c,point d){ double s1=cross(b-a,c-a); double s2=cross(b-a,d-a); //在调用这个函数之前要确保s1 s2不相等 return point(c.x*s2-d.x*s1,c.y*s2-d.y*s1)/(s2-s1); } //线段的表示:直接用直线的结构 typedef Line Segment; //点和线段的关系 //先用叉积看是不是共线,然后点积看p与v的两个端点是不是形成的钝角 bool point_on_seg(point p,Line v){ return sgn(cross(p-v.p1,v.p2-v.p1))==0&&sgn(dot(p-v.p1,p-v.p2))<=0; } //点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值 double dis_point_seg(point p,Segment v){ if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2)); return dis_point_line(point(p,v)); } //判断两个线段是否相交 //如果一条线段的两端在另一条线段的两侧,那么两个端点与另一线段产生的两个叉积正负相反 bool cross_segment(point a,point b,point c,point d){ double c1=cross(b-a,c-a),c2=cross(b-a,d-a); double d1=cross(d-c,a-c),d2=cross(d-c,b-c); return sgn(c1)*sgn(c2)<0&&sgn(d1)*sgn(d2)<0; //1:相交 0:不想交 } //两个线段的交点:先判断两条线段是否相交,如果相交转化为两条直线求交点
多边形
1.点在三角形的内部还是外部 2.求多边形的面积 3.求多边形的重心
//多边形 //1.判断点在多边形的内部:转角法,以P为起点引起一条水平线,检查与多边形每条边的相交情况,统计P穿过这些边的次数 //c=cross(p-j,i-j) u=i.y-p.y v=j.y-p.y int point_in_polygon(point pt,point *p,int n){ for(int i=0;i<n;i++){ if(p[i]==pt) return 3; //3.点在多边形的顶点上 } for(int i=0;i<n;i++){ Line v=Line(p[i],p[(i+1)%n]); if(point_on_seg(pt,v)) return 2; //2.点在多边形的边上 } int num=0; for(int i=0;i<n;i++){ int j=(i+1)%n; int c=sgn(cross(pt-p[j],p[i]-p[j])); int u=sgn(p[i].y-pt.y); int v=sgn(p[j].y-pt.y); if(c>0&&u<0&&v>=0) num++; if(c<0&&u>=0&&v<0) num--; } return num!=0; //1:点在内部 0:点在外部 } //求多边形的面积 //以原点为中心点划分三角形,然后求多边形的面积 double polygon_area(point *p,int n){ double area=0; for(int i=0;i<n;i++) area+=cross(p[i],p[(i+1)%n]); return area/2; } //求多边形的重心 //将多边形进行三角剖分,算出每个三角形的重心,三角形的重心是3点坐标的平均值,然后对每个三角形的邮箱面积加权平均??? point polygon_center(point *p,int n){ point ans(0,0); if(polygon_area(p,n)==0) return ans; for(int i=0;i<n;i++){ ans=ans+(p[i]+p[(i+1)%n])*cross(p[i],p[(i+1)%n]); } return ans/polygon_area(p,n)/6; }
hdu 1115 求N(3<N<1000000)边形的重心
#include<iostream> #include<cmath> using namespace std; const double pi=acos(-1.0); const double eps=1e-8; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} point operator / (double k){return point(x/k,y/k);} bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;} bool operator < (point B){ return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);} }; typedef point Vector; //点积和叉积 //点积是A*B=|A||B|cosC C是A,B之间的夹角 //A*B=A.x*B.x+A.y*B.y double dot(Vector A,Vector B) { return A.x*B.x+A.y*B.y; } double cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; } //求多边形的面积 //以原点为中心点划分三角形,然后求多边形的面积 double polygon_area(point *p,int n){ double area=0; for(int i=0;i<n;i++) area+=cross(p[i],p[(i+1)%n]); return area/2; } point polygon_center(point *p,int n){ point ans(0,0); if(polygon_area(p,n)==0) return ans; for(int i=0;i<n;i++){ ans=ans+(p[i]+p[(i+1)%n])*cross(p[i],p[(i+1)%n]); } return ans/polygon_area(p,n)/6; } int main(){ int t,n,i; point center; point p[100000]; scanf("%d",&t); while(t--){ scanf("%d",&n); for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y); center=polygon_center(p,n); printf("%.2f %.2f\n",center.x,center.y); } return 0; }
凸包
搞懂概念:给定一些点,求所有能把这些点包含在内的面积最小的多边形
主要是用点求凸包的算法:Graham算法O(nlong2n)、jarvis步进法O(nh),h是凸包上的顶点数
andrew算法(算法复杂度:nlogn):从做左下扫描到右上得到下凸包,从右到左扫描得到上凸包
hdu 1392 求凸包的周长
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> using namespace std; //hdu 1392 求凸包的周长 const double pi=acos(-1.0); const double eps=1e-8; const int maxn=104; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} point operator / (double k){return point(x/k,y/k);} bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;} bool operator < (point B){ return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);} }; typedef point Vector; //点积和叉积 //点积是A*B=|A||B|cosC C是A,B之间的夹角 //A*B=A.x*B.x+A.y*B.y double dot(Vector A,Vector B) { return A.x*B.x+A.y*B.y; } double cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; } double Distance(point a,point b){ return hypot(a.x-b.x,a.y-b.y); } int convex_hull(point *p,int n,point *ch){ sort(p,p+n); n=unique(p,p+n)-p; int v=0; //求下凸包,如果p[i]是右拐弯的,那么这个点不在凸包上, 往回退 for(int i=0;i<n;i++){ while(v>1&&sgn(cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--; ch[v++]=p[i]; } int j=v; //求上凸包 for(int i=n-2;i>=0;i--){ while(v>j&&sgn(cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--; ch[v++]=p[i]; } if(n>1) v--; return v; } int main(){ int n; point p[maxn],ch[maxn]; while(scanf("%d",&n)&&n){ for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y); int v=convex_hull(p,n,ch); double ans=0; if(v==1) ans=0; else if(v==2) ans=Distance(ch[0],ch[1]); else{ for(int i=0;i<v;i++) ans+=Distance(ch[i],ch[(i+1)%v]); } printf("%.2f\n",ans); } return 0; }
Graham扫描法:
凸包--Graham扫描法_小蒟蒻yyb的博客-CSDN博客
首先找到最靠近左下的那个点,这个点一定在凸包上(不难理解吧。。。画个图就知道了)
以这个点为极点,其他点按照极角排序
然后按照顺序依次访问所有点,判断可行性
struct Node { int x,y; }p[MAX],S[MAX];//p储存节点的位置,S是凸包的栈 inline bool cmp(Node a,Node b)//比较函数,对点的极角进行排序 { double A=atan2((a.y-p[1].y),(a.x-p[1].x)); double B=atan2((b.y-p[1].y),(b.x-p[1].x)); if(A!=B)return A<B; else return a.x<b.x; //这里注意一下,如果极角相同,优先放x坐标更小的点 } long long Cross(Node a,Node b,Node c)//计算叉积 { return 1LL*(b.x-a.x)*(c.y-a.y)-1LL*(b.y-a.y)*(c.x-a.x); } void Get()//求出凸包 { p[0]=(Node){INF,INF};int k; for(int i=1;i<=n;++i)//找到最靠近左下的点 if(p[0].y>p[i].y||(p[0].y==p[i].y&&p[i].x<p[0].x)) {p[0]=p[i];k=i;} swap(p[k],p[1]); sort(&p[2],&p[n+1],cmp);//对于剩余点按照极角进行排序 S[0]=p[1],S[1]=p[2];top=1;//提前在栈中放入节点 for(int i=3;i<=n;)//枚举其他节点 { if(top&&Cross(S[top-1],p[i],S[top])>=0) top--;//如果当前栈顶不是凸包上的节点则弹出 else S[++top]=p[i++];//加入凸包的栈中 } //底下这个玩意用来输出凸包上点的坐标 //for(int i=0;i<=top;++i) // printf("(%d,%d)\n",S[i].x,S[i].y); }
最近点对
分治法O(nlog2n)
合并的时候,如果两个点在两个集合里,就不能简单合并,如果集合s1中最小距离是d1,集合s2中最小距离是d2,那么在两个集合中间的点找距离他们
小于d1,d2的,记录在tmp_p[]中,但不能用暴力法直接列出点集tmp_p[]中的所有点对,否则会超时,按照y坐标对tmp_p[]中的点排序,再加上剪枝
hdu 1007 求最近点对的一半
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> //最近点对 using namespace std; const double pi=acos(-1.0); const double eps=1e-8; const int maxn=100010; const double INF=1e20; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; }; double Distance(point a,point b){ return hypot(a.x-b.x,a.y-b.y); } bool cmpxy(point a,point b){ return sgn(a.x-b.x)<0||sgn(a.x-b.x)==0&&sgn(a.y-b.y)<0; } bool cmpy(point a,point b){ return sgn(a.y-b.y)<0; } point p[maxn],tmp_p[maxn]; double closest_pair(int left,int right){ double dis=INF; if(left==right) return dis; if(left+1==right) return Distance(p[left],p[right]); int mid=(left+right)/2; double d1=closest_pair(left,mid); double d2=closest_pair(mid+1,right); dis=min(d1,d2); int k=0; for(int i=left;i<=right;i++){ if(fabs(p[mid].x-p[i].x)<=dis){ tmp_p[k++]=p[i]; } } sort(tmp_p,tmp_p+k,cmpy); for(int i=0;i<k;i++){ for(int j=i+1;j<k;j++){ if(tmp_p[j].y-tmp_p[i].y>=dis) break; //剪枝 dis=min(dis,Distance(tmp_p[i],tmp_p[j])); } } return dis; } int main(){ int n; while(scanf("%d",&n)&&n){ for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y); sort(p,p+n,cmpxy); printf("%.2f\n",closest_pair(0,n-1)/2); } return 0; }
旋转卡壳
是一种思想,可以用来求出凸包上最远点对距离
凸包的直径——旋转卡壳 - 小蒟蒻yyb - 博客园 (cnblogs.com)
poj会有编译错误,为啥
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> using namespace std; const int maxn=50010; const int INF=0x3f3f3f3f; struct point{ int x,y; }p[maxn],p0,s[maxn]; int n,top,T; bool cmp(point a,point b){ double aa=atan2(a.y-p0.y,a.x-p0.x); double bb=atan2(b.y-p0.y,b.x-p0.x); if(aa!=bb) return aa<bb; else return a.x<b.x; } long long cross(int x1,int y1,int x2,int y2){ return (1LL*x1*y2-1LL*x2*y1); } long long cross2(point a,point b,point c){ //以a为中心点 return cross((b.x-a.x),(b.y-a.y),(c.x-a.x),(c.y-a.y)); } void convex_hull(){ //寻找凸包 p0=(point){INF,INF}; int k=0; for(int i=0;i<n;i++){ //寻找最下方的点 if(p0.y>p[i].y||(p0.y==p[i].y)&&(p0.x>p[i].x)) p0=p[i],k=i; } swap(p[k],p[0]); sort(&p[1],&p[n],cmp); //按照极角对下方的点进行排序 s[0]=p[0]; s[1]=p[1]; top=1; //栈顶 for(int i=2;i<n;){ //求出凸包 if(top&&cross2(s[top-1],p[i],s[top])>=0) top--; else s[++top]=p[i++]; //s里面放的是凸包边上的点 } } long long distan(point a,point b){ return 1LL*(a.x-b.x)*(a.x-b.x)+1LL*(a.y-b.y)*(a.y-b.y); } long long getmax(){ //求出直径 long long re=0; if(top==1) //仅有两个点 return distan(s[0],s[1]); s[++top]=s[0]; //把第一个点放到最后 int j=2; for(int i=0;i<top;i++){ //枚举边 while(cross2(s[i],s[i+1],s[j])<cross2(s[i],s[i+1],s[j+1])) j=(j+1)%top; re=max(re,max(distan(s[i],s[j]),distan(s[i+1],s[j]))); } return re; } int main(){ scanf("%d",&n); for(int i=0;i<n;i++) scanf("%lld %lld",&p[i].x,&p[i].y); long long ans=INF,ss; convex_hull(); ans=getmax(); printf("%lld\n",ans); return 0; }
半平面交
给出一堆半平面,然后求出凸包先
算法:增量法O(n^2)
这个算法:O(nlong2n)
主要是增加新板平面的时候根据情况进行队列的增删改查
应用:hdu 2297 "RUN"(真看不出来是半平面,没学过的话)
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> using namespace std; const int maxn=50010; const double INF=1e12; const double pi=acos(-1.0); const double eps=1e-8; int sgn(double x){ if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} }; typedef point Vector; double cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; } struct Line{ point p; Vector v; double ang; Line(){} Line(point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);} bool operator < (Line &L){ //用于极角排序 return ang<L.ang; } }; bool onleft(Line L,point p){ return sgn(cross(L.v,p-L.p))>0; //点p在L的左边:即点p在L的外面 } point cross_point(Line a,Line b){ //两条直线的交点 Vector u=a.p-b.p; double t=cross(b.v,u)/cross(a.v,b.v); return a.p+a.v*t; } vector<point> HPI(vector<Line> L){ //求半平面交,返回凸多边形 int n=L.size(); sort(L.begin(),L.end()); //将所有半平面按照极角排序 int first,last; vector<point> p(n); //两个相交半平面的交点 vector<Line> q(n); //双端队列 vector<point> ans; //半平面形成的凸包 q[first=last=0]=L[0]; for(int i=0;i<n;i++){ //情况1:删除尾部的半平面 while(first<last&&!onleft(L[i],p[last-1])) last--; //情况2:删除首部的半平面 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]=cross_point(q[last-1],q[last]); } //情况3:删除队列尾部的无用半平面:尾部的交点在第一条线段外面 while(first<last&&!onleft(q[first],p[last-1])) last--; if(last-first<=1) return ans; p[last]=cross_point(q[last],q[first]); //计算队列尾首部的交点 for(int i=first;i<=last;i++) ans.push_back(p[i]); return ans; //返回凸多边形的交点 } int main(){ int t,n; cin>>t; while(t--){ cin>>n; vector<Line> L; L.push_back(Line(point(0,0),Vector(0,-1))); L.push_back(Line(point(0,INF),Vector(-1,0))); //添加反向y轴,y极大的向左的直线 while(n--){ double a,b; scanf("%lf %lf",&a,&b); L.push_back(Line(point(0,a),Vector(1,b))); } vector<point> ans=HPI(L); printf("%d\n",ans.size()-2); //去掉认为加上的两个点 } return 0; }
圆
圆的定义、点和圆的关系、直线和圆的关系、线段和圆的关系、直线和圆的交点(会用到前面的点、线知识)
#include<bits/stdc++.h> using namespace std; const double pi=acos(-1.0); const double eps=1e-8; //圆:圆的定义、直线和圆的关系、线段和圆的关系、直线和圆的交点(会用到前面的点、线知识) int sgn(double x){ if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} }; double Distance(point a,point b){ //return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //或者用库函数hypot() return hypot(a.x-b.x,a.y-b.y); } struct Line{ point p1,p2; Line(){} //两个点 Line(point p1,point p2):p1(p1),p2(p2){} //一个点,一个倾斜角angle在0,pi之间 Line(point p,double angle){ p1=p; if(sgn(angle-pi/2)==0) {p2=(p1+point(0,1));} else{p2=(p1+point(1,tan(angle)));} } Line(double a,double b,double c){ if(sgn(a)==0){ p1=point(0,-c/b); p2=point(1,-c/b); } if(sgn(b)==0){ p1=point(-c/a,0); p2=point(-c/a,1); } else{ p1=point(0,-c/b); p2=point(1,((-c-a)/b)); } } }; double dis_point_line(point p,Line v){ return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2); } //点在直线上的投影 point point_line_proj(point p,Line v){ double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1); return v.p1+(v.p2-v.p1)*k; } //点关于直线的对称点 point point_line_symmetry(point p,Line v){ point q=point_line_proj(p,v); return point(2*q.x-p.x,2*q.y-p.y); } //线段的表示:直接用直线的结构 typedef Line Segment; //点和线段的关系 //先用叉积看是不是共线,然后点积看p与v的两个端点是不是形成的钝角 bool point_on_seg(point p,Line v){ return sgn(cross(p-v.p1,v.p2-v.p1))==0&&sgn(dot(p-v.p1,p-v.p2))<=0; } //点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值 double dis_point_seg(point p,Segment v){ if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2)); return dis_point_line(point(p,v)); } //圆的定义 struct circle{ point c; double r; circle(){} circle(point c,double r):c(r),r(r){} circle(double a,double b,double _r){ c=point(a,b);r=_r; } }; //点和圆的关系:根据点到圆心的距离判断 int point_circle_relation(point p,circle c){ double dst=Distance(p,c.c); if(sgn(dst-c.r)<0) return 0; //在园内 if(sgn(dst-c.r)==0) return 1; //在圆上 else return 2; //点在圆外 } //直线和圆的关系:根据圆心到直线的距离判断 int line_circle_relation(Line v,circle c){ double dst=dis_point_line(c.c,v); if(sgn(dst-c.r)<0) return 0; //直线和圆相交 if(sgn(dst-c.r)==0) return 1; //直线和圆相切 return 2; //在圆外 } //线段和圆的关系 int seg_circle_relation(Segment v,circle c){ double dst=dis_point_seg(c.c,v); if(sgn(dst-c.r)<0) return 0; //线段在圆内 if(sgn(dst-c.r)==0) return 1; //线段和圆相切 return 2; //线段在圆外 } //直线和圆的交点 pa和pb是交点,返回的是交点的个数 int line_cross_circle(Line v,circle c,point &pa,point &pb){ if(line_circle_relation(v,c)==2) return 0 ; //无交点 point q=point_line_proj(c.c,v); //圆心在直线上的投影点 double d=dis_point_line(c.c,v); //圆心到直线的距离 double k=sqrt(c.r*c.r-d*d); if(sgn(k)==0){ pa=q;pb=q;return 1; //相切的情况 } point n=(v.p2-v.p1)/len(v.p2-v.p1); //单位向量 pa=q+n*k; pb=q-n*k; return 2; //两个交点 }
圆的模板的应用:hdu 5572
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> //hdu 5572 圆的模板应用 using namespace std; const double pi=acos(-1.0); const double eps=1e-8; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} point operator / (double k){return point(x/k,y/k);} }; typedef point Vector; double dot(Vector A,Vector B) { return A.x*B.x+A.y*B.y; } double len(Vector A){ return sqrt(dot(A,A)); //或者求长度的平方 } double len2(Vector A){ return dot(A,A); } double cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; } //两点之间的距离 double Distance(point a,point b){ return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //或者用库函数hypot() //return hypot(a.x-b.x,a.y-b.y); } struct Line{ point p1,p2; Line(){} //两个点 Line(point p1,point p2):p1(p1),p2(p2){} }; typedef Line Segment; //点和直线的关系 int point_line_relation(point p,Line v){ int c=sgn(cross(p-v.p1,v.p2-v.p1)); if(c<0) return 1; //1 p在v的左边 if(c>0) return 2; //2 p在v的右边 return 0; //p在v上 } //点和直线的距离 p v(p1,p2) //用叉积求p,p1,p2构成的平行四边形面积,然后除以平行四边形的底边长,就得到了高 double dis_point_line(point p,Line v){ return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2); } //点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值 double dis_point_seg(point p,Segment v){ if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2)); return dis_point_line(p,v); } //点在直线上的投影 point point_line_proj(point p,Line v){ double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1); return v.p1+(v.p2-v.p1)*k; } //点关于直线的对称点 point point_line_symmetry(point p,Line v){ point q=point_line_proj(p,v); return point(2*q.x-p.x,2*q.y-p.y); } //圆的定义 struct circle{ point c; double r; circle(){} circle(point c,double r):c(c),r(r){} circle(double a,double b,double _r){ c=point(a,b);r=_r; } }; //直线和圆的关系:根据圆心到直线的距离判断 int line_circle_relation(Line v,circle c){ double dst=dis_point_line(c.c,v); if(sgn(dst-c.r)<0) return 0; //直线和圆相交 if(sgn(dst-c.r)==0) return 1; //直线和圆相切 return 2; //在圆外 } //线段和圆的关系 int seg_circle_relation(Segment v,circle c){ double dst=dis_point_seg(c.c,v); if(sgn(dst-c.r)<0) return 0; //线段在圆内 if(sgn(dst-c.r)==0) return 1; //线段和圆相切 return 2; //线段在圆外 } //直线和圆的交点 pa和pb是交点,返回的是交点的个数 int line_cross_circle(Line v,circle c,point &pa,point &pb){ if(line_circle_relation(v,c)==2) return 0 ; //无交点 point q=point_line_proj(c.c,v); //圆心在直线上的投影点 double d=dis_point_line(c.c,v); //圆心到直线的距离 double k=sqrt(c.r*c.r-d*d); if(sgn(k)==0){ pa=q;pb=q;return 1; //相切的情况 } point n=(v.p2-v.p1)/len(v.p2-v.p1); //单位向量 pa=q+n*k; pb=q-n*k; return 2; //两个交点 } int main(){ int t; scanf("%d",&t); for(int cas=1;cas<=t;cas++){ circle o; point a,b,v; scanf("%lf %lf %lf",&o.c.x,&o.c.y,&o.r); scanf("%lf %lf %lf %lf",&a.x,&a.y,&v.x,&v.y); scanf("%lf %lf",&b.x,&b.y); Line l(a,a+v); Line t(a,b); //情况1:直线和圆不想交,而且直线经过点 if(point_line_relation(b,l)==0&&seg_circle_relation(t,o)>=1&&sgn(cross(b-a,v))==0) printf("Case #%d: Yes\n",cas); else{ point pa,pb; //情况2:直线和圆相切,不经过点 if(line_cross_circle(l,o,pa,pb)!=2) printf("Case #%d: No\n",cas); else{ //情况3:直线和圆相交 point cut; //直线和圆的碰撞点 if(Distance(pa,a)>Distance(pb,a)) cut=pb; else cut=pa; Line mid(cut,o.c); //圆心到碰撞点的线 point en=point_line_symmetry(a,mid); //镜像点 Line light(cut,en); //反射线 if(Distance(light.p2,b)>Distance(light.p1,b)) swap(light.p1,light.p2); if(sgn(cross(light.p2-light.p1,point(b.x-cut.x,b.y-cut.y)))==0) printf("Case #%d: Yes\n",cas); else printf("Case #%d: No\n",cas) ; } } } return 0; }
最小圆覆盖:就是找个点能够覆盖平面上所有的点
两种方法:几何算法(增量法)和模拟退火(很慢)
hdu 3007 裸题
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> //hdu 5572 圆的模板应用 using namespace std; const double pi=acos(-1.0); const int maxn=505; const double eps=1e-8; int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (point B){return point(x+B.x,y+B.y);} point operator - (point B){return point(x-B.x,y-B.y);} point operator * (double k){return point(x*k,y*k);} point operator / (double k){return point(x/k,y/k);} }; //两点之间的距离 double Distance(point a,point b){ return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //或者用库函数hypot() //return hypot(a.x-b.x,a.y-b.y); } //求三角形abc的外接圆的圆心 point circle_center(const point a,const point b,const point c){ point center; 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; center.x=a.x+(c1*b2-c2*b1)/d; center.y=a.y+(a1*c2-a2*c1)/d; return center; } //求最小圆覆盖,返回圆心c,半径r void min_cover_circle(point *p,int n,point &c,double &r){ //随机函数,打乱所有点,很重要,降低算法复杂度 random_shuffle(p,p+n); c=p[0];r=0; //初始值 for(int i=1;i<n;i++){ if(sgn(Distance(p[i],c)-r)>0){ //点i在圆的外部,需要拓展 c=p[i];r=0; //这里就直接更新圆(相当于直接进队首) for(int j=0;j<i;j++){ //重新检查前面所有的点 if(sgn(Distance(p[j],c)-r)>0) { //两点定圆 c.x=(p[j].x+p[i].x)/2; c.y=(p[j].y+p[i].y)/2; r=Distance(p[j],c); for(int k=0;k<j;k++){ //两点不能定圆那就三点定圆 if(sgn(Distance(p[k],c)-r)>0){ c=circle_center(p[i],p[j],p[k]); r=Distance(p[i],c); } } } } } } } int main(){ int n; point p[maxn]; point c;double r; while(~scanf("%d",&n)&&n){ for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y); min_cover_circle(p,n,c,r); printf("%.2f %.2f %.2f\n",c.x,c.y,r); } return 0; }
模拟退火
//模拟退火算法 :算法复杂度很高 void min_cover_circle2(point *p,int n,point &c,double &r){ double T=100.0; //初识温度 double delta=0.98 //降温系数 c=p[0]; int pos; while(T>eps){ //eps是终止温度 pos=0;r=0; //初识:p[0]是圆心,半径是0 for(int i=0;i<=n-1;i++){ //找距离圆心最远的点 if(Distance(c,p[i])>r){ r=Distance(c,p[i]); //距圆心最远的点肯定在圆周上 pos=i; } } c.x+=(p[pos].x-c.x)/r*T; //逼近最后的解 c.y+=(p[pos].y-c.y)/r*T; T*=delta; } }
三维几何
除了一些基础的点、向量、直线、线段的表示以外,还有三维点积、叉积、三角形的面积、平面(3个点)、平面法向量、直线和平面的关系
四面体的有向体积
最小球覆盖、三维凸包(懒得写了)
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<stack> #include<cstdio> #include<queue> #include<map> #include<vector> #include<set> //三维几何 int sgn(double x){ //判断x是不是等于0 if(fabs(x)<eps) return 0; else return x<0? -1:1; } struct point3{ double x,y,z; point3(){} point3(double x,double y,double z) x(x),y(y),z(z){} point3 operator + (point3 b){return point3(x+b.x,y+b.y,z+b.z);} point3 operator - (point3 b){return point3(x-b.x,y-b.y,z-b.z);} point3 operator * (double b){return point3(x*b,y*b,z*b);} point3 operator / (point3 b){return point3(x/k,y/k,z/k);} bool operator == (point3 b){ return sgn(x-b.x)==0&&sgn(y-b.y)==0&&sgn(z-b.z)==0; } }; typedef point3 Vector3; //三维向量 double Distance(Vector3 a,Vector3 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)); } //线和线段 struct Line3{ point3 p1,p2; Line3(){} Line3(point3 p1,point3 p2):p1(p1),p2(p2){} }; typedef Line3 Segment3; //三维点积 double dot3(Vector3 a,Vector3 b){ return a.x*b.x+a.y*b.y+a.z*b.z; } //求向量的长度 double len3(Vector3 a){ return sqrt(dot3(a,a)); } //或者是求长度的平方 double len32(Vector3 a){ return dot3(a,a); } //求向量A,B的夹角大小 double angle3(Vector3 a,Vector3 b){ return acos(dot3(a,b)/len3(a,a)/len3(b,b)); } //三维叉积:是一个向量,可以把三维叉积A,B的叉积看成垂直于A和B的向量 Vector3 cross3(Vector3 a,Vector3 b){ return point3(a.y*b.z-a.z-b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y-b.x); } //三角形面积:有向面积,先求三维叉积,然后取叉积的长度值 double area32(point3 a,point3 b,point3 c){ return len3(cross3(b-a,c-a)); } // 判断点p在不在三角形ABC内:在内的话进行三角剖分后的面积相加为原本ABC的面积 dcmp(area32(p,a,b)+area32(p,b,c)+area32(p,a,c),area32(a,b,c))==0 //平面:用三个点可以确定一个平面 struct plane{ point3 p1,p2,p3; plane(){} plane(point3 p1,point3 p2,point3 p3):p1(p1),p2(p2),p3(p3){} }; //平面法向量 point3 pvec(point3 a,point3 b,point3 c){ return cross3(b-a,c-a); } //或者这样也行 point3 pvec(plane x){ return cross(x.p2-x.p1,x-p3-x.p1); } //直线和平面的关系 int line_corss_plane(Line3 u,plane f,point3 &p){ point3 v=pvec(f); double x=dot3(v,u.p2-f.p1); double y=dot3(v,u.p1-f.p1); double d=x-y; if(sgn(x)==0&&sgn(y)==0) return -1; //直线在平面上 if(sgn(d)==0) return 0; //平行 p=((u.p1*x)-(u.p2*y))/d; //相交的点 return 1; } //四面体的有向体积 double volume4(point3 a,point3 b,point3 c,point3 d){ return dot3(cross(b-a,c-a),d-a); } //最小球覆盖 //三维凸包:增量法、懒得写了