Computational Geometry Template (一)

 
#include    
#include    
  
#define eps 1e-8   
#define zero(x) (((x)>0 (x):-(x))eps;  
}  
  
//判两点在线段异侧,点在线段上返回0   
bool opposite_side(point p1,point p2,line l)  
{  
    return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps;  
}  
  
//判两直线平行   
bool parallel(line u,line v)  
{  
    return zero((u.a.x-u.b.x)*(v.a.y-v.b.y)-(v.a.x-v.b.x)*(u.a.y-u.b.y));  
}  
  
//判两直线垂直   
bool perpendicular(line u,line v)  
{  
    return zero((u.a.x-u.b.x)*(v.a.x-v.b.x)+(u.a.y-u.b.y)*(v.a.y-v.b.y));  
}  
  
//判两线段相交,包括端点和部分重合   
bool intersect_in(line u,line v)  
{  
    if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))  
        return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);  
    return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);  
}  
  
//判两线段相交,不包括端点和部分重合   
bool intersect_ex(line u,line v)  
{  
    return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);  
}  
  
//计算两直线交点,注意事先判断直线是否平行!   
//线段交点请另外判线段相交(同时还是要判断是否平行!)   
point intersection(line u,line v)  
{  
    point ret=u.a;  
    double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))  
             /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));  
    ret.x+=(u.b.x-u.a.x)*t;  
    ret.y+=(u.b.y-u.a.y)*t;  
    return ret;  
}  
point intersection(point u1,point u2,point v1,point v2)  
{  
    point ret=u1;  
    double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
    ret.x+=(u2.x-u1.x)*t;  
    ret.y+=(u2.y-u1.y)*t;  
    return ret;  
}  
//点到直线上的最近点   
point ptoline(point p,line l)  
{  
    point t=p;  
    t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;  
    return intersection(p,t,l.a,l.b);  
}  
  
//点到直线距离   
double disptoline(point p,line l)  
{  
    return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b);  
}  
  
//点到线段上的最近点   
point ptoseg(point p,line l)  
{  
    point t=p;  
    t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;  
    if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps)  
        return distance(p,l.a)eps)  
        return distance(p,l.a)





=pi+pi)  
        dlng-=pi+pi;  
    if (dlng>pi)  
        dlng=pi+pi-dlng;  
    lat1*=pi/180,lat2*=pi/180;  
    return acos(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2));  
}  
//计算距离,r 为球半径   
double line_dist(double r,double lng1,double lat1,double lng2,double lat2)  
{  
    double dlng=fabs(lng1-lng2)*pi/180;  
    while (dlng>=pi+pi)  
        dlng-=pi+pi;  
    if (dlng>pi)  
        dlng=pi+pi-dlng;  
    lat1*=pi/180,lat2*=pi/180;  
    return r*sqrt(2-2*(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2)));  
}  
//计算球面距离,r 为球半径   
inline double sphere_dist(double r,double lng1,double lat1,double lng2,double lat2)  
{  
    return r*angle(lng1,lat1,lng2,lat2);  
}  
//外心   
point circumcenter(point a,point b,point c)  
{  
    line u,v;  
    u.a.x=(a.x+b.x)/2;  
    u.a.y=(a.y+b.y)/2;  
    u.b.x=u.a.x-a.y+b.y;  
    u.b.y=u.a.y+a.x-b.x;  
    v.a.x=(a.x+c.x)/2;  
    v.a.y=(a.y+c.y)/2;  
    v.b.x=v.a.x-a.y+c.y;  
    v.b.y=v.a.y+a.x-c.x;  
    return intersection(u,v);  
}  
//内心   
point incenter(point a,point b,point c)  
{  
    line u,v;  
    double m,n;  
    u.a=a;  
    m=atan2(b.y-a.y,b.x-a.x);  
    n=atan2(c.y-a.y,c.x-a.x);  
    u.b.x=u.a.x+cos((m+n)/2);  
    u.b.y=u.a.y+sin((m+n)/2);  
    v.a=b;  
    m=atan2(a.y-b.y,a.x-b.x);  
    n=atan2(c.y-b.y,c.x-b.x);  
    v.b.x=v.a.x+cos((m+n)/2);  
    v.b.y=v.a.y+sin((m+n)/2);  
    return intersection(u,v);  
}  
//垂心   
point perpencenter(point a,point b,point c)  
{  
    line u,v;  
    u.a=c;  
    u.b.x=u.a.x-a.y+b.y;  
    u.b.y=u.a.y+a.x-b.x;  
    v.a=b;  
    v.b.x=v.a.x-a.y+c.y;  
    v.b.y=v.a.y+a.x-c.x;  
    return intersection(u,v);  
}  
//重心   
//到三角形三顶点距离的平方和最小的点   
//三角形内到三边距离之积最大的点   
point barycenter(point a,point b,point c)  
{  
    line u,v;  
    u.a.x=(a.x+b.x)/2;  
    u.a.y=(a.y+b.y)/2;  
    u.b=c;  
    v.a.x=(a.x+c.x)/2;  
    v.a.y=(a.y+c.y)/2;  
    v.b=b;  
    return intersection(u,v);  
}  
//费马点   
//到三角形三顶点距离之和最小的点   
point fermentpoint(point a,point b,point c)  
{  
    point u,v;  
    double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y);  
    int i,j,k;  
    u.x=(a.x+b.x+c.x)/3;  
    u.y=(a.y+b.y+c.y)/3;  
    while (step>1e-10)  
    {  
        for (k=0; k<10; step/=2,k++)  
        {  
            for (i=-1; i<=1; i++)  
            {  
                for (j=-1; j<=1; j++)  
                {  
                    v.x=u.x+step*i;  
                    v.y=u.y+step*j;  
                    if(distance(u,a)+distance(u,b)+distance(u,c)>distance(v,a)+distance(v,b)+distance(v,c))  
                    {  
                        u=v;  
                    }  
                }  
            }  
        }  
    }  
    return u;  
}  
  
  
//矢量差 U - V   
point3 subt(point3 u,point3 v)  
{  
    point3 ret;  
    ret.x=u.x-v.x;  
    ret.y=u.y-v.y;  
    ret.z=u.z-v.z;  
    return ret;  
}  
//取平面法向量   
point3 pvec(plane3 s)  
{  
    return xmult(subt(s.a,s.b),subt(s.b,s.c));  
}  
point3 pvec(point3 s1,point3 s2,point3 s3)  
{  
    return xmult(subt(s1,s2),subt(s2,s3));  
}  
//两点距离,单参数取向量大小   
double distance(point3 p1,point3 p2)  
{  
    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));  
}  
  
  
///三维///   
//向量大小   
double vlen(point3 p)  
{  
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);  
}  
  
//判三点共线   
bool dots_inline(point3 p1,point3 p2,point3 p3)  
{  
    return vlen(xmult(subt(p1,p2),subt(p2,p3)))eps&&vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps;  
}  
  
//判两点在线段同侧,点在线段上返回0,不共面无意义   
bool same_side(point3 p1,point3 p2,line3 l)  
{  
    return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),sub