# [计算几何]【学习笔记】

哼 都写完了网页崩溃全没了 生气了 哼 本来已经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;
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
}
}

//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

1.精度问题

const double eps=1e-8;
inline int sgn(double x){
if(abs(x)<eps) return 0;
else return x<0?-1:1;
}

2.向量 点

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
}
}

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;
}

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];
}

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);
}
ABC

8.扫描线：见题目

posted @ 2017-02-01 23:25  Candy?  阅读(1981)  评论(2编辑  收藏