计算几何模板

【整理于2016年】

计算几何模板
const double eps = 1e-8;
const double pi = acos(-1.0);
using namespace std;
//采用eps的精度判断大/小于零
int epssgn(double x)
{
    if (fabs(x)<eps) return 0;
    else return x<0?-1:1;
}
//二维点类
struct Point
{
    double x,y;
    Point(double a=0,double b=0){x=a;y=b;}
};
typedef Point Vector;//矢量可坐标化表示，故用同一个结构体
//二维直线类，一般方程ax+by+c=0
struct Line
{
    double a,b,c,angle;
    Point p1,p2;
    Line(Point s,Point e)
    {
        a=s.y-e.y;
        b=e.x-s.x;
        c=s.x*e.y-e.x*s.y;
        angle=atan2(e.y-s.y,e.x-s.x);
        p1=s;p2=e;
    }
    Line(){}
};
//二维线段类
struct Segment
{
    Point s,e;
    Segment(Point a,Point b){s=a;e=b;}
    Segment(double x1,double y1,double x2,double y2)
    {
        s=Point(x1,y1);
        e=Point(x2,y2);
    }
    Segment(){}
};
//向量的加减及数乘
Vector operator + (Point a,Point b)
{
    return Vector(a.x+b.x,a.y+b.y);
}
Vector operator - (Point a,Point b)
{
    return Vector(a.x-b.x,a.y-b.y);
}
Vector operator * (Point a,double k)
{
    return Vector(a.x*k,a.y*k);
}
Vector operator / (Point a,double k)
{
    return Vector(a.x/k,a.y/k);
}
//求向量的模（长度）
double len(Vector a)
{
    return sqrt(a.x*a.x+a.y*a.y);
}
//得到sp-op和ep-op的叉积
//>0时ep在opsp的逆时针方向，<0时顺时针，=0时共线
double Cross(Point sp, Point ep, Point op)
{
    return (sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y);
}
//两向量求叉积，求三角形面积需要除以2
double Cross(Vector a,Vector b)
{
    return a.x*b.y-b.x*a.y;
}
//两向量求点积
double Dot(Vector a,Vector b)
{
    return a.x*b.x+a.y*b.y;
}
//求最大最小值
double max(double a,double b)
{
    if (a<b) return b; else return a;
}
double min(double a,double b)
{
    if (a>b) return b; else return a;
}
//三点共线
bool Inline(Point a,Point b,Point c)
{
    return !sgn(Cross(b,c,a));
}
//点a在线段bc上
bool PointOnSeg(const Point &a, const Point &b, const Point &c)
{
    return Inline(b,c,a) && Dot(a-c,a-b)<eps;
}
//求两点之间的直线距离
double dis(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
//判断两直线是否平行
int Parallel(Line l1,Line l2)
{
    return (fabs(l1.a*l2.b-l2.a*l1.b)<eps);
}
//判断两直线是否相等
int LineEqual(Line l1,Line l2)
{
    if (!Parallel(l1,l2)) return 0;
    else return (fabs(l1.a*l2.c-l2.a*l1.c)<eps && fabs(l1.b*l2.c-l2.b*l1.c)<eps);
}
//求点C到线段AB的距离 
double PointToSegDist(Point A,Point B,Point C)
{
    if (dis(A,B)<eps) return dis(B,C);
    if (epssgn(Dot(B-A,C-A))<0) return dis(A,C);
    if (epssgn(Dot(A-B,C-B))<0) return dis(B,C);
    return fabs(Cross(B-A,C-A))/dis(B,A);
}
//求线段两端AB到另一线段CD的距离
double TwoSegMinDist(Point A,Point B,Point C,Point D)
{
    return min(min(PointToSegDist(A,B,C),PointToSegDist(A,B,D)),
               min(PointToSegDist(C,D,A),PointToSegDist(C,D,B)));
}
Point SymPoint(Point p,Line l) //求二维平面上点p关于直线p1p2的对称点
{
    Point result;
    double a=l.p2.x-l.p1.x;
    double b=l.p2.y-l.p1.y;
    double t=((p.x-l.p1.x)*a+(p.y-l.p1.y)*b)/(a*a+b*b);
    result.x=2*l.p1.x+2*a*t-p.x;
    result.y=2*l.p1.y+2*b*t-p.y;
    return result;
}
//直接求两线段ab和cd的交点（确定相交的情况下）
Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d)
{
    double u = Cross(b,c,a), v = Cross(a,d,b);
    return Point((c.x * v + d.x * u)/(u+v), (c.y * v + d.y * u)/(u+v));
}
//判断线段s1e1与s2e2是否相交（含端点）
//不含端点的话将下面的<=改成<
int IsSegmentIntersect(Point s1, Point e1, Point s2, Point e2)
{
    if( min(s1.x,e1.x)<= max(s2.x,e2.x) &&
        min(s1.y,e1.y)<= max(s2.y,e2.y) &&
        min(s2.x,e2.x)<= max(s1.x,e1.x) &&
        min(s2.y,e2.y)<= max(s1.y,e1.y) &&
        Cross(s2,e2,s1)*Cross(s2,e2,e1)<=0 &&
        Cross(s1,e1,s2)*Cross(s1,e1,e2)<=0)
    return 1;
    return 0;
}
//知道直线上两点p1p2，判断直线与线段se是否相交，含端点
int IsLineIntersectSegment(Point p1,Point p2,Point s,Point e)
{
    if (Cross(p1,p2,s)*Cross(p1,p2,e)>eps) return 0;
    else return 1;
}
int IsLineIntersectSegment(Line l1,Point s,Point e)
{
    if (Cross(l1.p1,l1.p2,s)*Cross(l1.p1,l1.p2,e)>eps) return 0;
    else return 1;
}
//求两条直线l1和l2的交点
Point GetIntersect(Line l1, Line l2) 
{
    Point res;
    res.x=(l1.b*l2.c-l2.b*l1.c)/(l1.a*l2.b-l2.a*l1.b);
    res.y=(l1.c*l2.a-l2.c*l1.a)/(l1.a*l2.b-l2.a*l1.b);
    return res;
}
//求质量均匀分布的多边形的重心
Point BaryCenter(Point *p,int n)
{
    Point res(0,0);
    double s=0.0,t;
    int i;
    for (i=1;i<n-1;i++)
    {
        t=Cross(p[i]-p[0],p[i+1]-p[0])/2;   //分割成三角形，算面积
        s+=t;
        res.x+=(p[0].x+p[i].x+p[i+1].x)*t;  //将面积作为重量放在三角形的重心上
        res.y+=(p[0].y+p[i].y+p[i+1].y)*t;  //质量*坐标，三角形重心所需的除以3放在后面
    }
    res.x/=(3*s);
    res.y/=(3*s);
    return res;
}
//求多边形面积，要求点集按逆时针顺序
double ConvexPolygonArea(Point *p,int n)
{
    int i;
    double area=0;
    for (i=1;i<n-1;i++) area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return area/2;
}
//圆及相关
struct circle
{
    Point o;
    double r;
};
double dis(circle a, circle b)
{
    return sqrt((a.o.x - b.o.x) * (a.o.x - b.o.x) + (a.o.y - b.o.y) * (a.o.y - b.o.y));
}
//两圆相交面积：CircleIntersectArea
double CIA(circle a, circle b)
{
    double d = dis(a,b);
    if(d >= a.r + b.r) return 0;
    if(d <= fabs(a.r - b.r))
    {
        double r = a.r < b.r ? a.r : b.r;
        return pi*r*r;
    }
    double ang1 = acos((a.r * a.r + d * d - b.r * b.r)/2. / a.r / d );
    double ang2 = acos((b.r * b.r + d * d - a.r * a.r)/2. / b.r / d );
    double ret = ang1 * a.r * a.r + ang2 * b.r * b.r - d * a.r * sin(ang1);
    return ret;
}
//两圆环相交面积，不需其他模板只需一个公式：
//考虑圆环A和B，相交面积一致性地表达为
// S = A大B大 - A小B大 - A大B小 + A小B小
//所以用四遍CIA函数即可
//一般不用考虑距离判断相交，这个公式有一般性
//求以原点为圆心，过a、b，半径为r的优弧扇形面积（即圆心角小于180°）
//如果是以其他点为圆心，事先做一下坐标平移然后再调用
double SectorArea(Point a,Point b,double r)
{
    double theta=atan2(a.y,a.x)-atan2(b.y,b.x);
    while (theta<=0) theta+=2*pi;
    while (theta>=2*pi) theta-=2*pi;
    theta=min(theta,2*pi-theta);
    return r*r*theta/2;
}
//另一个版本：求以r点为圆心，R为半径，两端顶点为a,b的优弧扇形面积（不用特意平移处理到原点形式了）
double sectorarea(const Point &r, const Point &a, const Point &b,double R)
{
    double A2 = Dot(r-a,r-a), B2 = Dot(r-b,r-b), C2 = Dot(a-b,a-b);
    return R * R * acos((A2 + B2 - C2)*0.5 / sqrt(A2) / sqrt(B2) ) * 0.5;
}
// 线段与圆相交（判定）
//形参分别为：线段端点a，线段端点b，圆心r，圆半径R，相交点p1,p2
//返回值rtos,目测是r到线段的最小距离
double LineCrossCircle(const Point &a, const Point &b, const Point &r, double R, Point &p1, Point &p2)
{
    Point fp = LineCross(r, Point(r.x + a.y - b.y, r.y + b.x -a.x), a, b);
    double rtol = dis(r,fp);
    double rtos = PointOnSeg(fp,a,b)? rtol : min(dis(a,r),dis(b,r));
    double atob = dis(a,b);
    double fptoe = sqrt(R * R - rtol * rtol)/atob;
    if(rtos > R -eps) return rtos;
    p1 = fp + (a-b)*fptoe;
    p2 = fp + (b-a)*fptoe;
    return rtos;
}
//求以原点为圆心半径为r的园，与原点、a、b组成的三角形的重叠面积
// Circle & Triangle Intersect Area.
// circle : center(0,0),radius=r;
// triangle : a,b,and(0,0);
// maybe you need transform the coordinates to use this function.
double CTIA(Point a,Point b,double r)
{
    Point p[2];
    int num=0;
    int ina=(epssgn(len(a)-r)<0);
    int inb=(epssgn(len(b)-r)<0);
    if (ina)
    {
        if (inb) return fabs(Cross(a,b))/2.0; //情形1：三角形完全在圆内：直接求三角形面积
        else  //情形2：一个点在圆内，一个点在圆外：分割成一个完全在圆内的三角形和一个扇形
        {
            CircleCrossLine(a,b,Point(0,0),r,p,num);
            return SectorArea(b,p[0],r)+fabs(Cross(a,p[0]))/2.0;
        }
    }
    else
    {
        CircleCrossLine(a,b,Point(0,0),r,p,num);
        //情形2：一个点在圆内，一个点在圆外：分割成一个完全在圆内的三角形和一个扇形
        if (inb) return SectorArea(p[0],a,r)+fabs(Cross(p[0],b))/2.0;
        else
        {
            //情形4：两个点都在圆外