几何算法模板

1、几何算法/*COMPUTATIONALGEOMETRYROUTINES*WRITTENBY:LIUYu(C)2003*/叉乘/两个点的距离/点到直线距离/返回直线Ax+By+C=0的系数/线段/圆/两个圆的公共面积/矩形/根据下标返回多边形的边/两个矩形的公共面积/多边形,逆时针或顺时针给出x,y/多边形顶点/多边形的边/多边形的周长/判断点是否在线段上/判断两条线断是否相交，端点重合算相交/判断两条线断是否平行/判断两条直线断是否相交/直线相交的交点/判断是否简单多边形/求多边形面积/判断是否在多边形上/判断是否在多边形内部/点阵的凸包，返回一个多边形/ 最近点对的距离#include#in

2、clude#include#include#includeusingnamespacestd;typedefdoubleTYPE; /把double 定义为TYPE#defineAbs(x)(x)0)?(x):(-(x) /用Abs()这个宏定义一个绝对值函数#defineSgn(x)(x)(b)?(a):(b) /取大值#defineMin(a,b)(a)a)and(o-b)/叉乘TYPECross(constPOINT&a,constPOINT&b,constPOINT&o)return(a.x-o.x)*(b.y-o.y)-(b.x-o.x)*(a.y-o.y); /叉积模板/plana

3、rpointsdistance/两个点的距离TYPEDistance(constPOINT&a,constPOINT&b)returnSqrt(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z);structLINEPOINTa;POINTb;LINE();LINE(POINT_a_,POINT_b_):a(_a_),b(_b_); /直线由两点决定/点到直线距离doublePointToLine(POINTp0,POINTp1,POINTp2,POINT&cp)doubled=Distance(p1,p2);doubles=C

4、ross(p1,p2,p0)/d;cp.x=p0.x+s*(p2.y-p1.y)/d;cp.y=p0.y-s*(p2.x-p1.x)/d;returnAbs(s);/返回直线Ax+By+C=0的系数voidCoefficient(constLINE&L,TYPE&A,TYPE&B,TYPE&C)A=L.b.y-L.a.y;B=L.a.x-L.b.x;C=L.b.x*L.a.y-L.a.x*L.b.y;voidCoefficient(constPOINT&p,constTYPEa,TYPE&A,TYPE&B,TYPE&C)A=Cos(a);B=Sin(a);C=-(p.y*B+p.x*A);/线

5、段structSEGPOINTa;POINTb;SEG();SEG(POINT_a_,POINT_b_):a(_a_),b(_b_);/圆structCIRCLETYPEx;TYPEy;TYPEr;CIRCLE()CIRCLE(TYPE_x_,TYPE_y_,TYPE_r_):x(_x_),y(_y_),r(_r_); /圆由圆心和半径组成POINTCenter(constCIRCLE&circle)returnPOINT(circle.x,circle.y); /返回的是一个POINT类型的对象，他这里没有直接设出一个具体的对象，而是直接用类名来实现，也是可以的，而且他这里是把对象传递进来，

