几何算法库.doc

const double INF=1e100;
const double ZERO=1e-6;
const double PI=2*asin(1.0);

struct XYpoint{ //(x,y)
double x;
double y;
};

struct XYline{ // Ax+By+C=0;
double A;
double B;
double C;
};

struct XYsegment{
XYpoint a,b;
};

struct XYround{ // 圆
XYpoint center;
double r;
};

inline double min(double a,double b)
{
return a<b?a:b;
}
inline double max(double a,double b)
{
return a>b?a:b;
}

/********************************************/
/* 两点确定一条直线 */
/********************************************/
XYline makeLine(double x1,double y1,double x2,double y2)
{
XYline line;
line.A=(y2-y1);
line.B=(x1-x2);
line.C=y1*(x2-x1)+x1*(y1-y2);
return line;
}

/********************************************/
/* 两点间距离 */
/********************************************/
inline double dis_PP(double x1,double y1,double x2,double y2)
{
return sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
}

/********************************************/
/* 点直线距离 */
/********************************************/
double dis_PL(double x,double y,XYline line)
{
return fabs(line.A*x + line.B*y + line.C) / sqrt(line.A*line.A + line.B*line.B);
}

/********************************************/
/* 海伦公式求三角形面积 */
/********************************************/
double triangleArea(double x1,double y1,double x2,double y2,double x3,double y3)
{
double a=sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
double b=sqrt( (x1-x3)*(x1-x3) + (y1-y3)*(y1-y3) );
double c=sqrt( (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) );
double s=(a+b+c)/2;
return sqrt( s * (s-a) * (s-b) * (s-c) );
}
inline double tArea(XYpoint p1,XYpoint p2,XYpoint p3)
{
return fabs(0.5*((p3.y-p1.y)*(p2.x-p1.x)-(p2.y-p1.y)*(p3.x-p1.x)));
}

/********************************************/
/* 判断直线相交 */
/* 1 - 有交点 0 - 无交点 -1 - 重合 */
/********************************************/
int inter_LL( XYline line1 , XYline line2 , XYpoint &point )
{
if(line1.A*line2.B == line1.B*line2.A)
{
if(line1.A*line2.C == line1.C*line2.A && line1.B*line2.C == line1.C*line2.B)
return -1;// 重合
else
return 0;// 平行
}
else
{
point.x= -1 * (line1.C*line2.B-line2.C*line1.B) / (line1.A*line2.B-line2.A*line1.B);
point.y= -1 * (line1.C*line2.A-line2.C*line1.A) / (line1.B*line2.A-line2.B*line1.A);
// 交点
return 1;
}
}

/********************************************/
/* 两直线夹角 （弧度） */
/********************************************/
double angle_LL(XYline line1 , XYline line2)
{
double t;
if( (line1.A*line2.A + line1.B*line2.B)==0 )
t=INF;
else
t=fabs( (line1.A*line2.B - line2.A*line1.B) / (line1.A*line2.A + line1.B*line2.B) );
return atan(t);
}

/********************************************/
/* 过某点和直线垂直的直线 */
/********************************************/
XYline ver_PL( XYpoint p , XYline line )
{
XYline temp;
temp.A=line.B;
temp.B=-1 * line.A;
temp.C=-temp.A*p.x-temp.B*p.y;
return temp;
}

/********************************************/
/* 多边形面积 ( v[0]=v[n] ) */
/* n */
/* a=0.5*sigma{ X[i-1]*Y-X*Y[i-1] } */
/* i=1 */
/********************************************/
double polygonArea( int vcount , XYpoint ver[] )
{
int i;
double s=0;
XYpoint t;
if( vcount<3 )
return 0;
t=ver[vcount];
ver[vcount]=ver[0];
for(i=1;i<=vcount;i++)
s+=ver[i-1].x*ver.y - ver.x*ver[i-1].y ;
ver[vcount]=t;
return s/2 ;
}

/********************************************/
/* 返回(P1-P0)*(P2-P0)的叉积。 */
/* 若为正，则<P0,P1>在<P0,P2>的顺时针方向 */
/* 若为零，则<P0,P1><P0,P2>共线； */
/* 若为负，则<P0,P1>在<P0,P2>的在逆时针方向 */
/* 这个函数可以确定两条线段在交点处的转向 */
/* 比如确定p0p1和p1p2在p1处是左转还是右转， */
/* 只要求 (p2-p0)*(p1-p0)， */
/* 若<0则左转，>0则右转，=0则共线 */
/********************************************/
double multiply( XYpoint p1 , XYpoint p2 , XYpoint p0 )
{
/*
double x0=p0.x;
double x1=p1.x;
double y0=p0.y;
double y0=p1.y;
double r;
r=((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0));
if( r<0 ) return -1;
if( r==0 ) return 0;
if( r>0 ) return 1;
*/
return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));

}

Top

回复人： galois_godel() ( ) 信誉：108 2005-7-12 21:46:16 得分: 0

/********************************************/
/* 判断直线段相交 */
/********************************************/
bool inter_SS( XYsegment u , XYsegment v )
{
return( (max(u.a.x,u.b.x)>=min(v.a.x,v.b.x))&&
(max(v.a.x,v.b.x)>=min(u.a.x,u.b.x))&&
(max(u.a.y,u.b.y)>=min(v.a.y,v.b.y))&&
(max(v.a.y,v.b.y)>=min(u.a.y,u.b.y))&&
(multiply(v.a,u.b,u.a)*multiply(u.b,v.b,u.a)>=0)&&
(multiply(u.a,v.b,v.a)*multiply(v.b,u.b,v.a)>=0));
}

