二维计算几何模板
本模板主要来自刘汝佳的《算法竞赛入门经典》,部分缺省的地方我进行了补齐,且对每个函数作出了解释
/*****************************************************************
点和向量的声明,以及加减,数乘比较等基础操作
*****************************************************************/
const double pi=acos(-1.0);
struct Point{
double x,y;
Point(double x=0,double y=0):x(x),y(y){
}
};//点的结构体
typedef Point Vector;//同时也是向量
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);}//向量相减 ,得到向量ba
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< (const Point &a,const Point &b)
{
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}//方便对点进行排序,依据的排序规则为水平坐标小的就小一些,当水平坐标相等时就纵坐标小的小
const double eps=1e-10;//可以适当调整为1e-8
int dcmp(double x)
{
if(fabs(x)<eps) return 0;
if(x<0) return -1;
return 1;
}
bool operator==(const Point &a,const Point &b)
{
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;//对向量(点)的相等进行重载
}
/*********************************************
向量之间的一些基本运算
*********************************************/
//向量与向量的计算
double Dot(Vector a,Vector b)
{
return a.x*b.x+a.y*b.y;//向量的点乘
}
double Length(Vector a)
{
return sqrt(Dot(a,a));//计算向量a的长度
}
double Angle(Vector a,Vector b)
{
return acos(Dot(a,b) / Length(a) / Length(b));//计算向量a,b的夹角
}
double Cross(Vector a,Vector b)
{
return a.x*b.y-a.y*b.x;
//计算a×b,切记,叉乘有前后之分
//b的在a的逆时针方向时为正
}
double Area2(Point a,Point b,Point c)
{
return Cross(b-a,c-a);
//计算a,b,c三点构成的三角形的面积的两倍
//注意:ab(由a指向b的向量),ac(a指向c的向量)。
//ac在ab的逆时针方向时,面积为正,反之为负
}
Vector Rotate(Vector a,double rad)
{
return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
//rad是弧度,旋转方向为逆时针
}
Vector Normal(Vector a)
{
double L=Length(a);
return Vector(-a.y/L,a.x/L);
//求出a向量的单位法向量
//逆时针旋转90度的法向量
}
//还有一种用虚数表现点的向量的方法,慢一些
/********************************************************
#include<complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;
double Dot(Vector a,Vector b)
{
return real(conj(a)*b);
}
double Cross(Vector a,Vector b)
{
return imag(conj(a)*b);
}
Vector Rotate(Vector a,double rad)
{
return a*(exp(Point(0,rad)));
}
//conj返回该复数的共轭复数
//imag返回复数的虚部
//real返回复数的实部
************************************************************************************/
/***********************************************************************************
以下为点和直线
************************************************************************************/
struct Line{
Point p;
Vector v;
Line(double px,double py,double ex,double ey):p(px,py),v(ex-px,ey-py){}
Line(double px,double py,Point a):
p(px,py),v(a.x-px,a.y-py){}
Line(Point a,Point b):
p(a.x,b.x),v(b.x-a.x,b.y-a.y){}
Point point(int x)//已知一点的横坐标,求其点
{
double t=(x-p.x)/v.x;
return Point(x,p.y+t*v.y);
}
};
Point GetLineIntersection(Point p,Vector v,Point q,Vector w)
{
Vector u=p-q;
double t=Cross(w,u)/Cross(v,w);
return p+v*t;
//返回直线p+tv和q+tw的交点(v,w分别是两条线的方向向量,p和q分别是直线上的一点)
//需要注意的是Cross(v,w)不能共线
}
double Distancetoline(Point p,Point a,Point b)
{
Vector v1=b-a,v2=p-a;
return fabs(Cross(v1,v2))/Length(v1);
//求得p到直线ab的距离
}
double Distancetosegment(Point p,Point a,Point b)
{
if(a==b) return Length(p-a);
Vector v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
if(dcmp(Dot(v1,v3))>0) return Length(v3);
return fabs(Cross(v1,v2))/Length(v1);
//求得p到直线ab的最小距离
}
Point GetLineProjection(Point p,Point a,Point b)
{
Vector v=b-a;
return a+v*(Dot(v,p-a) /Dot(v,v));
//p在直线ab上的投影点
}
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)
{
double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),
c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-a1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
//判断a1a2与b1b2是否相交(交点不在a1,a2,b1,b2上)
//如果c1c2表示两直线共线
}
//当c1c2不同时为零时,还有一种一条线段在;另一条直线上的端点的情况,需要用下函数判断某点是否在直线上
bool OnSegment(Point p,Point a1,Point a2)
{
return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
//用叉乘判断时候在这条直线上,用点乘判断是否在线段中
//判断p点是否在线段a1a2上
}
/***************************************
多边形面积公式
***************************************/
double PolygonArea(Point* p,int n)
{
double area=0;
for(int i=1;i<n-1;i++)
{
area+=Cross(p[i]-p[0],p[i+1]-p[0]);
}
return area/2;
//注意此方法的得出的面积有向,在不能保证点的分布是顺时针还是啊逆时针的时候需要加上绝对值
//p是点的数组,n是多边形的边数,也是点数
}
/**********************************************************
与圆和球有关的计算问题
**********************************************************/
struct Circle{
Point c;//圆心
double r;//半径
Circle(Point c,double r):c(c),r(r){}
Point point(double rad)
{
return Point(c.x+cos(rad)*r,c.y+sin(rad)*r);// 圆上rad角处点的坐标
}
};
int getLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point>& sol)
{
double a=L.v.x,b=L.p.x-C.c.x,c=L.v.y,d=L.p.y-C.c.y;
double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-C.r*C.r;
double delta=f*f-4*e*g;
if(dcmp(delta)<0) return 0; //相离
if(dcmp(delta)==0) //相切
{
t1=t2=-f/(2*e);//t为交点的横坐标
sol.push_back(L.point(t1));
return 1;
}
//相交
t1=(-f-sqrt(delta))/(2*e);sol.push_back(L.point(t1));
t2=(-f+sqrt(delta))/(2*e);sol.push_back(L.point(t2));
return 2;
}
//向量的极角
double angle(Vector v) {return atan2(v.y,v.x);}
int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol)
{
double d=Length(c1.c-c2.c);
if(dcmp(d)==0)
{
if(dcmp(c1.r-c2.r)==0) return -1;//两圆相交
return 0;//当共圆心时,如半径大小不相同则必定无交点
}
if(dcmp(c1.r+c2.r-d)<0) return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0) return 0;
double a=angle(c2.c-c1.c);
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));
//c1c2到c1p1的角(p1是交点)
Point p1=c1.point(a-da),p2=c1.point(a+da);
sol.push_back(p1);
if(p1==p2) return 1;
sol.push_back(p2);
return 2;
//本函数用于求两圆之间的相交情况
}
int getTangents(Point p,Circle c,Vector* v)
{
Vector u=c.c-p;
double dist=Length(u);
if(dist<c.r) return 0;//点在圆内不可能出现切线
if(dcmp(dist-c.r)==0) {
v[0]=Rotate(u,pi/2);
return 1;
//点在圆上恰有一条切线
}
double ang=asin(c.r/dist);
v[0]=Rotate(u,-ang);
v[1]=Rotate(u,ang);
return 2;//两条切线
//本函数用于求一个点和圆所能产生的切线
}
int getTangents(Circle A,Circle B,Point* a,Point* b)
{
int cnt=0;
if(A.r<B.r)
{
swap(A,B);
swap(a,b);
}
int d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);//圆心距离的平方
int rdiff=A.r-B.r;
int rsum=A.r+B.r;
if(d2<rdiff*rdiff) return 0;
//两圆的圆心距离小于半径差
//内含
double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);//圆心连线的极角
if(d2==0&&A.r==B.r) return -1;//两圆重合
if(d2==rdiff*rdiff) {
a[cnt]=A.point(base);b[cnt]=B.point(base);
cnt++;
return 1;
}
//外公切线
double ang=acos((A.r-B.r)/sqrt(d2));
a[cnt]=A.point(base+ang);b[cnt]=B.point(base+ang);cnt++;
a[cnt]=A.point(base-ang);b[cnt]=B.point(base-ang);cnt++;
if(d2==rsum*rsum){
a[cnt]=A.point(base);b[cnt]=B.point(base+pi);cnt++;
}
else if(d2>rsum*rsum){//相离
double ang=acos((A.r+B.r)/sqrt(d2));
a[cnt]=A.point(base+ang);b[cnt]=B.point(base+ang);cnt++;
a[cnt]=A.point(base-ang);b[cnt]=B.point(base-ang);cnt++;
}
return cnt;
}
/**********************************************************
二维几何常用算法
**********************************************************/
typedef vector<Point> Polygon;
int isPointInPolygon(Point p,Polygon poly)
{
int wn=0;
int n=poly.size();
for(int i=0;i<n;i++)
{
if(isPointOnSegment(p,poly[i],poly[(i+1%n)])) return -1;//边界
int k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
int d1=dcmp(poly[i].y-p.y);
int d2=dcmp( poly[(i+1)%n].y-p.y);
if(k>0&&d1<=0&&d2>0) wn++;
if(k<0&&d2<=0&&d1>0) wn--;
}
if(wn!=0) return 1;//内部
return 0;//外部
//本函数用于判断点是否在多边形内部
}
int ConvexHull(Point* p,int n,Point* ch)
{
//输入点数组p,个数n,输出点数组ch,函数返回凸包的顶点数量
//输入点钟不能有重复点
//函数执行完成后p的顺序将会被破坏
// 如果在凸包的边界上不希望出现输入的点的话把两个<=改成<
//精度高时采用dcmp比较
sort(p,p+n);
int m=0;
for(int i=0;i<n;i++)
{
while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0) m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--)
{
while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0) m--;
ch[m++]=p[i];
}
if(n>1) m--;
return m;
}
(未完待续)