6、返回的是他的圆心/圆的面积TYPEArea(constCIRCLE&circle)returnPi*circle.r*circle.r;/两个圆的公共面积TYPECommonArea(constCIRCLE&A,constCIRCLE&B)TYPEarea=0.0;constCIRCLE&M=(A.rB.r)?A:B;constCIRCLE&N=(A.rB.r)?B:A; /按照圆半径的大小来把圆区分开来TYPED=Distance(Center(M),Center(N);if(DM.r-N.r) /判断出两个圆之间的关系是相交的TYPEcosM=(M.r*M.r+D*D-N.r*N.r)/(

7、2.0*M.r*D);TYPEcosN=(N.r*N.r+D*D-M.r*M.r)/(2.0*N.r*D);TYPEalpha=2.0*ArcCos(cosM);TYPEbeta=2.0*ArcCos(cosN);TYPETM=0.5*M.r*M.r*Sin(alpha);TYPETN=0.5*N.r*N.r*Sin(beta);TYPEFM=(alpha/360.0)*Area(M);TYPEFN=(beta/360.0)*Area(N);area=FM+FN-TM-TN;/相交圆的公共部分的面积就是通过一个扇形的面积减去三角形积，然后把这两个部分相加得到elseif(D=M.r-N.r)

8、/如果两圆是内含的关系，那么公共圆的面积就是一个圆的面积area=Area(N);returnarea;/矩形/矩形的线段/2/ -b/ | |/ 3| |1/ a-/0structRECTPOINTa;/左下点POINTb;/右上点RECT();RECT(constPOINT&_a_,constPOINT&_b_)a=_a_; b=_b_;/根据下标返回多边形的边（没有看懂是干什么的，或者说是具体怎么操作的）SEGEdge(constRECT&rect,intidx) /SEG是线段SEGedge;while(idx0)idx+=4;switch(idx%4)case0:edge.a=rec

9、t.a;edge.b=POINT(rect.b.x,rect.a.y);break;case1:edge.a=POINT(rect.b.x,rect.a.y);edge.b=rect.b;break;case2:edge.a=rect.b;edge.b=POINT(rect.a.x,rect.b.y);break;case3:edge.a=POINT(rect.a.x,rect.b.y);edge.b=rect.a;break;default:break;returnedge;TYPEArea(constRECT&rect)return(rect.b.x-rect.a.x)*(rect.b.y

10、-rect.a.y);/两个矩形的公共面积TYPECommonArea(constRECT&A,constRECT&B)TYPEarea=0.0;POINTLL(Max(A.a.x,B.a.x),Max(A.a.y,B.a.y); POINTUR(Min(A.b.x,B.b.x),Min(A.b.y,B.b.y); /这里很妙，可以直接把相交部分的左下方的点与右上方的点的坐标求出来if(LL.x=UR.x)&(LL.y=UR.y) /判断是否两个矩形是否相交area=Area(RECT(LL,UR);returnarea;/多边形,逆时针或顺时针给出x,ystructPOLYintn;/n个点

11、TYPE*x;/x,y为点的指针，首尾必须重合TYPE*y;POLY():n(0),x(NULL),y(NULL);POLY(int_n_,constTYPE*_x_,constTYPE*_y_)n=_n_;x=newTYPEn+1;memcpy(x,_x_,n*sizeof(TYPE);xn=_x_0;y=newTYPEn+1;memcpy(y,_y_,n*sizeof(TYPE);yn=_y_0;/多边形顶点POINTVertex(constPOLY&poly,intidx)idx%=poly.n; /idx应该指的是点的个数returnPOINT(poly.xidx,poly.yidx)

12、; /由于POLY里面的x,y是指针，那么就相当于是关于多边形顶点的坐标的一个数组，用来记录下所有多边形顶点的坐标信息/多边形的边SEGEdge(constPOLY&poly,intidx)idx%=poly.n;returnSEG(POINT(poly.xidx,poly.yidx),POINT(poly.xidx+1,poly.yidx+1);/多边形的周长TYPEPerimeter(constPOLY&poly)TYPEp=0.0;for(inti=0;ipoly.n;i+)p=p+Distance(Vertex(poly,i),Vertex(poly,i+1);returnp;bool

13、IsEqual(TYPEa,TYPEb)return(Abs(a-b)Epsilon); /基础的用于判断两个数是否相等boolIsEqual(constPOINT&a,constPOINT&b)return(IsEqual(a.x,b.x)&IsEqual(a.y,b.y); /重载一下用于判断两个点是否相等boolIsEqual(constLINE&A,constLINE&B)TYPEA1,B1,C1;TYPEA2,B2,C2;Coefficient(A,A1,B1,C1); /把直线A的系数返回给A1,B1,C1Coefficient(B,A2,B2,C2);returnIsEqual(

14、A1*B2,A2*B1)&IsEqual(A1*C2,A2*C1)&IsEqual(B1*C2,B2*C1); /判断两条直线是否相等/判断点是否在线段上boolIsOnSeg(constSEG&seg,constPOINT&p)return(IsEqual(p,seg.a)|IsEqual(p,seg.b)|(p.x-seg.a.x)*(p.x-seg.b.x)0|(p.y-seg.a.y)*(p.y-seg.b.y)=0)&(Cross(u.a,v.b,v.a)*Cross(v.b,u.b,v.a)=0)&(Max(u.a.x,u.b.x)=Min(v.a.x,v.b.x)&(Max(v.

15、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); /后面的四个比较用于限定线段相交时的一些点坐标之间的关系，只有这样限定后才可以保证线段相交，否则如果没有这四个条件的话，只能保证的是直线的相交。/判断两条线断是否平行boolIsParallel(constLINE&A,constLINE&B)TYPEA1,B1,C1;TYPEA2,B2,C2;Coefficient(A,A1,B1,C1);Coefficient(B,A2,B2,C2);retur

16、n(A1*B2=A2*B1)&(A1*C2!=A2*C1)|(B1*C2!=B2*C1); /只要系数成比例就行/判断两条直线断是否相交boolIsIntersect(constLINE&A,constLINE&B)return!IsParallel(A,B); /不平行就相交/直线相交的交点POINTIntersection(constLINE&A,constLINE&B)TYPEA1,B1,C1;TYPEA2,B2,C2;Coefficient(A,A1,B1,C1);Coefficient(B,A2,B2,C2);POINTI(0,0);I.x=-(B2*C1-B1*C2)/(A1*B2

17、-A2*B1);I.y=(A2*C1-A1*C2)/(A1*B2-A2*B1);returnI;/判断矩形是否在圆内boolIsInCircle(constCIRCLE&circle,constRECT&rect)return(circle.x-circle.r=rect.a.x)&(circle.x+circle.r=rect.a.y)&(circle.y+circle.r=rect.b.y);/判断是否简单多边形boolIsSimple(constPOLY&poly)if(poly.n3)returnfalse;SEGL1,L2;for(inti=0;ipoly.n-1;i+)L1=Edg

18、e(poly,i); /用Edge取出多边形的边for(intj=i+1;jpoly.n;j+) /用二重循环来对多边形上的边进行一一判断L2=Edge(poly,j);if(j=i+1) /相邻的两条边进行比较if(IsOnSeg(L1,L2.b)|IsOnSeg(L2,L1.a) returnfalse;/如果第二条直线的右边的点在第一条直线上或者第一条直线的坐边的点缀第二条直线上，说明这个多边形是凹的，不是简单的elseif(j=poly.n-i-1) /同上面的一样，不过这里取出的是左边的相邻直线if(IsOnSeg(L1,L2.a)|IsOnSeg(L2,L1.b) returnfa

19、lse;elseif(IsIntersect(L1,L2) returnfalse; /不相邻的直线但是相交，那么说明非简单/forj/forireturntrue;/求多边形面积TYPEArea(constPOLY&poly)if(poly.n3)returnTYPE(0); /如果边数小于3说明非多边形，没有面积可言doubles=poly.y0*(poly.xpoly.n-1-poly.x1);for(inti=1;ipoly.n;i+)s+=poly.yi*(poly.xi-1-poly.x(i+1)%poly.n);returns/2; /把多边形分割成三角形的和，不过这里的面积计算

20、的方法没有看懂/判断点是否在多边形上boolIsOnPoly(constPOLY&poly,constPOINT&p)for(inti=0;ipoly.n;i+)if(IsOnSeg(Edge(poly,i),p) returntrue;returnfalse;/判断点是否在多边形内部boolIsInPoly(constPOLY&poly,constPOINT&p)SEGL(p,POINT(Infinity,p.y); /Infinity=1e+10，L这条直线相当于是过p点且平行于x轴的直线intcount=0;for(inti=0;iS.b.y)?(S.a):(S.b); /q取这条多边形

21、边上y坐标大的点if(IsOnSeg(L,q)+count;elseif(!IsOnSeg(L,S.a)&!IsOnSeg(L,S.b)&IsIntersect(S,L)+count;return(count%2!=0);/这里用count这个计数器的奇偶来判断点是否在多边形内没有看懂/点阵的凸包，返回一个多边形，Graham-scan算法POLYConvexHull(constPOINT*set,intn)/不适用于点少于三个的情况POINT*points=newPOINTn;memcpy(points,set,n*sizeof(POINT);TYPE*X=newTYPEn;TYPE*Y=n

22、ewTYPEn;inti,j,k=0,top=2;for(i=1;in;i+)if(pointsi.ypointsk.y)|(pointsi.y=pointsk.y)&(pointsi.xpointsk.x)k=i; /找到p0点，一般取y坐标最小，如果y相同则取x坐标最小，即最左边的点std:swap(points0,pointsk); for(i=1;in-1;i+)k=i;for(j=i+1;j0)|(Cross(pointsj,pointsk,points0)=0)&(Distance(points0,pointsj)Distance(points0,pointsk)k=j;std:swap(pointsi,pointsk);X0=points0.x;Y0=points0.y;X1=points1.x;Y1=points1.y;X2=points2.x;Y2=points2.y;for(i=3;i=0)top-;+top;Xtop=pointsi.x;Ytop=pointsi.y;deletepoints;POLYpoly(