计算几何基础-代码部分

讲解部分看这篇博客

• 先来一个namespace geo，存下计算几何的基本类型和基本运算

//#pragma GCC optimize(2)
#include<cstdio>
#include<cmath>
#include<vector>
#include<algorithm>
using namespace std;
inline int read(){
	register int x,ch,f=1;while((ch=getchar())<48||57<ch)if(ch=='-')f=-1;
	x=ch^48;while(47<(ch=getchar())&&ch<58)x=x*10+(ch^48);
	return x*f;
}
namespace geo{
	#define op operator
	#define il inline 
	const double eps=1e-10;
	inline int cmp(double x){
		if(fabs(x)<eps)return 0;
		else return x<0?-1:1;
	}
	class vec {
		public:
		double x,y;
		vec(double x=0,double y=0):x(x),y(y) { }
		vec op+(const vec&b)const{
			return vec{x+b.x,y+b.y};
		} 
		vec op-(const vec&b)const{
			return vec{x-b.x,y-b.y};
		} 
		vec op*(double p)const{
			return vec{x*p,y*p};
		} 
		double op*(const vec& b)const{
			return x*b.x+y*b.y;
		}
		double len()const{
			return sqrt(*this**this);
		}
		double op^(const vec& b)const{
			return x*b.y-y*b.x;
		}
		vec op/(double p)const{
			return vec{x/p,y/p};
		} 
		bool op<(const vec&b)const{
			return x<b.x|| x==b.x&&y<b.y;
		}
		bool op==(const vec&b)const{
			return !cmp(x-b.x)&&!cmp(y-b.y)==0;
		}
		
	};
	typedef vec point;
	typedef vec Point;
	inline double len(const vec&x){
		return x.len();
	}
	inline int sqr(const vec& x){
		return x.x*x.x+x.y*x.y;
	}
	inline double area2(const point&a,const point&b,const point&c){
		return (b-a)^(c-a);
	}
	typedef struct Line{
		Point P;
		vec v;
		double ang;
		Line() { }
		Line(const Point&p,const vec&v): P(p),v(v){ ang=atan2(v.y,v.x);}
		bool operator<(const Line&L)const{
			return ang<L.ang;
		}
	}line;
	inline Point LineIntersection(point P,vec v,point Q,vec w){
		vec u=P-Q;
		double t=(w^u)/(v^w);
		return P+v*t;
	}
	inline Point Intersection(Line l1,Line l2){
		vec u=l1.P-l2.P;
		double t=(l2.v^u)/(l1.v^l2.v);
		return l1.P+l1.v*t;
	}
	inline bool OnSegment(point p,point a1,point a2){
		return cmp(area2(p,a1,a2))==0 && cmp((a1-p)*(a2-p))<0;
	}
	inline bool OnLeft(Line l,point P){
		return (l.v)*(P-l.P)>0;
	}
	typedef vector<Point> Polygon;
	inline double Area(const Polygon& p){
		double area=0;
		for(int i=1;i<p.size()-1;++i)
			area+=area2(p[0],p[i],p[i+1]);
		return area/2;
	}
	#undef op
	#undef il
};
using namespace geo;

• 求凸包的Andrew算法，凸包的点集存在ch[]里，返回值是点的个数。
  对效率要求不高则可以返回Polygon&会造成悬垂引用Polygon类型（或者inline+O2直接返回Polygon类型，编译器会作RVO优化）。

int Andrew(Point*p,int n,Point*ch){
	sort(p,p+n);
	int m=0;
	for(int i=0;i<n;++i){
		while(m>1&&area2(ch[m-2],ch[m-1],p[i])<=0)--m;
		ch[m++]=p[i];
	}
	int k=m;
	for(int i=n-2;~i;--i){
		while(m>k&&area2(ch[m-2],ch[m-1],p[i])<=0)--m;
		ch[m++]=p[i];
	}
	return m-(n>1);
}

• 旋转卡壳求凸包直径（凸包中任意两点的最大距离）
  这里为了精确，返回的是最大距离的平方。根据实际需要开方即可。

int diameter(point* p, int n){
	int ans=0;
	for(int i=0,j=1;i<n;++i){
		const point& A=p[i],B=p[(i+1)%n];
		while(cmp(area2(A,B,p[(j+1)%n])-area2(A,B,p[j]))>0)j=(j+1)%n;
			ans=max(ans,max(sqr(A-p[j]),sqr(B-p[j])));
	}
	return ans;
}