计算几何_模板

#define eps 1e-10

#define infinity 1e20

struct Lpoint

{

    double x,y;

Lpoint(double x=0,double y=0):x(x),y(y){}

};

typedef Lpoint Ldir;

/*struct Ldir

{

double dx,dy;

};*/

struct Llineseg

{

    Lpoint a,b;

};

struct Llinevec{

Lpoint s;

Ldir v;

Llinevec(Lpoint s=Lpoint(),Ldir v=Lpoint()):s(s),v(v){}

};

int dcmp( double x )    //控制精度

{

    if ( fabs(x) < eps ) return 0;

    else return x < 0 ? -1 : 1;

}

//点线法

typedef Lpoint 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 );

}

Vector operator*( Vector A, double p )      //向量数乘

{

    return Vector( A.x * p, A.y * p );

}

Vector operator/( Vector A, double p )      //向量数除

{

    return Vector( A.x / p, A.y / p );

}

bool operator<(const Lpoint& A,const Lpoint& B )   //两点比较

{

    return dcmp( A.x - B.x ) < 0 || ( dcmp( A.x - B.x ) == 0 && dcmp( A.y - B.y ) < 0 );

}

bool operator==( const Lpoint& a, const Lpoint& 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 ) );

}

double Angle( Vector A, Vector B )    //向量夹角

{

    return acos( Dot(A, B) / Length(A) / Length(B) );

}

double Cross( Vector A, Vector B )   //向量叉积

{

    return A.x * B.y - A.y * B.x;

}

Vector Rotate( Vector A, double rad )    //向量旋转

{

    return Vector( A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad) );

}

Vector Normal( Vector A )    //向量单位法向量

{

    double L = Length(A);

    return Vector( -A.y / L, A.x / L );

}

Lpoint GetLineIntersection( Lpoint P, Vector v, Lpoint Q, Vector w )   //两直线交点

{

    Vector u = P - Q;

    double t = Cross( w, u ) / Cross( v, w );

    return P + v * t;

}

double DistanceToLine( Lpoint P, Lpoint A, Lpoint B )    //点到直线的距离

{

    Vector v1 = B - A, v2 = P - A;

    return fabs( Cross( v1, v2 ) ) / Length(v1);

}

double DistanceToSegment( Lpoint P, Lpoint A, Lpoint 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);

    else if ( dcmp( Dot(v1, v3) ) > 0 ) return Length(v3);

    else return fabs( Cross( v1, v2 ) ) / Length(v1);

}

Lpoint GetLineProjection( Lpoint P, Lpoint A, Lpoint B )    // 点在直线上的投影

{

    Vector v = B - A;

    return A + v*( Dot(v, P - A) / Dot( v, v ) );

}

bool SegmentProperIntersection( Lpoint a1, Lpoint a2, Lpoint b1, Lpoint 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 - b1 );

    return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;

}

bool OnSegment( Lpoint p, Lpoint a1, Lpoint a2 )   //点在线段上，不包含端点

{

    return dcmp( Cross(a1 - p, a2 - p) ) == 0 && dcmp( Dot( a1 - p, a2 - p ) ) < 0;

}

//点点法

double points_distance(Lpoint a,Lpoint b){

double mid=(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);

return sqrt(mid);

} 

double xmulti(Lpoint p1,Lpoint p2,Lpoint p0)

{

    return ((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));

}

Lpoint intersection(Llineseg u,Llineseg v){

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

}

//外心   

Lpoint circumcenter(Lpoint a,Lpoint b,Lpoint c){  

    double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2.0;  

    double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2.0;  

    double d = a1 * b2 - a2 * b1;

Lpoint ret;

ret.x=a.x + (c1*b2 - c2*b1)/d;

ret.y=a.y + (a1*c2 - a2*c1)/d;

    return ret;  

}  

//内心 

Lpoint incenter(Lpoint a,Lpoint b,Lpoint c){ 

Llineseg 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); 

}  

//垂心 

Lpoint perpencenter(Lpoint a,Lpoint b,Lpoint c){ 

Llineseg 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); 

} 

//重心 

//到三角形三顶点距离的平方和最小的点 

//三角形内到三边距离之积最大的点 

Lpoint barycenter(Lpoint a,Lpoint b,Lpoint c){ 

Llineseg 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); 

} 

//费马点 

//到三角形三顶点距离之和最小的点 

Lpoint fermentpoint(Lpoint a,Lpoint b,Lpoint c){ 

Lpoint 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 (points_distance(u,a)+points_distance(u,b)

+points_distance(u,c)>points_distance(v,a)

+points_distance(v,b)+points_distance(v,c)) 

u=v; 

} 

return u; 

} 

int ponls(Llineseg l,Lpoint p)

{

    return ((xmulti(l.b,p,l.a)==0)&&

(((p.x-l.a.x)*(p.x-l.b.x)<0)||

((p.y-l.a.y)*(p.y-l.b.y)<0)));

}

double mx(double t1,double t2)

{

    return t1>t2?t1:t2;

}

double mn(double t1,double t2)

{

    return t1<t2?t1:t2;

}

int lsinterls(Llineseg u,Llineseg v)

{

    return ((mx(u.a.x,u.b.x)>=mn(v.a.x,v.b.x))&&

(mx(v.a.x,v.b.x)>=mn(u.a.x,u.b.x))&&

(mx(u.a.y,u.b.y)>=mn(v.a.y,v.b.y))&&

(mx(v.a.y,v.b.y)>=mn(u.a.y,u.b.y))&&

(xmulti(v.a,u.b,u.a)*xmulti(u.b,v.b,u.a)>=0)&&

(xmulti(u.a,v.b,v.a)*xmulti(v.b,u.b,v.a)>=0)

);

}

int Equal_Point(Lpoint p1,Lpoint p2){

return ((fabs(p1.x-p2.x)<eps)&&(fabs(p1.y-p2.y)<eps));

}

int lsinterls_A(Llineseg u,Llineseg v){

return ((lsinterls(u,v))&&(!Equal_Point(u.a,v.a))

&&(!Equal_Point(u.b,v.b))&&(!Equal_Point(u.b,v.a))

&&(!Equal_Point(u.b,v.b)));

}

int pinplg(int vcount,Lpoint polygon[],Lpoint q){

int c=0,i,n;

Llineseg l1,l2;

l1.a=q;l1.b=q;l1.b.x=infinity;n=vcount;//

for(i=0;i<vcount;i++){

l2.a=polygon[i];

l2.b=polygon[(i+1)%n];

if((lsinterls_A(l1,l2))||

((ponls(l1,polygon[(i+1)%n]))&&

((!ponls(l1,polygon[(i+2)%n]))&&

(xmulti(polygon[i],polygon[(i+1)%n],l1.a)*

xmulti(polygon[(i+1)%n],polygon[(i+2)%n],l1.a)>0)||

(ponls(l1,polygon[(i+2)%n]))&&

(xmulti(polygon[(i)%n],polygon[(i+2)%n],l1.a)*

xmulti(polygon[(i+2)%n],polygon[(i+3)%n],l1.a)>0))))

c++;

}

return (c%2!=0);

}

//