/********************************************/
/* 判断两点是否共点 */
/********************************************/
bool equal_PP(XYpoint a,XYpoint b)
{
if( fabs(a.x-b.x)<ZERO && fabs(a.y-b.y)<ZERO )
return true;
else
return false;
}

/********************************************/
/* 判断点p是否在线段l上 */
/********************************************/
bool online( XYsegment l , XYpoint p )
{
if( equal_PP(p,l.a) || equal_PP(p,l.b) )
return true;
return( (multiply(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 )) );
}

/********************************************/
/* 判断点是否在多边形内 */
/********************************************/
bool insidePolygon( int vcount , XYpoint ver[] , XYpoint point )
{
int i,c=0;
XYpoint t;
XYsegment line,edge;
t=ver[vcount];
ver[vcount]=ver[0];
line.a=line.b=point;
line.b.x=INF;
for(i=0;i<vcount;i++)
{
edge.a=ver; edge.b=ver[i+1];
if( online(edge,point) )
return true;
if( ver.y==ver[i+1].y )
continue;
if( min(ver.y,ver.y)!=point.y && inter_SS(line,edge) )
c++;
}
ver[vcount]=t;
return c%2 ;
}

/********************************************/
/* 得到线段的中点 */
/********************************************/
XYpoint middle( XYpoint p1 , XYpoint p2 )
{
XYpoint point;
point.x= (p1.x+p2.x)/2;
point.y= (p1.y+p2.y)/2;
return point;
}

/********************************************/
/* 判断线段是否在多边形内 */
/* 1 - 在内部 0 - 在外部 -1 - 相交 */
/********************************************/
int segmentPolygon( int vcount , XYpoint ver[] , XYsegment seg )
{
int i;
XYpoint point,t;
XYsegment edge;
t=ver[vcount];
ver[vcount]=ver[0];
for(i=0;i<vcount;i++)
{
edge.a=ver; edge.b=ver[i+1];
if( inter_SS(seg,edge) )//&& !online(edge,seg.a) && !online(edge,seg.b) )
return -1;
}
return insidePolygon( vcount , ver , middle(seg.a,seg.b) ) ;
}

/********************************************/
/* 计算线段（p0->p1)的幅角 [0,2*PI) */
/********************************************/
double angle_PP( XYpoint p0 , XYpoint p1 )
{
double x,y;
x=p1.x-p0.x;
y=p1.y-p0.y;
if( equal_PP( p0 , p1 ) )
return -1;
if( x==0 ) {
if( y>0 )
return PI/2;
else
return PI+PI/2;
}
if( x>0 && y>=0 ) return atan(y/x);
if( x>0 && y<0 ) return 2*PI+atan(y/x);
if( x<0 && y>=0 ) return PI+atan(y/x);
if( x<0 && y<0 ) return PI+atan(y/x);
}

/********************************************/
/* 卷包裹法寻找凸包 */
/********************************************/
int convexHull( int pointNum , XYpoint pointSet[] , XYpoint convex[] )
{
int i;
int verNum;
int lastPoint,nowPoint;
double lastAngle,nowAngle,tempAngle,t1,t2;
// 找到 最右 最下的点
for(i=0,lastPoint=0;i<pointNum;i++)
{
if( pointSet.y<pointSet[lastPoint].y )
lastPoint=i;
else if( pointSet.y==pointSet[lastPoint].y &&
pointSet.x<pointSet[lastPoint].x )
lastPoint=i;
}
// 寻找凸包
lastAngle=PI;
convex[0]=pointSet[lastPoint];
verNum=1;
do {
nowAngle=lastAngle-2*PI;
for(i=0;i<pointNum;i++)
{
tempAngle=angle_PP( convex[verNum-1] , pointSet );
if( tempAngle!=-1 )
{
t1=lastAngle-tempAngle; t2=lastAngle-nowAngle;
if(t1<0) t1+=2*PI;
if(t2<0) t2+=2*PI;
if( t1<t2 ) {
nowPoint=i; nowAngle=tempAngle;
}
else if( t1==t2 ) {
if( dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet.x,pointSet.y)
<dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet[nowPoint].x,pointSet[nowPoint].y) )
nowPoint=i; nowAngle=tempAngle;
}
}
}
lastAngle=nowAngle;
convex[verNum]=pointSet[nowPoint];
verNum++;
}
while( !equal_PP(convex[0],convex[verNum-1]) ) ;
return verNum-1;
}

/********************************************/
/* 两点确定一个圆 */
/********************************************/
XYround makeRound2(XYpoint p1,XYpoint p2)
{
XYround r;
r.center=middle(p1,p2);
r.r=dis_PP(p1.x,p1.y,r.center.x,r.center.y);
return r;
}

/********************************************/
/* 三点确定一个圆 */
/********************************************/
XYround makeRound3(XYpoint p1,XYpoint p2,XYpoint p3)
{
XYround r;
XYline l1,l2;
if( fabs( (p1.x-p2.x)/(p1.y-p2.y)-(p1.x-p3.x)/(p1.y-p3.y) )<ZERO )
{
r.r=INF;
return r;
}
l1=ver_PL( middle(p1,p2), makeLine(p1.x,p1.y,p2.x,p2.y) );
l2=ver_PL( middle(p1,p3), makeLine(p1.x,p1.y,p3.x,p3.y) );
inter_LL(l1,l2,r.center);
r.r=dis_PP(p1.x,p1.y,r.center.x,r.center.y);
return r;
}

/********************************************/
/* 判断点是否在圆内 */
/********************************************/
inline bool insideRound( XYround r,XYpoint p )
{
return ( dis_PP(p.x,p.y,r.center.x,r.center.y)<=r.r );
}