學計算幾何專題已經兩年了(還不是因為剛剛過年了.....)
哼 都寫完了網頁崩潰全沒了 生氣了 哼 本來已經12點了超想睡覺的現在又要重寫一遍 哼 我不睡覺了
先扔一個完整版方便復制
//By Candy? #include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> using namespace std; typedef long long ll; const int N=1e4,INF=1e9; inline int read(){ char c=getchar();int x=0,f=1; while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();} while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();} return x*f; } const double eps=1e-8; const double pi=acos(-1); inline int sgn(double x){ if(abs(x)<eps) return 0; else return x<0?-1:1; } struct Vector{ double x,y; Vector(double a=0,double b=0):x(a),y(b){} bool operator <(const Vector &a)const{ return sgn(x-a.x)<0||(sgn(x-a.x)==0&&sgn(y-a.y)<0); } void print(){printf("%lf %lf\n",x,y);} }; typedef Vector Point; 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 b){return Vector(a.x*b,a.y*b);} Vector operator /(Vector a,double b){return Vector(a.x/b,a.y/b);} bool operator ==(Vector a,Vector b){return sgn(a.x-b.x)==0&&sgn(a.y-b.y)==0;} 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 Len(Vector a){return sqrt(Dot(a,a));} double Len2(Vector a){return Dot(a,a);} double Angle(Vector a,Vector b){ return acos(Dot(a,b)/Len(a)/Len(b)); } Vector Normal(Vector a){ return Vector(-a.y,a.x);//counterClockwise } Vector Rotate(Vector a,double rad){ return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); } //Line:use p and v //struct Line{ // Point p; // Vector v; // Line(){} // Line(Point p,Vector v):p(p),v(v){} // bool operator <(const Line a)const{ // return sgn(Cross(v,a.v))>=0; // } //}; //bool OnLeft(Line l,Point p){ // return sgn(Cross(l.v,p-l.p))>=0; //} //Point LI(Line a,Line b){ // Vector v=a.p-b.p; // double t=Cross(b.v,v)/Cross(a.v,b.v); // return a.p+a.v*t; //} //Line:use s and t struct Line{ Point s,t; Line(){} Line(Point a,Point b):s(a),t(b){} bool operator <(Line a)const{ return sgn(Cross(t-s,a.t-a.s))>=0; } }; double DisTL(Point p,Point a,Point b){ Vector v1=b-a,v2=p-a; return abs(Cross(v1,v2)/Len(v1)); } double DisTS(Point p,Point a,Point b){ if(a==b) return Len(p-a); Vector v1=b-a,v2=p-a,v3=p-b; if(sgn(Dot(v1,v2))<0) return Len(v2); else if(sgn(Dot(v1,v3))>0) return Len(v3); else return abs(Cross(v1,v2)/Len(v1)); } bool OnLeft(Line l,Point p){ return sgn(Cross(l.t-l.s,p-l.s))>=0; } bool OnSeg(Point p,Point a,Point b){ return DisTL(p,a,b)==0&&sgn(Dot(p-a,p-b)<0&&!(p==a))&&!(p==b); } Point LI(Line a,Line b){ Vector v=a.s-b.s,v1=a.t-a.s,v2=b.t-b.s; double t=Cross(v2,v)/Cross(v1,v2); return a.s+v1*t; } bool isLSI(Line l1,Line l2){ Vector v=l1.t-l1.s,u=l2.s-l1.s,w=l2.t-l1.s; return sgn(Cross(v,u))!=sgn(Cross(v,w)); } bool isSSI(Line l1,Line l2){ return isLSI(l1,l2)&&isLSI(l1,l2); } //---chong he //bool isSSI(Line l1,Line l2){ // Vector v1=l1.t-l1.s,v2=l2.t-l2.s; // if(sgn(Cross(v1,v2))==0){ // int flag=0; // Vector u=l2.s-l1.s,w=l2.t-l1.s; // if(sgn(Dot(u,w))<0) flag=1; // u=l2.s-l1.t,w=l2.t-l1.t; // if(sgn(Dot(u,w))<0) flag=1; // return flag; // } // else return isLSI(l1,l2)&&isLSI(l2,l1); //} Point Circumcenter(Point a,Point b,Point c){ Point p=(a+b)/2,q=(a+c)/2; Vector v=Normal(b-a),u=Normal(c-a); if(sgn(Cross(v,u))==0){ if(sgn(Len(a-b)+Len(b-c)-Len(a-c))==0) return (a+c)/2; if(sgn(Len(a-b)+Len(a-c)-Len(b-c))==0) return (b+c)/2; if(sgn(Len(a-c)+Len(b-c)-Len(a-b))==0) return (a+b)/2; } return LI(Line(p,p+v),Line(q,q+u)); } Point Barycenter(Point a,Point b,Point c){ return (a+b+c)/3; } bool cmpPolar(Point a,Point b){ return sgn(Cross(a,b))>0; } int PointInPolygon(Point p,Point poly[],int n){ int wn=0; for(int i=1;i<=n;i++){ if(sgn(DisTS(p,poly[i],poly[i%n+1]))==0) return -1; int k=sgn(Cross(poly[i%n+1]-poly[i],p-poly[i])), d1=sgn(poly[i].y-p.y),d2=sgn(poly[i%n+1].y-p.y); if(k>0&&d1<=0&&d2>0) wn++; if(k<0&&d2<=0&&d1>0) wn--; } return (bool)wn; } double PolygonArea(Point p[],int n){ double s=0; for(int i=2;i<n;i++) s+=Cross(p[i]-p[1],p[i+1]-p[1]); return abs(s/2); } bool isConvex(Point poly[],int n){ int last=0,now=0; for(int i=1;i<=n;i++){ now=sgn(Cross(poly[i%n+1]-poly[i],poly[(i+1)%n+1]-poly[i%n+1])); if(last==0||now==0||now*last>0) last=now; else return false; } return true; } int ConvexHull(Point p[],int n,Point ch[]){//cannot handle repeat point sort(p+1,p+1+n); int m=0; for(int i=1;i<=n;i++){ while(m>1&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--; ch[++m]=p[i]; } int k=m; for(int i=n-1;i>=1;i--){//n-1 while(m>k&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--; ch[++m]=p[i]; } if(n>1) m--;//the first rpoint return m; } double RotatingCalipers(Point p[],int n){ if(n==1) return 0; if(n==2) return Len(p[1]-p[2]); int now=1; double ans=0; p[n+1]=p[1]; for(int i=1;i<=n;i++){ while(sgn(DisTL(p[now],p[i],p[i+1])-DisTL(p[now+1],p[i],p[i+1]))<=0) now=now%n+1; ans=max(ans,Len(p[now]-p[i])); ans=max(ans,Len(p[now]-p[i+1])); } return ans; } void iniPolygon(Point p[],int &n,double inf){ n=0; p[++n]=Point(-inf,-inf); p[++n]=Point(inf,-inf); p[++n]=Point(inf,inf); p[++n]=Point(-inf,inf); } Point t[N];int tn; void CutPolygon(Point p[],int &n,Point a,Point b){//get the left of a->b tn=0; Point c,d,e; for(int i=1;i<=n;i++){ c=p[i],d=p[i%n+1]; if(sgn(Cross(b-a,c-a))>=0) t[++tn]=c; if(isLSI(Line(a,b),Line(c,d))){ e=LI(Line(a,b),Line(c,d)); t[++tn]=e; } } n=tn;for(int i=1;i<=n;i++)p[i]=t[i]; } // iniPolygon(q,m,INF); // for(int i=1;i<=n;i++) CutPolygon(q,m,p[i%n+1],p[i]); double minCircleCover(Point p[],int n,Point &c){ random_shuffle(p+1,p+1+n); c=p[1]; double r=0; for(int i=2;i<=n;i++) if(sgn(Len(c-p[i])-r)>0){ c=p[i],r=0; for(int j=1;j<i;j++) if(sgn(Len(c-p[j])-r)>0){ c=(p[i]+p[j])/2,r=Len(c-p[i]); for(int k=1;k<j;k++) if(sgn(Len(c-p[k])-r)>0){ c=Circumcenter(p[i],p[j],p[k]); r=Len(c-p[i]); } } } return r; } //jie xi ji he double a,b,c; double f(Point p){return a*p.x+b*p.y+c;} Point abcLI(Line l){ double u=abs(f(l.s)),v=abs(f(l.t)); return Point(l.s.x*v+l.t.x*u,l.s.y*v+l.t.y*u)/(u+v); } double F(double x){return x;} //function inline double cal(double l,double r,double fl,double fr,double fm){ return (fl+fr+4*fm)*(r-l)/6; } double Simpson(double l,double r,double now,double fl,double fr,double fm){ double mid=(l+r)/2,flm=F((l+mid)/2),frm=F((mid+r)/2); double p=cal(l,mid,fl,fm,flm),q=cal(mid,r,fm,fr,frm); if(sgn(now-p-q)==0) return now; else return Simpson(l,mid,p,fl,fm,flm)+Simpson(mid,r,q,fm,fr,frm); } //sao miao xian
並不想多說什么 ljw的課件已經很好了
1.精度問題
const double eps=1e-8; inline int sgn(double x){ if(abs(x)<eps) return 0; else return x<0?-1:1; }
2.向量 點
向量+-*／ 比較
點積 cosΘ
叉積 sinΘ aXb>0 b在a的左面(逆時針方向)
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 Len(Vector a){return sqrt(Dot(a,a));} double Len2(Vector a){return Dot(a,a);} double Angle(Vector a,Vector b){ return acos(Dot(a,b)/Len(a)/Len(b)); }
法向量 (-y,x) 逆時針90
旋轉 想一個與坐標軸平行的直角然后一推就可以了
Vector Normal(Vector a){ return Vector(-a.y,a.x);//counterClockwise } Vector Rotate(Vector a,double rad){ return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); }
3.直線 線段
感覺用兩點比較方便
struct Line{ Point s,t; Line(){} Line(Point a,Point b):s(a),t(b){} bool operator <(Line a)const{ return sgn(Cross(t-s,a.t-a.s))>=0; } };
點到直線距離：面積法
線段最短距離：用點積特判是點到端點到距離的情況
點線位置關系 點在線段上
double DisTL(Point p,Point a,Point b){ Vector v1=b-a,v2=p-a; return abs(Cross(v1,v2)/Len(v1)); } double DisTS(Point p,Point a,Point b){ if(a==b) return Len(p-a); Vector v1=b-a,v2=p-a,v3=p-b; if(sgn(Dot(v1,v2))<0) return Len(v2); else if(sgn(Dot(v1,v3))>0) return Len(v3); else return abs(Cross(v1,v2)/Len(v1)); } bool OnLeft(Line l,Point p){ return sgn(Cross(l.t-l.s,p-l.s))>=0; } bool OnSeg(Point p,Point a,Point b){ return DisTL(p,a,b)==0&&sgn(Dot(p-a,p-b)<0&&!(p==a))&&!(p==b); }
直線相交交點：面積求比例，正負畫個圖看一下是否符合就行了
直線與線段相交：端點叉積不同
線段與線段相交：兩遍直線與線段(有時候特判重合)
Point LI(Line a,Line b){ Vector v=a.s-b.s,v1=a.t-a.s,v2=b.t-b.s; double t=Cross(v2,v)/Cross(v1,v2); return a.s+v1*t; } bool isLSI(Line l1,Line l2){ Vector v=l1.t-l1.s,u=l2.s-l1.s,w=l2.t-l1.s; return sgn(Cross(v,u))!=sgn(Cross(v,w)); } bool isSSI(Line l1,Line l2){ return isLSI(l1,l2)&&isLSI(l1,l2); }
4.三角形
外心：兩個中垂線交點，特判三點共線
重心：三個點平均
Point Circumcenter(Point a,Point b,Point c){ Point p=(a+b)/2,q=(a+c)/2; Vector v=Normal(b-a),u=Normal(c-a); if(sgn(Cross(v,u))==0){ if(sgn(Len(a-b)+Len(b-c)-Len(a-c))==0) return (a+c)/2; if(sgn(Len(a-b)+Len(a-c)-Len(b-c))==0) return (b+c)/2; if(sgn(Len(a-c)+Len(b-c)-Len(a-b))==0) return (a+b)/2; } return LI(Line(p,p+v),Line(q,q+u)); } Point Barycenter(Point a,Point b,Point c){ return (a+b+c)/3; }
5.多邊形
bool cmpPolar(Point a,Point b){ return sgn(Cross(a,b))>0; }
double PolygonArea(Point p[],int n){ double s=0; for(int i=2;i<n;i++) s+=Cross(p[i]-p[1],p[i+1]-p[1]); return abs(s/2); }
選一個基准點，不停叉積
判斷點在多邊形內：射線法
- 引一條向右的水平射線
- 注意判斷點在多邊形上的情況
- 逆時針穿過邊界，轉數+1，順時針則-1
- 注意射線穿過頂點的情況 下算上不算
int PointInPolygon(Point p,Point poly[],int n){ int wn=0; for(int i=1;i<=n;i++){ if(sgn(DisTS(p,poly[i],poly[i%n+1]))==0) return -1; int k=sgn(Cross(poly[i%n+1]-poly[i],p-poly[i])), d1=sgn(poly[i].y-p.y),d2=sgn(poly[i%n+1].y-p.y); if(k>0&&d1<=0&&d2>0) wn++; if(k<0&&d2<=0&&d1>0) wn--; } return (bool)wn; }
凸包：
判斷凸包：相鄰邊叉積同號
求凸包：Graham掃描法 我不重寫了！細節見代碼
bool isConvex(Point poly[],int n){ int last=0,now=0; for(int i=1;i<=n;i++){ now=sgn(Cross(poly[i%n+1]-poly[i],poly[(i+1)%n+1]-poly[i%n+1])); if(last==0||now==0||now*last>0) last=now; else return false; } return true; } int ConvexHull(Point p[],int n,Point ch[]){//cannot handle repeat point sort(p+1,p+1+n); int m=0; for(int i=1;i<=n;i++){ while(m>1&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--; ch[++m]=p[i]; } int k=m; for(int i=n-1;i>=1;i--){//n-1 while(m>k&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--; ch[++m]=p[i]; } if(n>1) m--;//the first rpoint return m; }
旋轉卡殼
- 被凸包上被一對平行直線卡住的點對叫對踵點
- 答案一定出在一對對踵點上
- 嘗試用通過旋轉一對平行直線，枚舉所有對踵點
- 按逆時針順序枚舉凸包上一條邊，則凸包上到該邊所在直線最遠的點也單調逆時針旋轉，相當於用一條直線卡對面一個頂點
double RotatingCalipers(Point p[],int n){ if(n==1) return 0; if(n==2) return Len(p[1]-p[2]); int now=1; double ans=0; p[n+1]=p[1]; for(int i=1;i<=n;i++){ while(sgn(DisTL(p[now],p[i],p[i+1])-DisTL(p[now+1],p[i],p[i+1]))<=0) now=now%n+1; ans=max(ans,Len(p[now]-p[i])); ans=max(ans,Len(p[now]-p[i+1])); } return ans; }
注意特判；變種很多 見題目
半平面交
- 初始化時加上一個范圍巨大的“框”
- 每次拿一個新的半平面切割原先的凸集
- 保留在新加直線左邊的點，刪除右邊的，有向直線與多邊形相交產生的新的點加入到新多邊形內
void iniPolygon(Point p[],int &n,double inf){ n=0; p[++n]=Point(-inf,-inf); p[++n]=Point(inf,-inf); p[++n]=Point(inf,inf); p[++n]=Point(-inf,inf); } Point t[N];int tn; void CutPolygon(Point p[],int &n,Point a,Point b){//get the left of a->b tn=0; Point c,d,e; for(int i=1;i<=n;i++){ c=p[i],d=p[i%n+1]; if(sgn(Cross(b-a,c-a))>=0) t[++tn]=c; if(isLSI(Line(a,b),Line(c,d))){ e=LI(Line(a,b),Line(c,d)); t[++tn]=e; } } n=tn;for(int i=1;i<=n;i++)p[i]=t[i]; }
最小圓覆蓋：三點確定圓 隨即增量法 O(n)
double minCircleCover(Point p[],int n,Point &c){ random_shuffle(p+1,p+1+n); c=p[1]; double r=0; for(int i=2;i<=n;i++) if(sgn(Len(c-p[i])-r)>0){ c=p[i],r=0; for(int j=1;j<i;j++) if(sgn(Len(c-p[j])-r)>0){ c=(p[i]+p[j])/2,r=Len(c-p[i]); for(int k=1;k<j;k++) if(sgn(Len(c-p[k])-r)>0){ c=Circumcenter(p[i],p[j],p[k]); r=Len(c-p[i]); } } } return r; }
6.辛普森積分
用二次函數擬合
我也是看過托馬斯微積分的....
double F(double x){return x;} //function inline double cal(double l,double r,double fl,double fr,double fm){ return (fl+fr+4*fm)*(r-l)/6; } double Simpson(double l,double r,double now,double fl,double fr,double fm){ double mid=(l+r)/2,flm=F((l+mid)/2),frm=F((mid+r)/2); double p=cal(l,mid,fl,fm,flm),q=cal(mid,r,fm,fr,frm); if(sgn(now-p-q)==0) return now; else return Simpson(l,mid,p,fl,fm,flm)+Simpson(mid,r,q,fm,fr,frm); }
7.解析幾何
求直線與直線方程交點的神奇方法
double a,b,c; double f(Point p){return a*p.x+b*p.y+c;} Point abcLI(Line l){ double u=abs(f(l.s)),v=abs(f(l.t)); return Point(l.s.x*v+l.t.x*u,l.s.y*v+l.t.y*u)/(u+v); }
8.掃描線：見題目